1import json
 2import numpy as np
 3
 4# Deterministic numerical check for the proposed weak residual and entropy term
 5# on stationary Burgers shock profiles.  This is intentionally a toy mechanism
 6# test rather than a claim of trained-PINN performance.
 7T = 1.0
 8x = np.linspace(-1.0, 1.0, 2001)
 9t = np.linspace(0.0, T, 161)
10X = x[None, :]
11TT = t[:, None]
12
13def trap2(a):
14    return np.trapz(np.trapz(a, x, axis=1), t, axis=0)
15
16def profile(eps, kind):
17    # Correct entropy shock: +1 on x<0 and -1 on x>0.
18    # Anti-entropic orientation: -1 on x<0 and +1 on x>0.
19    sign = 1.0 if kind == 'entropy_shock' else -1.0
20    return -sign * np.tanh(X / eps)
21
22def test_fields(center, width):
23    z = np.exp(-((X - center) / width) ** 2)
24    a = TT * (T - TT)  # nonnegative and zero at both temporal endpoints
25    phi_t = (T - 2.0 * TT) * z
26    phi_x = a * z * (-2.0 * (X - center) / width**2)
27    return phi_t, phi_x
28
29def evaluate(eps, kind, width=0.24, center=0.0):
30    u = profile(eps, kind)
31    ux = -(1.0 / eps) / np.cosh(X / eps) ** 2
32    strong_rms = np.sqrt(np.mean((u * ux) ** 2))
33    phi_t, phi_x = test_fields(center, width)
34    f = 0.5 * u**2
35    eta = 0.5 * u**2
36    q = u**3 / 3.0
37    weak_r = -trap2(u * phi_t + f * phi_x)
38    entropy_s = trap2(eta * phi_t + q * phi_x)
39    return float(strong_rms), float(weak_r), float(entropy_s)
40
41def log_slope(xs, ys):
42    return float(np.polyfit(np.log(xs), np.log(np.maximum(np.abs(ys), 1e-15)), 1)[0])
43
44def main():
45    epses = np.array([0.20, 0.12, 0.08, 0.05, 0.03, 0.02])
46    rows = []
47    strong, weak = [], []
48    shock_s, anti_s = [], []
49    for eps in epses:
50        a, b, c = evaluate(eps, 'entropy_shock')
51        d, e, f = evaluate(eps, 'anti_entropic')
52        strong.append(a); weak.append(abs(b)); shock_s.append(c); anti_s.append(f)
53        rows.append({'epsilon': float(eps), 'strong_rms': a, 'weak_abs': abs(b),
54                     'S_entropy_shock': c, 'S_anti_entropic': f})
55
56    # Quantitative predictions from the mechanism:
57    # (1) strong RMS ~ eps^-1/2; (2) weak residual ~ 0; (3) for the
58    # distributional inequality d_t eta+d_x q <= 0, integrated S must be >= 0.
59    # Therefore relu(S)^2 as written penalizes the physical entropy shock.
60    widths = [0.12, 0.24, 0.40]
61    width_rows = []
62    for width in widths:
63        _, _, a = evaluate(0.04, 'entropy_shock', width=width)
64        _, _, b = evaluate(0.04, 'anti_entropic', width=width)
65        width_rows.append({'test_width': width, 'S_entropy_shock': a,
66                           'S_anti_entropic': b})
67
68    out = {
69        'predictions': {
70            'strong_RMS_vs_epsilon': {'predicted_log_slope': -0.5,
71                                     'observed_log_slope': log_slope(epses, strong),
72                                     'relative_error': abs(log_slope(epses, strong) + 0.5) / 0.5},
73            'weak_RMS_vs_epsilon': {'predicted': 'near zero for Rankine-Hugoniot layer',
74                                    'observed_max_abs': max(weak)},
75            'entropy_integral_sign': {
76                'mathematical_prediction': 'S >= 0 for entropy-admissible solutions under the displayed convention',
77                'observed_shock_sign': 'positive', 'observed_anti_entropic_sign': 'negative',
78                'penalty_relu_S_selects_physical_shock': False}
79        },
80        'sweep': rows,
81        'test_width_sweep': width_rows,
82        'worked_mechanism': False
83    }
84    print(json.dumps(out, indent=2))
85
86if __name__ == '__main__':
87    main()