vB2 - FRCMFD-v2: TEST 0A — SOLITON RELAXATION (EXPANDED DOMAIN) revised
# ============================================================
# FRCMFD-v2: TEST 0A — SOLITON RELAXATION (EXPANDED DOMAIN)
# ============================================================
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 — SOLITON RELAXATION (EXPANDED DOMAIN)")
print("="*80)
# =============================================================================
# GRID SETUP — EXPANDED AXIAL DOMAIN
# =============================================================================
# Axial domain expanded for tail isolation
z_min = -80.0
z_max = 80.0
nz = 400
z_grid = np.linspace(z_min, z_max, nz)
dz = z_grid[1] - z_grid[0]
# Radial domain (match your existing setup)
nr = 199
r_max = 20.0
r_grid = np.linspace(0.0, r_max, nr)
dr = r_grid[1] - r_grid[0]
# =============================================================================
# MODEL PARAMETERS (MATCH TEST 2)
# =============================================================================
v = 1.0
mu = -1.0
lam = 1.0
kappa = 1.0
m = 0
S_max = 1.0
Psi_sat = 1.0
# =============================================================================
# BUILD OPERATORS (SAME STRUCTURE AS TEST 2)
# =============================================================================
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_operators_periodic(nz, dz):
main_l = np.full(nz, -2.0 / dz**2)
off_l = np.full(nz - 1, 1.0 / dz**2)
L_z = sp.diags([off_l, main_l, off_l], [-1, 0, 1], format="lil")
L_z[0, -1] = 1.0 / dz**2
L_z[-1, 0] = 1.0 / dz**2
off_d = np.full(nz - 1, 1.0 / (2.0 * dz))
D_z = sp.diags([-off_d, off_d], [-1, 1], format="lil")
D_z[0, -1] = -1.0 / (2.0 * dz)
D_z[-1, 0] = 1.0 / (2.0 * dz)
return L_z.tocsr(), D_z.tocsr()
L_r, W_r = build_radial_operator(r_grid, dr)
L_z, D_z = build_axial_operators_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) + sp.kron(L_z, I_r)
W_2D = sp.kron(sp.diags(np.ones(nz)*dz), W_r)
dV = W_2D.diagonal() * 2 * np.pi
r_mesh_2d = np.tile(r_grid, nz)
# =============================================================================
# ENERGY FUNCTIONAL (MATCH TEST 2)
# =============================================================================
def compute_energy(Psi):
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 float((kin_grad + pot_mass + pot_nonlinear + pot_tension + pot_centrifugal).real)
# =============================================================================
# FRCMFD-v2 FULL PHYSICS RELAXATION ENGINE
# =============================================================================
Psi = np.exp(-((z_grid.reshape(-1,1))**2 + (r_grid.reshape(1,-1))**2) / 20.0).ravel()
dt_relax = 0.01
n_steps = 15000
print("\n[Relaxing Soliton on Expanded Domain via FRCMFD-v2 Operators...]")
for step in range(n_steps):
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)
FRCMFD_residual = -(term_kin + term_mass + term_nonlinear + term_tension + term_centrifugal)
Psi = Psi + dt_relax * FRCMFD_residual
if step % 10 == 0:
norm = np.sqrt(np.sum(np.abs(Psi)**2 * dV))
if norm > 1e-12:
Psi = Psi * (10.0 / norm)
if step % 1500 == 0:
E = compute_energy(Psi)
print(f" step={step:5d}, E={E:.6e}, Peak Amplitude={np.max(np.abs(Psi)):.4f}")
# =============================================================================
# SAVE SOLITON
# =============================================================================
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
filename = f"test_0A_soliton_expanded_{timestamp}.npz"
np.savez(
filename,
Psi_soliton=Psi,
r_grid=r_grid,
z_grid=z_grid,
dr=dr,
dz=dz,
v=v,
mu=mu,
lam=lam,
kappa=kappa,
m=m,
S_max=S_max,
Psi_sat=Psi_sat
)
print(f"\n✓ Soliton saved to: {filename}")
print("="*80)