57 lines
1.7 KiB
Python
57 lines
1.7 KiB
Python
"""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))
|