PerturbedEquilibrium - BUGFIX! - Regularize stored xi_s in kinetic Clebsch displacements

examples/Solovev_kinetic_calculated_example/forcing.dat

@@ -0,0 +1,5 @@
+# Forcing data for perturbed equilibrium calculations
+# normalization: normal_field_T
+# Format: n  m  amplitude_real  amplitude_imag
+# Single mode test case: n=1, m=2, amplitude=1e-4 T
+1  2  1e-4  0.0

examples/Solovev_kinetic_calculated_example/gpec.toml

@@ -2,6 +2,8 @@
 # Exercises the real KineticForces (NTV) physics path at full strength; kinetic regression
 # fixture (et[1]≈1.894-1.525i). The equilibrium is generated analytically from the embedded
 # [SOL_INPUT] section, with n·T tied to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
+# Also the only deck that runs [PerturbedEquilibrium] on a kinetic solve, so it is what covers
+# the kinetic branch of compute_clebsch_displacements.
 
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
@@ -25,6 +27,19 @@ dw = 0.5                                # Triangularity parameter for wall shape
 tw = 0.05                               # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true                   # Equal arc length distribution of nodes on wall
 
+[ForcingTerms]
+forcing_data_file = "forcing.dat"       # Manual mode table (n, m, amplitude_real, amplitude_imag)
+forcing_data_format = "ascii"           # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
+
+[PerturbedEquilibrium]
+fixed_boundary = false                  # Use fixed boundary conditions
+output_eigenmodes = true                # Output eigenmode fields as b-fields
+compute_response = true                 # Compute plasma response to forcing
+compute_singular_coupling = true        # Compute singular layer coupling metrics
+verbose = false                         # Enable verbose logging
+write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
+reg_spot = 0.05                         # Regularization width for singular surfaces (0 = disabled)
+
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes

regression-harness/cases/solovev_kinetic_calculated.toml

@@ -111,6 +111,29 @@ extract = "value"
 label = "npert"
 noise_threshold = 0
 
+# Perturbed equilibrium — guards the kinetic branch of compute_clebsch_displacements.
+# This is the only case that runs PE on a kinetic solve; without it that branch is untested.
+[quantities.pe_clebsch_alpha_norm]
+h5path = "PerturbedEquilibrium/Response/xi_clebsch_alpha"
+type = "complex_matrix"
+extract = "norm"
+label = "||xi_clebsch_alpha||"
+noise_threshold = 1e-12
+
+[quantities.pe_clebsch_psi1_norm]
+h5path = "PerturbedEquilibrium/Response/dxi_clebsch_psidpsi"
+type = "complex_matrix"
+extract = "norm"
+label = "||dxi_clebsch_psi/dpsi||"
+noise_threshold = 1e-12
+
+[quantities.pe_b_theta_reg_norm]
+h5path = "PerturbedEquilibrium/Response/Jb_theta_reg"
+type = "complex_matrix"
+extract = "norm"
+label = "||Jb_theta_reg||"
+noise_threshold = 1e-12
+
 # Runtime (special: not from H5)
 [quantities.runtime]
 h5path = ""

src/PerturbedEquilibrium/FieldReconstruction.jl

@@ -19,7 +19,8 @@ where χ₁ = 2π·Ψ₀ [Park Phys. Plasmas 14, 052110 (2007) eq. 8-10].
 Clebsch displacement components for PENTRC (matches Fortran gpout_xclebsch):
     ξ^ψ         = xsp_mn    (unregularized)
     ∂ξ^ψ/∂ψ    = xmp1_mn   (regularized: xsp1 * singfac²/(singfac² + reg_spot²))
-    ξ^α         = xms_mn    (regularized: -A⁻¹(B·xmp1 + C·xsp), divided by χ₁ in output)
+    ξ^α         = xms_mn    (regularized: -A⁻¹(B·xmp1 + C·xsp) ideal, or ξ_s scaled by the same
+                             singfac factor when kinetic; divided by χ₁ in output)
 
 Contravariant displacement from Jacobian convolution (matches Fortran gpeq_contra):
     ξ^ψ·J(m) = Σ_dm jmat(dm) · xsp(m+dm)
@@ -334,8 +335,9 @@ Matches Fortran gpeq_sol regularization + gpout_xclebsch output convention:
 
 When reg_spot=0, clebsch_psi1 = xi_psi1 and clebsch_alpha = xi_s/χ₁ (no regularization).
 
-The regularized xms is computed as -A⁻¹(B·xmp1 + C·xsp) matching Fortran gpeq_sol,
-where A, B, C are the stability matrices evaluated at each ψ from `mats`.
+Ideal runs re-solve the regularized xms = -A⁻¹(B·xmp1 + C·xsp) from `mats.ideal`. Kinetic runs
+instead scale the stored ξ_s by the same singfac factor and invert nothing — the kinetic A is
+non-Hermitian and never re-inverted here.
 """
 function compute_clebsch_displacements(
     xi_psi_modes::Matrix{ComplexF64},
@@ -364,8 +366,20 @@ function compute_clebsch_displacements(
         return clebsch_psi, clebsch_psi1, clebsch_alpha
     end
 
-    # A/B/C of the active model, matching what the ODE integrated.
-    active_mats = mats.kinetic === nothing ? mats.ideal : mats.kinetic
+    # Kinetic runs regularize the stored ξ_s directly
+    if mats.kinetic !== nothing
+        hint = Ref(1)
+        for ipsi in 1:npsi
+            q = equil.profiles.q_spline(psi_grid[ipsi]; hint=hint)
+            for ipert in 1:mpert
+                singfac = (mlow + ipert - 1) - nn * q
+                reg_factor = singfac^2 / (singfac^2 + reg_spot^2)
+                clebsch_psi1[ipsi, ipert] = xi_psi1_modes[ipsi, ipert] * reg_factor
+                clebsch_alpha[ipsi, ipert] = xi_s_modes[ipsi, ipert] * reg_factor / chi1
+            end
+        end
+        return clebsch_psi, clebsch_psi1, clebsch_alpha
+    end
 
     # Per-thread workspaces: matrix ops and spline hints are not safe to share across threads.
     # Size by maxthreadid() and index by threadid() under :static scheduling (GPEC convention).
@@ -398,21 +412,16 @@ function compute_clebsch_displacements(
             xmp1_vec[ipert] = clebsch_psi1[ipsi, ipert]
         end
 
-        # Compute regularized xms = -A⁻¹(B·xmp1 + C·xsp) (matches Fortran gpeq_sol)
-        # Evaluate stability matrices at this psi
-        active_mats.A_spline(view(amat, :), psi_norm; hint=hint)
-        active_mats.B_spline(view(bmat, :), psi_norm; hint=hint)
-        active_mats.C_spline(view(cmat_buf, :), psi_norm; hint=hint)
+        # Compute regularized xms = -A⁻¹(B·xmp1 + C·xsp)
+        mats.ideal.A_spline(view(amat, :), psi_norm; hint=hint)
+        mats.ideal.B_spline(view(bmat, :), psi_norm; hint=hint)
+        mats.ideal.C_spline(view(cmat_buf, :), psi_norm; hint=hint)
 
         # xms = -(A\B)*xmp1 - (A\C)*xsp
         xsp_vec = view(xi_psi_modes, ipsi, :)
         mul!(xms_vec, bmat, xmp1_vec)                     # xms = B*xmp1
         mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im)  # xms += C*xsp
-        # cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by
-        # active_mats.A_spline each surface), avoiding a fresh factorization allocation per surface.
-        # NOTE: this assumes the ideal A (positive-definite Newcomb kinetic-energy form). The
-        # kinetic A is non-Hermitian and needs an LU, as compute_node_xi_s! does — see the
-        # `active_mats` binding above.
+        # factorize in-place to avoid allocations
         amat_fact = cholesky!(Hermitian(amat, :L))
         ldiv!(amat_fact, xms_vec)                          # xms = A\(B*xmp1 + C*xsp)
         xms_vec .*= -1                                     # xms = -A\(B*xmp1 + C*xsp)
@@ -985,10 +994,12 @@ function _apply_rzphi_transform(
     # Per-thread scratch (the immutable `ft` functor and `geom` are shared read-only): θ-space
     # transform inputs/outputs (length mtheta) and mode-space forward-DFT outputs (length mpert),
     # so the DFTs run in place with no per-surface allocation.
-    bufs = [(R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta),
-             psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta),
-             Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert))
-            for _ in 1:Threads.maxthreadid()]
+    bufs = [
+        (R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta),
+            psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta),
+            Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert))
+        for _ in 1:Threads.maxthreadid()
+    ]
 
     Threads.@threads :static for ipsi in 1:npsi
         buf = bufs[Threads.threadid()]