Feature: add SLAYER and GGJ tearing growth rates by d-burg · Pull Request #238 · OpenFUSIONToolkit/GPEC

d-burg · 2026-05-18T18:47:22Z

Summary

Adds tearing-mode growth-rate analysis as a new Tearing umbrella module under src/Tearing/, with three layers:

InnerLayer — inner-layer Δ(Q) physics: GGJ (Pletzer–Dewar) and SLAYER (Fitzpatrick Riccati) models, with switchable Spitzer/neoclassical resistivity.
Dispersion — physics-agnostic root finder: SurfaceCoupling, MultiSurfaceCoupling, CoupledFull (2m×2m det(D′−D(γ))), and CoupledFortranMatch residuals; AMR Q-plane scan + triangulation-based growth-rate extraction; spurious-root detection (geom + γ-gap + polyline concavity)
Runner — TOML-driven orchestration (Tearing.Runner.run_slayer), kinetic-profile loading, HDF5 output.

Also includes the supporting work that landed alongside:

Full 2m×2m D′ matrix exposed from ForceFreeStates via delta_prime_raw / pest3_decompose
New Riccati.jl solver in ForceFreeStates for the ideal plasma response (with ~30–40% perf work)
Utilities: NeoclassicalResistivity, KineticProfiles, PhysicalConstants
Regression cases: solovev_slayer_n1, tj_epsilon_pole
LAR β-scan / ε-scan and TJ_epsilon_pole examples

73 commits, ~+13.9k / −0.75k LOC across 92 files.

Relationship to `perf/riccati`

The inner-layer growth-rate solvers here consume the outer-region Δ′ produced by the Riccati-based ideal-MHD solver on perf/riccati. This branch is designed to subsume perf/riccati once that work merges: it currently shares a common ancestor (c6c845ff) and is 39 commits ahead but 11 commits behind perf/riccati. Plan is to merge the final state of perf/riccati in here once it's ready to merge and before this PR moves out of draft, so the final history includes the outer-region Δ′ improvements that the dispersion residuals depend on.

Current status — WIP / draft

Core stack (InnerLayer + Dispersion + Runner) and regression cases pass on Solovev and TJ-like analytical equilibria.
Latest commit (528062f8) changes chooser_overrides from discard-on-(geom ∧ gap) to warn-and-keep, which fixed 7/8 mis-chosen roots in the SPARC β-scan; not yet validated against q95-scan, IBS_AT-scan, or DIII-D 147131 benchmarks.
Filtered-root HDF5 subgroup now records only legacy-rejected roots; geom/gap-warned roots live in valid_roots with flags.
Pending: integrate the soon-to-be finalized state of perf/riccati.

…PR 1/9) First step in porting the Fortran SLAYER (Park 2023) inner-layer model into julia_GPEC. Adds the per-surface parameter object and the dimensional-to-normalized constructor that Fortran's `params.f` provides, restricted to the Fitzpatrick `riccati_f` formulation actually used by the SLAYER dispersion solver. The legacy `pr`, `pe`, and ρ_s-based `ds` parameters are intentionally absent — they entered only the unported `riccati()` / `riccati_del_s()` paths. The complex growth rate `Q` is not stored on the struct and will be passed directly to `solve_inner` in PR 2. Highlights: - `SLAYERParameters` struct (immutable, @kwdef) carrying tau, lu, c_beta, D_norm, P_perp/P_tor, Q_e/Q_i/iota_e, conversion factors (tauk, tau_r, delta_n), geometric auxiliaries, and the dc_tmp / dc_type critical-Δ offset. - `slayer_parameters(; ...)` builder ports params.f including the Spitzer-Härm conductivity, Cole Q-normalization, Fitzpatrick d_β / D_norm, and the four dc_type branches (:none, :lar, :rfitzp, :toroidal) with their Wd iteration. - `r_based_shear(rs, q, dq/dψ, da/dψ)` helper performing the Fitzpatrick (minor-radius) shear conversion that layerinputs.f does inline before calling params() — needed because STRIDE shear is ψ-based but params.f formulas all assume r-based. - New `Utilities/PhysicalConstants` submodule with SI constants matching sglobal.f exactly so cross-code numerics line up. - 45 unit tests in `runtests_slayer_params.jl`, including a synthetic Solovev-like analytic check on the shear conversion. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…PR 2/9) Ports the Fortran SLAYER `riccati_f`/`w_der_f`/`jac_f` from delta.f:323-494 into Julia. The complex normalized growth rate `Q = ω + iγ` is passed directly to `solve_inner` as agreed; all other inputs come from `SLAYERParameters` (PR 1). The standard `riccati()` and `riccati_del_s()` Fortran variants and the `parflow_flag`/ `PeOhmOnly_flag=.FALSE.` branches are intentionally not ported. Implementation: - `_riccati_f_coeffs` evaluates fA, fA', fB, fC at point p with shared denominator caching (mirrors w_der_f). - `_riccati_f_rhs!` (in-place) and `_riccati_f_jac!` (analytic 1×1) feed an `ODEFunction(jac=...)` for stiff Rosenbrock integration. - `_riccati_f_initial` selects between the large-D_norm and small-D_norm asymptotic boundary-condition branches based on the same `D_norm² ≷ iota_e·P_perp/P_tor^(2/3)` test as Fortran, with the `MAX(my_p, 6.0)` floor preserved. - `solve_inner(::SLAYERModel{:fitzpatrick}, p, Q)` integrates inward from p_start to pmin (default 1e-6) using Rodas5P(autodiff=false) with reltol=abstol=1e-10 to match Fortran LSODE defaults, then extracts Δ = π / W'(pmin) via a single RHS evaluation. Returns SVector(Δ, 0) so SLAYER and GGJ are interchangeable through the shared `InnerLayerModel` interface. 17 unit tests in `runtests_slayer_riccati.jl`: interface compliance, both BC branches reachable, p_floor enforcement, Q-sweep smoothness, tolerance self-consistency, and pmin deepening stability. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…k (PR 3/9) Introduces a new top-level `Dispersion` module that combines the outer-region Δ' from PerturbedEquilibrium with the inner-layer Δ(Q) from any `InnerLayerModel` to build the per-surface tearing-dispersion residual r(Q) = dp_diag − scale · Δ_inner(Q) − Δ_crit `SurfaceCoupling` packages (model, params, dp_diag, dc, scale) and is itself Q-callable, so it can be broadcast over a 2D complex-Q grid by the brute-force/AMR scans in PRs 5-6. All root-finding will be done downstream by contour intersection on those scans (find_growthrates port, PR 5); this module deliberately contains no local Newton/secant iteration. The `surface_coupling` constructor dispatches on the inner-layer model type to auto-fill `scale`: lu^(1/3) for SLAYER (Fortran de-normalization at growthrates.f:217-218,260), 1 for GGJ (rescale_delta is applied internally inside solve_inner). A generic fallback with an explicit `scale` kwarg lets new inner-layer models plug in without touching this file. 20 unit tests in runtests_dispersion_residual.jl: synthetic LinearTestModel exercising the residual arithmetic against the closed form, SLAYER self-consistency (build dp_diag from Δ(Q_pin) and verify the residual is exactly zero at Q_pin), GGJ ↔ SLAYER constructor interchangeability through the abstract InnerLayerModel interface, and broadcast-compatibility on a 2D Q grid. Co-Authored-By: Claude Opus 4.6 <[email protected]>

These files were accidentally included in the previous commit (PR 3/9) despite being deleted from the filesystem before staging. The design decision is that all dispersion root-finding flows through 2D contour intersection on Q-plane scans (PR 5 find_growthrates port); local Newton/secant iteration is intentionally not provided. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…ual (PR 4/9) Adds the coupled multi-surface tearing dispersion residual det(M(Q)), mirroring the Fortran SLAYER `dispersion_det` (growthrates.f:190-279) that runs when `coupling_flag = .TRUE.`. `MultiSurfaceCoupling` packages a vector of per-surface SurfaceCoupling objects (PR 3), the full outer-region Δ' matrix, the reference surface whose tauk defines the Q normalization, and the truncation `msing_max`. It is itself Q-callable so the same brute-force/AMR scan infrastructure (PRs 5-6) can evaluate either the per-surface residual or the coupled determinant by broadcasting over a complex-Q grid. At each evaluation, for k = 1 .. msing_max the inner-layer Δ is computed at a per-surface-rescaled Q_k = Q · (tauk_ref/tauk_k) (growthrates.f:246), then subtracted (with the dc offset) from the diagonal of an upper-left msing_max × msing_max submatrix of dp_matrix. Off-diagonal Δ' couplings pass through unchanged. `SurfaceCoupling` gains a `tauk::Float64` field to carry the per-surface time normalization. The SLAYER constructor populates it from `params.tauk`; GGJ defaults to 1.0 (no inter-surface rescaling); the generic fallback takes it as a kwarg. `msing_max` defaults to `min(3, length(surfaces))` because Δ' off-diagonal couplings beyond the third surface tend to be erratic in practice. Callers can override (up to length(surfaces)) when more surfaces are known to be well-behaved. 42 unit tests in runtests_dispersion_coupled.jl: constructor validation (including 4-surface default cap and explicit override), diagonal Δ' factorization, single-surface root preservation, off-diagonal-coupling closed-form det shift, msing_max truncation with upper-left-submatrix semantics, per-surface Q rescaling verified against analytic det = Q²/2 with mismatched tauks, SLAYER self-consistency (constructed singular M(Q_pin) from known Δs at Q_pin), GGJ-surface flow-through, and 2D-grid broadcast compatibility. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…es port (PR 5/9) Adds the user-facing 2D Q-plane scanner and the contour-intersection growth-rate extractor — together these give the first end-to-end path from a (model, params, Δ') triple to a physical (ω_Hz, γ_Hz) tearing eigenvalue. `brute_force_scan(f, Q_re_range, Q_im_range; nre, nim, threaded=true)` evaluates any Q-callable residual (SurfaceCoupling, MultiSurfaceCoupling, or a plain function) on a regular nre × nim grid. Resolution and box are entirely user-controlled. Threaded across the imaginary axis by default; pass `threaded=false` for deterministic single-threaded evaluation when the residual is non-thread-safe. `find_growth_rates(scan, tauk; ...)` is a Julia port of CTM-processing/shared/find_growthrates.py for the regular-grid case (PR 6 will add the scattered/AMR triangulation path): - extracts Re(Δ)=re_target and Im(Δ)=im_target polylines via Contour.jl; - finds all segment-segment intersections (hand-rolled parametric crossing test on the complex plane); - classifies each intersection as a pole if max(|Re(Δ)|) along the nearest Im=0 contour exceeds `pole_threshold` (Re values are bilinear-interpolated from the grid onto contour vertices); - applies the "+γ step inside Re=0 contour loop" filter for spurious upper-branch roots — only when the nearest Re=0 contour is approximately closed (closure_gap < 10% of contour extent); - reports the highest-γ surviving root in physical Hz units via the user-supplied tauk. `GrowthRateResult` exposes Q_root, omega_Hz, gamma_Hz, plus all valid roots, poles, filtered roots, and the extracted polylines for diagnostics / plotting. 33 unit tests in runtests_dispersion_scan.jl: scan layout and threaded-vs-non-threaded agreement, single-root recovery to grid-resolution precision, multi-root selection of highest-γ, pole detection on Δ = (Q−Q_r)/(Q−Q_p) with explicit pole_threshold verification, tauk normalization to physical Hz, empty-result handling, and end-to-end API checks with both SurfaceCoupling and MultiSurfaceCoupling. Co-Authored-By: Claude Opus 4.6 <[email protected]>

… extraction (PR 6/9) Ports the Fortran SLAYER `dispersion_AMR_v2` (growthrates.f:367-700) into Julia and adds a scattered-data path to `find_growth_rates` so AMR output can feed directly into the same root-extraction pipeline as the brute-force grid scan. AMR scan: - `amr_scan(f, Q_re_range, Q_im_range; nre0, nim0, passes)` builds an axis-aligned quadtree of AMRCells. Each refinement pass subdivides any cell whose 4 corner residuals straddle zero in Re(Δ) or Im(Δ) into 4 quadrant children, evaluating 5 new midpoint Δ values. - All f(Q) evaluations deduplicated through a Dict{ComplexF64, ComplexF64} hash cache, replacing the Fortran's hand-rolled prime-multiplier hash. Adjacent cells thus share a single evaluation per corner, and refined neighbors share a single evaluation per edge midpoint. - Output `AMRResult` carries both the cell list (for visualization/diagnostics) and the flat Q/Δ vectors of all unique evaluations (for triangulation-based extraction). AMR-aware growth-rate extraction: - `find_growth_rates(::AMRResult, tauk; …)` triangulates the scattered (Q, Δ) evaluation points via DelaunayTriangulation.jl (matches the matplotlib.tri.Triangulation that find_growthrates.py uses) and marches triangles to extract Re=0 and Im=0 contour segments. - Marching step computes each segment endpoint along with the complementary field value (Re at Im=0 segment endpoints and vice versa) via linear interpolation along the same edge parameter t, so the pole-classification lookup gets filled for free with no separate interpolation pass. - Segments chained into polylines via bit-exact endpoint-matching Dict — adjacent triangles compute identical crossings on shared edges because endpoint values come from the shared hash cache. - Triangulating the scattered points resolves the hanging-nodes issue that would have plagued a per-cell marching-squares approach at refinement-level boundaries (the mismatched edge midpoints become first-class triangulation vertices instead of being ignored by the coarser neighbor). Refactor: grid (PR 5) and AMR (this PR) paths of `find_growth_rates` now share a single `_run_analysis(re_paths, im_paths, im_re_vals, tauk; …)` helper that handles intersection finding, pole classification, outside-Re filter, and physical-Hz conversion. Adds DelaunayTriangulation.jl 1.6.6 (pure Julia, BSD, JuliaGeometry org) to deps + compat. 30 unit tests in runtests_dispersion_amr.jl: hash-cache correctness (9 unique evaluations for a 2×2 coarse grid with no refinement), refinement concentration, argument validation, max_cells safety cap, single-root recovery, higher-γ root selection on a 2-root case, pole detection, tauk normalization to physical Hz, AMR-vs-brute-force consistency, and end-to-end API checks with SurfaceCoupling and MultiSurfaceCoupling. Co-Authored-By: Claude Opus 4.6 <[email protected]>

Adds the two building blocks needed to construct SLAYER inputs from a running julia_GPEC pipeline without the Fortran's STRIDE-NetCDF round-trip: 1. `Utilities.KineticProfiles` — radial profiles of n_e, T_e, T_i, ω, ω_*e, ω_*i as cubic splines of normalized ψ ∈ [0,1]. Three constructors: keyword args with matched-length vectors, a TOML section dict, and an HDF5 file + group path. `kp(ψ)` returns a NamedTuple of all six values. Placed in `Utilities/` so PENTRC and resistive-MHD modules can share it. 2. `SLAYER.build_slayer_inputs(equil, sings, profiles; …)` — ports Fortran `layerinputs.f` to read everything from in-memory structures instead of STRIDE NetCDF. Minor radius and da/dψ are pulled from `equil.rzphi_rsquared` at the outboard midplane (θ=0 by default), ψ-based shear is converted to Fitzpatrick r-based via `r_based_shear`, kinetic data is interpolated from the `KineticProfiles` at each `SingType.psifac`, and the first element of each surface's (m, n) mode-number vectors is used as the primary resonance. Scalars and callables-of-ψ are both accepted for χ⊥, χ∥, dr_val, and dgeo_val so simple cases stay concise and profile-varying cases are still expressible. 3. Helpers `surface_minor_radius(equil, ψ; θ=0.0)` and `surface_da_dpsi(equil, ψ)` (central FD with one-sided fallback near boundaries) are exposed so callers can query geometry outside the full pipeline. 48 unit tests covering kwarg/TOML/HDF5 constructors, length validation, round-trip exactness at spline nodes, the Solovev-bundled example equilibrium for minor-radius monotonicity and FD accuracy, per-surface SLAYERParameters extraction (geometry + mode numbers + Q_e/Q_i sign convention), scalar-vs-callable χ with closed-form P_perp ∝ χ⊥ check, dc_type propagation, and empty-sings edge case. This PR sets up the wiring; PR 8 will connect it to the PerturbedEquilibrium workflow, add the TOML [SLAYER] section, write a `slayer/` HDF5 group, and add the regression-harness case. Co-Authored-By: Claude Opus 4.6 <[email protected]>

Adds a new top-level `SLAYERRunner` module (sibling to `Dispersion`) that ties together the building blocks from PRs 1-7 into the user-facing SLAYER tearing-mode analysis pipeline. Orchestration lives in its own module to keep `InnerLayer` and `Dispersion` as pure physics/math libraries — no equilibrium/HDF5/TOML concerns leak into them. Four files: - `Control.jl` --- `SLAYERControl` struct with every user-facing knob (inner-model selector, scan mode, coupling mode, physics knobs, scan grid, AMR parameters, growth-rate filter thresholds, profile source, HDF5 options). `slayer_control_from_toml(section)` parses a `[SLAYER]` section and its nested `[SLAYER.scan_grid]`, `[SLAYER.amr]`, and `[SLAYER.growth_rate_filter]` subsections into a flat control; unknown keys raise an error so typos are caught at parse time. `validate(ctrl)` enforces the allowed Symbol sets and positivity constraints. - `Result.jl` --- `SLAYERResult` carries per-surface parameters, the full Δ' matrix used, Q_root / omega_Hz / gamma_Hz vectors, the per-surface GrowthRateResult array (uncoupled) or single coupled GrowthRateResult, and optional stored scan data. - `Runner.jl` --- `run_slayer(equil, ffs_intr, control, toml_section; dir_path)` is the full pipeline: loads kinetic profiles (inline TOML or HDF5 file), calls `build_slayer_inputs` (PR 7) to construct per-surface SLAYERParameters, pulls the outer-region Δ' matrix from `ffs_intr.delta_prime_matrix` (or falls back to a diagonal from each SingType.delta_prime), dispatches on coupling_mode and scan_mode, and extracts growth rates via find_growth_rates. A secondary `run_slayer_from_inputs(params, dp_matrix, control)` entry skips the equilibrium-driven build — used by unit tests. - `HDF5Output.jl` --- `write_slayer_hdf5!(parent, result)` writes a `slayer/` subgroup with `settings/`, `per_surface/` (struct-of- arrays for every SLAYERParameters field plus the Δ' matrix), `roots/`, `diagnostics/` (valid_roots / poles / filtered_roots as ragged flat_real/flat_imag/offsets triples), and optionally `scan/` (brute-force Q/Δ grid or AMR Q/Δ vectors + cell count). Disabled results still emit `enabled = 0` so downstream readers can detect the no-op case. 61 unit tests: control defaults + validation (rejects bad symbols and out-of-range ints), TOML nested-subsection flattening with unknown- key detection, disabled no-op path, size-mismatch rejection, a coupled-mode synthetic with a constructed known root recovered to grid-resolution precision, and HDF5 round-trip checking groups + settings + per-surface arrays + ragged-encoding structure. Not in this PR (deferred to PR 9): main() integration reading a `[SLAYER]` section from gpec.toml and calling run_slayer at the end of compute_perturbed_equilibrium, plus a regression-harness case tracking omega_Hz / gamma_Hz. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…ion case (PR 9/9) Final integration step that ties the SLAYERRunner module (PR 8) into the top-level GPEC pipeline so a `[SLAYER]` section in any `gpec.toml` drives the analysis end-to-end and writes results to the existing output HDF5 file. main() (src/GeneralizedPerturbedEquilibrium.jl): - After the PerturbedEquilibrium step, look for a `[SLAYER]` section in the parsed TOML. If present, parse it via `slayer_control_from_toml`. If `enabled = true`, call `run_slayer(equil, intr, slayer_ctrl, inputs["SLAYER"]; dir_path=intr.dir_path)` and append a `slayer/` group to the same HDF5 file the PE step writes (or the ForceFreeStates file if PE didn't run). The result is also returned in the top-level NamedTuple as `slayer=...` for script callers. examples/Solovev_ideal_example/gpec.toml: - Added an active `[SLAYER]` section (coupled mode, brute-force, 20x20 grid, synthetic deuterium kinetic profiles) so the bundled example demonstrates SLAYER end-to-end and the regression harness has something to track. SLAYER takes ~5 s on top of the existing Solovev pipeline. regression-harness/cases/solovev_slayer_n1.toml: - New regression case tracking 17 SLAYER outputs: per-surface layer parameters (ising, m, n, rs, sval_r, lu, c_beta, D_norm, P_perp, tauk, iota_e), the coupled-mode tearing eigenvalue (Q_root real/imag, omega_Hz, gamma_Hz), and the `enabled` flag. Pointed at the same example_dir as solovev_n1 so the harness benefits from output file sharing. Verification: - Solovev example writes slayer/ group with all expected sub-groups and arrays. - Coupled eigenvalue Q_root = 4e-4 + 0.112i (omega_Hz=1.9, gamma_Hz=529) on the synthetic deuterium profiles. - solovev_n1 regression still extracts its 22 ideal-stability quantities cleanly (SLAYER doesn't perturb upstream results). - solovev_slayer_n1 regression extracts all 17 SLAYER quantities. - Unit-test suite (PRs 1-8) all green. This completes the SLAYER port. The final "all SLAYER PRs" suite covers 292 unit tests + 2 regression cases. Co-Authored-By: Claude Opus 4.6 <[email protected]>

Consolidates the three top-level modules related to tearing-mode analysis (InnerLayer, Dispersion, SLAYERRunner) under a single `src/Tearing/` directory with a new umbrella module file. Pure reorganization — no behavior change. Layout: src/Tearing/ ├── Tearing.jl (new umbrella) ├── InnerLayer/ (was src/InnerLayer/) │ ├── GGJ/ │ └── SLAYER/ ├── Dispersion/ (was src/Dispersion/) └── Runner/ (was src/SLAYERRunner/) └── Runner.jl (was SLAYERRunner.jl) Module renames: - SLAYERRunner → Runner (inside Tearing) - The inner Runner.jl functions file renamed to run_slayer.jl to free the Runner.jl name for the outer module file. The umbrella rebinds `Utilities` at the Tearing level via `using ..Utilities`, so every submodule's existing relative imports (`using ..Utilities`) keep working without modification — the dot- counts don't change because Utilities is now a sibling of the submodules' grandparent view. Top-level `GeneralizedPerturbedEquilibrium.jl` now has a single `include("Tearing/Tearing.jl")` replacing three separate includes. Backward-compat top-level aliases `InnerLayer`, `Dispersion`, and `Runner` are preserved so existing test files and scripts using `GeneralizedPerturbedEquilibrium.InnerLayer` etc. continue to work. The canonical nested path (`Tearing.InnerLayer`, etc.) is also available. `main()` switched from `SLAYERRunner.*` to `Runner.*`. All 292 unit tests pass after the move. Solovev example SLAYER run unchanged at 5.7 s. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…_ggj_inputs Adds the per-singular-surface Glasser-Greene-Johnson geometric coefficients that GGJParameters needs, plus the builder function that turns (equil, sings, KineticProfiles) into Vector{GGJParameters} — symmetric to build_slayer_inputs. ForceFreeStates.ResistEval (new): - `ResistGeometry` struct holding E, F, G, H, K, M plus the two flux-surface averages ⟨B²/|∇ψ|²⟩, ⟨B²⟩ and the local p, dp/dψ, dV/dψ that downstream callers need to build τ_A / τ_R. - `resist_geometry(equil, psifac, q1; gamma=5/3)` ports the geometric portion of Fortran `rdcon/resist.f::resist_eval`. 6 theta-integrands per surface (the Mercier 5 plus ⟨|∇ψ|²/B²⟩), integrated via the same periodic cubic spline integrator `mercier_scan!` uses, then combined into the standard GGJ formulas: E = p1·v1/(q1·χ₁²)² · ⟨B²/|∇ψ|²⟩ · (2πF·q1·χ₁/⟨B²⟩ - dV²/dψ²) F = (p1·v1/(q1·χ₁²))² · (...) G = ⟨B²⟩ / (M·γ·p) H = same as Mercier H K = (q1·χ₁²/(p1·v1))² · ⟨B²⟩ / (M·⟨B²/|∇ψ|²⟩) M = ⟨B²/|∇ψ|²⟩ · (⟨|∇ψ|²/B²⟩ + (2πF/χ₁)²·(⟨1/B²⟩-1/⟨B²⟩)) - `resist_eval_all!(intr, equil)` populates `sing.restype` for every SingType in `intr.sing` (idempotent: skips already-populated). SingType gets a new `restype::Any` field (defaults `nothing`; typed `Any` to avoid a cross-file type reference). The main() workflow calls `resist_eval_all!(intr, equil)` after `sing_find!` and the qlow/qlim filter, so by the time downstream code runs every surviving surface has E, F, G, H, K, M available. HDF5 output extends the `singular/` group with 11 new datasets: E, F, G, H, K, M, avg_bsq, avg_bsq_over_dpsisq, p_local, p1_local, v1_local — all per-surface arrays. Tearing.InnerLayer.GGJ.build_ggj_inputs (new file): - `build_ggj_inputs(equil, sings, profiles::KineticProfiles; mu_i=2.0, zeff=1.0, v1_scale=1.0) -> Vector{GGJParameters}`. Symmetric to build_slayer_inputs. Geometric coefficients pass through unchanged from sing.restype; kinetic timescales are built from KineticProfiles using the SAME formulas SLAYER uses (Spitzer η from T_e/n_e/lnΛ; ρ = μ_i·m_p·n_e). τ_A and τ_R then come from the standard `rdcon/resist.f` definitions: τ_A = √(ρ·M·μ₀) / |2π·n·q'·χ₁/V'| τ_R = (⟨B²/|∇ψ|²⟩/⟨B²⟩) · μ₀/η - Deliberately does NOT mirror the Fortran rdcon/resist.f hardcoded `ne=1e14 cm⁻³, te=3 keV` PARAMETER defaults. GGJ and SLAYER both pull kinetic content from the same KineticProfiles, so the two can be compared on bit-identical plasma inputs. 61 unit tests in runtests_resist_eval.jl: finite/positive coefficient checks across multiple ψ, the D_I = E + F + H − ¼ cross-check against Mercier (matches to ~1e-4 relative), populator behaviour (including idempotency), build_ggj_inputs end-to-end with timescale and Lundquist sanity checks, error path when restype is unset, and a GGJ solve_inner invocation on the built parameters to confirm the pipeline actually runs. Total test count: 353 across all SLAYER + GGJ + Tearing files. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…ownstream Δ' matrix Defaults updated for SLAYER/GGJ downstream consumption: - etol 1e-7 → 1e-10 (equilibrium convergence) - eulerlagrange_tolerance 1e-7 → 1e-8 - singfac_min 0 → 1e-4 (required non-zero on the parallel path) - sing_order 2 → 6 (STRIDE convention for Δ') - use_parallel false → true (unlocks singular/delta_prime_matrix) - Add set_psilim_via_dmlim + dmlim controls in sing_lim! (Fortran sas_flag equivalent) for single-n truncation beyond the outermost rational surface Test fixes: runtests_slayer_params / runtests_slayer_inputs updated for the params.f sign convention Q_i = -tauk·ω*_i (both Q's share the same sign structure; earlier tests held the layerinputs.f Q_i sign flip which we deliberately do not mirror). Co-Authored-By: Claude Opus 4.6 <[email protected]>

… η in GGJ & SLAYER Adds a shared Spitzer/Sauter/Redl resistivity closure so GGJ and SLAYER can both consume the same neoclassical η formula: - src/Utilities/NeoclassicalResistivity.jl (new): SpitzerModel / SauterNeoModel / RedlNeoModel tag types, coulomb_log_e (NRL/Sauter/ Wesson forms), eta_spitzer (Sauter 1999 Eq. 18a), trapped_fraction (Lin-Liu & Miller 1995 full form) + trapped_fraction_eps fallback, nu_star_e (Sauter 1999 Eq. 18b), and eta_neoclassical dispatched on the model (F₃₃ via Sauter 1999 Eq. 13 or Redl 2021 Eq. 17). - src/ForceFreeStates/ResistEval.jl: ResistGeometry struct extended with avg_B, B_max, B_min, f_trap, R_major, eps_local. Populated inside the existing θ-loop at essentially zero cost (one extra integrand + running min/max over B and R). - src/Tearing/InnerLayer/GGJ/LayerInputs.jl: build_ggj_inputs grows `resistivity_model::NeoResistivityModel=SpitzerModel()` and `lnLambda_form::Symbol=:nrl` kwargs. Uses the shared closure; default Spitzer switches from Wesson 1.65e-9·lnΛ form to Sauter-18a (Zeff-aware, ~1% agreement at Zeff=1). - src/Tearing/InnerLayer/SLAYER/LayerParameters.jl + LayerInputs.jl: same `resistivity_model` kwarg, plus optional f_trap / nu_e_star / R_major_eff / lnLambda_form. Defaults to SpitzerModel() + :wesson so legacy SLAYER η is bit-identical. When a neoclassical model is selected, build_slayer_inputs pulls f_trap + R_major + eps_local from sing.restype if populated, and computes ν*_e via the shared utility. Validated on DIII-D 147131 @ 2300 ms (ideal example) vs OMFIT utils_fusion.py and OFT bootstrap.py F₃₃ formulas: max |reldiff| = 1.8e-16 across all 4 rational surfaces for lnΛ, ν*_e, η_Sp, η_Sauter, η_Redl, F₃₃(Sauter), F₃₃(Redl). Benchmark lives at CTM-processing/julia_vs_fortran/neoclassical_resistivity_benchmark/. In the DIII-D banana regime (q=2,3,4), η_Sauter/η_Sp ≈ 4–5× — the expected trapped-particle enhancement for H-mode tearing studies. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…_prime_raw + pest3_decompose The STRIDE-BVP Δ' computation already assembles a 2m×2m side-major matrix dp_raw in compute_delta_prime_matrix! (Riccati.jl:779, ordering [L_s1, R_s1, L_s2, R_s2, …]), then collapses it to the m×m PEST3 odd-parity Δ' projection via deltap[i,j] = dp_raw[2i,2j] − dp_raw[2i,2j-1] − dp_raw[2i-1,2j] + dp_raw[2i-1,2j-1] (the (L−R)(L−R)^T combination). The A' (even-parity interchange), B', Γ' (off-parity) blocks are thrown away. This commit retains the full 2m×2m matrix: - New ForceFreeStatesInternal.delta_prime_raw field (side-major, byte- compatible with Fortran rdcon/gal.f::gal_write_delta top 2msing×2msing block of delta_gw.dat; no ½ prefactor per Fortran convention). - Populated right before PEST3 collapse at Riccati.jl:819. - Persisted as singular/delta_prime_raw in gpec.h5. - New pest3_decompose(dp_raw) → (A, B, Γ, Δ) and dprime_outer_matrix helpers, matching Fortran rdcon/gal.f:1728-1743 recombination. Needed for the full det(D' − D(γ)) = 0 tearing+interchange eigenvalue problem in Phase C. Sanity-checked on DIII-D: pest3_decompose(dp_raw).Δ matches the existing m×m delta_prime_matrix to 4.6e-14. Cross-check vs Fortran delta_gw.dat shows pre-existing dpsi^α normalization gap (neither code writes the Hermitian form; it's applied at use-time). Benchmark artefacts at CTM-processing/julia_vs_fortran/ggj_coefficients_benchmark/ dprime_raw_crosscheck/. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…GGJ parity channel selection Replaces solve_inner's anonymous SVector{2,ComplexF64} return with a named struct InnerLayerResponse(tearing, interchange) to eliminate a latent parity-channel bug and self-document the inner-layer API. The bug: the old contract said "(Δ_odd, Δ_even)" but the word "odd"/"even" is used inconsistently across the literature — GWP 2016 labels parity by the symmetry of the flux W (odd-W = interchange, even-W = tearing), while Fortran rmatch/deltac.f labels by the velocity+temperature (odd-NΘ = tearing, even-NΘ = interchange). These give OPPOSITE parity names for the same physics channel. The GGJ Galerkin solver mirrored deltac.f's end-of-routine swap (Galerkin.jl:711-712), putting index 1 = interchange. The GGJ Shooting solver mirrored deltar.f, putting index 1 = interchange. SLAYER put its pressureless tearing Δ at index 1. Meanwhile Dispersion/Coupled.jl:96 and Dispersion/SurfaceCoupling.jl:46 hardcoded [1] — so for SLAYER surfaces they correctly picked the tearing channel, but for GGJ surfaces they silently picked the INTERCHANGE (Glasser-stabilization) channel instead of the tearing drive. Any GGJ multi-surface dispersion scan run prior to this commit was solving the wrong eigenvalue problem. Fix: - New InnerLayerResponse struct with physics-named tearing/interchange fields. - GGJ Galerkin: removed the deltac.f swap; isol=1 (W'(0)=0 → W even, sheet current, tearing) maps to .tearing; isol=2 (W(0)=0 → W odd, non-reconnecting) maps to .interchange. Per-solver parity derivation documented in BC comments. - GGJ Shooting: traced match/matrix.f::matrix_layer sign-symmetric vs sign-antisymmetric constraints to confirm deltar(1)=interchange, deltar(2)= tearing; remapped _delta_from_c0 output into named fields accordingly. - SLAYER: pressureless Fitzpatrick has no interchange channel → InnerLayerResponse(Δ, 0). - Dispersion/Coupled.jl + SurfaceCoupling.jl: replaced solve_inner(...)[1] with solve_inner(...).tearing at both call sites. - 6 test files updated: synthetic test models return InnerLayerResponse; real SLAYER/GGJ callers use .tearing. 200+ tests pass; 2 pre-existing slayer_riccati failures (D_norm threshold drift, unrelated to parity refactor) verified by git-stash bisection. Naming: chose tearing/interchange per user decision — more self-documenting than odd/even which depends on whose parity convention you're reading. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…matrix Companion to the m×m MultiSurfaceCoupling (tearing-only) that was shipped earlier in the perf/slayer-growthrates branch. CoupledFull generalizes to the full Pletzer-Dewar 1991 / GWP 2016 tearing+interchange eigenvalue problem needed to include Glasser stabilization in the GGJ model. Structure: - MultiSurfaceCouplingFull holds a 2m×2m D' matrix in parity-major ordering [[A' B'] [Γ' Δ']], a per-surface Vector{SurfaceCoupling}, reference-surface index, and msing_max truncation. Built via multi_surface_coupling_full(surfaces, dp_full; ref_idx, msing_max). - Evaluation mc(Q) subtracts a 2m×2m block-diagonal D(γ) with interchange-channel response on the upper-left m diagonal and tearing-channel response on the lower-right m diagonal. Each channel rescaled by per-surface tauk_ref/tauk_k and sc.scale; sc.dc critical offset subtracted from the tearing channel only. Tests (20): constructor validation, pressureless SLAYER-like reduction to det(A')·det(Δ'−Δ_t) via block-diagonal outer, Schur-complement identity for the full coupling case, Q-rescaling via tauk ratios, interchange-channel physical activation, dprime_outer_matrix round-trip against pest3_decompose, msing_max truncation preserves parity-block structure. Paired with a Julia↔Fortran inner-layer GGJ Galerkin benchmark (at CTM-processing/julia_vs_fortran/inner_layer_benchmark/) that runs rmatch's deltac_run qscan on the DIII-D resistive example and the matching Julia solve_inner(GGJModel(:galerkin), ...) at identical (E,F,G,H,K,M,τ_A,τ_R,v1) inputs and Q grid. The benchmark finds a uniform 2.10× factor Julia/Fortran across BOTH channels and ALL Q (not a pole/convergence artifact) — to be investigated as a follow-up; the eigenvalue problem topology is insensitive to this uniform factor so the CoupledFull machinery is usable as-is for root finding via contour-intersection. Co-Authored-By: Claude Opus 4.6 <[email protected]>

…matrix Adds MultiSurfaceCouplingFortran — a literal Julia port of Fortran rmatch/match.f::match_delta (fulldomain=0 branch). This is the full Pletzer-Dewar 4m×4m tearing+interchange coupled dispersion matrix, with the inner-layer amplitudes d^j_± kept as explicit DOFs alongside the outer-region amplitudes C^j_{L,R}, coupled by the ±1 matching identity C^j_L = d^j_+ − d^j_- C^j_R = −(d^j_+ + d^j_-) Motivation: the naive 2m×2m form det(D' − diag(Δ_int, Δ_tear)) = 0 (shipped earlier as CoupledFull) is structurally incorrect because D' lives in the (L,R) side-major basis while the inner-layer output (Δ_tearing, Δ_interchange) lives in the (+,-) parity basis. The two cannot be subtracted directly without an explicit basis transform (Wang-Glasser-Brennan-Liu-Park 2020, Phys. Plasmas 27, 122503, Eq. 11a-11d). Fortran rmatch avoids the transform by keeping both sets of amplitudes alive in a 4m-DOF linear system. This commit mirrors that choice. Validation on DIII-D resistive example (n=1, msing=4): - Julia 4m×4m |det| ∈ [4.6e31, 3.5e39] vs Fortran rmatch [4.0e32, 6.3e36] — same order of magnitude in the same regions. - Same dipolar pole structure at origin, same green/magenta contour sign-change network in both codes. Julia shows some extra contour noise in the lower half-plane consistent with the known uniform 2.10× inner-layer factor + STRIDE-BVP vs Galerkin outer-solve drift (both documented in CTM-processing/julia_vs_fortran/ inner_layer_benchmark/FINDINGS.md). CoupledFull (2m×2m) stays untouched — it remains exported for reference and its 20 tests still pass, but its determinant values should not be used for physical root finding. Use multi_surface_coupling_fortran for that. The patched Fortran rmatch (match_detgrid subroutine added for apples-to-apples grid scans) lives in ../GPEC/rmatch/match.f in the user's local tree; not part of this commit. 26 new unit tests in runtests_dispersion_coupled_fortran.jl covering constructor validation, 1-surface 4x4 hand-verified determinant, 2-surface Fortran-assembly equivalence, Q rotation shift, scale factor, msing_max truncation, pressureless (SLAYER-like) smoke test, GGJ-like m=3 smoke test. Co-Authored-By: Claude Opus 4.6 <[email protected]>

Adds an `inner_kwargs::NamedTuple` field to `MultiSurfaceCouplingFortran` so callers can forward Galerkin grid-tuning parameters (pfac, xfac, nx, nq) to `solve_inner` at every Q evaluation. Matches the Fortran rmatch `&DELTAC_LIST` namelist convention and enables apples-to-apples Julia↔ Fortran dispersion comparisons. Added test verifies the kwarg reaches solve_inner. All 31 existing CoupledFortranMatch tests continue to pass. Context: investigation of the apparent 2.091× Julia↔Fortran discrepancy on DIII-D GGJ inner-layer output revealed it was a **benchmark configuration error**, not a code bug. Fortran rmatch rescales τ_R by η_rdcon/η_user at match.f:212-213 (a deliberate optimization for the η-scan workflow — lets users rerun rmatch at different resistivity without redoing rdcon). When our Julia benchmark drivers fed the raw τ_R from delta_gw.dat into GGJParameters, they were comparing Julia at the "rdcon resistivity" to Fortran at the rmatch.in resistivity. Fix: set rmatch.in::eta to match the value baked into delta_gw.dat. With matched eta, Julia↔Fortran agree to 0.4% across all Q and both channels, with clean 4m×4m determinant agreement in the detgrid benchmark (192×192 narrow-box scan, |det| ranges overlap to < 0.5%). Benchmark updates (in CTM-processing sibling repo, untracked): - run_fortran_deltac_qscan.py + run_fortran_detgrid.py: eta forced to match delta_gw.dat (5.089e-9) - compare_detgrid.py: SLAYER-convention axes (growth on y, rotation on x) and 3-panel layout (Fortran 4m×4m, Julia 4m×4m, Julia m×m — dropped the CoupledFull 2m×2m since it was shown to be structurally wrong). - FINDINGS.md: full write-up of the eta-rescale root cause. Co-Authored-By: Claude Opus 4.6 <[email protected]>

Overhaul of `build_slayer_inputs` + `solve_inner(::SLAYERModel{:fitzpatrick})` so that Julia and Fortran SLAYER produce identical coupled-dispersion det(Q) scans at every plot-frame Q, on the same (geqdsk, kinetic file, slayer.in namelist) inputs. Verified by quantitative 4-hypothesis test at TJ ε=0.001 and β=0.1 benchmark cases: hypothesis median Re median Im J(Q) ~ F(Q) identity +1.01 +1.02 <- eps J(Q) ~ F(Q) identity +0.99 +1.01 <- beta (the three reflection hypotheses all give off-axis ratios) Before this patch the eps_0.001 ratio was (+1.10, -0.98) — a clean Im-axis reflection in Riccati p-space that produced a visually "flipped-about-ω=0" magenta (Im det=0) contour despite all normalized SLAYER parameters (τ_k, S, D_norm, P_perp, P_tor, Q_e, Q_i, d_beta) matching Fortran to <1%. ### `LayerInputs.jl::build_slayer_inputs` Four new kwargs + internal ω_*e/ω_*i computation (port of Fortran `slayer/layerinputs.f:456-459`): * `bt` now also supports a scalar override in addition to a callable or `nothing` (F-spline default). * `R0 = nothing` override magnetic-axis R; default `equil.ro`. Lets the benchmark driver pass the geqdsk RMAXIS literal so both codes use the same reference axis. * `rs_method = :midplane` keeps original θ=0 outboard-midplane chord behaviour by default; `:fsa` activates a θ-mean of √rzphi_rsquared that matches Fortran STRIDE's `issurfint` / `a_surf` flux-surface-averaged minor radius. * `z_i = 1.0` ion charge for the diamagnetic formula; hardcoded to 1 for main D ion in Fortran `layerinputs.f:399`. * `compute_omega_star = true` when `true`, per-surface ω_*e / ω_*i are re-derived from cubic-spline derivatives of (n_e, T_e, T_i) carried in `profiles`, using χ₁ = 2π·equil.psio and the formulae ω_*e = (2π/χ₁)·(T_e·dn_e/dψ / n_e + dT_e/dψ) ω_*i = -(2π/(z_i·χ₁))·(T_i·dn_e/dψ / n_e + dT_i/dψ) (the main-ion density is taken equal to n_e by quasi-neutrality, matching the gpeckf staging convention and Fortran's kin%f(1) after read_kin). Fortran's elementary-charge `e` cancels when T_e, T_i are in eV and dT/dψ is scaled by e, giving the equivalent form above. Setting `compute_omega_star=false` preserves the legacy behaviour where `profiles.omega_e` and `profiles.omega_i` are used as-is (for backward compatibility). ### `Riccati.jl::solve_inner(::SLAYERModel{:fitzpatrick})` Replaced `Q_c = ComplexF64(Q)` (raw pass-through) with the Wick- rotation+conjugate: Q_c = im * conj(ComplexF64(Q)) Fortran `slayer/growthrates.f:337,340` applies `g_tmp = q_in * ifac` with `ifac = (0, +1)` (from `sglobal.f:105`). The algebraically natural Julia port would be `Q_c = Q * im`, but empirically that gives `Julia_det(Q) = Fortran_det(-Q)` (180° rotation), and `Q_c = Q * (-im)` gives `Julia_det(Q) = Fortran_det(-conj(Q))` (Im-axis reflection). The form `im * conj(Q)` substitutes into Julia's Riccati so that `-conj(Q_c) = im·Q` — matching Fortran's internal `g_tmp` — and yields identity. Root cause of the residual Im-axis reflection in Julia's Riccati (suspected: branch selector in `_riccati_f_initial` large-D vs small-D regime, or in the asymptotic `W_bound` sign convention) is not yet identified and is tracked in `~/Desktop/plasma/CTM-processing/CONVENTIONS.md` §4 TODO. Once found, `Q_c = Q * im` should be restored to match Fortran's `ifac` literally. ### Upstream fixes that unblocked this Prior attempts to resolve Julia↔Fortran SLAYER disagreement stalled on three issues that this patch exposes and resolves cleanly: 1. `equil.config.b0exp` (which the benchmark driver was passing as `bt`) is a TOML normalization constant (default 1.0, user- set 2.0), **not** the geqdsk BCENTR. With `bt` now acceptable as a scalar kwarg, the benchmark driver feeds the geqdsk BCENTR literal directly; τ_k J/F ratio went from 5.12× (ε=0.001) / 21.5× (β=0.1) to 1.0009 / 1.0070. 2. `equil.ro` is the GS solver-found axis R, not the geqdsk RMAXIS header value. The new `R0` kwarg lets the driver pass the literal so both codes use the same axis reference. 3. Julia's `surface_minor_radius(..., theta=0)` is outboard- midplane only, not flux-surface-averaged. Fortran STRIDE's `a_surf` IS flux-surface-averaged. The new `rs_method=:fsa` aligns the conventions. After these three plus the Wick-rotation+conjugate, all SLAYER normalized params agree sub-percent across both test cases and the coupled-dispersion panels are pixel-level identical between Julia and Fortran. Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

…erkin scratch buffers Two performance-motivated changes that came out of the julia_vs_fortran benchmark work. Both preserve numerical output exactly (no behaviour change beyond thread-scheduling nondeterminism in the residual evaluations, and even that is serialised before cache insertion so the final result set is deterministic). ### `ContourSearchAMR.jl::amr_scan` Added `parallel = Threads.nthreads() > 1` kwarg and a bulk-eval helper `_bulk_eval_into_cache!` that: * partitions the set of Q-values needed this phase into already-cached vs new (keeps uniqueness), * evaluates all new points via `Threads.@threads` when `parallel=true` and more than one Julia thread is available, * pushes the results into the shared `Dict{ComplexF64,ComplexF64}` cache serially afterwards so no Dict data races occur. Used in both the initial nre0 × nim0 coarse-grid phase and in each refinement pass. The per-call evaluation of `f` (typically a `MultiSurfaceCoupling` or `MultiSurfaceCouplingFortran` closure) is thread-safe because each invocation constructs its own per-surface solver state — the only shared mutable state is the cache, which the helper handles serially. Deterministic output regardless of thread count. On the 100×100 + 4-pass benchmark scan this cut Julia SLAYER AMR from ~60s to ~15s on an Apple M2 Max (8 threads). ### `GGJ/Galerkin.jl::GalerkinWorkspace` + `_assemble_and_solve!` Added five preallocated scratch buffers to `GalerkinWorkspace` (`cell_mat_buf`, `cell_mat_ext_buf`, `cell_rhs_ext_buf`, `ab_buf`, `rhs_buf`) sized to the max case (`np+1=4`) used at any cell type, and re-use them via `fill!(buf, 0)` inside the per-cell loop. Previously each cell called `zeros(ComplexF64, ...)` which accumulated thousands of MiB of allocations over a full dispersion scan. Same numerical output; the cell-matrix sub-slices are explicitly zeroed before use and smaller cells (e.g. `CT_EXT` with `cell.np=1`) rely on the remaining buffer elements staying zero from the previous `fill!` call. Measured on the TJ ε=0.001 benchmark (nx=256, cutoff=20, tol_res=1e-7, msing=2): Galerkin det evaluation dropped from ~4.2 MiB allocs / call to ~30 kiB / call, with a corresponding 20-25% wall-time reduction in the GGJ AMR scan. Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

…ar coupled residual `MultiSurfaceCouplingFortran` (aka the 4m×4m Pletzer-Dewar tearing+ interchange dispersion matrix, port of Fortran `rmatch/match.f::match_delta` fulldomain=0 branch) was adding `+ sc.dc` to BOTH the inner-layer interchange and tearing Δ channels before assembling the coupled matching block: # CoupledFortranMatch.jl, before: delta1 = resp.interchange * sc.scale + sc.dc # WRONG delta2 = resp.tearing * sc.scale + sc.dc # WRONG The code comment claimed this was "per the Fortran convention (χ_parallel shift that acts on the outer diagonal before matching)." That is NOT in Fortran — `match.f:508-519` assembles the fulldomain=0 block directly from the raw `delta1 = deltar(ising, 1)` / `delta2 = deltar(ising, 2)` with no Δ_crit offset anywhere: ! Fortran match.f (fulldomain=0): delta1 = deltar(ising, 1) delta2 = deltar(ising, 2) mat(idx3, idx3) = -delta1 mat(idx3, idx4) = delta2 mat(idx4, idx3) = -delta1 mat(idx4, idx4) = -delta2 The Δ_crit proxy represents a slab-layer χ_parallel-matching correction and is meaningful only for tearing-only models like SLAYER (which drops the interchange channel and needs a proxy for the missing Glasser/ Mercier stabilization). GGJ's 4m×4m Pletzer-Dewar matching already includes the interchange channel explicitly (`resp.interchange`), so adding `sc.dc` double-counts that physics. ### Fix 1. `CoupledFortranMatch.jl:179-180`: drop `+ sc.dc` on both channels. delta1 / delta2 are now the raw inner-layer outputs, matching match.f:508-519 bit-for-bit. 2. `SurfaceCoupling.jl`: remove the `dc::Real=0.0` kwarg from `surface_coupling(model::GGJModel, ...)`. The SLAYER and generic overloads still accept it — SLAYER genuinely needs it for its slab-layer Δ_crit subtraction. The `SurfaceCoupling.dc` struct field is hard-wired to 0 for GGJ callers, making the API reflect the physics. ### Tests - `test/runtests_dispersion_coupled.jl`: 42 / 42 pass - `test/runtests_dispersion_residual.jl`: 20 / 20 pass (Both test files construct `surface_coupling(GGJModel, ...)` with positional args only — no call sites broken.) ### Impact For the julia_vs_fortran benchmark, this is a no-op when the driver was already passing `dc=0.0` for GGJ (the safe default we settled on earlier in the session). The fix prevents the footgun of anyone else accidentally passing a nonzero `dc` to a GGJ coupling and getting physically wrong results. Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

… and v1 GGJ: - LayerInputs.jl: changed `v1 = 1.0` placeholder to `v1 = rg.v1_local / equil.params.volume`. This is the dV/dψ normalization that `rescale_delta` consumes as `v1^(2*p1)` to convert raw Galerkin Δ to outer-region matching units. Matches Fortran resist.f:144 (`sing%restype%v1 = v1/volume`) and match.f:1078 (`deltar = deltar * sfac**(2*p1/3) * v1**(2*p1)`). Previously, on realistic shaped equilibria where v1_local/volume != 1, Julia's GGJ Δ disagreed with Fortran by `(v1_local/volume)^(2*p1)`. Analytical TJ/Solovev cases hid the bug because v1_local/volume happens to hover near unity there. SLAYER: - LayerInputs.jl: changed `dr_val = 0.0` default to `dr_val = nothing`. When `nothing` is passed, build_slayer_inputs auto-derives the per-surface resistive interchange index `D_R = E + F + H²` from `sing.restype` (already populated by `resist_eval_all!`). Without this, the slayer_panels benchmark driver was reading a scalar dr_val=-0.1 from a Fortran namelist and applying it uniformly to every surface, producing dc_tmp values that didn't match Fortran's per-surface STRIDE-derived values. With `nothing` default, dc_type in {:lar, :rfitzp, :toroidal} now produces a non-zero per-surface dc_tmp without manual configuration. dgeo_val behaves analogously but errors clearly if dc_type=:toroidal is requested without an explicit value (auto-derive needs ⟨|∇ψ|²⟩ FSA which isn't yet exposed in ResistGeometry — TODO). NOTE on Fortran/STRIDE divergence: Julia uses D_R correctly per Connor-Hastie-Helander 2015 (PPCF 57 065001) Eq. 59. Fortran STRIDE has a one-character bug in stride_netcdf.f:100 — `dr_rationals(i) = locstab%f(1)/respsi` uses index 1 (= D_I, the Mercier criterion) instead of index 2 (= D_R, the resistive interchange). Julia and Fortran will therefore disagree on dc_tmp magnitude by ~D_I/D_R per surface (~3-4× on DIII-D) until that upstream Fortran bug is fixed. The disagreement is documented at the build_slayer_inputs docstring. Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

…rowth_rates Adds `pole_threshold_adaptive::Bool = false` to SLAYERControl. When true, `run_slayer_from_inputs` overrides `control.pole_threshold` per scan with `|mean(Δ)|` over the dispersion-residual array before calling `find_growth_rates`. Backward-compatible (default false uses the literal `pole_threshold`). Justification: the hardcoded default `pole_threshold=10.0` is too restrictive when |Δ| spans 8+ orders of magnitude (typical for SLAYER coupled-dispersion scans). All intersections then get classified as poles and zero roots are returned. The adaptive recipe — empirically matching the Python `10·median(|Δ|)` heuristic and the omfit `|mean(Deltas_AMR)|` recipe — yields the correct root identification on the DIIID benchmark and TJ βₚ scan cases (verified at βₚ=0.1 coupled_rfitzp: 6 roots / 8 poles vs 0 roots with the static threshold). Plumbing changes: - Control.jl: new field + docstring - HDF5Output.jl: written to /slayer/settings/pole_threshold_adaptive - run_slayer.jl: `_pole_threshold_for(scan)` closure dispatches per-scan - Runner.jl: import Statistics.mean

@threads

… default 1 (serial) eliminates DIII-D 147131 thread-race The parallel BVP path in `parallel_eulerlagrange_integration` was always invoking `Threads.@threads :static` over the FM chunks, ignoring the `parallel_threads` field on `ForceFreeStatesControl`. On numerically delicate equilibria (e.g. DIII-D 147131 at βₚ ≈ 0.07) this exposed a sub-tolerance nondeterminism: chunk crossings whose post-jump matrices depend on the order of independent FP operations across threads, producing intermittently divergent FM matrices and intermittent BVP failures. The algorithm is correct; the wall-time interleaving of parallel chunks was perturbing it within tolerance. Fix: * `Riccati.jl`: branch on `bvp_threads = clamp(parallel_threads, 1, julia_nthreads)`. `bvp_threads == 1` runs the chunks serially on the calling thread (race-free, bit-deterministic). Otherwise, the existing `:static` parallel path is used. * `ForceFreeStatesStructs.jl`: document `parallel_threads` semantics, default `1`, and the cost (~14% slower than 2-thread on DIII-D 147131 reference). Verified: with `parallel_threads = 1` (default) and `JULIA_NUM_THREADS = 2`, the DIII-D 147131 βₚ=0.07 reference Δ' diagonal matches CONVENTIONS.md §6 exactly: q=2: +7.92 - 0.03i q=3: -5.24 - 0.30i q=4: -40.20 + 209.91i q=5: +126.6 - 169.24i in 54.5 s wall (single 4-singular-surface coupled BVP). No regressions on TJ. Production scans should keep the default; users with robust equilibria and strict wall-time budgets can opt in to `parallel_threads > 1` knowing the trade-off. Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

…BVP speedup; bit-identical Δ' in 15-trial DIII-D 147131 sweep) Empirical reliability sweep on DIII-D 147131 βₚ≈0.07 (5 trials at each of parallel_threads ∈ {1, 2, 4}, JULIA_NUM_THREADS=4, post-JIT, single Julia session) showed: parallel_threads | wall (avg, single 4-singular-surface coupled BVP) -----------------|------------------------------------------------- 1 (serial) | 9.25 s — bit-deterministic by construction 2 | 7.37 s — bit-identical Δ' in all 5 trials (+20.3%) 4 | 7.51 s — bit-identical Δ' in all 5 trials (+18.9%) Δ′ diagonals were bit-identical across all 15 trials and matched the §6 reference values exactly. Speedup saturates at 2 threads — the BVP has ~10 FM chunks, so 2 threads is enough to amortize them; 4 adds scheduling overhead with no benefit on this BVP. Bumping default to 2 captures the ~20% wall-time win on production scans. The serial path remains available (`parallel_threads = 1`) as a deterministic fallback if the historical intermittent race re-manifests on a delicate equilibrium. Documentation in `ForceFreeStatesControl` docstring updated to record the trade-off and the empirical reliability data. Use `parallel_threads = 1` (NOT `use_parallel = false`) if a parallel run ever diverges — `use_parallel = false` produces silently wrong Δ' values (see CONVENTIONS.md §7). Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>

…aster) The Fitzpatrick `riccati_f` ODE is a 1-equation system. The prior code modeled `W` as a 1-element `Vector{ComplexF64}` with an in-place RHS (`_riccati_f_rhs!(dW, W, params, x)`); every Rosenbrock stage allocated fresh `dW` intermediates. Converting `W` to a `ComplexF64` scalar with an out-of-place RHS removes those per-stage heap allocations and lets stage updates stay on the stack. Per-call benchmark (1000 calls, Rodas5P, identical inputs): vector form: 1.62 ms / call scalar form: 0.96 ms / call (41% faster) Signature changes: _riccati_f_rhs!(dW, W, params, x) -> nothing --> _riccati_f_rhs(W::Number, params, x) -> ComplexF64 _riccati_f_jac!(J, W, params, x) -> nothing --> _riccati_f_jac(W::Number, params, x) -> ComplexF64 solve_inner ODE state: u0 = ComplexF64[W_bound]; ODEFunction{true}(...) --> u0 = ComplexF64(W_bound); ODEFunction{false}(...) Solver-agnostic. Rodas5P stays the default. The change works equally well under any OrdinaryDiffEq stiff solver (Rosenbrock / SDIRK / BDF) since they all support scalar `u0` via the out-of-place form. Validation (against the temporary baseline at SLAYER_coupling_paper/ regression_temporary/, 88 TJ records frozen pre-change): TJ uncoupled_2over1_rfitzp at βₚ=0.001 γ baseline = +4.0552247503e+00 kHz γ scalar = +4.0551819762e+00 kHz relative drift = 1.05e-5 (within solver-replacement noise) TJ coupled_rfitzp at βₚ=0.07 (exercises full BVP path) γ baseline = -8.1071602485e-03 kHz γ scalar = -8.1071881463e-03 kHz relative drift = 3.44e-6 n_valid_roots = 26, n_poles = 27 (exact match to baseline topology) check_regression.py --dry --scope tj : 88/88 pass (5e-4 abs/rel tolerance on integrator outputs, exact match on topology fields). Production wall-time on the coupled-BVP case: baseline (vector form): ~14 min (slowest of 4 parallel cases per βₚ) scalar form: ~10 min (~29% reduction) In contrast to the prior KenCarp4 solver-swap attempt (commit 5a9026a8, reverted as 2b1e1b0f), which looked like a 38% per-call win in synthetic tests but came out 17% SLOWER in production, this change shows consistent gains from per-call benchmark through to full production scan. The reason the wins translate cleanly: the scalar form makes the existing solver faster without changing its convergence path or step-control behaviour, so production characteristics scale linearly from the micro-benchmark. The companion KenCarp4 swap stays deferred (tracked in todos) until we have direct production-side per-Q timing instrumentation to understand the bench/production discrepancy. Test infrastructure also committed: profiling/profile_slayer_amr.jl CPU + alloc profile harness profiling/test_riccati_solver_convergence.jl 7-solver convergence sweep

…tring Empirical finding from Phase 2.5 of the AMR speedup work: sub-percent floating-point differences between ODE solvers cascade through the AMR's zero-crossing flagging and produce structurally different cell trees, not just numerically-noisy Δ values. Concrete observation on TJ coupled_rfitzp at βₚ=0.07 under the scalar ODE form (commit b17e0b43): Solver SLAYER wall γ valid_roots poles Rodas5P ~10 min -8.107e-3 kHz 26 27 KenCarp4 ~9 min -8.107e-3 kHz 43 34 KenCarp4 is per-call faster (consistent with the convergence-test results), but its slightly different Δ at AMR cell corners flips many "refine" / "no-refine" decisions and lands on a substantially different final cell list. The most-unstable root (γ) agrees to 2.1e-5 relative, but the inventory of secondary roots and poles differs by ~17 / ~7. Implication: solver swaps are NOT pure per-call optimizations. Future attempts need to be validated against the topology fields (`n_valid_roots`, `n_poles`), not just γ. The temporary regression harness at SLAYER_coupling_paper/regression_temporary/check_regression.py already treats these as exact-match fields, which correctly gates solver swaps. The 92-record baseline serves as a topology fingerprint.

…30% additional per-call speedup) The Fitzpatrick `riccati_f` ODE coefficients fA, fA', fB, fC use parameters (Q, Q_e, Q_i, P_perp, P_tor, D_norm, iota_e) that are CONSTANT across the integration. The prior code recomputed `Q*(Q+iQi)`, `Q+iQe`, `D²·iota_e⁻¹` etc. at every RHS evaluation — tens of thousands of redundant multiplications per `solve_inner` call. This commit lifts the x-independent quantities into a `_RiccatiConsts` struct built once per `solve_inner` call: Q_plus_iQe constant part of denom = (Q + iQe + x²) A = Q · (Q + iQi) fB constant term B = (Q + iQi)·(P_perp + P_tor) fB · x² coefficient C = P_perp · P_tor fB · x⁴ coefficient E = P_perp + (Q + iQi)·D² fC · x² coefficient G = P_tor · D² / iota_e fC · x⁴ coefficient The hot RHS (`_riccati_f_rhs`) and Jacobian (`_riccati_f_jac`) now access only the bundled constants and `x`, doing ~3 muls + 1 division per call instead of ~10 muls + 2 divisions. Per-call benchmark (1000 calls, Rodas5P, identical inputs): prior (scalar form, post b17e0b43): 0.96 ms / call precompute (this commit): 0.67 ms / call (-30% per call) cumulative vs vector-form baseline: 1.62 → 0.67 ms (-59%, 2.42× faster) Validation against the temporary baseline at SLAYER_coupling_paper/ regression_temporary/: TJ coupled_rfitzp at βₚ=0.07 (full BVP path) γ baseline = -8.1071602485e-03 kHz γ precompute = -8.1071881463e-03 kHz relative drift = 3.44e-6 (same as scalar-only Phase 2.3 baseline) n_valid_roots = 26, n_poles = 27 (exact match to baseline topology) check_regression.py --dry --scope tj : 88/88 pass Production wall on TJ coupled_rfitzp at βₚ=0.07: vector-form baseline: ~14 min scalar form (Phase 2.3): ~10 min scalar + precompute: ~9 min (~36% cumulative reduction) The active SLAYER step alone is now ~41% faster than baseline. Production wall scales sub-linearly because main() / find_growth_rates / file-write overheads remain unchanged. Implementation note — algebraic simplification rejected: A natural further optimization is `fA' = 1 − 2·fA` (algebraic identity: (denom − 2p²)/denom = 1 − 2·(p²/denom) = 1 − 2·fA). It saves one complex division per call. However, when tested, the integrator's adaptive stepping near marginal stability compounded ULP-level differences in fA' across thousands of steps, producing ~3e-3 relative γ drift versus this form's 3e-6. The drift was within the regression's abs-tolerance gate but still a real precision regression. Reverted — kept the explicit `(denom − 2·p²)/denom` form, which preserves bit-identical Δ at warm benchmark points vs the scalar-form baseline.

@warn

…tion kwargs Two additive kwargs to support convergence-vs-resolution studies and graceful behaviour when the cell-count safety rail is hit: snapshot_callback::Union{Nothing,Function} = nothing If provided, called at the end of each AMR pass (and once for the initial grid, pass=0) with arguments (pass::Int, cells::Vector{AMRCell}, cache::Dict{ComplexF64,ComplexF64}). The callback receives live references; copy if persistence is needed. Used by convergence studies to extract intermediate γ at each pass count from a SINGLE AMR run (avoids re-running for every target pass). max_cells_action::Symbol = :error :error (default, prior behaviour) raises when length(cells) > max_cells. :warn_truncate logs a @warn, stops further refinement in the current pass, and exits the outer pass loop — leaving a usable AMRResult with the partial cell tree. Useful for resolution-sweep studies that deliberately push max_cells to bound runtime. Backward compatibility: defaults preserve the exact prior behaviour. Validated via regression rerun of TJ coupled_rfitzp at βₚ=0.07 (88/88 pass, γ + topology bit-identical to pre-change baseline).

d-burg · 2026-06-17T20:07:16Z

Branch update since 2026-05-18

Commit summary

f50b4bd (06-16) Removed all Fortran-source citations from the Tearing backend; cite papers/textbooks only. Renamed Fortran-named identifiers: CoupledFortranMatch.jl -> CoupledFullMatch.jl, MultiSurfaceCouplingFortran -> MultiSurfaceCouplingFull, multi_surface_coupling_fortran -> multi_surface_coupling_full (+ test file). Dropped one-off profiling/ scripts from VCS.
6154cf5 (06-15) Verified the riccati_f dispersion solver and del_s layer-width diagnostic term-by-term against the reference layer implementation and Fitzpatrick's published layer theory.
1e1c3e3 (06-15) Regression harness: NaN-safe JSON caching/diff so SLAYER no-root cases (NaN Q_root) cache and compare correctly.
369911c (06-15) Neoclassical resistivity now the default (tau_R = mu0*rs^2/eta, Sauter F_33), with :spitzer / :spitzer_harm toggles; derived the Wick-rotation/sign convention; scoped GGJ as a Delta-diagnostic (parity channels via ggj_inner_deltas, future-work warning on contour-matching gamma, no crash); loud Delta'=0 warning; try/catch around the SLAYER stage; @autodocs for Tearing.Dispersion/Runner; full HDF5 settings provenance.
335f08e (06-12) Robust SLAYER gamma extraction (isolate-then-polish, validity gate, adaptive pole threshold) + pre-merge audit fixes.
a12e25f (06-08) New resistive-layer-thickness (del_s Riccati) diagnostic.
6f2a76e (05-29) Pre-merge audit: Coupled* consolidation, multi-n test conditioning, SLAYER robustness.
9609432 (05-27) Unmasked pre-existing SLAYER + multi-n fullruns test failures after the perf/riccati merge.

Net feature state

SLAYER tearing growth rates: equilibrium + Delta' -> per-surface two-fluid layer Riccati -> complex-plane scan (brute-force or AMR) -> contour-intersection root with polish/validity gating. Uncoupled and coupled (reduced m x m and full 4m x 4m Pletzer-Dewar) paths.
del_s resistive-layer-width diagnostic per rational surface.
Neoclassical resistivity (Sauter/Redl) default with Spitzer / Spitzer-Härm fallbacks, configurable in [SLAYER].
GGJ inner layer available as a two-parity Delta diagnostic; contour-matching gamma extraction marked future work.
TOML-driven runner, HDF5 output with full settings provenance, regression cases (solovev_slayer_n1, diiid_slayer_n1).

Reviewer notes

Changes the default SLAYER resistivity: tau_R-derived outputs (Lundquist S, P_perp, D_norm, tauk, growth rates) move vs prior baselines; stable/unstable surface pattern is unchanged (DIII-D still gives an unstable 2/1).
Renames the exported multi_surface_coupling_fortran / MultiSurfaceCouplingFortran symbols.
Tests green (dispersion, SLAYER, resist-eval, fullruns); regression harness confirms only resistivity-dependent quantities shift.
Deferred: ~300 Fortran citations remain in non-Tearing modules (pre-existing develop code), to be cleaned on a separate branch.

Copilot

Pull request overview

Adds a new tearing-mode growth-rate analysis stack to GPEC by introducing a Tearing umbrella module (InnerLayer + Dispersion + Runner) and wiring it into the main pipeline, along with supporting utilities (kinetic profiles, resistivity, constants) and ForceFreeStates outputs needed to provide Δ′ inputs.

Changes:

Introduce src/Tearing/ with SLAYER (Riccati-based) and GGJ inner-layer models, dispersion residuals/scans (brute-force + AMR), and a TOML-driven runner with HDF5 output.
Extend ForceFreeStates to expose/store raw and projected Δ′ matrices and to compute GGJ geometric coefficients per singular surface (ResistEval).
Add extensive unit tests, new example configurations, and regression-harness cases for SLAYER workflows.

Reviewed changes

Copilot reviewed 60 out of 62 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
test/runtests.jl	Adds new tearing-related test includes to the suite.
test/runtests_slayer_runner.jl	End-to-end runner/control/HDF5 structure tests for SLAYER.
test/runtests_slayer_riccati.jl	Validates SLAYER Riccati solver behavior and diagnostics.
test/runtests_slayer_params.jl	Tests SLAYER parameter construction and resistivity closures.
test/runtests_slayer_inputs.jl	Builds SLAYER inputs from equilibrium + kinetic profiles.
test/runtests_resist_eval.jl	Tests GGJ geometric coefficient evaluation and GGJ input builder.
test/runtests_kinetic_profiles.jl	Tests `KineticProfiles` constructors and HDF5/TOML loading.
test/runtests_dispersion_scan.jl	Tests brute-force scan and contour-based growth-rate extraction.
test/runtests_dispersion_residual.jl	Tests per-surface residual arithmetic (`SurfaceCoupling`).
test/runtests_dispersion_polish.jl	Tests local root polishing and validity gating behavior.
test/runtests_dispersion_amr.jl	Tests AMR scan path and triangulation-based extraction.
src/Utilities/Utilities.jl	Registers/exports new Utilities submodules (constants, profiles, resistivity).
src/Utilities/PhysicalConstants.jl	Adds SI constants aligned with Fortran conventions.
src/Utilities/KineticProfiles.jl	Adds spline-based kinetic profile container + TOML/HDF5 loaders.
src/Tearing/Tearing.jl	New umbrella module binding InnerLayer/Dispersion/Runner.
src/Tearing/Runner/Runner.jl	Runner module wiring and exports for the SLAYER pipeline.
src/Tearing/Runner/Result.jl	Defines `SLAYERResult` and empty-result helper.
src/Tearing/Runner/HDF5Output.jl	Writes SLAYER results into an HDF5 `slayer/` group.
src/Tearing/Runner/Control.jl	Defines `SLAYERControl`, TOML parsing/flattening, validation.
src/Tearing/InnerLayer/SLAYER/SLAYER.jl	SLAYER inner-layer module entrypoint and exports.
src/Tearing/InnerLayer/SLAYER/Riccati.jl	Implements Fitzpatrick Riccati Δ(Q) solver with cached constants + Jacobian.
src/Tearing/InnerLayer/SLAYER/LayerThickness.jl	Adds `del_s` Riccati layer-thickness diagnostic (`LayerWidths`).
src/Tearing/InnerLayer/InnerLayerInterface.jl	Defines shared inner-layer interfaces/types (`InnerLayerResponse`, etc.).
src/Tearing/InnerLayer/InnerLayer.jl	InnerLayer module now includes GGJ + SLAYER and exports both.
src/Tearing/InnerLayer/GGJ/Reference.jl	Adjusts GGJ reference fixture construction style.
src/Tearing/InnerLayer/GGJ/LayerInputs.jl	Builds physical GGJParameters from equilibrium+sings+profiles.
src/Tearing/InnerLayer/GGJ/GGJParameters.jl	GGJParameters now subtype `InnerLayerParameters`; doc cleanup.
src/Tearing/InnerLayer/GGJ/GGJ.jl	Exposes `build_ggj_inputs` and resistivity-model types.
src/Tearing/Dispersion/SurfaceCoupling.jl	Introduces callable per-surface dispersion residual container.
src/Tearing/Dispersion/Dispersion.jl	New dispersion module that ties residuals/scans/extraction together.
src/Tearing/Dispersion/CoupledFullMatch.jl	Implements full 4m×4m Pletzer–Dewar matching determinant.
src/Tearing/Dispersion/Coupled.jl	Implements reduced m×m coupled tearing-only determinant residual.
src/Tearing/Dispersion/BruteForceScan.jl	Adds regular-grid Q-plane scan with optional threaded evaluation.
src/InnerLayer/SLAYER/Slayer.jl	Removes old placeholder SLAYER file (migration to `src/Tearing/`).
src/InnerLayer/InnerLayerInterface.jl	Removes old inner-layer interface file (migration to `src/Tearing/`).
src/GeneralizedPerturbedEquilibrium.jl	Wires in `Tearing` module + backward-compatible top-level aliases; runs SLAYER stage and writes outputs.
src/ForceFreeStates/Riccati.jl	Persists raw Δ′ matrix and adds `pest3_decompose` helper.
src/ForceFreeStates/ResistEval.jl	Adds per-surface GGJ geometric coefficient evaluation and populator.
src/ForceFreeStates/ForceFreeStatesStructs.jl	Adds `restype` to `SingType` and `delta_prime_raw` to internal state.
src/ForceFreeStates/ForceFreeStates.jl	Includes new `ResistEval.jl` in module assembly.
src/Equilibrium/DirectEquilibrium.jl	Makes separatrix crossing search robust to non-monotone ψ outside LCFS via pre-scan bracketing.
regression-harness/src/reporter.jl	Allows NaN JSON parsing; formatting/indent fixes.
regression-harness/src/extractor.jl	Emits/compares JSON arrays with NaNs preserved; treats NaN==NaN in diffs.
regression-harness/cases/solovev_slayer_n1.toml	Adds SLAYER regression case spec (Solovev analytic).
regression-harness/cases/diiid_slayer_n1.toml	Adds SLAYER regression case spec (DIII-D-like example).
Project.toml	Adds `DelaunayTriangulation` dependency/compat for AMR triangulation extraction.
examples/Solovev_ideal_example/gpec.toml	Adds an inline SLAYER config block + synthetic profiles.
examples/DIIID-like_SLAYER_example/gpec.toml	Adds a new example config for SLAYER on a DIII-D-like equilibrium.
docs/src/inner_layer.md	Updates docs page to describe `Tearing` module and adds autodocs sections.
docs/make.jl	Renames nav entry from “Inner Layer” to “Tearing”.
.gitignore	Ignores scratch/profiling and Claude agent working dirs.

logan-nc · 2026-06-18T15:47:43Z

@d-burg it looks like there are surprisingly few technical review requests above given the size of this PR. Nice!

Is this ready for @matt-pharr to give a human review?

Please coordinate with @matt-pharr on how to use the existing read_kinetic_profiles function from Equilibrium and use that to set all the local rational surface quantities. Note that the old listing of things in the namelist IO should not make it to Julia GPEC - it does not play well with multi-n, which we support in Julia. Add columns to the input ascii if you need. Lets have one, consistent interface for resistive and kinetic physics.

…56 grid Co-Authored-By: Claude Opus 4.8 <[email protected]>

Co-Authored-By: Claude Opus 4.8 <[email protected]>

…ptures q=6) Co-Authored-By: Claude Opus 4.8 <[email protected]>

Co-Authored-By: Claude Opus 4.8 <[email protected]>

Add read_kinetic_file (extension dispatch) and write_kinetic_h5 for a GPEC HDF5 kinetic schema (psi, n_i, n_e, T_i, T_e, omega_E; optional omega_tor, chi_e, chi_phi). Legacy 6-column ASCII (.gpeckf/.kin) still reads via the same path; the NTV loader is refactored onto the dispatcher with bit-identical results. Adds round-trip / ASCII-vs-HDF5 equivalence tests and docs; .gitignore tracks kinetic .h5 inputs. Co-Authored-By: Claude Opus 4.8 <[email protected]>

…brium, analytic ExB, HDF5 profiles Edge-current-matched equilibrium (q=6 retained at psihigh=0.995; fixed mpsi to avoid edge under-resolution, see #302). Kinetic profiles carry an analytic ExB rotation (omega_E = omega_tor - omega_dia) and synthetic transport diffusivities chi_e/chi_phi, shipped in the new .h5 (gpec.toml points at it; .gpeckf retained for back-compat). Relax NTV psi-quadrature tolerances as an interim mitigation (see #303). Re-pin the parallel-integration regression (et[1], Delta-prime diagonal) for the updated equilibrium. Co-Authored-By: Claude Opus 4.8 <[email protected]>

… torque Drop leftover println/printf debug in tpsi! that flooded stdout during the psi-quadrature (see #303). Co-Authored-By: Claude Opus 4.8 <[email protected]>

Drop the unverified 'parallel and standard ODE paths agree to rel_diff=0' claim from the DIIID-like testset (it only runs the parallel path); note the cross-path agreement is verified on the Solovev testset. Addresses PR review. Co-Authored-By: Claude Opus 4.8 <[email protected]>

… guard, separatrix) Four fixes from the Copilot review on PR #238: - LayerParameters.jl: iota_e = Q_e/(Q_e - Q_i) is singular when the electron/ion diamagnetic frequencies are degenerate (Q_e == Q_i). The old guard set iota_e = 0, which then produced Inf/NaN in the downstream Riccati coefficients (1/iota_e). Now throw a clear ArgumentError naming the cause. - LayerParameters.jl: fix the SLAYERParameters doc table, which listed Q_i = +tauk*omega_i; the code and sign-convention section use -tauk*omega_i. - run_slayer.jl: raise a clear error when eltype(params) does not match control.inner_model (previously an opaque MethodError in _build_surface_coupling); and warn that coupling_mode=:coupled with a GGJ model uses the reduced m x m tearing-only determinant, dropping the GGJ interchange (Glasser) stabilization — use multi_surface_coupling_full for coupled GGJ. - DirectEquilibrium.jl: the separatrix pre-scan only bracketed on f_prev*f_curr < 0, missing an exact psi=0 at a sampled point; handle the exact-zero case by returning that R directly. Behavior unchanged on the example cases (DIII-D SLAYER still gives the unstable 2/1 at +256 Hz); SLAYER test suite passes. Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

…into feature/tearing-growthrates

…le; plumb χ⊥/χ_φ Per the PR review request to have one consistent kinetic-profile interface for resistive and kinetic physics, SLAYER now reads kinetic profiles through the shared Equilibrium.read_kinetic_file (the standardized KineticProfileData object) instead of its own inline [SLAYER.profiles] TOML block. - run_slayer reads control.profile_file (HDF5 GPEC kinetic schema or ASCII) via read_kinetic_file; builds the spline KineticProfiles plus χ⊥(ψ)/χ_φ(ψ) callables from the file's chi_e / chi_phi and threads them into build_slayer_inputs (previously chi_perp/chi_tor were scalar-1.0 defaults). The scalar control.chi_perp/chi_tor remain as fallbacks when the file carries no χ (e.g. ASCII tables); a warning fires in that case. - Remove the inline-profile path: drop SLAYERControl.profile_source and the [SLAYER.profiles] handling, and remove Utilities.kinetic_profiles_from_toml / kinetic_profiles_from_h5 (the namelist-style listing) and their exports. The KineticProfiles interpolant type stays as the internal spline container. - Examples retargeted to the standardized files: the DIII-D SLAYER example uses the shipped TkMkr_D3Dlike_Hmode_kinetic.h5 (real n_e/T_e/T_i/omega_E + χ⊥/χ_φ); the Solovev example uses a generated synthetic solovev_kinetic.h5 (flat n_e, declining T, χ⊥=χ_φ=1) preserving its sanity-check behavior. - Update the KineticProfiles and slayer_control_from_toml tests accordingly. Verified: full suite green; χ⊥≠χ_φ from the DIII-D file now gives P_perp≠P_tor per surface (e.g. 2/1 P_perp=50.3, P_tor=34.7), with the 2/1 unstable. Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

…o inner surfaces - SLAYER χ loading: a chi_e/chi_phi dataset that is absent OR all-zero now falls back to the scalar control.chi_perp/chi_tor (χ must be positive; χ=0 ⇒ τ_⊥→∞). This lets a kinetic file keep the chi_e/chi_phi keys set to zero to defer to the scalar values, rather than omitting them. - Regression harness: add a `first_N` extract mode (real values of the leading N vector elements). Use it in diiid_slayer_n1 to golden-pin only the inner three rational surfaces (2/1, 3/1, 4/1) for Q_root/ω/γ/no_root — the Δ'/γ contour search is numerically unreliable on the outermost surfaces (5/1, 6/1, 7/1 near the edge), so those are not tracked. Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

The Solovev analytic equilibrium does not yield meaningful tearing physics through SLAYER — its Δ' / inner-layer dispersion produces no usable growth- rate root (the coupled determinant returns NaN), so the solovev_slayer_n1 regression case only pinned a no-root placeholder plus structural params. SLAYER examples and regression now use the realistic DIII-D-like geqdsk equilibrium exclusively (examples/DIIID-like_SLAYER_example, diiid_slayer_n1). - Remove regression-harness/cases/solovev_slayer_n1.toml - Remove the [SLAYER] block from the Solovev example (now ideal-only) with a note pointing to the DIII-D SLAYER example - Remove the generated examples/Solovev_ideal_example/solovev_kinetic.h5 Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

Resolve conflicts from develop's local-stability refactor (Mercier.jl removed, Bal.jl -> Ballooning.jl with D_I from det(d0bar), new RootAreaWeighted/Kinetic/ CoordinateInvariant modules) against the tearing-growthrates stack. Conflict resolutions: - examples/DIIID-like_ideal_example/gpec.toml: keep the tearing-branch grid config (grid_type=ldp, mpsi=256, psihigh=0.995, explicit power_*=0) that the diiid_n1 regression is pinned to; graft develop's fuller jac_type comment. - DirectEquilibrium.jl: take develop's Fortran-faithful Newton separatrix (dr=-psi/psi_r with restart shifts); the auto-merged call sites already use its 2-arg signature. - GeneralizedPerturbedEquilibrium.jl: keep resist_* imports, take develop's ballooning_alpha_boundary(s), drop the now-deleted mercier_scan!; adopt develop's "Root-area-weighted" wording (matches rootA_* fields). Fix breakages from develop's API changes (not flagged as textual conflicts): - runtests_resist_eval.jl / runtests_slayer_inputs.jl: develop's analytic setup_equilibrium now requires a SolovevConfig from [SOL_INPUT] (two-arg form). - runtests_resist_eval.jl: the deleted mercier_scan! reference D_I is rewired onto the ballooning det(d0bar) route (prepare_ballooning_coefficients(...).di). That is an independent discretization from the surface-average E,F,H integrals, so the cross-check tolerance is relaxed 1e-3 -> 8e-2 (measured ~2-5% agreement). Full test suite green (52 groups, 0 failures). Regression harness pending: the separatrix change is expected to shift diiid_n1 equilibrium quantities and may need a re-pin. Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

…reeStates flag refactor The develop merge replaced ForceFreeStatesControl's bal_flag/mer_flag with a single local_stability_flag (unified Mercier + ballooning) and removed delta_mband. The SLAYER example gpec.toml — a feature-branch file the merge did not touch — still set the old keys, so its [ForceFreeStates] section raised a keyword-argument error when run against the merged code (caught only by the end-to-end regression, not the unit suite). - bal_flag=false + mer_flag=true -> local_stability_flag=true (SLAYER needs the Mercier D_I/D_R local-stability profiles) - drop delta_mband (no longer a ForceFreeStatesControl field) - replace the "PENTRC" legacy-name reference in the kinetic_source comment with "kinetic NTV model" per the example-annotation convention Verified: diiid_slayer_n1 runs clean against the merged tree (6 surfaces, roots for 2/1,3/1,4/1, enabled=1). Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

d-burg · 2026-06-26T23:50:32Z

@matt-pharr I believe this is ready for human review!

Copilot's four inline comments are all addressed in code and marked resolved.

Per Nik's request, SLAYER no longer has its own namelist-style profile listing. Kinetic profiles (n_e, T_e, T_i, omega_E, and chi_perp/chi_phi) are read through the shared Equilibrium.read_kinetic_file / KineticProfileData interface.

develop is now merged in and the conflicts are resolved.

Verification:

full unit suite is green
regression harness (pre-merge -> merged) clean: diiid_n1 shows ~1% drift in beta/li/volume from develop's new separatrix (expected); Solovev cases unchanged; diiid_slayer_n1 runs clean (2/1, 3/1, 4/1 roots).

Do we have anything to coordinate on structuring the inner-layer matching for RPEC?

SLAYERControl.dr_val/dgeo_val defaulted to 0.0, which the layer-input builder treats as an explicit value, so the D_R = E+F+H² auto-derivation (reached only when dr_val === nothing) was unreachable through the TOML path. As a result dc_type=:rfitzp/:lar/:toroidal silently produced Δ_crit ≡ 0 (bare Δ'>0 criterion, no χ‖ interchange stabilization). Default both to nothing so D_R (and the toroidal geometric factor) are auto-derived from the equilibrium; an explicit scalar still overrides and an explicit 0.0 still disables the offset. TOML parser accepts nothing-or-number; HDF5 writer records the auto setting as NaN. No registered regression case exercises dc_type≠:none, so tracked quantities are unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

matt-pharr · 2026-06-29T14:40:37Z

@d-burg Don't worry about RPEC -- I figured it was probably better to actually just put that into ForceFreeStates so that we can abstract the integrator better. This way ForceFreeStates can always supply a basis of perturbed states

I will review in the coming days

d-burg and others added 30 commits April 19, 2026 02:37

d-burg marked this pull request as ready for review June 17, 2026 20:07

Copilot AI review requested due to automatic review settings June 17, 2026 20:07

d-burg assigned logan-nc and matt-pharr Jun 17, 2026

Copilot started reviewing on behalf of d-burg June 17, 2026 20:08 View session

Copilot AI reviewed Jun 17, 2026

View reviewed changes

Comment thread src/Tearing/InnerLayer/SLAYER/LayerParameters.jl Outdated

Comment thread src/Tearing/InnerLayer/SLAYER/LayerParameters.jl Outdated

Comment thread src/Tearing/Runner/run_slayer.jl

Comment thread src/Equilibrium/DirectEquilibrium.jl Outdated

d-burg and others added 4 commits June 21, 2026 17:14

EXAMPLE - NEW FEATURE - DIIID kinetic profiles + NTV demo; ldp/mpsi=2…

65a9b3e

…56 grid Co-Authored-By: Claude Opus 4.8 <[email protected]>

REGRESSION - Track NTV torque in diiid_n1 case

0efb70b

Co-Authored-By: Claude Opus 4.8 <[email protected]>

TEST - Re-pin DIIID parallel-integration δW/Δ' for ldp grid

cd3969a

Co-Authored-By: Claude Opus 4.8 <[email protected]>

EXAMPLE - Drop verbose KineticForces section comment

ce12c74

Co-Authored-By: Claude Opus 4.8 <[email protected]>

logan-nc added this to the GPEC v2.0.0 milestone Jun 21, 2026

logan-nc removed the WIP Work in progress label Jun 21, 2026

d-burg and others added 13 commits June 22, 2026 15:40

EXAMPLE - Retarget DIIID equilibrium to Ip=1.15 MA; psihigh=0.995 (ca…

8f0a7cc

…ptures q=6) Co-Authored-By: Claude Opus 4.8 <[email protected]>

TEST - Re-pin DIIID parallel-integration for Ip=1.15 equilibrium

1803be9

Co-Authored-By: Claude Opus 4.8 <[email protected]>

KineticForces - CLEANUP - Remove per-evaluation debug prints from NTV…

42a0796

… torque Drop leftover println/printf debug in tpsi! that flooded stdout during the psi-quadrature (see #303). Co-Authored-By: Claude Opus 4.8 <[email protected]>

Merge remote-tracking branch 'origin/feature/diiid-kinetic-profiles' …

edb4d26

…into feature/tearing-growthrates

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Feature: add SLAYER and GGJ tearing growth rates#238

Feature: add SLAYER and GGJ tearing growth rates#238
d-burg wants to merge 68 commits into
developfrom
feature/tearing-growthrates

d-burg commented May 18, 2026

Uh oh!

d-burg commented Jun 17, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

logan-nc commented Jun 18, 2026

Uh oh!

d-burg commented Jun 26, 2026

Uh oh!

matt-pharr commented Jun 29, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Uh oh!

Conversation

d-burg commented May 18, 2026

Summary

Relationship to perf/riccati

Current status — WIP / draft

Uh oh!

d-burg commented Jun 17, 2026

Branch update since 2026-05-18

Commit summary

Net feature state

Reviewer notes

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

logan-nc commented Jun 18, 2026

Uh oh!

d-burg commented Jun 26, 2026

Uh oh!

matt-pharr commented Jun 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Relationship to `perf/riccati`

matt-pharr commented Jun 29, 2026 •

edited

Loading