FRCMFD-v2: TEST 0A — TOROIDAL SOLITON GENERATOR (m=1)
"""
FRCMFD-v2: TEST 0A — TOROIDAL SOLITON GENERATOR (m=1)
======================================================
Generates a stationary toroidal soliton with m=1 (non-zero angular winding).
The centrifugal term -v²m²/r² Ψ is ACTIVE in this configuration.
Output:
test_0A_toroidal_m1_.npz
"""
import numpy as np
import scipy.sparse as sp
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')
print("="*80)
print("FRCMFD-v2: TEST 0A — TOROIDAL SOLITON GENERATOR (m=1)")
print("="*80)
# =============================================================================
# PARAMETERS (Toroidal configuration: m=1)
# =============================================================================
# Grid parameters
r_max = 40.0
z_max = 80.0 # Extended domain for toroidal shape
nr_full = 200
nz = 400
dr = r_max / nr_full
dz = z_max / nz
# Radial grid (exclude r=0)
r_grid_full = np.linspace(dr/2, r_max - dr/2, nr_full)
r_grid = r_grid_full[1:]
nr = len(r_grid)
# Axial grid (symmetric about 0)
z_grid = np.linspace(-z_max/2, z_max/2, nz)
# Physical parameters (toroidal regime)
v_substrate = 1.0 # Substrate speed
mu = -1.0 # Linear restoring (negative mass)
lam = 0.4 # Nonlinear focusing (toroidal regime)
kappa = 0.2 # Tension coupling (toroidal regime)
m = 1 # Winding mode (TOROIDAL — active centrifugal term)
S_max = 2.0 # Maximum tension
Psi_sat = 0.8 # Saturation threshold
# Relaxation parameters
dtau = 0.005
tau_max = 400.0
n_steps = int(tau_max / dtau)
n_save = 200
n_save_print = n_save * 5
print(f"\n[Grid] nr={nr}, nz={nz}, DOF={nr*nz:,}")
print(f" Lr = {r_grid[-1]:.2f}, Lz = {z_grid[-1]-z_grid[0]:.2f}")
print(f" dr={dr:.4f}, dz={dz:.4f}")
print(f"\n[Parameters] (TOROIDAL m=1)")
print(f" v={v_substrate}, mu={mu}, lam={lam}, kappa={kappa}")
print(f" m={m} (centrifugal term ACTIVE: -v²m²/r² Ψ)")
print(f" S_max={S_max}, Psi_sat={Psi_sat}")
print(f" dtau={dtau}, tau_max={tau_max}, steps={n_steps:,}")
# =============================================================================
# BUILD OPERATORS
# =============================================================================
print("\n[Building 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[0] = dr * 0.5
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_periodic(nz, dz):
main_diag = np.ones(nz) * (-2.0 / dz**2)
off_diag = np.ones(nz - 1) / dz**2
L_z = sp.diags([off_diag, main_diag, off_diag], [-1, 0, 1], format='lil')
L_z[0, nz-1] = 1.0 / dz**2
L_z[nz-1, 0] = 1.0 / dz**2
W_z = sp.diags(np.ones(nz) * dz, format='csr')
return L_z.tocsr(), W_z
L_r, W_r = build_radial_operator(r_grid, dr)
L_z, W_z = build_axial_operator_periodic(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"✓ L_2D shape: {L_2D.shape}, nnz={L_2D.nnz:,}")
# =============================================================================
# INITIAL GUESS (Toroidal ring)
# =============================================================================
print("\n[Generating initial guess (toroidal ring)...]")
R0 = 6.0
sigma = 2.0
Psi_2d = np.zeros((nz, nr), dtype=complex)
for j, z in enumerate(z_grid):
for i, r in enumerate(r_grid):
rho = np.sqrt((r - R0)**2 + z**2)
Psi_2d[j, i] = np.exp(-rho**2 / (2 * sigma**2))
Psi_2d = Psi_2d / np.max(np.abs(Psi_2d))
# Toroidal phase winding
theta = np.arctan2(z_grid[:, None], r_grid[None, :])
Psi_2d = Psi_2d * np.exp(1j * m * theta)
Psi = Psi_2d.ravel()
print(f"✓ Initial max|Ψ| = {np.max(np.abs(Psi)):.4f}")
print(f"✓ Toroidal winding m={m} applied")
# =============================================================================
# ENERGY AND GRADIENT
# =============================================================================
def compute_energy(Psi):
psi_sq = np.abs(Psi)**2
kin_grad = -0.5 * v_substrate**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_substrate**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 compute_gradient(Psi):
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_substrate**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_substrate**2 * m**2 * Psi / (r_mesh_2d**2 + 1e-12)
return term_kin + term_mass + term_nonlinear + term_tension + term_centrifugal
def residual_norm(grad):
return np.sqrt(np.sum(np.abs(grad)**2 * dV))
# =============================================================================
# RELAXATION
# =============================================================================
print("\n" + "="*80)
print("RUNNING IMAGINARY-TIME RELAXATION")
print("∂Ψ/∂τ = -δE/δΨ*")
print("="*80)
times = []
energies = []
max_amps = []
center_amps = []
residuals = []
Psi_current = Psi.copy()
for step in range(n_steps):
grad = compute_gradient(Psi_current)
Psi_current = Psi_current - dtau * grad
if step % n_save == 0:
tau = step * dtau
times.append(tau)
E = compute_energy(Psi_current)
energies.append(E)
max_amp = np.max(np.abs(Psi_current))
max_amps.append(max_amp)
Psi_2d = Psi_current.reshape((nz, nr))
mid_z = nz // 2
r_idx = np.argmin(np.abs(r_grid - R0))
center_amp = np.abs(Psi_2d[mid_z, r_idx])
center_amps.append(center_amp)
resid = residual_norm(grad)
residuals.append(resid)
if step % n_save_print == 0 and step > 0:
print(f" τ={tau:6.2f}, E={E:.6e}, max|Ψ|={max_amp:.4f}, resid={resid:.2e}, center|Ψ|={center_amp:.4f}")
if step > n_save * 5 and len(residuals) > 5:
if residuals[-1] < 1e-8:
print(f"\n✓ Convergence reached at τ={step*dtau:.2f}")
break
print(f"\n✓ Relaxation complete after {len(energies)} saved steps")
print(f" Final max|Ψ| = {max_amps[-1]:.4f}")
print(f" Final center|Ψ| = {center_amps[-1]:.4f}")
print(f" Final energy = {energies[-1]:.6e}")
print(f" Final residual = {residuals[-1]:.2e}")
Psi_soliton = Psi_current
# =============================================================================
# SAVE
# =============================================================================
print("\n[Verifying and saving toroidal soliton...]")
expected_size = nr * nz
actual_size = len(Psi_soliton)
print(f" Expected flattened size: {expected_size}")
print(f" Actual size: {actual_size}")
assert actual_size == expected_size, f"Shape mismatch: {actual_size} vs {expected_size}"
try:
test_reshape = Psi_soliton.reshape((nz, nr))
print(f" Reshape test: {test_reshape.shape} = (nz={nz}, nr={nr}) ✅")
except Exception as e:
raise ValueError(f"Reshape failed: {e}")
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
soliton_file = f"test_0A_toroidal_m1_{timestamp}.npz"
np.savez(
soliton_file,
Psi_soliton=Psi_soliton,
r_grid=r_grid,
z_grid=z_grid,
dr=dr,
dz=dz,
v=v_substrate,
mu=mu,
lam=lam,
kappa=kappa,
m=m,
S_max=S_max,
Psi_sat=Psi_sat,
)
print(f"\n✅ TOROIDAL SOLITON SAVED: {soliton_file}")
print(f" m={m} (centrifugal term ACTIVE)")
print(f" Peak amplitude: {np.max(np.abs(Psi_soliton)):.4f}")
print(f" Energy: {energies[-1]:.6e}")
print(f" Residual: {residuals[-1]:.2e}")
print("\n" + "="*80)
print("TEST 0A (TOROIDAL m=1) COMPLETE")
print("="*80)
print(f"\nUse this soliton in Test 2 with:")
print(f"soliton_file = \"{soliton_file}\"")
print("v_boost = 0.30 * v_substrate")
print("k_boost = v_boost / v_substrate")
print("="*80)
