Feature: add SLAYER and GGJ tearing growth rates

.gitignore

@@ -8,6 +8,8 @@
 *.i
 *.out
 *.h5
+# kinetic-profile inputs are tracked despite the blanket *.h5 (output) ignore
+!examples/**/*kinetic*.h5
 .vscode
 *.DS_Store
 *.pdf
@@ -26,3 +28,9 @@ pestotv
 regression-harness/*.sqlite
 regression-harness/*.sqlite-shm
 regression-harness/*.sqlite-wal
+scratch/
+.claude/worktrees/
+.claude/agent-memory/
+
+# Local profiling scratch (one-off study scripts, not part of the package)
+profiling/

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 [deps]
 AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -36,6 +37,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 AdaptiveArrayPools = "0.3.5"
 Contour = "0.6.3"
+DelaunayTriangulation = "1.6.6"
 DelimitedFiles = "1.9.1"
 DiffEqCallbacks = "4.9.0"
 Documenter = "1.14.1"

docs/make.jl

@@ -35,7 +35,7 @@ makedocs(;
             "KineticForces" => "kinetic_forces.md",
             "Forcing Terms" => "forcing_terms.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md",
-            "Inner Layer" => "inner_layer.md",
+            "Tearing" => "inner_layer.md",
             "Analysis" => "analysis.md",
             "Utilities" => "utilities.md"
        ],

docs/src/inner_layer.md

@@ -1,9 +1,10 @@
-# Inner Layer Module
+# Tearing Module
 
-The InnerLayer module provides abstract scaffolding for resistive inner-layer
-models used in matched asymptotic expansions for resistive MHD stability
-analysis. It currently includes the GGJ (Glasser–Greene–Johnson) shooting
-method for computing the inner-layer response.
+The `Tearing` module groups the resistive tearing-mode analysis stack:
+`InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and
+SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and
+contour-intersection root extraction), and `Runner` (user-facing TOML
+configuration, profile loading, and HDF5 output).
 
 ## InnerLayer
 
@@ -16,3 +17,21 @@ Modules = [GeneralizedPerturbedEquilibrium.InnerLayer]
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ]
 ```
+
+## SLAYER
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER]
+```
+
+## Dispersion
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Dispersion]
+```
+
+## Runner
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Runner]
+```

docs/src/kinetic_forces.md

@@ -4,6 +4,27 @@ Kinetic torque and energy calculations for perturbed equilibria.
 Implements neoclassical toroidal viscosity (NTV) from the PENTRC formulation
 ([Logan & Park, 2013](citations.md); [Logan, 2015](citations.md)).
 
+## Kinetic profile file formats
+
+Kinetic profiles are read by `read_kinetic_file`, which dispatches on the file
+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
+    carries a `units` attribute; the root carries `schema_version` and
+    `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
+    diffusivities m²/s. Write files with `write_kinetic_h5`.
+  - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table
+    `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`,
+when present, are carried for the resistive-layer (SLAYER) analysis and ignored
+here.
+
 ## Profile Scaling Knobs
 
 Seven scaling factors are available on `KineticForcesControl` to modify kinetic

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -0,0 +1,97 @@
+# DIII-D-like SLAYER tearing-mode growth-rate example.
+# Reuses the equilibrium and H-mode kinetic profiles of the sibling
+# DIIID-like_ideal_example (geqdsk referenced by relative path; not
+# duplicated). Runs equilibrium + ForceFreeStates + SLAYER and skips
+# PerturbedEquilibrium (ForceFreeStates.force_termination = true).
+
+[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)
+power_bp = 0                   # Poloidal field power exponent for Jacobian
+power_b = 0                    # Toroidal field power exponent for Jacobian
+power_r = 0                    # Major radius power exponent for Jacobian
+grid_type = "log_asymptotic"   # Radial grid packing type
+psilow = 1e-4                  # Lower limit of normalized flux coordinate
+psihigh = 0.9995               # Upper limit of normalized flux coordinate
+mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
+psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+mtheta = 256                   # Number of poloidal grid points
+newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+etol = 1e-10                   # Error tolerance for equilibrium solver
+force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
+
+[Wall]
+shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.2415                     # Distance from plasma (conformal) or shape parameter
+aw = 0.05                      # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                       # Elongation parameter for wall shapes
+cw = 0                         # Offset of wall center from major radius
+dw = 0.5                       # Triangularity parameter for wall shapes
+tw = 0.05                      # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true          # Equal arc length distribution of nodes on wall
+
+[ForceFreeStates]
+force_termination = true       # Run FFS + SLAYER, skip PerturbedEquilibrium
+local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
+mat_flag = true                # Construct coefficient matrices for diagnostic purposes
+ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
+vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
+
+psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
+qlow = 1.02                    # Integration initiated at q determined by min(q0, qlow)...
+qhigh = 1e3                    # Integration terminated at q limit determined by min(qa, qhigh)...
+sing_start = 0                 # Start integration at the sing_start'th rational from the axis (psilow)
+
+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
+mthvac = 512                   # Number of points used in splines over poloidal angle at the plasma-vacuum interface
+thmax0 = 1                     # Linear multiplier on the automatic choice of theta integration bounds
+
+kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
+kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
+eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations
+save_interval = 3              # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
+singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e4                    # Maximum fraction of solutions allowed before re-normalized
+
+# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details)
+use_parallel          = true   # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix)
+parallel_threads      = 1      # serial/bit-deterministic BVP — keeps the regression Δ' (and hence γ) reproducible
+populate_dense_xi     = false  # No PerturbedEquilibrium here; the dense EL pass has no consumer (auto-disabled under force_termination anyway). SLAYER needs only delta_prime_matrix from the parallel BVP.
+truncate_at_dW_peak   = false  # Edge-dW scan stays diagnostic; integration domain set by qhigh / psihigh / dmlim
+set_psilim_via_dmlim  = true   # TRUE for diverted geqdsks — q → ∞ at separatrix, so dmlim truncation avoids the δW kink instability at negligible domain cost
+dmlim                 = 0.2    # Truncate integration at (last_rational_q + dmlim) / n
+
+[SLAYER]
+# SLAYER tearing-mode growth-rate analysis on the DIII-D-like H-mode equilibrium.
+# Runs in the force_termination path (PE skipped). Kinetic profiles (incl. the
+# χ⊥/χ_φ heat/momentum diffusivities) are read from the standardized HDF5
+# kinetic file shared with the kinetic/NTV physics. Per-surface (uncoupled)
+# analysis: the headline is the unstable 2/1; the validity gate drops surfaces
+# with no real root (e.g. the 5/1, whose Δ' BVP yields a huge uncancellable Δ').
+enabled                 = true
+inner_model             = "slayer_fitzpatrick"
+scan_mode               = "amr"
+coupling_mode           = "uncoupled"
+dc_type                 = "none"
+mu_i                    = 2.0
+zeff                    = 1.0
+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(|Δ|)
+filter_above_poles      = true
+filter_outside_re       = true
+polish_roots            = true       # refine each root to the true zero; enables the validity gate
+store_scan              = false      # per-surface scans omitted from HDF5 (uncoupled has one per surface)
+# Standardized kinetic-profile file (GPEC HDF5 kinetic schema), shared with the
+# sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥), chi_phi (χ_φ).
+profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5"
+
+[SLAYER.scan_grid]
+Q_re_range = [-2.0, 2.0]
+Q_im_range = [-0.5, 3.0]
+nre        = 41
+nim        = 31