ForceFreeStates - BUGFIX! - ⚠️ Keep both crossings of a rational surface on reverse-shear profiles

examples/LAR_reversed_shear_example/gpec.toml

@@ -0,0 +1,44 @@
+# Large-aspect-ratio equilibrium with reversed shear: n=1 ideal stability.
+# A hollow current profile makes q fall from 2.40 on the axis to a minimum of 1.82 at ψ_N ≈ 0.37
+# and rise to 2.84 at the edge, so q = 2 is crossed twice, with q′ < 0 on the inner crossing and
+# q′ > 0 on the outer one. Regression fixture for keeping both crossings as separate surfaces; the
+# equilibrium is generated analytically from the embedded [LAR_INPUT] section, with no wall.
+
+[Equilibrium]
+eq_type = "lar"                         # Type of the input 2D equilibrium file
+jac_type = "hamada"                     # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
+grid_type = "core_edge_packed"          # Radial grid packing type
+psilow = 0.01                           # Lower limit of normalized poloidal flux
+psihigh = 0.99                          # Upper limit of normalized poloidal flux
+mpsi = 128                              # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
+mtheta = 256                            # Number of poloidal grid points
+
+[Wall]
+shape = "nowall"                        # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+
+[ForceFreeStates]
+vac_flag = true                         # Compute plasma, vacuum, and total energies for free-boundary modes
+qlow = 1.02                             # Integration initiated at q determined by min(q0, qlow)
+qhigh = 1e3                             # Integration terminated at q limit determined by min(qa, qhigh)
+nn_low = 1                              # Smallest toroidal mode number to include
+nn_high = 1                             # Largest toroidal mode number to include
+delta_mlow = 8                          # Expands lower bound of Fourier harmonics
+delta_mhigh = 8                         # Expands upper bound of Fourier harmonics
+mtheta_vacuum = 960                     # Number of points used in splines over poloidal angle at the plasma-vacuum interface
+eulerlagrange_tolerance = 1e-9          # Relative tolerance for ODE integration of Euler-Lagrange equations
+singfac_min = 1e-4                      # Fractional distance from rational q at which ideal jump enforced
+integrator = "forward"                  # Integration formalism: "forward" (dense ξ), "riccati" (chunked propagator BVP, singular-surface Δ′), or "galerkin" (outer-region singular Galerkin Δ′)
+truncate_at_rational_offset = false     # Truncate integration at (last_rational_q + rational_offset_fraction)/n; false for analytic equilibria
+
+# Large-aspect-ratio analytic equilibrium parameters (eq_type = "lar"); see LargeAspectRatioConfig in src/Equilibrium.
+[LAR_INPUT]
+r0 = 10.0                               # Major radius of the plasma [m]
+a = 1.0                                 # Minor radius of the plasma [m]
+beta0 = 1e-3                            # Beta on axis (normalized pressure)
+q0 = 2.4                                # Safety factor on axis
+p_pres = 2.0                            # Exponent of the pressure profile p00·(1 − (r/a)²)^p_pres
+p_sig = 2.0                             # Exponent of the current-related profile
+sigma_type = "hollow"                   # Current profile: "default", "wesson" or "hollow" (off-axis peak, reversed shear)
+hollow = 6.0                            # Off-axis current enhancement of the "hollow" profile
+mtheta = 128                            # Number of grid points in the poloidal direction
+ma = 128                                # Number of grid points in the radial direction

regression-harness/cases/lar_reversed_shear_n1.toml

@@ -0,0 +1,201 @@
+# Regression case: large-aspect-ratio reversed-shear n=1 ideal stability.
+# q falls from 2.40 on the axis to 1.82 and rises to 2.84 at the edge, so q = 2 is crossed twice;
+# the two crossings must stay separate singular surfaces with opposite-sign q′. Tracks the
+# free-boundary energies, the singular surfaces and the integration for that fixture. Each
+# [quantities.*] block names an HDF5 path, how to extract it, and its noise floor.
+[case]
+name = "lar_reversed_shear_n1"
+description = "Large-aspect-ratio reversed-shear equilibrium (hollow current), n=1, ideal stability"
+example_dir = "examples/LAR_reversed_shear_example"
+
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
+[quantities.et_real]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "total energy Re(et[1])"
+noise_threshold = 1e-10
+order = 10
+
+[quantities.et_imag]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "imag_first"
+label = "total energy Im(et[1])"
+noise_threshold = 1e-10
+order = 11
+
+[quantities.ep_real]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "plasma energy Re(ep[1])"
+noise_threshold = 1e-10
+order = 12
+
+[quantities.ev_real]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "vacuum energy Re(ev[1])"
+noise_threshold = 1e-10
+order = 13
+
+[quantities.vacuum_eigenvalue]
+h5path = "ForceFreeStates/FreeBoundaryStability/vacuum_eigenvalue"
+type = "real_scalar"
+extract = "value"
+label = "vacuum matrix min eigenvalue"
+noise_threshold = 1e-10
+order = 14
+
+# Energies — full eigenvalue arrays
+[quantities.ep_all]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies"
+type = "complex_vector"
+extract = "all_complex"
+label = "plasma energy (all)"
+noise_threshold = 1e-10
+order = 20
+
+[quantities.ev_all]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies"
+type = "complex_vector"
+extract = "all_complex"
+label = "vacuum energy (all)"
+noise_threshold = 1e-10
+order = 21
+
+[quantities.et_all]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "all_complex"
+label = "total energy (all)"
+noise_threshold = 1e-10
+order = 22
+
+# Integration
+[quantities.nstep]
+h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep"
+type = "int_scalar"
+extract = "value"
+label = "ODE steps (saved)"
+noise_threshold = 0
+order = 30
+
+[quantities.nstep_total]
+h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep_total"
+type = "int_scalar"
+extract = "value"
+label = "ODE steps (total)"
+noise_threshold = 0
+order = 31
+
+# Equilibrium
+[quantities.q0]
+h5path = "Equilibrium/q_axis"
+type = "real_scalar"
+extract = "value"
+label = "q0"
+noise_threshold = 1e-12
+order = 40
+
+[quantities.q95]
+h5path = "Equilibrium/q_95"
+type = "real_scalar"
+extract = "value"
+label = "q95"
+noise_threshold = 1e-12
+order = 41
+
+[quantities.betat]
+h5path = "Equilibrium/beta_t"
+type = "real_scalar"
+extract = "value"
+label = "beta_t"
+noise_threshold = 1e-12
+order = 42
+
+[quantities.betan]
+h5path = "Equilibrium/beta_N"
+type = "real_scalar"
+extract = "value"
+label = "beta_n"
+noise_threshold = 1e-12
+order = 43
+
+# Singular surfaces
+[quantities.msing]
+h5path = "SingularSurfaces/rational_count"
+type = "int_scalar"
+extract = "value"
+label = "# singular surfaces"
+noise_threshold = 0
+order = 50
+
+[quantities.sing_psi]
+h5path = "SingularSurfaces/rational_psi"
+type = "real_vector"
+extract = "all_real"
+label = "singular psi locations"
+noise_threshold = 1e-8
+order = 51
+
+[quantities.sing_q]
+h5path = "SingularSurfaces/rational_q"
+type = "real_vector"
+extract = "all_real"
+label = "singular q values"
+noise_threshold = 1e-8
+order = 52
+
+[quantities.sing_dqdpsi]
+h5path = "SingularSurfaces/dqdpsi"
+type = "real_vector"
+extract = "all_real"
+label = "singular dq/dpsi (opposite signs on the two q = 2 crossings)"
+noise_threshold = 1e-8
+order = 53
+
+# Mode numbers
+[quantities.mpert]
+h5path = "Info/mpert"
+type = "int_scalar"
+extract = "value"
+label = "mpert"
+noise_threshold = 0
+order = 60
+
+[quantities.npert]
+h5path = "Info/npert"
+type = "int_scalar"
+extract = "value"
+label = "npert"
+noise_threshold = 0
+order = 61
+
+# Profiles (checksums for fast comparison of large arrays)
+[quantities.q_profile]
+h5path = "Equilibrium/Profiles/q"
+type = "real_vector"
+extract = "checksum"
+label = "q profile (checksum)"
+noise_threshold = 0
+order = 70
+
+[quantities.pressure_profile]
+h5path = "Equilibrium/Profiles/mu0p"
+type = "real_vector"
+extract = "checksum"
+label = "pressure profile (checksum)"
+noise_threshold = 0
+order = 71
+
+# Runtime (special: not from H5)
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999

src/Equilibrium/Grid/GridRefinement.jl

@@ -469,21 +469,34 @@ surface (no knot-on-surface jump) and *stable across* grid refinement — a narr
 wide neighbors would be consistent but noisy. Non-bracket knots within `w` of the center, or
 within `min_spacing` of either new bracket knot, are dropped; endpoints always win; bracket knots
 outside the open domain are dropped, and knots closer than `collapse_atol` are collapsed to keep
-the grid strictly increasing.
+the grid strictly increasing. Each half-width is capped at half the distance to the neighboring
+centers, so adjacent brackets never interleave and every center keeps its own interval, and the
+stencil part of it also at half the distance to the nearer domain edge, so a near-tangent surface
+(q′ → 0, a huge stencil) cannot clear the grid. A cap that cuts a bracket below its stencil floor
+warns, since Δ′ there is sampled on a narrower stencil than intended.
 """
 function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}, min_half_widths::Vector{Float64}, min_spacing::Float64; collapse_atol::Float64=1e-7)
     isempty(centers) && return copy(grid)
     length(centers) == length(min_half_widths) || error("bracket_mandatory_nodes: centers and min_half_widths length mismatch")
     lo, hi = grid[1], grid[end]
-    order = sortperm(centers)
-    kept = Tuple{Float64,Bool}[(g, false) for g in grid]  # (value, is_bracket_knot)
-    last_c = -Inf
-    for idx in order
-        c, hw = centers[idx], min_half_widths[idx]
+    # In-domain centers in increasing order; centers closer than collapse_atol share one bracket.
+    cs, hws = Float64[], Float64[]
+    for idx in sortperm(centers)
+        c = centers[idx]
         (lo < c < hi) || continue
-        c - last_c < collapse_atol && continue  # duplicate rational (same q via several (m,n))
+        (isempty(cs) || c - cs[end] >= collapse_atol) || continue
+        push!(cs, c)
+        push!(hws, min_half_widths[idx])
+    end
+    kept = Tuple{Float64,Bool}[(g, false) for g in grid]  # (value, is_bracket_knot)
+    capped = Tuple{Float64,Float64,Float64}[]  # (center, capped half-width, stencil floor)
+    for (i, c) in enumerate(cs)
         k = clamp(searchsortedlast(grid, c), 1, length(grid) - 1)
-        w = max(hw, 0.5 * (grid[k+1] - grid[k]))  # blend to half the local spacing
+        cap_nb = 0.5 * min(i > 1 ? c - cs[i-1] : Inf, i < length(cs) ? cs[i+1] - c : Inf)
+        cap_edge = 0.5 * min(c - lo, hi - c)
+        # Stencil floor capped by neighbors and edges; blend to half the local spacing capped by neighbors only.
+        w = max(min(hws[i], cap_nb, cap_edge), min(0.5 * (grid[k+1] - grid[k]), cap_nb))
+        w < hws[i] && push!(capped, (c, w, hws[i]))
         left, right = c - w, c + w
         # Drop non-bracket, non-endpoint knots inside the zone or crowding a new bracket knot.
         filter!(kept) do t
@@ -492,8 +505,10 @@ function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}
         end
         left > lo && push!(kept, (left, true))
         right < hi && push!(kept, (right, true))
-        last_c = c
     end
+    isempty(capped) || @warn "$(length(capped)) rational surface(s) have grid bracket half-widths capped below their Δ′ stencil " *
+                             "(neighboring surface or domain edge too close, or a near-tangent surface with q′ ≈ 0); Δ′ there is less accurate. " *
+                             "(ψ, half-width, stencil): $(join(capped, ", "))"
     sort!(kept; by=first)
     merged = Float64[]
     for (v, _) in kept

src/Equilibrium/Sources/Analytic/LargeAspectRatio.jl

@@ -11,7 +11,8 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi
   - `q0`: The safety factor on axis.
   - `p_pres`: The exponent for the pressure profile, defined as `p00 * (1 - (r / a)^2)^p_pres`.
   - `p_sig`: The exponent that determines the shape of the current-related function profile.
-  - `sigma_type`: The type of sigma profile, can be "default" or "wesson". If "wesson", the sigma profile is defined as `sigma0 * (1 - (r / a)^2)^p_sig`.
+  - `sigma_type`: The type of sigma profile, can be "default", "wesson" or "hollow". If "wesson", the sigma profile is defined as `sigma0 * (1 - (r / a)^2)^p_sig`. If "hollow", it is `sigma0 * (1 + hollow * (r / a)^2) * (1 - (r / a)^2)^p_sig`, an off-axis current peak that makes q fall from `q0` to a minimum before rising (reversed shear) when `hollow` is large enough.
+  - `hollow`: The off-axis current enhancement of the "hollow" sigma profile (unused by the other profiles).
   - `mtheta`: The number of grid points in the poloidal direction.
   - `ma`: The number of grid points in the radial direction.
   - `zeroth`: If set to true, it neglects the Shafranov shift
@@ -26,6 +27,7 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi
     p_pres::Float64 = 2.0
     p_sig::Float64 = 1.0
     sigma_type::String = "default"
+    hollow::Float64 = 0.0
     mtheta::Int = 128
     ma::Int = 128
     zeroth::Bool = false
@@ -86,6 +88,14 @@ function lar_init_conditions(rmin::Float64, lar_input::LargeAspectRatioConfig)
     return r, y
 end
 
+# Current profile σ(x), x = r/a, with σ(0) = sigma0 for every `sigma_type`.
+function _lar_sigma(x::Float64, sigma0::Float64, lar_input::LargeAspectRatioConfig)
+    p_sig = lar_input.p_sig
+    lar_input.sigma_type == "wesson" && return sigma0 * (1 - x^2)^p_sig
+    lar_input.sigma_type == "hollow" && return sigma0 * (1 + lar_input.hollow * x^2) * (1 - x^2)^p_sig
+    return sigma0 / (1 + x^(2 * p_sig))^(1 + 1 / p_sig)
+end
+
 """
     lar_der(dy, r, y, lar_input)
 
@@ -111,19 +121,13 @@ function lar_der(dy::Vector{Float64}, r::Float64, y::Vector{Float64}, lar_input:
     beta0 = lar_input.beta0
     p00 = beta0 / 2.0
     p_pres = lar_input.p_pres
-    p_sig = lar_input.p_sig
 
-    sigma_type = lar_input.sigma_type
     x = r / lar_a
     xfac = 1 - x^2
     pp = -2 * p_pres * p00 * x * xfac^(p_pres - 1)  # p'
     sigma0 = 2.0 / (q0 * lar_r0)
 
-    sigma = if sigma_type == "wesson"
-        sigma0 * xfac^p_sig
-    else
-        sigma0 / (1 + x^(2 * p_sig))^(1 + 1 / p_sig)
-    end
+    sigma = _lar_sigma(x, sigma0, lar_input)
 
     # LAR equations have singularities at r=0; ensure we never evaluate there
     @assert r > 0.0 "LAR equations require r > 0 (called with r=$r)"
@@ -159,7 +163,6 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
     lar_r0 = lar_input.r0
     q0 = lar_input.q0
     beta0 = lar_input.beta0
-    sigma_type = lar_input.sigma_type
 
     p00 = beta0 / 2.0
     lar_input.p_pres = max(lar_input.p_pres, 1.001)
@@ -194,9 +197,7 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
         x = r / lar_a
         xfac = 1 - x^2
         pval = p00 * xfac^lar_input.p_pres
-        sigma = (sigma_type == "wesson") ?
-                sigma0 * xfac^lar_input.p_sig :
-                sigma0 / (1 + x^(2 * lar_input.p_sig))^(1 + 1 / lar_input.p_sig)
+        sigma = _lar_sigma(x, sigma0, lar_input)
         q = r^2 * y[2] / (lar_r0 * y[1])
         temp[i, 1] = r
         temp[i, 2] = y[1]

src/ForceFreeStates/CoreTypes.jl

@@ -70,7 +70,7 @@ A mutable struct holding internal state variables for stability calculations.
     kinsing_scan_cond::Vector{Float64} = Float64[]
     kinsing_scan_threshold::Float64 = 0.0
     psilim::Float64 = 0.0
-    psilow::Float64 = 0.0   # lower integration bound; raised above the axis by sing_min! when qlow > qmin (RDCON gal)
+    psilow::Float64 = 0.0   # lower integration bound; the outermost crossing of q = qlow (sing_min!), else the first grid node
     qlim::Float64 = 0.0
     q1lim::Float64 = 0.0
     debug_settings::DebugSettings = DebugSettings()
@@ -139,7 +139,7 @@ gpec.toml.
   - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil
   - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (the drift-kinetic matrices of the KineticForces model, evaluated on the `[KineticForces]` profiles; needs `integrator = "forward"` and is deck-driven through `main`)
   - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components.
-  - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil
+  - `qlow::Float64` - Integration starts at the outermost crossing of q = qlow, excluding the q < qlow region inside it (a core, or a reversed-shear ring); no effect when q ≥ qlow everywhere
   - `dW_edge_scan_start::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [dW_edge_scan_start, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below).
   - `truncate_at_dW_peak::Bool` - When `true` and `dW_edge_scan_start < psilim`, the edge-dW scan's peak is adopted as the physical plasma edge, reproducing the ode_record_edge heuristic of Fortran STRIDE. `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, fixups recorded beyond it are discarded, `q1lim` and the rational-surface list are re-derived at the new edge, and on the Riccati path the Δ' chunks are cut back to it. The result is the same physics as setting `psihigh` to the peak by hand. Δ' still depends on where in the edge band the peak falls, so set this deliberately when you mean to measure against the peak-defined edge (e.g. edge-mode regimes); leave at `false` (default) for the full domain at `qhigh` / `psihigh` / `rational_offset_fraction`.
   - `write_outputs_to_HDF5::Bool` - Write results to HDF5 format