Here are the exact functions from your successful Test‑0 run.
Here are the exact functions from your successful Test‑0 run.
✅ STEP 1 — Your exact Test‑0 functions
1. The Real‑Time Acceleration (from Test‑0B)
python
def acceleration_real(Psi, L_2D, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m):
"""
Real-time acceleration = -δE/δΨ* (canonical hyperbolic PDE)
This matches compute_gradient exactly.
"""
psi_sq = np.abs(Psi)**2
# Saturation and its derivative
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)
# Variational gradient components (IDENTICAL to compute_gradient)
term_kin = -v**2 * (L_2D @ Psi)
term_mass = mu * Psi
term_nonlinear = lam * psi_sq * Psi
term_tension_exact = 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_exact + term_centrifugal
# Canonical real-time acceleration = -gradient
return -gradient
2. The Energy Functional (from Test‑0A)
python
def compute_energy(Psi, L_2D, dV, r_mesh_2d, v, mu, lam, kappa, S_max, Psi_sat, m):
"""Compute energy functional E[Ψ]."""
psi_sq = np.abs(Psi)**2
# Kinetic energy (time derivative = 0 for stationary state)
kin_time = 0.0
# Gradient energy: ½v²∫|∇Ψ|² dV = -½v²∫ Ψ* L Ψ dV (since L is negative-definite)
kin_grad = -0.5 * v**2 * np.real(np.sum(np.conj(Psi) * (L_2D @ Psi) * dV))
# Potential energies
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)
# Centrifugal potential (positive definite barrier)
pot_centrifugal = 0.5 * v**2 * m**2 * np.sum(psi_sq / (r_mesh_2d**2 + 1e-12) * dV)
return (kin_time + kin_grad + pot_mass + pot_nonlinear + pot_tension + pot_centrifugal).real
3. The 2D Laplacian Operator (from Test‑0A)
python
def build_radial_operator(r_grid, dr):
"""Self-adjoint radial Laplacian: L_r = W_r⁻¹ M"""
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):
"""Self-adjoint axial Laplacian (symmetric stencil)"""
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
# Build 2D operator:
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)
4. The parameters that worked
python
v = 1.0
mu = -1.0
lam = 0.40
kappa = 0.20
m = 1
S_max = 2.0
Psi_sat = 0.8
dt_real = 0.001
📋 Summary of what you've provided
Function Purpose Status
acceleration_real() Real-time RHS (∂²Ψ/∂t² = -δE/δΨ*) ✅ Exact from Test‑0B
compute_energy() Energy functional ✅ Exact from Test‑0A
L_2D 2D self-adjoint Laplacian ✅ Exact from Test‑0A
Parameters v, mu, lam, kappa, S_max, Psi_sat, m ✅ Validated