ClaudeandClaude Opus 5.5
10069feadb
Modal magnitude R=hypot(a,b) with amplitude error and separate phase; free-intercept ln R rate fit with per-interval rates, log-fit residual, t* units and finite-positive check; spurious-pressure naming with raw rms; INVALID_CONTROL vs CONTROL_SCREEN_FAIL (refined results kept), mixed-source check, raw ratios to both controls; docs narrow the cubic-term scope, startup and placement limitations. oracle.py still byte-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
246 lines
10 KiB
Python
246 lines
10 KiB
Python
"""CPU-only oblique shear-wave reference and composite-leaf diagnostics.
|
|
|
|
Convention: integer mode MODE=(2,1) on the L x L periodic root domain, so the
|
|
wavevector k = (2*pi/L)*(2,1) points 26.565 deg from +x (mode (1,2) would be
|
|
63.435 deg). Velocity is polarized along E = (-1,2)/sqrt(5), perpendicular to k.
|
|
The frozen Taylor-Green oracle is reused for constants and geometry only.
|
|
"""
|
|
|
|
import math
|
|
|
|
import numpy as np
|
|
from oracle import CX, CY, W, U, TIMES, centers
|
|
|
|
MODE = (2, 1)
|
|
PHI = 0.0
|
|
E = np.array([-MODE[1], MODE[0]], dtype=float) / math.hypot(*MODE)
|
|
ANGLE_DEG = math.degrees(math.atan2(MODE[1], MODE[0]))
|
|
RE = 100 # Re = U L / nu, the same nu as Taylor-Green Re100 at equal N
|
|
FIT_FROM_TSTAR = 0.12 # late-time samples only; reduces startup sensitivity
|
|
INITIALIZATION = (
|
|
"physical FP64 equilibrium G at actual cell centers; rho=1 (exact p const)"
|
|
)
|
|
# Per-row screens. phase_drift_abs, residual_l2 and the run-level rate screen
|
|
# use eps(n) and are preregistered diagnostic choices, not error bounds.
|
|
LIMITS = dict(velocity_l2=0.05, spurious_pressure_scaled=0.02, mass_relative=1e-5)
|
|
# With equal absolute ceiling and floor (0.02) the comparative pressure budget
|
|
# adds no constraint unless the root control exceeds 0.01; raw ratios are kept.
|
|
FLOORS = dict(rate_error=0.002, residual_l2=0.002, spurious_pressure_scaled=0.02)
|
|
BASELINE_FACTOR = 2.0
|
|
|
|
|
|
def wavevector(length):
|
|
return 2 * np.pi / length * np.array(MODE, dtype=float)
|
|
|
|
|
|
def decay_rate(length, viscosity):
|
|
"""Exact amplitude decay rate per root step: nu |k|^2."""
|
|
k = wavevector(length)
|
|
return viscosity * float(k @ k)
|
|
|
|
|
|
def eps(n):
|
|
"""(|k| h_root)^2 with h_root = 1: dimensionless, coefficient one."""
|
|
k = wavevector(n)
|
|
return float(k @ k)
|
|
|
|
|
|
def exact(x, y, t, length, viscosity, amplitude=U):
|
|
k = wavevector(length)
|
|
a = amplitude * np.exp(-decay_rate(length, viscosity) * t)
|
|
s = np.sin(k[0] * x + k[1] * y + PHI)
|
|
p = np.zeros_like(s)
|
|
return a * E[0] * s, a * E[1] * s, p # kinematic pressure, constant
|
|
|
|
|
|
def equilibrium(rho, ux, uy):
|
|
# Same second-order D2Q9 equilibrium as oracle.initial.
|
|
cu = CX[:, None, None] * ux + CY[:, None, None] * uy
|
|
return (
|
|
W[:, None, None]
|
|
* rho
|
|
* (1 + 3 * cu + 4.5 * cu * cu - 1.5 * (ux * ux + uy * uy))
|
|
)
|
|
|
|
|
|
def initial(shape, level, length, viscosity):
|
|
x, y = centers(shape, level)
|
|
ux, uy, p = exact(x, y, 0.0, length, viscosity)
|
|
return equilibrium(1 + 3 * p, ux, uy)
|
|
|
|
|
|
def _leaves(arrays, masks):
|
|
accum = []
|
|
for level, (g, mask) in enumerate(zip(arrays, masks)):
|
|
if g.shape != (9, *mask.shape):
|
|
raise ValueError("population/mask shape mismatch")
|
|
values = g[:, mask].astype(float)
|
|
if not np.isfinite(values).all():
|
|
raise FloatingPointError("nonfinite leaf populations")
|
|
rho = values.sum(axis=0)
|
|
if (rho <= 0).any():
|
|
raise FloatingPointError("nonpositive leaf density")
|
|
x, y = centers(mask.shape, level)
|
|
area = np.full(rho.size, 0.25**level)
|
|
lev = np.full(rho.size, level)
|
|
accum.append((rho, CX @ values / rho, CY @ values / rho, x[mask], y[mask], area, lev))
|
|
return [np.concatenate(v) for v in zip(*accum)]
|
|
|
|
|
|
def metrics(arrays, masks, t, length, viscosity, initial_mass=None):
|
|
rho, ux, uy, x, y, w, lev = _leaves(arrays, masks)
|
|
volume = w.sum()
|
|
if not np.isclose(volume, length**2, rtol=0, atol=1e-10):
|
|
raise ValueError("leaf partition does not cover domain once")
|
|
k = wavevector(length)
|
|
theta = k[0] * x + k[1] * y + PHI
|
|
s, c = np.sin(theta), np.cos(theta)
|
|
v = ux * E[0] + uy * E[1]
|
|
# Weighted least squares on the actual composite quadrature (no
|
|
# continuous-orthogonality assumption).
|
|
gram = np.array([[w @ (s * s), w @ (s * c)], [w @ (s * c), w @ (c * c)]])
|
|
a_hat, b_hat = np.linalg.solve(gram, [w @ (v * s), w @ (v * c)])
|
|
ex, ey, _ = exact(x, y, t, length, viscosity)
|
|
exact_norm = w @ (ex * ex + ey * ey)
|
|
fit = a_hat * s + b_hat * c
|
|
longitudinal = (ux * k[0] + uy * k[1]) / math.sqrt(k @ k)
|
|
amplitude = U * math.exp(-decay_rate(length, viscosity) * t)
|
|
scale = amplitude * amplitude / 4 # continuum reference scale, not the TG norm
|
|
modal = math.hypot(a_hat, b_hat)
|
|
p = (rho - 1) / 3
|
|
p = p - w @ p / volume # one global gauge
|
|
mass = float(w @ rho)
|
|
return dict(
|
|
tstar=t * U / length,
|
|
mass=mass,
|
|
mass_relative=(
|
|
0.0 if initial_mass is None else abs(mass - initial_mass) / initial_mass
|
|
),
|
|
a_hat=float(a_hat),
|
|
b_hat=float(b_hat),
|
|
amplitude_exact=amplitude,
|
|
amplitude_ratio=float(a_hat / amplitude),
|
|
modal_amplitude=modal,
|
|
amplitude_error=modal / amplitude - 1,
|
|
phase_drift=float(math.atan2(b_hat, a_hat)),
|
|
phase_drift_abs=abs(float(math.atan2(b_hat, a_hat))),
|
|
velocity_l2=float(
|
|
np.sqrt(w @ ((ux - ex) ** 2 + (uy - ey) ** 2) / exact_norm)
|
|
),
|
|
residual_l2=float(
|
|
np.sqrt(w @ ((ux - fit * E[0]) ** 2 + (uy - fit * E[1]) ** 2) / exact_norm)
|
|
),
|
|
longitudinal_l2=float(np.sqrt(w @ (longitudinal**2) / exact_norm)),
|
|
spurious_pressure_rms=float(np.sqrt(w @ (p * p) / volume)),
|
|
spurious_pressure_scaled=float(np.sqrt(w @ (p * p) / volume) / scale),
|
|
spurious_pressure_max_scaled=float(np.max(np.abs(p)) / scale),
|
|
spurious_pressure_scaled_by_level=[
|
|
float(np.sqrt(w[lev == i] @ (p[lev == i] ** 2) / w[lev == i].sum()) / scale)
|
|
for i in range(len(masks))
|
|
],
|
|
rho_min=float(rho.min()),
|
|
rho_max=float(rho.max()),
|
|
max_speed=float(np.max(np.hypot(ux, uy))),
|
|
leaf_cells=int(len(rho)),
|
|
)
|
|
|
|
|
|
def fit_rate(rows, length, viscosity):
|
|
"""Free-intercept least squares of ln R, R = hypot(a_hat, b_hat), late samples.
|
|
|
|
The intercept absorbs only a constant amplitude prefactor once a single
|
|
exponential dominates; per-interval rates and the fit residual are kept so
|
|
a lingering transient stays visible. Rates are per root step of the record
|
|
(uniform 2N uses its own nu and doubled step counts).
|
|
"""
|
|
late = [r for r in rows if r["tstar"] >= FIT_FROM_TSTAR - 1e-12]
|
|
amp = np.array([math.hypot(r["a_hat"], r["b_hat"]) for r in late])
|
|
if len(late) < 3 or not (np.isfinite(amp).all() and (amp > 0).all()):
|
|
raise ValueError("need >= 3 late samples with finite positive R")
|
|
t = np.array([r["step"] for r in late], dtype=float)
|
|
(slope, intercept) = np.polyfit(t, np.log(amp), 1)
|
|
residual = np.log(amp) - (slope * t + intercept)
|
|
exact_rate = decay_rate(length, viscosity)
|
|
to_tstar = length / U # root steps per unit t*
|
|
return dict(
|
|
rate=float(-slope), exact_rate=exact_rate,
|
|
rate_tstar=float(-slope * to_tstar),
|
|
exact_rate_tstar=exact_rate * to_tstar, # = 5 (2 pi)^2 / Re
|
|
rate_error=abs(-slope / exact_rate - 1),
|
|
interval_rate_errors=[
|
|
float(-np.log(amp[i + 1] / amp[i]) / (t[i + 1] - t[i]) / exact_rate - 1)
|
|
for i in range(len(t) - 1)
|
|
],
|
|
log_fit_residual_max=float(np.max(np.abs(residual))),
|
|
)
|
|
|
|
|
|
def summarize(record):
|
|
rows = record["rows"]
|
|
n = record["n"]
|
|
out = fit_rate(rows, n, record["viscosity"])
|
|
out.update(
|
|
eps=eps(n),
|
|
complete=record["completed"] and len(rows) == len(TIMES),
|
|
residual_l2=max(r["residual_l2"] for r in rows),
|
|
phase_drift_abs=max(r["phase_drift_abs"] for r in rows),
|
|
spurious_pressure_scaled=max(r["spurious_pressure_scaled"] for r in rows),
|
|
spurious_pressure_rms=max(r["spurious_pressure_rms"] for r in rows),
|
|
amplitude_error_abs=max(abs(r["amplitude_error"]) for r in rows),
|
|
velocity_l2=max(r["velocity_l2"] for r in rows),
|
|
mass_relative=max(r["mass_relative"] for r in rows),
|
|
)
|
|
out["absolute_pass"] = out["complete"] and all(
|
|
out[key] <= out["eps"] for key in ("rate_error", "phase_drift_abs", "residual_l2")
|
|
) and all(out[key] <= limit for key, limit in LIMITS.items())
|
|
return out
|
|
|
|
|
|
def _usable(record):
|
|
try:
|
|
out = summarize(record)
|
|
except (KeyError, ValueError):
|
|
return None
|
|
finite = all(math.isfinite(v) for v in out.values() if isinstance(v, float))
|
|
return out if out["complete"] and finite else None
|
|
|
|
|
|
def assess(base_control, fine_control, refined):
|
|
"""Uniform root-spacing and finest-spacing controls versus refined runs.
|
|
|
|
INVALID_CONTROL: missing, incomplete, nonfinite or mixed-source controls.
|
|
CONTROL_SCREEN_FAIL: complete finite controls that fail a screen; refined
|
|
results are still computed and kept.
|
|
"""
|
|
if base_control is None or fine_control is None:
|
|
return dict(status="INVALID_CONTROL", reason="missing control")
|
|
n = base_control["n"]
|
|
if (base_control["depth"], fine_control["depth"], fine_control["n"]) != (1, 1, 2 * n):
|
|
raise ValueError("controls must be uniform N and uniform 2N")
|
|
if any(r["depth"] != 2 or r["n"] != n for r in refined):
|
|
raise ValueError("refined runs must be depth 2 at root N")
|
|
sources = {repr(r.get("source")) for r in [base_control, fine_control, *refined]}
|
|
base, fine = _usable(base_control), _usable(fine_control)
|
|
if base is None or fine is None or len(sources) != 1:
|
|
return dict(status="INVALID_CONTROL",
|
|
reason="incomplete/nonfinite control or mixed source")
|
|
report = dict(base_control=base, fine_control=fine, refined={}, status="PASS")
|
|
if not (base["absolute_pass"] and fine["absolute_pass"]):
|
|
report["status"] = "CONTROL_SCREEN_FAIL"
|
|
for r in refined:
|
|
s = summarize(r)
|
|
s["budget"] = {
|
|
key: max(BASELINE_FACTOR * base[key], floor) for key, floor in FLOORS.items()
|
|
}
|
|
s["budget_pass"] = all(s[key] <= s["budget"][key] for key in FLOORS)
|
|
s["control_ratios"] = {
|
|
name: {key: (s[key] / c[key] if c[key] > 0 else None)
|
|
for key in (*FLOORS, "spurious_pressure_rms")}
|
|
for name, c in (("root", base), ("finest", fine))
|
|
}
|
|
s["pass"] = s["absolute_pass"] and s["budget_pass"]
|
|
report["refined"][r["explosion"]] = s
|
|
if report["status"] == "PASS" and not s["pass"]:
|
|
report["status"] = "FAIL"
|
|
return report
|