ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap

docs/development/references.md

@@ -68,6 +68,19 @@ The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calc
   - Published: Physics of Plasmas **24**, 032505 (2017)
   - Describes: Self-consistent coupling with neoclassical effects
 
+## InnerLayer Module (SLAYER)
+
+The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are based on:
+
+- **Burgess et al. (2026)**: "Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer"
+  - Location: `docs/resources/2026-Burgess-Tearing Stability Prediction Combining Toroidal Calculations With a Two-Fluid Slab Layer.pdf`
+  - Describes: The two-fluid slab layer, its four regimes, and the `S^(1/3)` outer-inner matching
+
+- **Fitzpatrick (2025)**: "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation"
+  - Location: `docs/resources/2025-Fitzpatrick-Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation.pdf`
+  - Published: Nuclear Fusion (2025), doi 10.1088/1741-4326/ae4fdd (open access: arXiv:2511.07666)
+  - Describes: Diffusive-resistive layer width (Eq. 100) and the resistive-layer-overlap truncation `0 < Ψ < 1 − ε_c` (Sect. 5.9)
+
 ## Resistive MHD Stability Analysis (Future Work)
 
 GPEC will eventually implement resistive MHD stability analysis based on:

docs/src/citations.md

@@ -131,6 +131,29 @@ Combines the toroidal outer-region calculation with the two-fluid slab layer mod
 
 ---
 
+> R. Fitzpatrick, "Response of a magnetically diverted tokamak plasma to a resonant magnetic perturbation,"
+> *Nuclear Fusion* (2025).
+> DOI: [10.1088/1741-4326/ae4fdd](https://doi.org/10.1088/1741-4326/ae4fdd) (open access: arXiv:2511.07666)
+
+Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain truncation](@ref).
+
+---
+
+> R. Fitzpatrick, "Investigation of tearing mode stability near ideal stability boundaries via asymptotic matching techniques,"
+> *Physics of Plasmas* **32**, 062509 (2025).
+> DOI: [10.1063/5.0273568](https://doi.org/10.1063/5.0273568) (arXiv:2503.24184)
+
+The critical island width `W_d` and its short/long-mean-free-path parallel conductivity blend, used for the critical-Δ offset in `src/InnerLayer/SLAYER/LayerParameters.jl`.
+
+---
+
+> J. W. Connor, C. J. Ham, R. J. Hastie, and Y. Q. Liu, "The role of thermal conduction in tearing mode theory,"
+> *Plasma Physics and Controlled Fusion* **57**, 065001 (2015). arXiv:1410.7240
+
+The toroidal (Eq. 59) and large-aspect-ratio (Eq. 61) critical-Δ offsets from finite parallel thermal conduction, building on H. Lütjens, J.-F. Luciani, and X. Garbet, *Physics of Plasmas* **8**, 4267 (2001).
+
+---
+
 ## ErrorFields Module
 
 > N. C. Logan et al., "Robustness of the tokamak error field correction tolerance scaling,"

docs/src/kinetic_forces.md

@@ -12,8 +12,9 @@ extension:
   - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets
     at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`,
     `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional
-    `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``), and
-    `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``). Each dataset
+    `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``),
+    `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``), and `zeff`
+    (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
     diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`)
@@ -24,9 +25,12 @@ extension:
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
     compatibility. Header rows are skipped.
 
-The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`,
+The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`/`zeff`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored
-here.
+here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file
+has one, otherwise from `n_i/n_e` with the `[SLAYER]` impurity charge `zimp`, otherwise
+from the scalar `[SLAYER] zeff`. Supply `zeff` whenever the impurity content is not a
+single species fixed by `n_i/n_e`.
 
 ## Multi-ion runs
 

docs/src/stability.md

@@ -72,8 +72,8 @@ integrator = "forward"
 
 `riccati_eulerlagrange_integration` (the default) is our implementation of the STRIDE
 approach [Glasser 2018b], built on the dual Riccati reformulation [Glasser 2018a].  It
-decomposes the radial domain into
-independent chunks, integrates each chunk's fundamental-matrix (FM) propagator in parallel
+decomposes the radial domain into independent chunks, integrates each chunk's
+fundamental-matrix (FM) propagator in parallel
 using `Threads.@threads`, then multiplies the propagators in order and applies each
 singular-surface crossing serially.  It is the only driver that produces the inter-surface
 ``\Delta'`` matrix.  Because chunk endpoints are all it stores, `u_store` is sparse and stays
@@ -253,6 +253,104 @@ The Galerkin integrator computes the same quantity in the same PEST-3 convention
 publishes it on the same path, so downstream consumers (SLAYER among them) never branch
 on which formalism ran.
 
+## Integration domain truncation
+
+The Euler–Lagrange integration runs from the axis out to `psilim`, which is set from the
+target safety factor `qlim`:
+
+- `qlim` starts at `min(qmax, qhigh)`, where `qmax` is q at the equilibrium's `psihigh`.
+- With `set_psilim_via_dmlim = true` (the default, for diverted equilibria where q → ∞ at the
+  separatrix), `qlim` moves to `(last rational + dmlim) / n`, so the integration ends a fixed
+  fraction of a rational spacing past the last surface instead of close to one. The offset is
+  then stepped back by `1/n` until it no longer exceeds the bound it started from, so `dmlim`
+  never extends the domain past `qmax`, `qhigh`, or the layer-overlap cap below. A bound that
+  falls on a rational, or less than `dmlim / n` above one, therefore excludes it: with
+  `dmlim = 0.2`, both `qhigh = 4.0` and `qhigh = 4.1` give `qlim = 3.2`.
+- With `set_psilim_via_dmlim = false` (limited or analytic equilibria with finite edge q),
+  `qhigh` and `psihigh` set the domain directly.
+
+### Resistive-layer overlap
+
+Near the separatrix the rational surfaces crowd together faster than their resistive layers
+narrow, so eventually adjacent layers overlap. Past that point no surface keeps a
+well-separated inner region and the matched-asymptotic treatment is not defined; Sect. 5.9 of
+Fitzpatrick (2025), listed under [InnerLayer Module](@ref) in the citations, restricts the
+response calculation to the region inside it.
+
+GPEC locates that point with a scan of the rational surfaces of each toroidal mode number in
+the run; for multi-`n` runs the innermost overlap point, over all `n`, is the one used. Each surface gets a layer solve from the kinetic profiles, and its
+diffusive-resistive width (Fitzpatrick 2025, Eq. 100) is converted to a normalized-flux width
+so it can be compared with the spacing between neighbouring surfaces. The cut is placed at the
+last surface before two adjacent layers first touch.
+
+Surfaces beyond the equilibrium grid are located on the logarithmic divergence of q at the
+separatrix, taken from the equilibrium's own poloidal flux rather than fitted. Each X-point on
+the last closed flux surface contributes one term,
+
+```math
+q(\psi) = q_K + \sum_i A_i \ln\frac{\psi_{X_i} - \psi_K}{\psi_{X_i} - \psi},
+\qquad A_i = \frac{F}{2\pi R_{X_i}\sqrt{-\det H_i}},
+```
+
+where ``\psi_{X_i}`` is the flux at the saddle, ``H_i`` the Hessian of the poloidal flux there,
+``F = R B_\phi``, and ``q_K`` the traced q on the outermost grid surface ``\psi_K``. One term
+describes a single null and two a double null, connected or not. The X-points are found by
+following each minimum of the poloidal field along the outermost surface to the saddle it leads
+to, so nulls the closed surface does not pass (in a divertor leg, beside a coil) are never
+candidates. The law is used only if it reproduces the traced shear over the outer knots, and a
+plasma with no X-point on its boundary (a limited plasma) is not extrapolated. Divertors with a
+second null beside the primary one or on the leg (snowflake, x-point target) are not modelled.
+
+The X-point flux need not coincide with the equilibrium's boundary: a reconstruction that places
+its boundary slightly inside the separatrix has finite q at ``\psi = 1``, and rational surfaces
+between that boundary and the X-point lie in the real plasma but outside the equilibrium GPEC
+holds. The scan places only the first of them, to judge whether the outermost surface inside is
+clean, and scores it only when the kinetic profiles extend that far. `psihigh` remains the upper
+limit throughout: the overlap point can only narrow the domain, and when no overlap can be
+established the domain keeps its requested limit.
+
+The scan takes its kinetic profiles from the `[SLAYER]` `profile_file`, through the same control
+and loader the SLAYER analysis uses, so both see the same profiles, the same HDF5 group and the
+same χ⊥(ψ) — the layer widths cannot disagree between them. A deck without a `[SLAYER]`
+`profile_file` does not get the scan. Because the file is resolved from the run's TOML, the scan
+only runs on the `gpec.toml` entry point: the programmatic `solve(prob, alg)` path has no inputs to
+resolve it from, so the cap cannot be applied there and setting the flag warns instead. Its result is always written to
+`ForceFreeStates/LayerOverlap/`, but it only constrains the domain when
+`psilim_from_layer_overlap = true`. It then caps `qlim` before the `dmlim` step, so `dmlim`
+still selects the final surface from inside the cap. A cap beyond `psihigh` has no effect.
+
+The flag is off by default, so every shipped deck is unchanged by it. With the shipped
+DIII-D-like profiles the overlap point is ψ ≈ 0.9985: outside the ideal and Riccati decks'
+`psihigh = 0.995`, where the cap could not bind even if enabled, but inside the SLAYER deck's
+0.9995, where enabling it would exclude the q = 7 surface. `examples/DIIID-like_truncation_example`
+raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and q = 8, and the resulting
+domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen `psihigh = 0.995` of the
+ideal deck. Whenever `psilim` lands past the overlap point, with the flag on or off, GPEC warns;
+on the shipped decks that includes the SLAYER deck.
+
+**Which perpendicular diffusivity.** The two formulas that take `chi_perp` want different
+quantities. The layer width goes through `P_perp = τ_R/τ_E` with `τ_E = r²/χ_E`, and χ_E is the
+one-fluid perpendicular *energy* diffusivity, the gradient-weighted combination of the electron and
+ion channels, `χ_E = (n_e χ_e ∇T_e + n_i χ_i ∇T_i)/(n_e ∇T_e + n_i ∇T_i)`. The critical-Δ island
+width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` describes electron-temperature flattening against an electron parallel
+conductivity, so its χ_⊥ is the electron channel. GPEC currently feeds both from one profile, the
+kinetic file's `chi_e`, which in the shipped DIII-D-like data is the electron power-balance
+diffusivity: the right channel for the critical Δ, and a proxy for χ_E in the layer width. This is
+the same single-χ convention as Fitzpatrick's TJ code, which also assumes `T_e = T_i`.
+Separating the two is deferred: it needs an ion `chi_i` in the kinetic-profile schema, and the
+power-balance ion diffusivity available upstream is not yet usable (it comes back negative across
+most of the radius on the cases checked). Since both widths scale as `χ^{1/4}`, a factor-2 error in
+χ moves them by about 19%, but that can be enough to move the overlap cut by a surface.
+
+| Dataset | Contents |
+|---|---|
+| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius |
+| `delta_s_abs`, `width_delta_s` | Riccati layer thickness, in metres and as a Δψ width |
+| `width_visco`, `width_dr` | visco-resistive and diffusive-resistive (Eq. 100) Δψ widths; `width_dr` sets the cut |
+| `extrapolated` | 1 where the surface lies beyond the equilibrium grid |
+| `psilim_overlap`, `first_overlap_index` | the cut, and the first surface whose layer overlaps its neighbour |
+| `applied` | 1 when the cap was passed to the truncation (`dmlim` / `qhigh` may still cut deeper) |
+
 ## Configuration reference
 
 All `ForceFreeStates` options are set in the `[ForceFreeStates]` section of `gpec.toml`.

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -72,7 +72,8 @@ scan_mode               = "amr"      # Q-plane scan strategy: "amr" (adaptive re
 coupling_mode           = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant)
 dc_type                 = "none"     # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal"
 mu_i                    = 2.0        # Ion mass in proton-mass units (2.0 = deuterium)
-zeff                    = 1.0        # Effective charge
+zeff                    = 1.0        # fallback only; Z_eff(ψ) comes from the kinetic file (its zeff profile, else n_i/n_e with zimp)
+zimp                    = 6.0        # Impurity charge closing quasi-neutrality when Z_eff(ψ) is derived from n_i/n_e (6.0 = carbon)
 chi_perp                = 1.0       # fallback only; the kinetic file supplies χ⊥(ψ)
 chi_tor                 = 1.0       # fallback only; the kinetic file supplies χ_φ(ψ)
 pole_threshold_adaptive = true       # SLAYER |Δ| spans many orders; adapt to 10·median(|Δ|)

examples/DIIID-like_truncation_example/gpec.toml

@@ -0,0 +1,50 @@
+# DIII-D-like H-mode equilibrium — n=1 truncation-only case for the resistive-layer overlap cap.
+# Reuses the equilibrium and kinetic profiles of the sibling DIIID-like_ideal_example (referenced
+# by relative path; not duplicated). psihigh is pushed well past the layer-overlap point (psi ~ 0.9985,
+# where adjacent resistive layers first touch), so psilim_from_layer_overlap binds: the cap pulls qlim
+# to 7 and dmlim steps back to 6.2, excluding q=7 and q=8. Ideal free-boundary stability only, with no local
+# stability, forcing or perturbed equilibrium: the energies are kept so a regression case can pin
+# what the capped domain does to the physics, not just to the bookkeeping.
+
+[Equilibrium]
+eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
+eq_type = "efit"               # Type of the input 2D equilibrium file
+jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
+grid_type = "auto"             # Radial grid packing type ("auto" = two-pass measured-curvature grid)
+psilow = 1e-4                  # Lower limit of normalized poloidal flux
+psihigh = 0.9999               # Upper limit of normalized poloidal flux (past the psi ~ 0.9985 layer-overlap point, so the cap binds)
+mpsi = 0                       # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
+psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
+mtheta = 256                   # Number of poloidal grid points
+newq0 = 0                      # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation)
+etol = 1e-10                   # Error tolerance for equilibrium solver
+force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
+
+[ForceFreeStates]
+local_stability_flag = false   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
+vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
+force_termination = true       # Terminate after force-free states (skip perturbed equilibrium calculations)
+
+psiedge = 1.0                  # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
+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
+
+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
+
+# 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)
+set_psilim_via_dmlim      = true   # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix)
+dmlim                     = 0.2    # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true)
+psilim_from_layer_overlap = true   # Cap psilim where adjacent resistive layers overlap (the scan is recorded either way)
+
+[SLAYER]
+# The growth-rate analysis is not run here; the section is present because the layer-overlap scan
+# takes its kinetic profiles, and its chi_perp(psi), from the SLAYER profile_file.
+enabled = false                # Run the SLAYER tearing-mode analysis
+profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled