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)

Popular posts from this blog

THE GOLDEN BALLROOM/BUNKER

Conceptual Summary #2: (∂t2​S−c2∇2S+βS3)=σ(x,t)⋅FR​(C[Ψ])

ICE PROUDLY ANNOUNCES NEW “ELITE” TASK FORCE COMMANDER JEREMY DEWITTE