Gpd Numerical Convergence

<execution_context>

Called from $gpd-numerical-convergence command and referenced by $gpd-verify-work for numerical phases. </purpose>

<core_principle> A numerical result without demonstrated convergence and error bars is not a result -- it is a number. Numerical validation establishes that the number means something physical by proving it is independent of numerical artifacts (grid size, time step, basis truncation, Monte Carlo statistics) to within stated precision.

The validation hierarchy:

1. Does the code run without errors? (Necessary but trivially insufficient)
2. Does it reproduce known analytical results? (Benchmark validation)
3. Does it conserve what should be conserved? (Conservation law check)
4. Does it converge as discretization is refined? (Convergence test)
5. Is it stable under perturbation? (Stability analysis)
6. Are the error bars honest? (Error estimation) </core_principle>

Load project state and conventions before beginning validation:

• Run:

INIT=$(/home/qol/.gpd/venv/bin/python -m gpd.runtime_cli --runtime codex --config-dir ./.codex --install-scope local init phase-op --include state,config)
if [ $? -ne 0 ]; then
  echo "ERROR: gpd initialization failed: $INIT"
  # STOP — display the error to the user and do not proceed.
fi

• If init succeeds (non-empty JSON with state_exists: true): Extract convention_lock for unit system (needed to verify dimensional consistency of benchmarks). Extract active approximations and their validity ranges (informs which convergence tests are most critical). Extract intermediate_results for previously computed quantities to validate.
• If init fails or state_exists is false (standalone usage): Proceed with explicit specification of the computation to validate. The user must provide the unit system and computation details directly.

Convention context is important for numerical validation: unit system determines what "reasonable" values are, and approximation validity ranges determine the expected convergence behavior.

Convention verification (if project exists):

CONV_CHECK=$(/home/qol/.gpd/venv/bin/python -m gpd.runtime_cli --runtime codex --config-dir ./.codex --install-scope local --raw convention check 2>/dev/null)
if [ $? -ne 0 ]; then
  echo "WARNING: Convention verification failed — review before validating convergence"
  echo "$CONV_CHECK"
fi

For each computation, determine:

Classify each computation:

Before testing convergence of new results, verify the code reproduces known analytical solutions.

For each benchmark:

1. Identify the benchmark:

   • Exactly solvable special case (harmonic oscillator, free particle, Ising 2D, hydrogen atom)
   • Published numerical result with high precision
   • Analytical limit (weak coupling, high temperature, large N)

2. Run the computation at high resolution for the benchmark case.

3. Compare quantitatively:

   # Template benchmark test
   computed_value = run_computation(benchmark_params)
   known_value = analytical_result(benchmark_params)
   
   abs_error = abs(computed_value - known_value)
   rel_error = abs_error / abs(known_value) if known_value != 0 else abs_error
   
   print(f"Computed:  {computed_value:.12e}")
   print(f"Known:     {known_value:.12e}")
   print(f"Abs error: {abs_error:.2e}")
   print(f"Rel error: {rel_error:.2e}")
   
   # Assessment:
   # rel_error < machine_epsilon * condition_number: PASSED
   # rel_error < tolerance: PASSED (with stated precision)
   # rel_error > tolerance: FAILED (investigate)
   

4. Record result:

   Benchmark	Known Value	Computed	Rel. Error	Status
   {name}	{value}	{value}	{error}	PASS/FAIL

If ANY benchmark fails: Stop. Debug before proceeding. Wrong results at high resolution cannot be fixed by convergence testing. </step>

Every physical system has conserved quantities. Verify they are conserved numerically.

For each conservation law:

# Template conservation check
def check_conservation(trajectory, conserved_quantity_fn, times):
    """Check that a conserved quantity stays constant."""
    Q_values = [conserved_quantity_fn(trajectory[t]) for t in times]
    Q_initial = Q_values[0]

    drift = [abs(Q - Q_initial) / abs(Q_initial) for Q in Q_values]
    max_drift = max(drift)
    avg_drift = np.mean(drift)

    print(f"Conservation of {conserved_quantity_fn.__name__}:")
    print(f"  Initial value: {Q_initial:.12e}")
    print(f"  Max drift:     {max_drift:.2e}")
    print(f"  Avg drift:     {avg_drift:.2e}")

    # Assessment:
    # Symplectic integrator: drift bounded, no secular growth
    # Non-symplectic: linear drift acceptable if small enough
    # Any method: exponential drift = UNSTABLE
    return max_drift

Record:

For each numerical parameter, vary it systematically and measure convergence.

Design the convergence study

For parameter h (grid spacing, time step, etc.), use a geometric sequence:

• h, h/2, h/4, h/8, h/16 (factor-of-2 refinement, 5 levels)
• Or h, h/3, h/9, h/27 (factor-of-3, for higher-order methods)

Minimum: 3 levels (absolute minimum for estimating convergence order) Recommended: 5 levels (robust estimation with error on the error)

Run the convergence study

import numpy as np

def convergence_study(compute_fn, param_name, param_values, reference_value=None):
    """Systematic convergence test."""
    results = []
    for p in param_values:
        value = compute_fn(**{param_name: p})
        results.append({'param': p, 'value': value})

    # Compute successive differences
    for i in range(1, len(results)):
        delta = results[i]['value'] - results[i-1]['value']
        results[i]['delta'] = delta

    # Estimate convergence order from Richardson ratios
    for i in range(2, len(results)):
        r = results[i-1]['param'] / results[i]['param']  # refinement ratio
        if abs(results[i]['delta']) > 0 and abs(results[i-1]['delta']) > 0:
            ratio = abs(results[i]['delta']) / abs(results[i-1]['delta'])
            p_est = np.log(ratio) / np.log(1/r)
            results[i]['conv_order'] = p_est

    # Richardson extrapolation (using last two values)
    if len(results) >= 2 and 'conv_order' in results[-1]:
        p = results[-1]['conv_order']
        r = results[-2]['param'] / results[-1]['param']
        extrapolated = (r**p * results[-1]['value'] - results[-2]['value']) / (r**p - 1)
        error_estimate = abs(extrapolated - results[-1]['value'])
    else:
        extrapolated = results[-1]['value']
        error_estimate = abs(results[-1]['value'] - results[-2]['value']) if len(results) >= 2 else None

    return results, extrapolated, error_estimate

Assess convergence quality

Multi-parameter convergence

When multiple parameters exist (grid spacing AND time step AND basis size), test each independently while holding others at their finest values. Then verify that the combined refinement gives consistent results.

Cross-convergence check: The result should be independent of the ORDER in which parameters are refined. If refining h first gives a different limit than refining dt first, there is a coupling between discretization errors that must be understood.

Convergence Pitfall Detection

Standard convergence testing can miss these failure modes. Check for each explicitly.

Oscillatory integrals: Standard quadrature (Gauss, Simpson) fails for integrals of the form integral f(x) exp(iomegax) dx when omega is large. The error estimate may oscillate rather than decrease monotonically.

• Detection: error estimate oscillates or fails to decrease with refinement
• Fix: Use Filon's method, Levin's method, or stationary phase approximation. For Fourier transforms, use FFT with appropriate windowing.

Stiff ODEs: Explicit integrators (RK4, Euler) require absurdly small time steps for stiff systems (those with widely separated timescales). The solution may appear to converge but to the wrong answer.

• Detection: compare implicit (BDF, RADAU) vs explicit (RK4, RK45) methods. If they give different results at the same dt, the system is stiff and the explicit method is wrong.
• Fix: Use implicit methods (scipy.integrate.solve_ivp with method='Radau' or 'BDF'). Monitor the stiffness ratio (largest/smallest eigenvalue of Jacobian).

Critical slowing down: Near phase transitions in Monte Carlo simulations, the autocorrelation time diverges as tau ~ L^z with z ~ 2 for local algorithms. Standard convergence tests on individual configurations may appear fine while ensemble averages are severely undersampled.

• Detection: autocorrelation time grows with system size. If tau_int > N_samples/100, the simulation is too short.
• Fix: Use cluster algorithms (Swendsen-Wang, Wolff) near criticality, which have z ~ 0.2-0.5. Or use parallel tempering / multicanonical methods.

Catastrophic cancellation: When computing small differences of large numbers (e.g., E_binding = E_total - E_parts), the relative error of the difference can be much larger than the relative error of the individual terms.

• Detection: relative error of difference >> relative error of summands. Compute condition number: |a + b| / |a - b|. If >> 1, catastrophic cancellation is occurring.
• Fix: Reformulate the problem to avoid the subtraction (e.g., compute binding energy directly). Use higher precision arithmetic (mpmath, float128). Use compensated summation (Kahan algorithm).

Adaptive mesh artifacts: When using adaptive mesh refinement (AMR) or adaptive quadrature, standard convergence tests (compare N and 2N points) are inapplicable because the mesh is different at each resolution level.

• Detection: compare adaptive result with uniform mesh at the finest adaptive resolution. If they disagree significantly, the adaptive algorithm may be misallocating resolution.
• Fix: Verify that error indicators drive refinement correctly. Check that the solution is smooth on the final adaptive mesh (no under-resolved features at refinement boundaries). Run at multiple overall tolerance levels and verify convergence of the final result.

Perturbation stability

Does the result change significantly under small perturbations?

# Vary initial conditions slightly
for epsilon in [1e-2, 1e-4, 1e-6, 1e-8]:
    result_perturbed = compute(initial_conditions + epsilon * random_perturbation)
    sensitivity = abs(result_perturbed - result_unperturbed) / epsilon
    print(f"eps={epsilon:.0e}  sensitivity={sensitivity:.6e}")

# If sensitivity grows with decreasing epsilon: ill-conditioned
# If sensitivity is constant: well-conditioned (Lyapunov stable)
# If sensitivity is very large: chaotic regime (need ensemble statistics)

Floating-point precision stability

# Compare different precisions
result_f32 = compute(dtype=np.float32)
result_f64 = compute(dtype=np.float64)
# result_f128 = compute(dtype=np.float128)  # if available

rel_diff_32_64 = abs(result_f32 - result_f64) / abs(result_f64)
print(f"float32 vs float64: relative diff = {rel_diff_32_64:.2e}")

# If rel_diff > sqrt(machine_epsilon_f32): catastrophic cancellation suspected
# If rel_diff ~ machine_epsilon_f32: well-conditioned

Algorithm stability

For time-stepping algorithms:

• Check CFL condition: dt < C * dx (for explicit methods)
• Check that energy is bounded (no exponential growth)
• Check that solutions remain physical (positive densities, subluminal velocities)

For iterative algorithms:

• Check that residual decreases monotonically
• Check that convergence rate matches theoretical expectation
• Flag if iteration count exceeds expected bound </step>

Construct a complete error budget for each computed quantity:

Total error: Combine in quadrature (if independent) or use maximum (if correlated):

total_error_quadrature = np.sqrt(sum(e**2 for e in errors))  # independent
total_error_conservative = sum(abs(e) for e in errors)  # worst case

Dominant error: Identify which source dominates and whether further refinement is worthwhile. </step>

---