MONAD SERIES 11 — N* IDENTIFICATION & DIAGNOSTICS PIPELINE (FIXED) AND SOLVER
#!/usr/bin/env python3
"""
================================================================================
MONAD SERIES 11 — N* IDENTIFICATION & DIAGNOSTICS PIPELINE (FIXED)
================================================================================
FIXED: Initialization-order bug - diagnostics now called AFTER init_gaussian()
No NoneType errors in decompose_pi()
================================================================================
"""
import sys
import time as time_module
import json
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from pathlib import Path
from datetime import datetime
from typing import Dict, List, Optional, Tuple
import warnings
warnings.filterwarnings('ignore')
# ==============================================================================
# CRITICAL FIX: IMPORT SOLVER FROM DISK
# ==============================================================================
sys.path.insert(0, "/content")
from monad_series11_pure import MonadSolver11Pure
print("✅ MonadSolver11Pure imported successfully from /content/")
# ==============================================================================
# ENHANCED SOLVER WITH DIAGNOSTICS (FIXED INITIALIZATION ORDER)
# ==============================================================================
class MonadSolver11Diagnostics(MonadSolver11Pure):
"""
Extended Series 11 solver with additional diagnostics for N* identification.
FIXED: Removed premature _collect_diagnostics from __init__
Now called AFTER init_gaussian() via overridden method
"""
def __init__(self, config):
super().__init__(config)
# Initialize diagnostic history
self.diagnostic_history = {
'steps': [],
'H_total': [],
'H_kinetic': [],
'H_gradient': [],
'H_potential': [],
'metric_components': [],
'constitutive_residuals': [],
'field_stats': [],
'energy_balance': []
}
# Initialize H0 tracking (will be set after init_gaussian)
self.H0 = None
self.H0_components = None
self._diagnostics_initialized = False
def init_gaussian(self, amplitude=2.5, width=4.0):
"""
Override init_gaussian to call parent first, then collect diagnostics.
This fixes the NoneType bug by ensuring fields exist before diagnostics.
"""
# 1. Call parent to allocate Pxx, Pxy, Pyy arrays
super().init_gaussian(amplitude=amplitude, width=width)
# 2. Now that fields are initialized, collect step 0 invariants safely
self._collect_diagnostics(initial=True)
self._diagnostics_initialized = True
return self
def _collect_diagnostics(self, initial=False):
"""
Collect all diagnostics at current step.
FIXED: Added safety gate to prevent NoneType errors.
"""
# Safety gate: ensure fields exist before computing
if self.Pxx is None or self.Pxy is None or self.Pyy is None:
print("⚠️ WARNING: Fields not initialized, skipping diagnostics")
return None
# Hamiltonian components
H_comps = self._compute_hamiltonian_components()
# Store initial values if this is the first call
if initial or self.H0 is None:
self.H0 = H_comps['H_total']
self.H0_components = H_comps
# Reconstructed metric
metric = self._compute_reconstructed_metric()
# Constitutive residuals
residuals = self._compute_constitutive_residuals()
# Field statistics
field_stats = self._compute_field_stats()
# Energy balance
energy_balance = self._compute_energy_balance()
# Only store in history if not initial
if not initial:
self.diagnostic_history['steps'].append(self.step)
self.diagnostic_history['H_total'].append(H_comps['H_total'])
self.diagnostic_history['H_kinetic'].append(H_comps['H_kinetic'])
self.diagnostic_history['H_gradient'].append(H_comps['H_gradient'])
self.diagnostic_history['H_potential'].append(H_comps['H_potential'])
self.diagnostic_history['metric_components'].append(metric)
self.diagnostic_history['constitutive_residuals'].append(residuals)
self.diagnostic_history['field_stats'].append(field_stats)
self.diagnostic_history['energy_balance'].append(energy_balance)
return {
'H_components': H_comps,
'metric': metric,
'residuals': residuals,
'field_stats': field_stats,
'energy_balance': energy_balance,
'is_initial': initial
}
def get_diagnostics(self):
"""Public method to collect diagnostics at current step."""
if not self._diagnostics_initialized:
print("⚠️ WARNING: Diagnostics called before init_gaussian()")
return None
return self._collect_diagnostics(initial=False)
def _compute_hamiltonian_components(self):
"""Compute Hamiltonian components separately."""
dx = self.dx
# Get constitutive results
results = self._evaluate_monad_from_pi(
self.Pxx, self.Pxy, self.Pyy,
self.Pi_max, self.anchor,
self.kappa, self.eta,
self.beta, self.gamma,
self.m2, self.alpha, self.delta
)
S = results['S']
Psi = results['Psi']
Lam = results['Lam']
# Kinetic energy
H_kin = 0.5 * (self.Uxx**2 + self.Uxy**2 + self.Uyy**2)
# Gradient energy
gPxx_x, gPxx_y = self._cpu_grad_2d(self.Pxx, dx)
gPxy_x, gPxy_y = self._cpu_grad_2d(self.Pxy, dx)
gPyy_x, gPyy_y = self._cpu_grad_2d(self.Pyy, dx)
H_grad = 0.5 * self.c2 * (gPxx_x**2 + gPxx_y**2 +
gPxy_x**2 + gPxy_y**2 +
gPyy_x**2 + gPyy_y**2)
# Potential energy
ps = Psi**2
H_pot = (0.5 * self.beta * S**2 + 0.25 * self.gamma * S**4 +
0.5 * self.m2 * ps + 0.5 * self.alpha * Lam**2 + 0.25 * self.delta * Lam**4 +
self.kappa * S * ps + self.eta * ps * Lam)
# Total
H_total = H_kin + H_grad + H_pot
return {
'H_total': float(np.sum(H_total) * dx**2),
'H_kinetic': float(np.sum(H_kin) * dx**2),
'H_gradient': float(np.sum(H_grad) * dx**2),
'H_potential': float(np.sum(H_pot) * dx**2),
}
def _compute_reconstructed_metric(self):
"""Compute reconstructed metric g(Π) from primitive tensor."""
results = self._evaluate_monad_from_pi(
self.Pxx, self.Pxy, self.Pyy,
self.Pi_max, self.anchor,
self.kappa, self.eta,
self.beta, self.gamma,
self.m2, self.alpha, self.delta
)
S = results['S']
Psi = results['Psi']
Lam = results['Lam']
Psi_val = results['Psi_val']
metric = {
'g_tt': float(np.mean(Psi_val)),
'g_xx': float(np.mean(S)),
'g_yy': float(np.mean(Lam)),
'g_xy': float(np.mean(Psi)),
'g_tt_std': float(np.std(Psi_val)),
'g_xx_std': float(np.std(S)),
'g_yy_std': float(np.std(Lam)),
'g_xy_std': float(np.std(Psi)),
'g_det': float(np.mean(S * Lam - Psi**2)),
}
return metric
def _compute_constitutive_residuals(self):
"""Compute residuals of the constitutive map."""
results = self._evaluate_monad_from_pi(
self.Pxx, self.Pxy, self.Pyy,
self.Pi_max, self.anchor,
self.kappa, self.eta,
self.beta, self.gamma,
self.m2, self.alpha, self.delta
)
D_drift = results['D_drift']
return {
'D_drift_mean': float(np.mean(D_drift)),
'D_drift_std': float(np.std(D_drift)),
'D_drift_max': float(np.max(D_drift)),
'D_drift_min': float(np.min(D_drift)),
}
def _compute_field_stats(self):
"""Compute statistics of the primitive tensor field Π."""
return {
'Pxx_mean': float(np.mean(self.Pxx)),
'Pxx_std': float(np.std(self.Pxx)),
'Pxx_max': float(np.max(self.Pxx)),
'Pxx_min': float(np.min(self.Pxx)),
'Pxy_mean': float(np.mean(self.Pxy)),
'Pxy_std': float(np.std(self.Pxy)),
'Pxy_max': float(np.max(self.Pxy)),
'Pxy_min': float(np.min(self.Pxy)),
'Pyy_mean': float(np.mean(self.Pyy)),
'Pyy_std': float(np.std(self.Pyy)),
'Pyy_max': float(np.max(self.Pyy)),
'Pyy_min': float(np.min(self.Pyy)),
'Uxx_mean': float(np.mean(self.Uxx)),
'Uxx_std': float(np.std(self.Uxx)),
'Uxx_max': float(np.max(self.Uxx)),
'Uxx_min': float(np.min(self.Uxx)),
'Uxy_mean': float(np.mean(self.Uxy)),
'Uxy_std': float(np.std(self.Uxy)),
'Uxy_max': float(np.max(self.Uxy)),
'Uxy_min': float(np.min(self.Uxy)),
'Uyy_mean': float(np.mean(self.Uyy)),
'Uyy_std': float(np.std(self.Uyy)),
'Uyy_max': float(np.max(self.Uyy)),
'Uyy_min': float(np.min(self.Uyy)),
}
def _compute_energy_balance(self):
"""Compute energy balance components."""
comps = self._compute_hamiltonian_components()
if self.H0_components is not None:
dH_total = comps['H_total'] - self.H0_components['H_total']
dH_kinetic = comps['H_kinetic'] - self.H0_components['H_kinetic']
dH_gradient = comps['H_gradient'] - self.H0_components['H_gradient']
dH_potential = comps['H_potential'] - self.H0_components['H_potential']
else:
dH_total = 0.0
dH_kinetic = 0.0
dH_gradient = 0.0
dH_potential = 0.0
return {
'dH_total': float(dH_total),
'dH_kinetic': float(dH_kinetic),
'dH_gradient': float(dH_gradient),
'dH_potential': float(dH_potential),
'dH_total_rel': float(dH_total / max(abs(self.H0_components['H_total']), 1e-30)) if self.H0_components else 0.0,
}
def _evaluate_monad_from_pi(self, *args, **kwargs):
"""Wrapper for evaluate_monad_from_pi."""
from monad_series11_pure import evaluate_monad_from_pi
return evaluate_monad_from_pi(*args, **kwargs)
def _cpu_grad_2d(self, *args, **kwargs):
"""Wrapper for cpu_grad_2d."""
from monad_series11_pure import cpu_grad_2d
return cpu_grad_2d(*args, **kwargs)
def run_with_diagnostics(self, steps, diag_interval=1000, drift_threshold=0.05):
"""Run simulation with comprehensive diagnostics."""
# Ensure diagnostics are initialized
if not self._diagnostics_initialized:
raise RuntimeError("Must call init_gaussian() before run_with_diagnostics()")
# Memory gate for large N runs
if self.N >= 256:
print(f"⚠️ MEMORY GATE: N={self.N} - Large grid detected")
print(f" Estimated memory usage: ~{self.N**2 * 6 * 8 / 1024 / 1024:.1f} MB for fields")
print(f" Total cells: {self.N**2:,}")
print(f"\n{'='*60}")
print(f"SERIES 11 — DIAGNOSTICS SOLVER")
print(f"{'='*60}")
print(f"Grid: {self.N}×{self.N}, dt={self.dt:.6f}, κ={self.kappa:.3f}")
print(f"Diagnostics: every {diag_interval} steps")
print(f"{'-'*60}")
start = time_module.time()
# Initial diagnostics already collected in init_gaussian()
H0 = self.H0
print(f"Step {0:6d} | H_total = {H0:.8e}")
print(f" | H_kin = {self.H0_components['H_kinetic']:.8e}")
print(f" | H_grad = {self.H0_components['H_gradient']:.8e}")
print(f" | H_pot = {self.H0_components['H_potential']:.8e}")
for step in range(1, steps+1):
self.rk3_step()
if step % diag_interval == 0:
diag = self.get_diagnostics()
if diag is None:
print("⚠️ WARNING: Diagnostics returned None, skipping")
continue
H = diag['H_components']['H_total']
dr = abs(H - H0) / max(abs(H0), 1e-30)
print(f"Step {step:6d} | H_total={H:.8e} | ΔH/H={dr:.4e}")
print(f" | H_kin={diag['H_components']['H_kinetic']:.8e}")
print(f" | H_grad={diag['H_components']['H_gradient']:.8e}")
print(f" | H_pot={diag['H_components']['H_potential']:.8e}")
print(f" | g_det={diag['metric']['g_det']:.4e}")
print(f" | D_drift={diag['residuals']['D_drift_mean']:.4e}")
if np.isnan(H) or dr > drift_threshold:
print(f"❌ INSTABILITY at step {step}")
break
elapsed = time_module.time() - start
Hf = self.diagnostic_history['H_total'][-1] if self.diagnostic_history['H_total'] else 0
dr = abs(Hf - H0) / max(abs(H0), 1e-30)
print(f"{'-'*60}")
print(f"Completed {len(self.diagnostic_history['steps'])} diagnostic intervals in {elapsed:.2f}s")
print(f"{'='*60}")
return {
'steps': step,
'H0': float(H0),
'H_final': float(Hf),
'drift_rel': float(dr),
'elapsed': float(elapsed),
'diagnostic_history': self.diagnostic_history
}
# ==============================================================================
# AUTOMATED CHARTING FUNCTIONS
# ==============================================================================
def generate_charts(results, output_dir):
"""Generate automated charts from pipeline results."""
output_path = Path(output_dir)
output_path.mkdir(exist_ok=True)
# Extract data by N
data_by_N = {}
for r in results:
N = r['N']
if N not in data_by_N:
data_by_N[N] = []
data_by_N[N].append(r)
# Skip if no valid runs
valid_runs = [r for r in results if not r.get('failed', False)]
if not valid_runs:
print("⚠️ No valid runs to chart")
return
# 1. Drift vs N
fig, ax = plt.subplots(figsize=(10, 6))
Ns = sorted(data_by_N.keys())
drifts = []
errors = []
for N in Ns:
runs = data_by_N[N]
d = [r['drift_pct'] for r in runs if not r.get('failed', False)]
if d:
drifts.append(np.mean(d))
errors.append(np.std(d) if len(d) > 1 else 0)
else:
drifts.append(None)
errors.append(0)
valid_N = [Ns[i] for i in range(len(Ns)) if drifts[i] is not None]
valid_drifts = [d for d in drifts if d is not None]
valid_errors = [errors[i] for i in range(len(Ns)) if drifts[i] is not None]
if valid_N:
ax.errorbar(valid_N, valid_drifts, yerr=valid_errors, fmt='o-', capsize=5, markersize=8)
ax.set_xlabel('Grid Size N')
ax.set_ylabel('Energy Drift (%)')
ax.set_title('Hamiltonian Conservation vs Grid Resolution')
ax.grid(True, alpha=0.3)
ax.set_xscale('log')
ax.set_yscale('log')
fig.savefig(output_path / 'drift_vs_N.png', dpi=150, bbox_inches='tight')
plt.close()
# 2. Kinetic vs Gradient Energy
fig, ax = plt.subplots(figsize=(10, 6))
for N in Ns[:4]:
runs = data_by_N[N]
for r in runs:
if r.get('diagnostic_history'):
hist = r['diagnostic_history']
steps = hist.get('steps', [])
H_kin = hist.get('H_kinetic', [])
H_grad = hist.get('H_gradient', [])
if steps and H_kin and H_grad:
ax.plot(steps, H_kin, '--', label=f'N={N} (Kinetic)', alpha=0.7)
ax.plot(steps, H_grad, '-', label=f'N={N} (Gradient)', alpha=0.7)
ax.set_xlabel('Steps')
ax.set_ylabel('Energy')
ax.set_title('Kinetic vs Gradient Energy Evolution')
ax.legend()
ax.grid(True, alpha=0.3)
fig.savefig(output_path / 'kinetic_vs_gradient.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"✅ Charts saved to {output_path}/")
# ==============================================================================
# MARKDOWN REPORT GENERATOR
# ==============================================================================
def generate_markdown_report(results, timestamp):
"""Generate automated Markdown report."""
report_lines = []
report_lines.append("# MONAD SERIES 11 — N* IDENTIFICATION REPORT")
report_lines.append("")
report_lines.append(f"**Generated:** {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
report_lines.append("")
# Summary stats
total = len(results)
failed = len([r for r in results if r.get('failed', False)])
success = total - failed
report_lines.append("## 📊 Overall Statistics")
report_lines.append("")
report_lines.append("| Metric | Value |")
report_lines.append("|--------|-------|")
report_lines.append(f"| Total runs | {total} |")
report_lines.append(f"| Successful | {success} ({success/total*100:.1f}%) |")
report_lines.append(f"| Failed | {failed} ({failed/total*100:.1f}%) |")
report_lines.append("")
# Results by N
report_lines.append("## 📐 N* Sweep Results")
report_lines.append("")
report_lines.append("| N | Steps | Drift % | Status |")
report_lines.append("|---|-------|---------|--------|")
for r in results:
if 'N_star' in r.get('label', ''):
status = "✅ OK" if not r['failed'] else f"❌ {r['failure_step']:,}"
report_lines.append(f"| {r['N']} | {r['steps']:,} | {r['drift_pct']:.4f} | {status} |")
report_lines.append("")
# Check if any valid runs
if success == 0:
report_lines.append("⚠️ **No successful runs completed.**")
report_lines.append("")
report_lines.append("All runs failed during initialization due to the NoneType bug.")
report_lines.append("The fix has been applied - re-run the pipeline.")
report_lines.append("")
report_lines.append("---")
report_lines.append("*Report generated automatically by Monad Series 11 N* Pipeline*")
return "\n".join(report_lines)
# ==============================================================================
# SINGLE RUN EXECUTOR
# ==============================================================================
def run_single_run(config, target_steps=25000, diag_interval=1000,
drift_threshold=0.05, label="", verbose=True):
"""Run a single configuration with diagnostics."""
if verbose:
print(f"\n▶️ Running: κ={config['kappa']:.3f} | N={config['N']} | dt={config['dt']:.1e} | steps={target_steps:,}")
if label:
print(f" Label: {label}")
print("-"*60)
start_time = time_module.perf_counter()
try:
# Create solver
solver = MonadSolver11Diagnostics(config)
# CRITICAL FIX: Initialize fields BEFORE diagnostics
solver.init_gaussian(amplitude=2.5)
# Now fields exist, diagnostics are collected in init_gaussian
# Run with diagnostics
result_metrics = solver.run_with_diagnostics(
steps=target_steps,
diag_interval=diag_interval,
drift_threshold=drift_threshold
)
elapsed = time_module.perf_counter() - start_time
H0 = result_metrics.get('H0', 0.0)
Hf = result_metrics.get('H_final', 0.0)
drift_rel = result_metrics.get('drift_rel', 0.0)
completed_steps = result_metrics.get('steps', target_steps)
diag_history = result_metrics.get('diagnostic_history', {})
physics_ok = drift_rel <= drift_threshold
result = {
'kappa': config['kappa'],
'N': config['N'],
'dt': config['dt'],
'steps': completed_steps,
'steps_target': target_steps,
'drift_rel': float(drift_rel),
'drift_pct': float(drift_rel * 100),
'H0': float(H0),
'H_final': float(Hf),
'elapsed': elapsed,
'physics_ok': physics_ok,
'failed': not physics_ok,
'failure_step': completed_steps if not physics_ok else None,
'failure_drift': drift_rel if not physics_ok else None,
'label': label,
'diagnostic_history': diag_history,
'timestamp': datetime.now().isoformat()
}
if verbose:
status = "❌ FAILED" if result['failed'] else "✅ SUCCESS"
print(f"{status} | Steps: {completed_steps:,} | Drift: {drift_rel*100:.4f}%")
print(f" Time: {elapsed:.2f}s")
return result
except Exception as e:
elapsed = time_module.perf_counter() - start_time
print(f"❌ Run failed with exception: {e}")
import traceback
traceback.print_exc()
return {
'kappa': config['kappa'],
'N': config['N'],
'dt': config['dt'],
'steps': 0,
'steps_target': target_steps,
'drift_rel': 1.0,
'drift_pct': 100.0,
'H0': 0,
'H_final': 0,
'elapsed': elapsed,
'physics_ok': False,
'failed': True,
'failure_step': 0,
'failure_drift': 1.0,
'label': label,
'diagnostic_history': {},
'error': str(e)
}
# ==============================================================================
# TEST SETS
# ==============================================================================
def run_n_star_sweep(base_config, all_results, verbose=True):
"""TEST 1: N* Sweep (N = [16, 32, 64, 128, 256])"""
print("\n" + "="*80)
print("🔬 TEST 1: N* Sweep")
print("="*80)
print("Goal: Find finite-resolution saturation point N*")
print("Config: κ=0.50, N=[16,32,64,128,256], 100K steps")
print("="*80)
N_values = [16, 32, 64, 128, 256]
results = []
for N in N_values:
config = base_config.copy()
config['N'] = N
config['dx'] = 0.4 * (64.0 / N)
config['kappa'] = 0.50
if N >= 256:
print(f"\n⚠️ Large N={N} detected - memory usage will be high")
print(f" Cells: {N**2:,}, Fields: ~{N**2 * 6 * 8 / 1024 / 1024:.1f} MB")
diag_interval = max(100, int(100000 / 50))
result = run_single_run(
config=config,
target_steps=100000,
diag_interval=diag_interval,
drift_threshold=0.05,
label=f"N_star_{N}",
verbose=verbose
)
all_results.append(result)
results.append(result)
print("\n📊 N* Sweep Results:")
print("-"*40)
print(f"{'N':>6} | {'Steps':>10} | {'Drift %':>10} | {'Status':>10}")
print("-"*40)
for r in results:
status = "✅ OK" if not r['failed'] else f"❌ {r['failure_step']:,}"
print(f"{r['N']:>6} | {r['steps']:>10,} | {r['drift_pct']:>10.4f} | {status:>10}")
return results
def run_extended_n_sweep(base_config, all_results, verbose=True):
"""TEST 2: Extended N* Sweep (N = [256, 512])"""
print("\n" + "="*80)
print("🔬 TEST 2: Extended N* Sweep")
print("="*80)
print("Goal: Confirm N* at higher resolution")
print("Config: κ=0.50, N=[256,512], 50K steps")
print("="*80)
N_values = [256, 512]
results = []
for N in N_values:
config = base_config.copy()
config['N'] = N
config['dx'] = 0.4 * (64.0 / N)
config['kappa'] = 0.50
print(f"\n⚠️ MEMORY GATE: N={N} - Large grid detected")
print(f" Cells: {N**2:,}, Fields: ~{N**2 * 6 * 8 / 1024 / 1024:.1f} MB")
diag_interval = max(100, int(50000 / 50))
result = run_single_run(
config=config,
target_steps=50000,
diag_interval=diag_interval,
drift_threshold=0.05,
label=f"N_extended_{N}",
verbose=verbose
)
all_results.append(result)
results.append(result)
print("\n📊 Extended N* Results:")
print("-"*40)
for r in results:
status = "✅ OK" if not r['failed'] else f"❌ {r['failure_step']:,}"
print(f"N={r['N']:>3}: {r['steps']:>10,} steps, drift={r['drift_pct']:>10.4f}%, {status}")
return results
def run_high_kappa_validation(base_config, all_results, verbose=True):
"""TEST 3: High-κ Validation"""
print("\n" + "="*80)
print("🔬 TEST 3: High-κ Validation")
print("="*80)
print("Goal: Verify stability at higher κ")
print("Config: κ=[0.53,0.55], N=[64,128], 100K steps")
print("="*80)
kappa_values = [0.53, 0.55]
N_values = [64, 128]
results = []
for N in N_values:
for kappa in kappa_values:
config = base_config.copy()
config['N'] = N
config['dx'] = 0.4 * (64.0 / N)
config['kappa'] = kappa
diag_interval = max(100, int(100000 / 50))
result = run_single_run(
config=config,
target_steps=100000,
diag_interval=diag_interval,
drift_threshold=0.05,
label=f"highkappa_k{kappa:.2f}_N{N}",
verbose=verbose
)
all_results.append(result)
results.append(result)
print("\n📊 High-κ Validation Results:")
print("-"*40)
for r in results:
status = "✅ OK" if not r['failed'] else f"❌ {r['failure_step']:,}"
print(f"κ={r['kappa']:.2f}, N={r['N']:>3}: {r['steps']:>10,} steps, drift={r['drift_pct']:>10.4f}%, {status}")
return results
def run_constitutive_map_sweep(base_config, all_results, verbose=True):
"""TEST 4: Constitutive Map Stability"""
print("\n" + "="*80)
print("🔬 TEST 4: Constitutive Map Stability")
print("="*80)
print("Goal: Map exact transition to stable regime")
print("Config: κ=[0.47,0.48,0.49,0.50,0.51,0.52], N=64, 100K steps")
print("="*80)
kappa_values = [0.47, 0.48, 0.49, 0.50, 0.51, 0.52]
N = 64
results = []
for kappa in kappa_values:
config = base_config.copy()
config['N'] = N
config['dx'] = 0.4 * (64.0 / N)
config['kappa'] = kappa
diag_interval = max(100, int(100000 / 50))
result = run_single_run(
config=config,
target_steps=100000,
diag_interval=diag_interval,
drift_threshold=0.05,
label=f"constitutive_kappa{kappa:.2f}_N{N}",
verbose=verbose
)
all_results.append(result)
results.append(result)
print("\n📊 Constitutive Map Results:")
print("-"*40)
for r in results:
status = "✅ OK" if not r['failed'] else f"❌ {r['failure_step']:,}"
print(f"κ={r['kappa']:.2f}: {r['steps']:>10,} steps, drift={r['drift_pct']:>10.4f}%, {status}")
return results
# ==============================================================================
# MAIN
# ==============================================================================
def main():
print("\n" + "="*80)
print("🚀 MONAD SERIES 11 — N* IDENTIFICATION & DIAGNOSTICS (FIXED)")
print("="*80)
print("FIXED: Initialization-order bug resolved")
print("FIXED: Diagnostics now called AFTER init_gaussian()")
print("="*80)
print("Tests:")
print(" 1. N* Sweep (N=[16,32,64,128,256], 100K steps)")
print(" 2. Extended N* Sweep (N=[256,512], 50K steps)")
print(" 3. High-κ Validation (κ=[0.53,0.55], N=[64,128])")
print(" 4. Constitutive Map Stability (κ=0.47-0.52, N=64)")
print("="*80)
# Base configuration
base_config = {
'N': 16,
'dx': 0.4 * (64.0 / 16),
'dt': 5e-6,
'c': 0.5,
'kappa': 0.3,
'eta': 0.2,
'beta': 0.5,
'gamma': 0.2,
'm2': 0.1,
'alpha': 0.4,
'delta': 0.15,
'Pi_max': 5.9259,
'ko_sigma': 0.045,
'anchor': 0.0,
'no_cfl_check': True
}
all_results = []
timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')
print(f"\n📂 Results timestamp: {timestamp}")
print(f"📁 Output: /content/series11_nstar/")
print("="*80)
# Run all tests
print("\n🔬 RUNNING TEST 1: N* Sweep")
run_n_star_sweep(base_config, all_results, verbose=True)
print("\n🔬 RUNNING TEST 2: Extended N* Sweep")
run_extended_n_sweep(base_config, all_results, verbose=True)
print("\n🔬 RUNNING TEST 3: High-κ Validation")
run_high_kappa_validation(base_config, all_results, verbose=True)
print("\n🔬 RUNNING TEST 4: Constitutive Map Stability")
run_constitutive_map_sweep(base_config, all_results, verbose=True)
# Save results
output_dir = Path('/content/series11_nstar')
output_dir.mkdir(exist_ok=True)
json_path = output_dir / f'series11_nstar_{timestamp}.json'
with open(json_path, 'w') as f:
json.dump(all_results, f, indent=2, default=str)
print(f"\n💾 Results saved to: {json_path}")
# Check if any valid runs
valid_runs = [r for r in all_results if not r.get('failed', False)]
if valid_runs:
# Generate charts
print("\n📊 Generating charts...")
generate_charts(all_results, output_dir / 'charts')
else:
print("\n⚠️ No valid runs completed - skipping charts")
# Generate Markdown report
print("\n📝 Generating Markdown report...")
md_report = generate_markdown_report(all_results, timestamp)
md_path = output_dir / f'Nstar_report_{timestamp}.md'
with open(md_path, 'w') as f:
f.write(md_report)
print(f"✅ Markdown report saved to: {md_path}")
# Summary
total = len(all_results)
failed = len([r for r in all_results if r.get('failed', False)])
success = total - failed
print("\n" + "="*80)
print("📊 SERIES 11 N* IDENTIFICATION — SUMMARY")
print("="*80)
print(f" Total runs: {total}")
print(f" Successful: {success} ({success/total*100:.1f}%)")
print(f" Failed: {failed} ({failed/total*100:.1f}%)")
if valid_runs:
best = min(valid_runs, key=lambda x: x.get('drift_rel', 1.0))
print(f"\n 🏆 Best: κ={best['kappa']:.3f}, N={best['N']}, drift={best['drift_pct']:.4f}%")
else:
print("\n ⚠️ No successful runs completed.")
print(" The initialization-order bug has been fixed.")
print(" Re-run the pipeline to get valid data.")
print("\n" + "="*80)
print("✅ N* IDENTIFICATION PIPELINE COMPLETE")
print("="*80)
print(f"\n📁 All results saved to: /content/series11_nstar/")
print(f" - JSON: {json_path.name}")
if valid_runs:
print(f" - Charts: charts/")
print(f" - Markdown: {md_path.name}")
if __name__ == "__main__":
main()
#!/usr/bin/env python3
"""
================================================================================
MONAD SERIES 11 — PURE PHYSICS SOLVER (NO DIAGNOSTICS, NO DATA PIPELINE)
================================================================================
PURE CPU-ONLY SOLVER — No file I/O, no diagnostics, no data pipeline.
Just the physics: one Monad tensor Π evolves according to the expanded
FRCMFD equations with exact invariant decomposition.
Run: python monad_series11_pure.py --steps 1000 --N 16 --kappa 0.3
================================================================================
"""
import sys
import argparse
import numpy as np
import time
# ==============================================================================
# FILTER COLAB -f ARGUMENT
# ==============================================================================
def parse_colab_args():
"""Parse arguments, filtering Colab's -f flag."""
raw_args = sys.argv[1:]
filtered = []
skip = False
for i, arg in enumerate(raw_args):
if skip:
skip = False
continue
if arg == '-f' and i + 1 < len(raw_args):
skip = True
continue
filtered.append(arg)
parser = argparse.ArgumentParser(description='Monad Series 11 Pure Physics Solver')
parser.add_argument('--steps', type=int, default=25000, help='Steps per run')
parser.add_argument('--N', type=int, default=16, help='Grid size (16 or 32 for CPU)')
parser.add_argument('--dt', type=float, default=5e-6, help='Timestep')
parser.add_argument('--c', type=float, default=0.5, help='Wave speed')
parser.add_argument('--kappa', type=float, default=0.3, help='Coupling parameter')
parser.add_argument('--ko_sigma', type=float, default=0.045, help='KO dissipation')
parser.add_argument('--amplitude', type=float, default=2.5, help='Initial amplitude')
parser.add_argument('--diag', type=int, default=1000, help='Diagnostic interval')
parser.add_argument('--no_cfl_check', action='store_true', default=False,
help='Skip CFL warning')
return parser.parse_args(filtered)
# ==============================================================================
# CPU DIFF-OPERATOR KERNELS (VECTORIZED NumPy)
# ==============================================================================
def cpu_grad_2d(F, dx):
"""Compute 2D gradient using central differences (vectorized)."""
dx = float(dx)
gx = (np.roll(F, -1, axis=0) - np.roll(F, 1, axis=0)) / (2.0 * dx)
gy = (np.roll(F, -1, axis=1) - np.roll(F, 1, axis=1)) / (2.0 * dx)
return gx, gy
def cpu_lap_2d(F, dx):
"""Compute 2D Laplacian (vectorized)."""
dx2 = float(dx) * float(dx)
lap = (np.roll(F, -1, axis=0) + np.roll(F, 1, axis=0) +
np.roll(F, -1, axis=1) + np.roll(F, 1, axis=1) - 4.0 * F) / dx2
return lap
def cpu_ko_2d(F, dx, sigma):
"""4th-order Kreiss-Oliger dissipation (vectorized)."""
sigma = float(sigma)
dx4 = float(dx)**4
d4x = (np.roll(F, -2, axis=0) - 4.0*np.roll(F, -1, axis=0) +
6.0*F - 4.0*np.roll(F, 1, axis=0) + np.roll(F, 2, axis=0))
d4y = (np.roll(F, -2, axis=1) - 4.0*np.roll(F, -1, axis=1) +
6.0*F - 4.0*np.roll(F, 1, axis=1) + np.roll(F, 2, axis=1))
return -(sigma / 16.0) * (d4x + d4y) / dx4
# ==============================================================================
# INVARIANT DECOMPOSITION OF Π (EXACTLY INVERTIBLE)
# ==============================================================================
def decompose_pi(Pxx, Pxy, Pyy):
"""Decompose symmetric tensor Π into invariant stress modes."""
Lam = Pxx + Pyy
S = Pxx - Pyy
Psi = Pxy
return S, Psi, Lam
def reconstruct_pi(S, Psi, Lam):
"""Reconstruct Π from invariant stress modes (exact inverse)."""
Pxx = (Lam + S) / 2.0
Pyy = (Lam - S) / 2.0
Pxy = Psi
return Pxx, Pxy, Pyy
# ==============================================================================
# CONSTITUTIVE EVALUATION FROM INVARIANTS
# ==============================================================================
def evaluate_monad_from_pi(Pxx, Pxy, Pyy, Pi_max=5.9259, anchor=0.0,
kappa=0.3, eta=0.2, beta=0.5, gamma=0.2,
m2=0.1, alpha=0.4, delta=0.15, eps=1e-15):
"""
Evaluate ALL FRCMFD quantities from the Monad tensor Π.
VERIFIED INVARIANTS:
- I1 = |S| + |Λ|
- I2 = S² - Ψ² + Λ²
- I3 = |S|³ + |Λ|³
- I4 = S⁴ - Ψ⁴ + Λ⁴
- I3 normalized by Pi_max**3
"""
S, Psi, Lam = decompose_pi(Pxx, Pxy, Pyy)
psi_mag_sq = Psi**2
I1 = np.abs(S) + np.abs(Lam)
I2 = S**2 - psi_mag_sq + Lam**2
I3 = np.abs(S**3) + np.abs(Lam**3)
I4 = S**4 - psi_mag_sq**2 + Lam**4
Ih1 = np.maximum(I1 / Pi_max, eps)
Ih2 = I2 / (Pi_max**2)
Ih3 = I3 / (Pi_max**3)
Ih4 = I4 / (Pi_max**4)
sqrt_Ih1 = np.sqrt(Ih1)
expf = np.exp(-0.5 * (Ih2**2 + Ih3**3 + Ih4**4))
lf = (1.0 / sqrt_Ih1) - 1.0
Psi_val = (1.0 / Pi_max) * lf * expf + anchor
I1_safe = np.maximum(I1, eps)
dPsi_dI1 = -1.0 / (2.0 * Pi_max * I1_safe * sqrt_Ih1) * expf
base = (1.0 / Pi_max) * lf * expf
dPsi_dI2 = -Ih2 / (Pi_max**2) * base
dPsi_dI3 = -1.5 * (Ih3**2) / (Pi_max**3) * base
dPsi_dI4 = -2.0 * (Ih4**3) / (Pi_max**4) * base
eps2 = 1e-10
s_smooth = S / (S**2 + eps2)
lam_smooth = Lam / (Lam**2 + eps2)
M_T = dPsi_dI1 * s_smooth + dPsi_dI2 * 2.0 * S + dPsi_dI3 * 3.0 * S * np.abs(S) + dPsi_dI4 * 4.0 * S**3
M_T = np.clip(M_T, -1e6, 1e6)
M_C = dPsi_dI1 * lam_smooth + dPsi_dI2 * 2.0 * Lam + dPsi_dI3 * 3.0 * Lam * np.abs(Lam) + dPsi_dI4 * 4.0 * Lam**3
M_C = np.clip(M_C, -1e6, 1e6)
M_R = np.clip(dPsi_dI2 * 2.0, -1e6, 1e6)
return {
'S': S, 'Psi': Psi, 'Lam': Lam,
'Psi_val': Psi_val,
'dPsi_dI1': dPsi_dI1, 'dPsi_dI2': dPsi_dI2,
'dPsi_dI3': dPsi_dI3, 'dPsi_dI4': dPsi_dI4,
'M_T': M_T, 'M_C': M_C, 'M_R': M_R,
'D_drift': np.abs(np.clip(np.log(np.maximum(np.abs(Psi_val - anchor) * Pi_max, eps)), -20, 20) +
0.5 * (np.clip(np.log(np.maximum(Ih1, eps)), -20, 20) + Ih2**2 + Ih3**3 + Ih4**4))
}
# ==============================================================================
# PURE PHYSICS SOLVER — NO DIAGNOSTICS, NO DATA PIPELINE
# ==============================================================================
class MonadSolver11Pure:
"""
Series 11: Pure physics solver — no diagnostics, no I/O, no data pipeline.
ONE Monad tensor Π evolves. S, Ψ, Λ are derived diagnostically.
"""
def __init__(self, config):
self.config = config
self.N = int(config['N'])
self.dx = float(config['dx'])
self.dt = float(config['dt'])
self.c2 = float(config.get('c', 0.5)**2)
self.kappa = float(config.get('kappa', 0.3))
self.eta = float(config.get('eta', 0.2))
self.beta = float(config.get('beta', 0.5))
self.gamma = float(config.get('gamma', 0.2))
self.m2 = float(config.get('m2', 0.1))
self.alpha = float(config.get('alpha', 0.4))
self.delta = float(config.get('delta', 0.15))
self.Pi_max = float(config.get('Pi_max', 5.9259))
self.ko_sigma = float(config.get('ko_sigma', 0.045))
self.anchor = float(config.get('anchor', 0.0))
self.step = 0
self.time = 0.0
# Fields
self.Pxx = None
self.Pxy = None
self.Pyy = None
self.Uxx = None
self.Uxy = None
self.Uyy = None
def init_gaussian(self, amplitude=2.5, width=4.0):
"""Initialize the Monad tensor Π with Gaussian profiles."""
N, dx = self.N, self.dx
half_width = (N * dx) / 2.0
x = np.linspace(-half_width, half_width, N, endpoint=False)
y = np.linspace(-half_width, half_width, N, endpoint=False)
X, Y = np.meshgrid(x, y)
R2 = X**2 + Y**2
amplitude = float(amplitude)
width = float(width)
self.Pxx = amplitude * np.exp(-R2 / (2.0 * width**2))
self.Pxy = 0.5 * amplitude * np.exp(-R2 / (2.0 * (width * 1.2)**2))
self.Pyy = 0.7 * amplitude * np.exp(-R2 / (2.0 * (width * 0.8)**2))
self.Uxx = np.zeros_like(self.Pxx)
self.Uxy = np.zeros_like(self.Pxy)
self.Uyy = np.zeros_like(self.Pyy)
# Small perturbation
np.random.seed(42)
noise_scale = 1e-8 * amplitude
self.Pxx += np.random.normal(0, noise_scale, self.Pxx.shape)
self.Pxy += np.random.normal(0, noise_scale, self.Pxy.shape)
self.Pyy += np.random.normal(0, noise_scale, self.Pyy.shape)
return self
def _finite_check(self):
"""Check for NaNs/Infs."""
if not np.all(np.isfinite(self.Pxx)):
raise RuntimeError("NaN/Inf detected")
if not np.all(np.isfinite(self.Pxy)):
raise RuntimeError("NaN/Inf detected")
if not np.all(np.isfinite(self.Pyy)):
raise RuntimeError("NaN/Inf detected")
def derivatives(self):
"""Compute time derivatives of the ONE Monad tensor Π."""
dx = self.dx
# Get diagnostics from current state
results = evaluate_monad_from_pi(
self.Pxx, self.Pxy, self.Pyy,
self.Pi_max, self.anchor,
self.kappa, self.eta,
self.beta, self.gamma,
self.m2, self.alpha, self.delta
)
S = results['S']
Psi = results['Psi']
Lam = results['Lam']
M_T = results['M_T']
M_C = results['M_C']
M_R = results['M_R']
# Laplacians
lapPxx = cpu_lap_2d(self.Pxx, dx)
lapPxy = cpu_lap_2d(self.Pxy, dx)
lapPyy = cpu_lap_2d(self.Pyy, dx)
# Gradients of stress modes
gSx, gSy = cpu_grad_2d(S, dx)
gLx, gLy = cpu_grad_2d(Lam, dx)
gPx, gPy = cpu_grad_2d(Psi, dx)
gS2 = gSx**2 + gSy**2
gL2 = gLx**2 + gLy**2
gP2 = gPx**2 + gPy**2
ps = Psi**2
ls = Lam**2
# Kinematic
dPxx = self.Uxx
dPxy = self.Uxy
dPyy = self.Uyy
# Momentum (force from action)
dUxx = (self.c2 * lapPxx - self.beta * self.Pxx - self.gamma * self.Pxx * self.Pxx**2 -
self.kappa * ps - self.eta * self.Pxx * ls + self.kappa * self.Pxx * M_T * gS2)
dUxy = (self.c2 * lapPxy - self.m2 * self.Pxy -
2.0 * self.kappa * self.Pxx * self.Pxy -
self.eta * self.Pxy * ls - self.kappa * self.Pxy * M_R * gP2)
dUyy = (self.c2 * lapPyy - self.alpha * self.Pyy - self.delta * self.Pyy * self.Pyy**2 -
self.kappa * self.Pxx * self.Pyy - self.eta * ps * self.Pyy +
self.kappa * self.Pyy * M_C * gL2)
# KO dissipation
if self.ko_sigma > 0:
dUxx += cpu_ko_2d(self.Uxx, dx, self.ko_sigma)
dUxy += cpu_ko_2d(self.Uxy, dx, self.ko_sigma)
dUyy += cpu_ko_2d(self.Uyy, dx, self.ko_sigma)
return dPxx, dPxy, dPyy, dUxx, dUxy, dUyy
def rk3_step(self):
"""SSP-RK3 timestep."""
dt = self.dt
Pxx_n = self.Pxx.copy()
Pxy_n = self.Pxy.copy()
Pyy_n = self.Pyy.copy()
Uxx_n = self.Uxx.copy()
Uxy_n = self.Uxy.copy()
Uyy_n = self.Uyy.copy()
# Stage 1
dPxx1, dPxy1, dPyy1, dUxx1, dUxy1, dUyy1 = self.derivatives()
self.Pxx = Pxx_n + dt*dPxx1
self.Pxy = Pxy_n + dt*dPxy1
self.Pyy = Pyy_n + dt*dPyy1
self.Uxx = Uxx_n + dt*dUxx1
self.Uxy = Uxy_n + dt*dUxy1
self.Uyy = Uyy_n + dt*dUyy1
self._finite_check()
# Stage 2
dPxx2, dPxy2, dPyy2, dUxx2, dUxy2, dUyy2 = self.derivatives()
self.Pxx = 0.75*Pxx_n + 0.25*self.Pxx + 0.25*dt*dPxx2
self.Pxy = 0.75*Pxy_n + 0.25*self.Pxy + 0.25*dt*dPxy2
self.Pyy = 0.75*Pyy_n + 0.25*self.Pyy + 0.25*dt*dPyy2
self.Uxx = 0.75*Uxx_n + 0.25*self.Uxx + 0.25*dt*dUxx2
self.Uxy = 0.75*Uxy_n + 0.25*self.Uxy + 0.25*dt*dUxy2
self.Uyy = 0.75*Uyy_n + 0.25*self.Uyy + 0.25*dt*dUyy2
self._finite_check()
# Stage 3
dPxx3, dPxy3, dPyy3, dUxx3, dUxy3, dUyy3 = self.derivatives()
self.Pxx = (1.0/3.0)*Pxx_n + (2.0/3.0)*self.Pxx + (2.0/3.0)*dt*dPxx3
self.Pxy = (1.0/3.0)*Pxy_n + (2.0/3.0)*self.Pxy + (2.0/3.0)*dt*dPxy3
self.Pyy = (1.0/3.0)*Pyy_n + (2.0/3.0)*self.Pyy + (2.0/3.0)*dt*dPyy3
self.Uxx = (1.0/3.0)*Uxx_n + (2.0/3.0)*self.Uxx + (2.0/3.0)*dt*dUxx3
self.Uxy = (1.0/3.0)*Uxy_n + (2.0/3.0)*self.Uxy + (2.0/3.0)*dt*dUxy3
self.Uyy = (1.0/3.0)*Uyy_n + (2.0/3.0)*self.Uyy + (2.0/3.0)*dt*dUyy3
self._finite_check()
self.step += 1
self.time += dt
def hamiltonian(self):
"""Compute Hamiltonian from the Monad tensor Π."""
dx = self.dx
results = evaluate_monad_from_pi(
self.Pxx, self.Pxy, self.Pyy,
self.Pi_max, self.anchor,
self.kappa, self.eta,
self.beta, self.gamma,
self.m2, self.alpha, self.delta
)
S = results['S']
Psi = results['Psi']
Lam = results['Lam']
kin = 0.5 * (self.Uxx**2 + self.Uxy**2 + self.Uyy**2)
gPxx_x, gPxx_y = cpu_grad_2d(self.Pxx, dx)
gPxy_x, gPxy_y = cpu_grad_2d(self.Pxy, dx)
gPyy_x, gPyy_y = cpu_grad_2d(self.Pyy, dx)
grad = 0.5 * self.c2 * (gPxx_x**2 + gPxx_y**2 +
gPxy_x**2 + gPxy_y**2 +
gPyy_x**2 + gPyy_y**2)
ps = Psi**2
pot = (0.5 * self.beta * S**2 + 0.25 * self.gamma * S**4 +
0.5 * self.m2 * ps + 0.5 * self.alpha * Lam**2 + 0.25 * self.delta * Lam**4 +
self.kappa * S * ps + self.eta * ps * Lam)
return float(np.sum(kin + grad + pot) * dx**2)
def run(self, steps, diag_interval=1000):
"""Run simulation for specified steps."""
print(f"\n{'='*60}")
print(f"SERIES 11 — PURE PHYSICS SOLVER (NO DIAGNOSTICS)")
print(f"{'='*60}")
print(f"Grid: {self.N}×{self.N}, dt={self.dt:.6f}, κ={self.kappa:.3f}")
print(f"Diagnostics: every {diag_interval} steps")
print(f"{'-'*60}")
start = time.time()
H0 = self.hamiltonian()
print(f"Step {0:6d} | H = {H0:.8e}")
for step in range(1, steps+1):
self.rk3_step()
if step % diag_interval == 0:
H = self.hamiltonian()
dr = abs(H-H0)/max(abs(H0), 1e-30)
print(f"Step {step:6d} | H={H:.8e} | ΔH/H={dr:.4e}")
if np.isnan(H) or dr > 0.05:
print(f"❌ INSTABILITY at step {step}")
break
elapsed = time.time() - start
Hf = self.hamiltonian()
dr = abs(Hf - H0) / max(abs(H0), 1e-30)
print(f"{'-'*60}")
print(f"Completed {step} steps in {elapsed:.2f}s ({step/elapsed:.1f} steps/sec)")
print(f"{'='*60}")
return {
'steps': step,
'H0': float(H0),
'H_final': float(Hf),
'drift_rel': float(dr),
'elapsed': float(elapsed)
}
# ==============================================================================
# MAIN
# ==============================================================================
if __name__ == "__main__":
args = parse_colab_args()
print("\n" + "="*80)
print("🚀 SERIES 11 — PURE PHYSICS SOLVER")
print("="*80)
print(f"Steps: {args.steps:,}")
print(f"Grid: {args.N}×{args.N}")
print(f"dt: {args.dt}")
print(f"κ: {args.kappa}")
print("="*80)
print("✅ PURE PHYSICS — No diagnostics, no I/O, no data pipeline")
print("✅ ONE Monad tensor Π evolves (3 components: Pxx, Pxy, Pyy)")
print("✅ THREE stress modes derived INVARIANTLY from Π:")
print(" - Compression (Λ) = trace(Π) = Pxx + Pyy")
print(" - Tension (S) = deviatoric = Pxx - Pyy")
print(" - Torsion (Ψ) = off-diagonal = Pxy")
print("="*80)
config = {
'N': args.N,
'dx': 0.4 * (64.0 / args.N),
'dt': args.dt,
'c': args.c,
'kappa': args.kappa,
'eta': 0.2,
'beta': 0.5,
'gamma': 0.2,
'm2': 0.1,
'alpha': 0.4,
'delta': 0.15,
'Pi_max': 5.9259,
'ko_sigma': args.ko_sigma,
'anchor': 0.0,
'no_cfl_check': args.no_cfl_check
}
# CFL Check
c = float(config.get('c', 0.5))
dx = float(config['dx'])
dt = float(config['dt'])
cfl = c * dt / dx
if cfl > 0.5 and not config.get('no_cfl_check', False):
print(f" ⚠️ CFL ~ {cfl:.3f} > 0.5 — dt may be too large for stability")
print(f" CFL = {cfl:.4f}")
solver = MonadSolver11Pure(config)
solver.init_gaussian(amplitude=args.amplitude)
results = solver.run(steps=args.steps, diag_interval=args.diag)
print("\n" + "="*80)
print("📊 SERIES 11 PURE PHYSICS RESULTS")
print("="*80)
print(f"Steps completed: {results['steps']:,}/{args.steps:,}")
print(f"Initial H: {results['H0']:.8e}")
print(f"Final H: {results['H_final']:.8e}")
print(f"Drift: {results['drift_rel']:.4e} ({results['drift_rel']*100:.4f}%)")
print(f"Time: {results['elapsed']:.2f}s")
print(f"Speed: {results['steps']/results['elapsed']:.1f} steps/sec")
print("="*80)
