"""Independent 1D continuum oracle; no LBM closure/header or GPU imports."""

import numpy as np

H, U, NU, DELTA = 32.0, 0.04, 0.02, 1.0
TIMES = tuple(round(v * H**2 / (np.pi**2 * NU)) for v in (0.1, 0.2, 0.4))


def solve(cs, n=256, times=TIMES, *, h=H, u0=U, nu=NU, delta=DELTA, cfl=0.4):
    """FV centered face flux plus SSPRK2, zero Dirichlet at 0,h.

    du/dt = d/dy[(nu + (cs*delta)^2*|du/dy|) du/dy].
    The oracle refinement dx=h/n MUST NOT replace the LES filter delta.
    """
    if not np.isfinite(cfl) or not 0 < cfl <= 0.4:
        raise ValueError("cfl must be finite and in (0, .4]")
    dx = h / n
    y = (np.arange(n) + 0.5) * dx
    u = u0 * np.sin(np.pi * y / h)
    a = (cs * delta) ** 2

    def gradients(v):
        g = np.empty(n + 1)
        g[1:-1] = np.diff(v) / dx
        g[0], g[-1] = 2 * v[0] / dx, -2 * v[-1] / dx
        return g

    def rhs(v):
        g = gradients(v)
        flux = (nu + a * np.abs(g)) * g
        return np.diff(flux) / dx

    t = 0.0
    out = {}
    for target in times:
        if target < t:
            raise ValueError("times must be nondecreasing and nonnegative")
        while t < target:
            # Nonlinear flux derivative is nu+2*a*abs(gradient).
            dt = min(
                target - t, cfl * dx * dx / (nu + 2 * a * np.max(np.abs(gradients(u))))
            )
            v = u + dt * rhs(u)
            u = 0.5 * u + 0.5 * (v + dt * rhs(v))
            t += dt
        out[int(target)] = u.copy()
    return y, out


def amplitude(y, u, h=H):
    s = np.sin(np.pi * y / h)
    return float(np.dot(u, s) / np.dot(s, s))


def relative(a, b):
    return float(np.linalg.norm(a - b) / np.linalg.norm(b))