1"""Invariant-domain reconstruction toy implementation for 2-D ideal-gas Euler states."""
 2import numpy as np
 3
 4GAMMA = 1.4
 5
 6def pressure(w, gamma=GAMMA):
 7    w = np.asarray(w, dtype=float)
 8    rho = w[..., 0]
 9    kinetic = 0.5 * np.sum(w[..., 1:3] ** 2, axis=-1) / rho
10    return (gamma - 1.0) * (w[..., 3] - kinetic)
11
12def admissible(w, rho_floor=1e-10, p_floor=1e-10):
13    w = np.asarray(w)
14    return bool(np.all(w[..., 0] >= rho_floor) and np.all(pressure(w) >= p_floor))
15
16def density_theta_bound(w0, wr, rho_floor=1e-10):
17    """Largest theta satisfying affine density floor, assuming w0 is admissible."""
18    r0, rr = float(w0[0]), float(wr[0])
19    if rr >= rho_floor:
20        return 1.0
21    return float(np.clip((r0-rho_floor)/(r0-rr), 0.0, 1.0))
22
23def limited_state(w0, wr, rho_floor=1e-10, p_floor=1e-10, iterations=60):
24    """Return state and largest tested theta on segment w0 + theta*(wr-w0)."""
25    w0, wr = np.asarray(w0, float), np.asarray(wr, float)
26    if not admissible(w0, rho_floor, p_floor):
27        raise ValueError("cell average must be admissible")
28    hi = density_theta_bound(w0, wr, rho_floor)
29    # Pressure positivity is checked by bisection; hi is already density-safe.
30    if admissible(w0 + hi*(wr-w0), rho_floor, p_floor):
31        return w0 + hi*(wr-w0), hi
32    lo = 0.0
33    for _ in range(iterations):
34        mid = 0.5*(lo+hi)
35        if admissible(w0 + mid*(wr-w0), rho_floor, p_floor):
36            lo = mid
37        else:
38            hi = mid
39    return w0 + lo*(wr-w0), lo
40
41def reconstruct_face(cell, raw, floor=1e-8):
42    return limited_state(cell, raw, floor, floor)
43
44def minmod(a, b):
45    return np.sign(a)*np.minimum(np.abs(a), np.abs(b)) if a*b > 0 else 0.0
46
47def scalar_reconstruction(cells, raw_faces):
48    """Toy scalar analogue: classical minmod vs hard invariant reconstruction."""
49    classical=[]; safe=[]; thetas=[]
50    for c, r in zip(cells, raw_faces):
51        # use density-like scalar component for a transparent baseline
52        slope=minmod(r[0]-c[0], c[0]-max(c[0]-0.2, 1e-5))
53        classical.append(c[0] + slope)
54        s,t=limited_state(c,r,1e-8,1e-8)
55        safe.append(s[0]); thetas.append(t)
56    return np.array(classical), np.array(safe), np.array(thetas)