diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index 8282d285f3..cd6d045ead 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/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) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 2f1345eca3..d744632b19 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/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) diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index e4cb5c3434..8abd245116 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/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) diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml index 22001feb1e..1131d95513 100644 --- a/examples/DIIID-like_truncation_example/gpec.toml +++ b/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) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index eded9dce8f..10a59ae27a 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/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. diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 024ad168a9..86e19fd488 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/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. diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 60a9fe8b6f..de462bb440 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/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 diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index 679ecd4a9a..066b108b1d 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/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. diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index fe8544c9af..7b63274bde 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/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 diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5703ae0a58..5d8565cce9 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/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. diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index ad4132dcc9..975f216eff 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/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