PerturbedEquilibrium - BUGFIX! - Converge the singular coupling with singfac_min = 1e-6 in the forward-integrator decks

examples/DIIID-like_error_field_example/gpec.toml

@@ -61,7 +61,7 @@ kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, o
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
 save_interval = 1              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
-singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6             # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3         # Column-norm threshold that triggers solution renormalization
 
 # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details)

examples/DIIID-like_ideal_example/gpec.toml

@@ -56,7 +56,7 @@ kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, o
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
 save_interval = 1              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
-singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6             # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3         # Column-norm threshold that triggers solution renormalization
 
 # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details)

examples/DIIID-like_ideal_multi_n_example/gpec.toml

@@ -58,7 +58,7 @@ kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, o
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
 save_interval = 1              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
-singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6             # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3         # Column-norm threshold that triggers solution renormalization
 
 # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details)

examples/DIIID-like_truncation_example/gpec.toml

@@ -34,7 +34,7 @@ delta_mlow = 8                 # Expands lower bound of Fourier harmonics
 delta_mhigh = 8                # Expands upper bound of Fourier harmonics
 
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6             # Fractional distance from rational q at which ideal jump enforced
 
 # Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details)
 integrator                = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)

examples/Solovev_ideal_example/gpec.toml

@@ -74,7 +74,7 @@ mtheta_vacuum = 960           # Number of points used in splines over poloidal a
 kinetic_source = "fixed"      # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6            # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3        # Column-norm threshold that triggers solution renormalization
 save_interval = 1             # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 

examples/Solovev_kinetic_NTV_example/gpec.toml

@@ -71,7 +71,7 @@ mtheta_vacuum = 960           # Number of points used in splines over poloidal a
 kinetic_source = "fixed"      # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6            # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3        # Column-norm threshold that triggers solution renormalization
 save_interval = 1             # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 

examples/Solovev_kinetic_calculated_example/gpec.toml

@@ -68,7 +68,7 @@ kinetic_source = "calculated" # Kinetic matrix source — exercises KineticForce
 kinetic_factor = 1.0          # Full-strength kinetic matrices (the "calculated" path is the real physics; no perturbation scaling)
 integrator = "forward"        # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6            # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3        # Column-norm threshold that triggers solution renormalization
 write_outputs_to_HDF5 = true  # Write outputs to the HDF5 file
 verbose = false               # Enable verbose logging

examples/Solovev_kinetic_multiion_example/gpec.toml

@@ -72,7 +72,7 @@ mtheta_vacuum = 960           # Number of points used in splines over poloidal a
 kinetic_source = "fixed"      # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6            # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3        # Column-norm threshold that triggers solution renormalization
 save_interval = 1             # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 

examples/a10_kinetic_example/gpec.toml

@@ -51,7 +51,7 @@ kinetic_source = "calculated"  # Kinetic matrix source: "fixed" test matrices, o
 kinetic_factor = 1.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode; 1.0 = full strength)
 integrator = "forward"         # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+singfac_min = 1e-6             # Fractional distance from rational q at which ideal jump enforced
 renorm_threshold = 1e3         # Column-norm threshold that triggers solution renormalization
 write_outputs_to_HDF5 = true   # Write outputs to the HDF5 file
 force_termination = true       # Stop after force-free states — the kinetic stability eigenvalues are the result of interest

src/ForceFreeStates/CoreTypes.jl

@@ -131,7 +131,7 @@ gpec.toml.
   - `renorm_threshold::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e3`. The trigger compares column-norm growth, not conditioning, so near-parallel columns can leave an error in the free-boundary energies that does not shrink with `eulerlagrange_tolerance`.
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
-  - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced
+  - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced. The default 1e-4 suits the Riccati Δ′, whose asymptotic matching gains a spurious imaginary part as it shrinks; `"forward"` runs feeding PerturbedEquilibrium need about 1e-6 for the singular coupling to converge at the outer rational surfaces.
   - `truncate_at_rational_offset::Bool` - Truncate the integration domain at `(last_rational_q + rational_offset_fraction) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`rational_offset_fraction = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `truncate_at_rational_offset = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `truncate_at_rational_offset = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + rational_offset_fraction / n" depends on which `n`); when `truncate_at_rational_offset = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`.
   - `psilim_from_layer_overlap::Bool` - Cap the integration domain where adjacent rational surfaces' resistive layers overlap, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9: past that point no surface retains a well-separated inner region. Applied as an upper bound on `qlim` at the top of `sing_lim!`, so `rational_offset_fraction` / `qhigh` still select the final surface from inside it. The scan runs whenever the run's `[SLAYER]` section names a readable `profile_file` — the same profiles and χ⊥(ψ) the SLAYER analysis uses — and is recorded under `ForceFreeStates/LayerOverlap/` either way; this flag only decides whether it constrains the domain. It resolves that file from the run's TOML, so the programmatic `solve(prob, alg)` entry point cannot run it and the flag is inert there (with a warning). Off on every shipped deck; with the shipped DIII-D-like profiles the overlap point (ψ ≈ 0.9985) lies outside the 0.995 decks' domain but inside the SLAYER deck's 0.9995. See the integration-domain truncation section of the stability documentation. Default `false`.
   - `rational_offset_fraction::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `truncate_at_rational_offset` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here.

src/PerturbedEquilibrium/SingularCoupling.jl

@@ -225,6 +225,15 @@ Build the singular-layer coupling matrices of `resp.ffs` onto `resp`, matching F
 resonant area-weighted field, current, (w/2)², penetrated field and Δ' per unit forcing mode, plus
 the rational-surface metadata and the edge-to-solution coefficients `C_coeffs`. The matrices are
 stored acting on b̃_x, with flux-space copies in `resp.coupling_flux` for applying a forcing.
+
+## Normalization
+
+  - Rows: the resonant and penetrated rows are a resonant flux divided by that surface's own area
+    `rational_area[row]`, so both are fields in T. `C_resonant_current`, `C_island_width_sq` and
+    `C_delta_prime` carry their own normalizations.
+  - Columns: each field-space coupling matrix is right-multiplied by the control-surface `flux_conform`, so it
+    acts on the root-area-weighted field b̃ in T and its singular values do not depend on `jac_type`.
+    Applied vectors use the flux-space copies in `resp.coupling_flux` on the run's own forcing flux.
 """
 function build_coupling_operators!(resp::PlasmaResponse; verbose::Bool=false)
     ffs = resp.ffs