From f1cd7b7b31bad3aa924db6e889f27a604d40ab69 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 18 Aug 2026 15:18:09 -0400 Subject: [PATCH] PE - BUGFIX! - Regularize stored xi_s in kinetic Clebsch displacements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit compute_clebsch_displacements factorized the active A with cholesky!(Hermitian(amat, :L)). In a kinetic run that A is the kinetic A (amat + kwmat[:,:,1] + ktmat[:,:,1]), which is non-Hermitian, so the upper triangle was silently discarded. The run completed without error and returned wrong regularized quantities. Mirror Fortran gpeq.f:117-123: under kin_flag GPEC scales the stored xi_s by singfac^2/(singfac^2 + reg_spot^2) and inverts nothing. It never holds a kinetic A here — idcon_matrix rebuilds A analytically from the metric tensors, ideal-only. Our stored xi_s is the right analogue; compute_node_xi_s! already builds it with the correct LU. The ideal branch is unchanged and now reads mats.ideal explicitly. Kinetic PE results move: xi_clebsch_alpha by 18%, and Jb_theta_reg/Jb_zeta_reg and Jxi_theta_reg/Jxi_zeta_reg by order their own magnitude. Unregularized quantities and reg_spot = 0 are unaffected. The branch had never executed in any example or regression case, so add [ForcingTerms]/[PerturbedEquilibrium] to Solovev_kinetic_calculated_example and track xi_clebsch_alpha, dxi_clebsch_psidpsi and Jb_theta_reg in its regression case. Co-Authored-By: Claude Opus 5 --- .../forcing.dat | 5 ++ .../gpec.toml | 15 ++++++ .../cases/solovev_kinetic_calculated.toml | 23 +++++++++ .../FieldReconstruction.jl | 49 ++++++++++++------- 4 files changed, 73 insertions(+), 19 deletions(-) create mode 100644 examples/Solovev_kinetic_calculated_example/forcing.dat diff --git a/examples/Solovev_kinetic_calculated_example/forcing.dat b/examples/Solovev_kinetic_calculated_example/forcing.dat new file mode 100644 index 000000000..42eb94703 --- /dev/null +++ b/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 diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 827582a7c..4d7fe2704 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/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 diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/solovev_kinetic_calculated.toml index 172416458..2ffe3cbe0 100644 --- a/regression-harness/cases/solovev_kinetic_calculated.toml +++ b/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 = "" diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 8c349ae1d..b4f66087b 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/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()]