AUTOMATED ZERO-DRIFT TRACKER — Secant Method Optimization for κ*

INSTALL MONAD SERIES 11.1 — PURE PHYSICS SOLVER FIRST - THEN RUN AUTOMATED ZERO-DRIFT TRACKER #!/usr/bin/env python3 """ ================================================================================ MONAD SERIES 11.1 — PURE PHYSICS SOLVER (CLARITY + SAFE REGULARIZATION) ================================================================================ This is a light, non-physical-change refactor of the original solver. Purpose: - Make the primitive / constitutive distinction explicit: - Psi_field: the primitive evolved field (Pxy) - Psi_const: the constitutive reconstruction Ψ(Ik(Π)) - Use safer sqrt-based regularizers for sign-like smoothing. - Improve finite-check messages for easier debugging. - Preserve all original PDEs, constitutive functional forms and dynamics. - Keep Π as the sole prognostic state (Version A behavior). - Add clarifying docstrings/comments where helpful. Run: python monad_series11.1_purephysics.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.1 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 # ============================================================================== # NOTE: This function computes diagnostic/constitutive quantities from Π. # It does NOT alter Π itself. The solver treats Π as primitive. # Return keys: # - Psi_field : the primitive Pxy passed in # - Psi_const : the constitutive reconstruction Ψ(Ik(Π)) # - dPsi_const_dI* : derivatives of Ψ_const wrt invariants (for modulatory terms) # ============================================================================== 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 Π. This returns both the primitive field (Psi_field == Pxy) and the constitutive reconstruction Psi_const = Ψ(Ik). The solver uses Psi_field for evolution and Psi_const for diagnostics / emergent geometry. """ S, Psi_field, Lam = decompose_pi(Pxx, Pxy, Pyy) psi_mag_sq = Psi_field**2 # Invariants 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) # Constitutive functional form (unchanged): expf = np.exp(-0.5 * (Ih2**2 + Ih3**3 + Ih4**4)) lf = (1.0 / sqrt_Ih1) - 1.0 Psi_const = (1.0 / Pi_max) * lf * expf + anchor # Derivatives of Psi_const wrt invariants (chain-rule components) 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 # Safer smooth sign approximations (use sqrt-based regularizer) eps2 = 1e-10 s_smooth = S / np.sqrt(S**2 + eps2) lam_smooth = Lam / np.sqrt(Lam**2 + eps2) # Compose modulatory multipliers using derivatives of Psi_const 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) # Constitutive diagnostic / drift measure (unchanged form) D_drift = np.abs( np.clip(np.log(np.maximum(np.abs(Psi_const - anchor) * Pi_max, eps)), -20, 20) + 0.5 * (np.clip(np.log(np.maximum(Ih1, eps)), -20, 20) + Ih2**2 + Ih3**3 + Ih4**4) ) return { 'S': S, 'Psi_field': Psi_field, 'Psi_const': Psi_const, 'Lam': Lam, 'dPsi_const_dI1': dPsi_dI1, 'dPsi_const_dI2': dPsi_dI2, 'dPsi_const_dI3': dPsi_dI3, 'dPsi_const_dI4': dPsi_dI4, 'M_T': M_T, 'M_C': M_C, 'M_R': M_R, 'D_drift': D_drift } # ============================================================================== # PURE PHYSICS SOLVER — NO DIAGNOSTICS, NO DATA PIPELINE # ============================================================================== class MonadSolver11Pure: """ Series 11.1: Pure physics solver — clarity fixes and safe regularization. KEY PRINCIPLE (Version A): - Π (Pxx, Pxy, Pyy) is the primitive dynamical variable and evolves. - Ψ_const = Ψ(Ik(Π)) is reconstructed diagnostically (emergent geometry). - By default, Ψ_const does NOT feed back into the time evolution. (Constitutive feedback is a well-defined separate extension.) """ 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 and report which field fails for easier debugging.""" if self.Pxx is None or self.Pxy is None or self.Pyy is None: raise RuntimeError("Fields are not initialized (Pxx/Pxy/Pyy are None)") if not np.all(np.isfinite(self.Pxx)): raise RuntimeError("NaN/Inf detected in Pxx") if not np.all(np.isfinite(self.Pxy)): raise RuntimeError("NaN/Inf detected in Pxy") if not np.all(np.isfinite(self.Pyy)): raise RuntimeError("NaN/Inf detected in Pyy") if not np.all(np.isfinite(self.Uxx)): raise RuntimeError("NaN/Inf detected in Uxx") if not np.all(np.isfinite(self.Uxy)): raise RuntimeError("NaN/Inf detected in Uxy") if not np.all(np.isfinite(self.Uyy)): raise RuntimeError("NaN/Inf detected in Uyy") def derivatives(self): """Compute time derivatives of the ONE Monad tensor Π. IMPORTANT: This uses the primitive field (Psi_field == Pxy) for dynamics. Constituive quantities (Psi_const) are used only where modulatory multipliers computed from derivatives of Psi_const are intentionally included. """ dx = self.dx # Diagnostics from current state (returns both primitive and constitutive quantities) 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_field = results['Psi_field'] # primitive Pxy 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 (use primitive field gradients for dynamics) gSx, gSy = cpu_grad_2d(S, dx) gLx, gLy = cpu_grad_2d(Lam, dx) gPx, gPy = cpu_grad_2d(Psi_field, dx) # primitive gradients gS2 = gSx**2 + gSy**2 gL2 = gLx**2 + gLy**2 gP2 = gPx**2 + gPy**2 ps_field = Psi_field**2 ls = Lam**2 # Kinematic dPxx = self.Uxx dPxy = self.Uxy dPyy = self.Uyy # Momentum (force from action) - uses primitive psi_field where originally designed dUxx = (self.c2 * lapPxx - self.beta * self.Pxx - self.gamma * self.Pxx * self.Pxx**2 - self.kappa * ps_field - 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_field * 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 Π. NOTE: For emergent metric / potential diagnostics we use Psi_const (the constitutive reconstruction). This makes the Hamiltonian reflect the emergent geometry produced by Ψ_const(Ik). """ 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_const = results['Psi_const'] # emergent constitutive reconstruction 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_const**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.1 — PURE PHYSICS SOLVER (CLARITY + SAFE REGULARIZATION)") 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 try: Hf = self.hamiltonian() except Exception: # If hamiltonian computation fails, fallback to last known H (H0) to avoid crash Hf = H0 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.1 — PURE PHYSICS SOLVER (CLARITY RELEASE)") 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.1 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) #!/usr/bin/env python3 """ AUTOMATED ZERO-DRIFT TRACKER — Secant Method Optimization for κ* Finds and verifies the perfect energy-conserving coupling coefficient using: - Run A: κ = 0.58 (lower bound) - Run B: κ = 0.61 (upper bound) - Secant root-finding for rapid convergence - Run C: κ = κ*_estimate (verification) - Polynomial refinement with all 3 points - Run D: κ = κ*_refined (200k decisive test) - Fully automated — no manual intervention required - Saves complete summary to JSON + TAR BASELINE ANCHOR: monad_series11.1_purephysics.py (UNCHANGED) """ import sys import json import time import numpy as np from pathlib import Path from datetime import datetime import tarfile # ============================================================================== # BASELINE ANCHOR - EXPLICIT PATH IMPORT (BYPASSES ALL MODULE CACHE ISSUES) # ============================================================================== import os import importlib.util file_path = "/content/monad_series11.1_purephysics.py" if not os.path.exists(file_path): # Try alternative names alt = "/content/monad_series11_pure.py" if os.path.exists(alt): file_path = alt else: raise FileNotFoundError("❌ Could not find monad_series11.1_purephysics.py") print(f"✅ Loading solver from: {file_path}") spec = importlib.util.spec_from_file_location("solver_module", file_path) solver_module = importlib.util.module_from_spec(spec) spec.loader.exec_module(solver_module) MonadSolver11Pure = solver_module.MonadSolver11Pure evaluate_monad_from_pi = solver_module.evaluate_monad_from_pi cpu_grad_2d = solver_module.cpu_grad_2d print("✅ Solver imported successfully") # ============================================================================== # CONFIGURATION # ============================================================================== BASELINE_CONFIG = { 'N': 64, 'dt': 1.25e-6, 'ko_sigma': 0.045, 'amplitude': 2.5, 'seed': 12345, 'diag_interval': 5000, } OUTDIR = Path("./zero_drift_bracketing") OUTDIR.mkdir(parents=True, exist_ok=True) # Known bracket from previous runs KAPPA_LOW = 0.55 KAPPA_HIGH = 0.65 # ============================================================================== # ENERGY COMPONENT DECOMPOSITION # ============================================================================== def compute_energy_components(solver): """Decompose Hamiltonian into components.""" dx = solver.dx results = evaluate_monad_from_pi( solver.Pxx, solver.Pxy, solver.Pyy, solver.Pi_max, solver.anchor, solver.kappa, solver.eta, solver.beta, solver.gamma, solver.m2, solver.alpha, solver.delta ) S = results.get('S', solver.Pxx - solver.Pyy) Psi_const = results.get('Psi_const', results.get('Psi_field', solver.Pxy)) Lam = results.get('Lam', solver.Pxx + solver.Pyy) # KINETIC H_kin = 0.5 * (solver.Uxx**2 + solver.Uxy**2 + solver.Uyy**2) H_kin_s = float(np.sum(H_kin) * dx**2) # GRADIENT gPxx_x, gPxx_y = cpu_grad_2d(solver.Pxx, dx) gPxy_x, gPxy_y = cpu_grad_2d(solver.Pxy, dx) gPyy_x, gPyy_y = cpu_grad_2d(solver.Pyy, dx) H_grad = 0.5 * solver.c2 * (gPxx_x**2 + gPxx_y**2 + gPxy_x**2 + gPxy_y**2 + gPyy_x**2 + gPyy_y**2) H_grad_s = float(np.sum(H_grad) * dx**2) # POTENTIAL ps = Psi_const**2 H_pot = (0.5 * solver.beta * S**2 + 0.25 * solver.gamma * S**4 + 0.5 * solver.m2 * ps + 0.5 * solver.alpha * Lam**2 + 0.25 * solver.delta * Lam**4 + solver.kappa * S * ps + solver.eta * ps * Lam) H_pot_s = float(np.sum(H_pot) * dx**2) H_total_s = H_kin_s + H_grad_s + H_pot_s return { 'H_total': H_total_s, 'H_kinetic': H_kin_s, 'H_gradient': H_grad_s, 'H_potential': H_pot_s, } # ============================================================================== # SINGLE RUN ENGINE # ============================================================================== def run_single(kappa, steps, label, outdir=OUTDIR): """Run a single simulation at given kappa for given steps.""" config = BASELINE_CONFIG.copy() diag_interval = min(config['diag_interval'], steps // 10) solver_config = { 'N': config['N'], 'dx': 0.4 * (64.0 / config['N']), 'dt': config['dt'], 'c': 0.5, 'kappa': 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': config['ko_sigma'], 'anchor': 0.0, 'no_cfl_check': True } print(f"\n{'='*60}") print(f"⚡ RUN: {label}") print(f" κ = {kappa:.5f} | Steps = {steps} | dt = {config['dt']:.2e}") print(f"{'='*60}") solver = MonadSolver11Pure(solver_config) solver.init_gaussian(amplitude=config['amplitude']) history = [] H0 = compute_energy_components(solver) history.append({'step': 0, **H0}) print(f"Step 0: H_total={H0['H_total']:.8e}") for step in range(1, steps + 1): solver.rk3_step() if step % diag_interval == 0 or step == steps: H = compute_energy_components(solver) history.append({'step': step, **H}) if step % (diag_interval * 5) == 0 or step == steps: dr = (H['H_total'] - H0['H_total']) / max(abs(H0['H_total']), 1e-30) print(f" Step {step:6d}: H={H['H_total']:.8e} | Drift={dr*100:+.6f}%") dH = history[-1]['H_total'] - history[0]['H_total'] rel_drift = dH / max(abs(history[0]['H_total']), 1e-30) # Save data timestamp = datetime.utcnow().strftime("%Y%m%d_%H%M%S") data_file = outdir / f"{label}_{timestamp}.json" with open(data_file, 'w') as f: json.dump(history, f, indent=2, default=float) print(f" Completed {steps} steps") print(f" ΔH = {dH:+.6e} ({rel_drift*100:+.6f}%)") print(f" Data saved to: {data_file}") return dH, history, rel_drift # ============================================================================== # MAIN AUTOMATED BRACKETING ENGINE # ============================================================================== def main(): print("="*80) print("🎯 ZERO-DRIFT AUTOMATED BRACKETING ENGINE") print("="*80) print(f"Baseline bracket: κ ∈ [{KAPPA_LOW}, {KAPPA_HIGH}]") print("Method: Secant root-finding + polynomial refinement") print("Sequence: Run A → Run B → Secant → Run C → Refine → Run D") print("="*80) results_db = [] # ====================================================================== # STEP 1: Run A — Lower bound (κ = 0.58) # ====================================================================== print("\n" + "="*80) print("STEP 1: Run A — Lower Bound (κ = 0.58)") print("="*80) k_A = 0.58 dH_A, hist_A, rel_A = run_single(k_A, 50000, "Run_A_LowerBound") results_db.append({ 'label': 'Run A (Lower Bound)', 'kappa': k_A, 'steps': 50000, 'dH': dH_A, 'rel_drift': rel_A, 'drift_pct': rel_A * 100 }) # ====================================================================== # STEP 2: Run B — Upper bound (κ = 0.61) # ====================================================================== print("\n" + "="*80) print("STEP 2: Run B — Upper Bound (κ = 0.61)") print("="*80) k_B = 0.61 dH_B, hist_B, rel_B = run_single(k_B, 50000, "Run_B_UpperBound") results_db.append({ 'label': 'Run B (Upper Bound)', 'kappa': k_B, 'steps': 50000, 'dH': dH_B, 'rel_drift': rel_B, 'drift_pct': rel_B * 100 }) # ====================================================================== # STEP 3: Secant method — Estimate κ* # ====================================================================== print("\n" + "="*80) print("STEP 3: Secant Root-Finding") print("="*80) # Check sign change if dH_A * dH_B >= 0: print(f"⚠️ WARNING: No sign change between κ={k_A} (ΔH={dH_A:+.6e}) and κ={k_B} (ΔH={dH_B:+.6e})") print(" Proceeding with secant estimate anyway...") # Secant method: k_C = k_A - dH_A * (k_B - k_A) / (dH_B - dH_A) slope = (dH_B - dH_A) / (k_B - k_A) k_secant = k_A - (dH_A / slope) # Clip to safe range k_secant = np.clip(k_secant, KAPPA_LOW, KAPPA_HIGH) print(f" Secant estimate: κ* ≈ {k_secant:.6f}") print(f" Based on: ΔH_A={dH_A:+.6e}, ΔH_B={dH_B:+.6e}") print(f" Slope: {slope:.6e}") # ====================================================================== # STEP 4: Run C — Secant estimate verification # ====================================================================== print("\n" + "="*80) print(f"STEP 4: Run C — Secant Estimate (κ = {k_secant:.6f})") print("="*80) k_C = k_secant dH_C, hist_C, rel_C = run_single(k_C, 50000, f"Run_C_SecantEstimate_{k_C:.6f}") results_db.append({ 'label': f'Run C (Secant Estimate)', 'kappa': k_C, 'steps': 50000, 'dH': dH_C, 'rel_drift': rel_C, 'drift_pct': rel_C * 100 }) # ====================================================================== # STEP 5: Polynomial refinement using all 3 points # ====================================================================== print("\n" + "="*80) print("STEP 5: Polynomial Refinement (3-point linear fit)") print("="*80) k_points = [k_A, k_B, k_C] dH_points = [dH_A, dH_B, dH_C] # Linear fit: dH(k) = a*k + b coeffs = np.polyfit(k_points, dH_points, 1) a, b = coeffs[0], coeffs[1] k_star = -b / a # Clip to safe range k_star = np.clip(k_star, KAPPA_LOW, KAPPA_HIGH) print(f" Linear fit: ΔH = {a:.6e} * κ + {b:.6e}") print(f" κ* (zero-intercept): {k_star:.6f}") # Calculate slope at κ* slope_star = a print(f" Slope at κ*: d(ΔH)/dκ = {slope_star:.6e}") # ====================================================================== # STEP 6: Run D — The Decisive 200k Test # ====================================================================== print("\n" + "="*80) print(f"STEP 6: Run D — Decisive 200k Verification (κ = {k_star:.6f})") print("="*80) dH_D, hist_D, rel_D = run_single(k_star, 200000, f"Run_D_Verification_κ{k_star:.6f}") results_db.append({ 'label': f'Run D (200k Verification)', 'kappa': k_star, 'steps': 200000, 'dH': dH_D, 'rel_drift': rel_D, 'drift_pct': rel_D * 100 }) # ====================================================================== # STEP 7: Summary Report # ====================================================================== print("\n" + "="*80) print("📋 FINAL SUMMARY LEDGER") print("="*80) print(f"{'Run':<25} | {'κ':>12} | {'Steps':>8} | {'ΔH':>14} | {'Drift (%)':>12}") print("-"*80) for r in results_db: print(f"{r['label']:<25} | {r['kappa']:12.6f} | {r['steps']:8d} | {r['dH']:14.6e} | {r['drift_pct']:12.5f}") print("-"*80) print(f"\n🎯 OPTIMAL κ* = {k_star:.6f}") print(f" Decisive 200k test drift: {rel_D*100:.5f}%") print(f" Slope at κ*: d(ΔH)/dκ = {slope_star:.6e}") # Interpret slope if abs(slope_star) > 1.0: print(" 📈 System is HIGHLY sensitive to κ (steep slope)") print(" → Fine-tuned balance point") elif abs(slope_star) > 0.1: print(" 📊 System is MODERATELY sensitive to κ") print(" → Moderate tuning required") else: print(" 📉 System is WEAKLY sensitive to κ (shallow slope)") print(" → Broad neutral basin (robust)") # Determine verdict if abs(rel_D) < 0.0001: verdict = "✅ EXCELLENT: Near-perfect energy conservation at κ*!" elif abs(rel_D) < 0.001: verdict = "✅ GOOD: Very small drift at κ*." elif abs(rel_D) < 0.005: verdict = "⚠️ MODERATE: Significant drift at κ*, but bounded." else: verdict = "⚠️ HIGH: Consider further refinement." print(f" Verdict: {verdict}") # ====================================================================== # STEP 8: Save Summary # ====================================================================== timestamp = datetime.utcnow().strftime("%Y%m%d_%H%M%S") summary_file = OUTDIR / f"bracketing_summary_{timestamp}.json" summary_data = { 'timestamp': timestamp, 'kappa_star': float(k_star), 'slope_at_kappa_star': float(slope_star), 'runs': results_db, 'verdict': verdict, 'method': 'Secant root-finding + linear refinement', 'baseline_bracket': [KAPPA_LOW, KAPPA_HIGH], 'config': BASELINE_CONFIG } with open(summary_file, 'w') as f: json.dump(summary_data, f, indent=2, default=float) print(f"\n📁 Summary saved to: {summary_file}") # ====================================================================== # STEP 9: Create TAR of results # ====================================================================== print("\n" + "="*80) print("📦 CREATING TAR ARCHIVE") print("="*80) tar_name = f"zero_drift_bracketing_{timestamp}.tar" tar_path = Path("/content") / tar_name with tarfile.open(tar_path, 'w') as tar: for f in OUTDIR.iterdir(): if f.is_file(): tar.add(f, arcname=f.name) print(f" Added: {f.name}") size_kb = tar_path.stat().st_size / 1024 print(f"\n📦 Created: {tar_name} ({size_kb:.1f} KB)") print("\n" + "="*80) print("✅ BRACKETING COMPLETE") print("="*80) print(f"📁 Results saved to: {OUTDIR.resolve()}") print(f"📦 TAR archive: {tar_path}") print("\nTO DOWNLOAD:") print(" 1. Click the folder icon on the left sidebar (Files)") print(" 2. Navigate to /content/") print(f" 3. Right-click on: {tar_name}") print(" 4. Select 'Download'") print("="*80) if __name__ == "__main__": main()