FRCMFD-v2: TEST 1 — VELOCITY SATURATION & DISPERSIVE PHASE WAKE
"""
FRCMFD-v2: TEST 1 — VELOCITY SATURATION & DISPERSIVE PHASE WAKE
================================================================
This script tests the core hypothesis: a moving toroidal soliton experiences
directional retension overload, creating an asymmetric phase wake.
INPUT: Stationary soliton from Test 0A (test_0A_soliton_20260523_192753.npz)
OUTPUT: Asymmetry ratio vs velocity, high-k power, energy conservation
METHOD:
1. Load validated stationary soliton
2. Apply velocity boost: Ψ_boosted = Ψ_soliton * exp(1j * k * z)
3. Initialize Ψ_dot = -v_soliton * ∂Ψ/∂z
4. Evolve with same real-time engine that preserved soliton perfectly
5. Measure forward/rear phase gradient asymmetry
6. Sweep velocities: 0.3v, 0.5v, 0.7v, 0.85v, 0.95v
EXPECTED SIGNATURE (if hypothesis correct):
✓ Asymmetry ratio > 1.2 at high velocity
✓ Asymmetry increases monotonically with velocity
✓ High-k power increases with velocity (ripples in wake)
✓ Energy conserved (< 1% drift)
"""
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from datetime import datetime
import json
import warnings
warnings.filterwarnings('ignore')
print("="*80)
print("FRCMFD-v2: TEST 1 — VELOCITY SATURATION & DISPERSIVE PHASE WAKE")
print("="*80)
# =============================================================================
# LOAD VALIDATED SOLITON FROM TEST 0A
# =============================================================================
print("\n[Loading stationary soliton from Test 0A...]")
# Explicit path to the most recent soliton
soliton_file = "/content/test_0A_soliton_20260523_192753.npz"
print(f"✓ Loading: {soliton_file}")
data = np.load(soliton_file)
Psi_soliton = data['Psi_soliton']
r_grid = data['r_grid']
z_grid = data['z_grid']
dr = float(data['dr'])
dz = float(data['dz'])
v = float(data['v'])
mu = float(data['mu'])
lam = float(data['lam'])
kappa = float(data['kappa'])
m = int(data['m'])
S_max = float(data['S_max'])
Psi_sat = float(data['Psi_sat'])
nr = len(r_grid)
nz = len(z_grid)
print(f"✓ Grid: nr={nr}, nz={nz} ({nr*nz:,} DOF)")
print(f"✓ Soliton amplitude: {np.max(np.abs(Psi_soliton)):.4f}")
print(f"✓ Parameters: v={v}, mu={mu}, lam={lam}, kappa={kappa}, m={m}, S_max={S_max}, Psi_sat={Psi_sat}")
# =============================================================================
# REBUILD OPERATORS (same as Test 0A)
# =============================================================================
print("\n[Rebuilding operators...]")
def build_radial_operator(r_grid, dr):
nr = len(r_grid)
r_face = np.zeros(nr + 1)
r_face[0] = r_grid[0] - dr/2
for i in range(1, nr + 1):
r_face[i] = r_grid[i-1] + dr/2
flux_right = r_face[1:] / dr
flux_left = r_face[:-1] / dr
main_diag = -(flux_left + flux_right)
lower_diag = flux_left[1:]
upper_diag = flux_right[:-1]
M = sp.diags([lower_diag, main_diag, upper_diag], [-1, 0, 1], format='csr')
w_r = r_grid * dr
W_r = sp.diags(w_r, format='csr')
W_r_inv = sp.diags(1.0 / w_r, format='csr')
return W_r_inv @ M, W_r
def build_axial_operator(nz, dz):
main_diag = np.ones(nz) * (-2.0 / dz**2)
upper = np.ones(nz - 1) / dz**2
lower = np.ones(nz - 1) / dz**2
L_z = sp.diags([lower, main_diag, upper], [-1, 0, 1], format='csr')
W_z = sp.diags(np.ones(nz) * dz, format='csr')
return L_z, W_z
L_r, W_r = build_radial_operator(r_grid, dr)
L_z, W_z = build_axial_operator(nz, dz)
I_r = sp.eye(nr, format='csr')
I_z = sp.eye(nz, format='csr')
L_2D = sp.kron(I_z, L_r, format='csr') + sp.kron(L_z, I_r, format='csr')
W_2D = sp.kron(W_z, W_r, format='csr')
dV = W_2D.diagonal() * 2 * np.pi
r_mesh_2d = np.tile(r_grid, nz)
print(f"✓ Operators rebuilt: L_2D {L_2D.shape}, nnz={L_2D.nnz:,}")
# =============================================================================
# ENERGY AND ACCELERATION FUNCTIONS (EXACT FROM TEST 0)
# =============================================================================
def compute_energy(Psi, L_2D, dV, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m):
psi_sq = np.abs(Psi)**2
kin_grad = -0.5 * v**2 * np.real(np.sum(np.conj(Psi) * (L_2D @ Psi) * dV))
pot_mass = -0.5 * mu * np.sum(psi_sq * dV)
pot_nonlinear = 0.25 * lam * np.sum(psi_sq * psi_sq * dV)
S = S_max * np.tanh(psi_sq / (Psi_sat**2))
pot_tension = 0.5 * kappa * np.sum(S * psi_sq * dV)
pot_centrifugal = 0.5 * v**2 * m**2 * np.sum(psi_sq / (r_mesh_2d**2 + 1e-12) * dV)
return (kin_grad + pot_mass + pot_nonlinear + pot_tension + pot_centrifugal).real
def acceleration_real(Psi, L_2D, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m):
psi_sq = np.abs(Psi)**2
S = S_max * np.tanh(psi_sq / (Psi_sat**2))
dS = (S_max / (Psi_sat**2)) * (1.0 / np.cosh(psi_sq / (Psi_sat**2))**2)
term_kin = -v**2 * (L_2D @ Psi)
term_mass = mu * Psi
term_nonlinear = lam * psi_sq * Psi
term_tension = kappa * (S + psi_sq * dS) * Psi
term_centrifugal = +v**2 * m**2 * Psi / (r_mesh_2d**2 + 1e-12)
gradient = term_kin + term_mass + term_nonlinear + term_tension + term_centrifugal
return -gradient
# =============================================================================
# BOOST FUNCTION
# =============================================================================
def boost_soliton(Psi_soliton, z_grid, v_soliton, v=1.0):
"""
Apply velocity boost to stationary soliton.
Ψ_boosted = Ψ_soliton * exp(i * k * z)
where k = v_soliton / v (from dispersion relation)
"""
nz = len(z_grid)
nr = len(Psi_soliton) // nz
z_mesh = np.tile(z_grid.reshape(-1, 1), (1, nr))
k = v_soliton / v
phase = np.exp(1j * k * z_mesh)
Psi_boosted = Psi_soliton.reshape((nz, nr)) * phase
return Psi_boosted.ravel()
def initialize_velocity(Psi_boosted, z_grid, v_soliton, dz):
"""
Initialize Psi_dot = -v_soliton * ∂Ψ/∂z
"""
nz = len(z_grid)
nr = len(Psi_boosted) // nz
Psi_2d = Psi_boosted.reshape((nz, nr))
dPsi_dz = np.gradient(Psi_2d, dz, axis=0)
Psi_dot = -v_soliton * dPsi_dz
return Psi_dot.ravel()
# =============================================================================
# PHASE WAKE DIAGNOSTIC
# =============================================================================
def measure_asymmetry(Psi, z_grid, r_grid, dz, nr_mid_ratio=0.5):
"""
Measure forward/rear phase gradient asymmetry.
Returns:
asymmetry_ratio: forward_gradient / rear_gradient
forward_gradient: mean phase gradient ahead of soliton
rear_gradient: mean phase gradient behind soliton
high_k_power: fraction of power in high-frequency modes
"""
nz = len(z_grid)
nr = len(r_grid)
Psi_2d = Psi.reshape((nz, nr))
# Mid-radius index (avoid edges)
mid_r_idx = int(nr * nr_mid_ratio)
# Phase gradient along z at mid-radius
phase = np.angle(Psi_2d[:, mid_r_idx])
phase_gradient = np.abs(np.gradient(phase, dz))
# Find soliton center (peak amplitude at mid-radius)
amplitude = np.abs(Psi_2d[:, mid_r_idx])
center_z_idx = np.argmax(amplitude)
# Boundary safety
margin = min(20, nz // 10)
if center_z_idx < margin or center_z_idx > nz - margin:
return 0.0, 0.0, 0.0, 0.0
# Forward face (ahead, z > center)
forward_slice = phase_gradient[center_z_idx:center_z_idx + margin]
forward_gradient = np.mean(forward_slice)
# Rear face (behind, z < center)
rear_slice = phase_gradient[center_z_idx - margin:center_z_idx]
rear_gradient = np.mean(rear_slice)
# Asymmetry ratio
asymmetry_ratio = forward_gradient / max(rear_gradient, 1e-10)
# High-k power (ripples in phase gradient)
fft = np.abs(np.fft.fft(phase_gradient))
high_k_power = np.sum(fft[len(fft)//4:]) / np.sum(fft + 1e-15)
return asymmetry_ratio, forward_gradient, rear_gradient, high_k_power
# =============================================================================
# RUN VELOCITY SWEEP
# =============================================================================
print("\n" + "="*80)
print("RUNNING VELOCITY SATURATION TESTS")
print("="*80)
# Velocity fractions to test (as fraction of substrate speed v)
velocity_fractions = [0.3, 0.5, 0.7, 0.85, 0.95]
results = {}
# Real-time parameters
dt = 0.001
t_max = 50.0
n_steps = int(t_max / dt)
n_save = 200
n_save_print = n_save * 5
for v_frac in velocity_fractions:
v_soliton = v_frac * v
print(f"\n[Test: v_soliton = {v_frac:.2f}v = {v_soliton:.3f}]")
# Boost the soliton
Psi = boost_soliton(Psi_soliton, z_grid, v_soliton, v)
Psi_dot = initialize_velocity(Psi, z_grid, v_soliton, dz)
# Storage
times = []
energies = []
max_amps = []
asymmetry_ratios = []
high_k_powers = []
for step in range(n_steps):
# Velocity Verlet (symplectic, 2nd order)
acc = acceleration_real(Psi, L_2D, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m)
Psi_dot_half = Psi_dot + 0.5 * dt * acc
Psi_new = Psi + dt * Psi_dot_half
acc_next = acceleration_real(Psi_new, L_2D, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m)
Psi_dot_new = Psi_dot_half + 0.5 * dt * acc_next
Psi = Psi_new
Psi_dot = Psi_dot_new
# Save diagnostics
if step % n_save == 0:
t = step * dt
times.append(t)
max_amps.append(np.max(np.abs(Psi)))
E = compute_energy(Psi, L_2D, dV, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m)
energies.append(E)
# Measure wake asymmetry (only after initial transient, t > 5)
if t > 5.0:
asym, fwd, rev, hk = measure_asymmetry(Psi, z_grid, r_grid, dz)
asymmetry_ratios.append(asym)
high_k_powers.append(hk)
if step % n_save_print == 0 and step > 0:
print(f" t={t:6.2f}, max|Ψ|={max_amps[-1]:.4f}, E={E:.6e}")
# Aggregate results
energy_drift = (energies[-1] - energies[0]) / abs(energies[0]) * 100
final_max_amp = max_amps[-1]
amplitude_preservation = final_max_amp / max_amps[0] if max_amps[0] > 0 else 0
# Average asymmetry over last 5 measurements
if len(asymmetry_ratios) >= 5:
avg_asymmetry = np.mean(asymmetry_ratios[-5:])
avg_high_k = np.mean(high_k_powers[-5:])
else:
avg_asymmetry = 0.0
avg_high_k = 0.0
results[v_frac] = {
'velocity': v_soliton,
'energy_drift': energy_drift,
'final_max_amp': final_max_amp,
'amplitude_preservation': amplitude_preservation,
'asymmetry_ratio': avg_asymmetry,
'high_k_power': avg_high_k,
}
print(f" ✓ Energy drift = {energy_drift:.4f}%")
print(f" ✓ Amplitude preservation = {amplitude_preservation*100:.1f}%")
print(f" ✓ Asymmetry ratio = {avg_asymmetry:.4f}")
print(f" ✓ High-k power = {avg_high_k:.4f}")
# =============================================================================
# ANALYSIS AND PLOTTING
# =============================================================================
print("\n" + "="*80)
print("ANALYSIS: TEST 1 RESULTS")
print("="*80)
velocities = list(results.keys())
asymmetry_vals = [results[v]['asymmetry_ratio'] for v in velocities]
highk_vals = [results[v]['high_k_power'] for v in velocities]
drift_vals = [results[v]['energy_drift'] for v in velocities]
preserve_vals = [results[v]['amplitude_preservation'] for v in velocities]
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# Plot 1: Asymmetry vs velocity
axes[0, 0].plot(velocities, asymmetry_vals, 'ro-', linewidth=2, markersize=8)
axes[0, 0].axhline(1.0, color='k', linestyle='--', alpha=0.3, label='Symmetric')
axes[0, 0].axhline(1.2, color='g', linestyle=':', alpha=0.5, label='Support threshold')
axes[0, 0].set_xlabel('Velocity (fraction of v)')
axes[0, 0].set_ylabel('Asymmetry Ratio (forward/rear)')
axes[0, 0].set_title('Phase Wake Asymmetry vs Velocity')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].legend()
# Plot 2: High-k power vs velocity
axes[0, 1].plot(velocities, highk_vals, 'bs-', linewidth=2, markersize=8)
axes[0, 1].set_xlabel('Velocity (fraction of v)')
axes[0, 1].set_ylabel('High-k Power (ripple fraction)')
axes[0, 1].set_title('Wake Ripple Signature vs Velocity')
axes[0, 1].grid(True, alpha=0.3)
# Plot 3: Energy conservation
axes[1, 0].plot(velocities, drift_vals, 'g^-', linewidth=2, markersize=8)
axes[1, 0].axhline(0.0, color='k', linestyle='--', alpha=0.3)
axes[1, 0].set_xlabel('Velocity (fraction of v)')
axes[1, 0].set_ylabel('Energy Drift (%)')
axes[1, 0].set_title('Energy Conservation vs Velocity')
axes[1, 0].grid(True, alpha=0.3)
# Plot 4: Amplitude preservation
axes[1, 1].plot(velocities, preserve_vals, 'm^-', linewidth=2, markersize=8)
axes[1, 1].axhline(0.8, color='r', linestyle='--', alpha=0.5, label='80% threshold')
axes[1, 1].axhline(1.0, color='k', linestyle='-', alpha=0.3)
axes[1, 1].set_xlabel('Velocity (fraction of v)')
axes[1, 1].set_ylabel('Amplitude Preservation')
axes[1, 1].set_title('Soliton Integrity vs Velocity')
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].legend()
plt.tight_layout()
plt.savefig('test_1_velocity_saturation.png', dpi=150)
print("\n✓ Plot saved: test_1_velocity_saturation.png")
# =============================================================================
# VERDICT
# =============================================================================
print("\n" + "="*80)
print("VERDICT")
print("="*80)
max_asym = max(asymmetry_vals)
trend_increasing = asymmetry_vals[-1] > asymmetry_vals[0] if len(asymmetry_vals) >= 2 else False
if max_asym > 1.2 and trend_increasing:
print("\n✅ POSITIVE EVIDENCE FOR VELOCITY SATURATION")
print(f" Asymmetry ratio increases with velocity ({asymmetry_vals[0]:.3f} → {asymmetry_vals[-1]:.3f})")
print(f" Maximum asymmetry = {max_asym:.3f} > 1.2 threshold")
print(" Phase wake shows forward/rear asymmetry consistent with directional saturation.")
verdict = "SUPPORTED"
elif max_asym < 1.1:
print("\n⚠️ INCONCLUSIVE EVIDENCE")
print(f" Asymmetry ratio small (max = {max_asym:.3f} < 1.1)")
print(" Phase wake is approximately symmetric or too weak to detect.")
print(" May need higher velocity or parameter adjustment.")
verdict = "INCONCLUSIVE"
else:
print("\n❌ NEGATIVE EVIDENCE FOR VELOCITY SATURATION")
print(f" Asymmetry ratio max = {max_asym:.3f}, trend increasing = {trend_increasing}")
print(" Phase wake lacks predicted structure or shows no clear velocity dependence.")
verdict = "CHALLENGED"
print(f"\n📊 VERDICT: {verdict}")
print(f"\nFinal results:")
for v_frac in velocities:
print(f" v={v_frac:.2f}v: asymmetry={results[v_frac]['asymmetry_ratio']:.4f}, "
f"high-k={results[v_frac]['high_k_power']:.4f}, "
f"energy drift={results[v_frac]['energy_drift']:.4f}%")
# =============================================================================
# SAVE RESULTS
# =============================================================================
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
test_name = f"test_1_velocity_saturation_{timestamp}"
results_json = {
'test_name': 'Test 1: Velocity Saturation & Dispersive Phase Wake',
'timestamp': timestamp,
'soliton_source': soliton_file,
'parameters': {
'v': v, 'mu': mu, 'lambda': lam, 'kappa': kappa,
'S_max': S_max, 'Psi_sat': Psi_sat, 'm': m,
'dt': dt, 't_max': t_max,
},
'velocity_fractions': velocities,
'results': {str(k): v for k, v in results.items()},
'verdict': verdict,
'summary': {
'max_asymmetry': float(max_asym),
'trend_increasing': trend_increasing,
}
}
with open(f"{test_name}_results.json", 'w') as f:
json.dump(results_json, f, indent=2, default=str)
print(f"\n✓ Results saved: {test_name}_results.json")
print("\n" + "="*80)
print("TEST 1 COMPLETE")
print("="*80)
