Merge pull request #238 from OpenFUSIONToolkit/feature/tearing-growthrates

Feature: add SLAYER and GGJ tearing growth rates

Author: d-burg

.gitignore

@@ -28,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

@@ -36,7 +36,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/developer_notes.md

@@ -20,6 +20,98 @@ Where CODE is the module name (EQUIL, ForceFreeStates, VAC, PERTURBED EQUILIBRIU
 
 The regression harness **must be run on every pull request before merging into `develop`**. It is the project's primary safeguard for tracking how numerical results evolve across changes, so it is only useful if every PR exercises it. When you open a PR, paste the regression report into the PR thread so reviewers can see what moved (and what did not). If your change touches a quantity that is not yet tracked, add a new regression case — or extend an existing one — in the same PR.
 
+### Open problem: pinning grid-sensitive Δ′ robustly
+
+The ideal-MHD Δ′ values pinned in `test/runtests_parallel_integration.jl`
+(`delta_prime_matrix` diagonal) and tracked by the harness are **single-point
+snapshots**: one value at one `psi_accuracy` on one grid. They exist to catch
+unintended changes, and are explicitly *not* converged Δ′. The extraction is
+intrinsically grid-sensitive — a `psi_accuracy` scan from 2e-3 to 2.5e-4 swings
+the DIII-D-like `dpm[1,1]` by roughly 50% (about 6.2 to 9.9), and switching the
+grid generator moves it again. Any such pin therefore encodes an arbitrary point
+on a varying curve, and re-pinning is required whenever the grid changes, which
+weakens it as a regression signal.
+
+A more defensible criterion is to pin the **plateau**: scan Δ′ across
+`psi_accuracy` and the integration-truncation controls, and take the value where
+the result is stationary — the mode of the resulting distribution rather than any
+single sample. Where a plateau exists it is a property of the physics rather than
+of the discretization, so it would survive grid changes and would be a genuine
+convergence statement.
+
+A `psi_accuracy` scan on the DIII-D-like SLAYER deck (2e-3 down to 3.125e-5, a
+64x range) shows what a plateau search is up against:
+
+| psi_accuracy | knots | implied | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) |
+|---|---|---|---|---|---|---|
+| 2.0e-3   | 172 | –    | 6.850 | -5.783 | -16.308 | 165.3 |
+| 1.0e-3   | 205 | –    | 6.387 | -5.873 | -16.867 | 154.3 |
+| 5.0e-4   | 252 | –    | 8.003 | -6.105 | -16.413 | 192.9 |
+| 2.5e-4   | 307 | 522  | 8.567 | -6.235 | -16.529 | 206.2 |
+| 1.25e-4  | 372 | 666  | 8.146 | -6.285 | -16.649 | 196.2 |
+| 6.25e-5  | 457 | 916  | 8.911 | -6.477 | -16.510 | 214.4 |
+| 3.125e-5 | 559 | 1226 | 9.514 | -6.260 | -16.485 | 228.7 |
+
+Three things follow. First, convergence is per-surface: the outer surface
+`dpm[3,3]` (q=4) *does* plateau under the auto grid — the last four points agree
+to about 1% — while `dpm[2,2]` is marginal and `dpm[1,1]` (q=2) never settles,
+still moving ~7% between the two tightest grids. A plateau detector must
+therefore report per-surface rather than pass/fail for the whole diagonal.
+(As the `ldp` scan below shows, q=2 is not inherently unconvergeable — it is the
+auto grid that prevents it from settling.)
+
+Second, the growth rate is linear in Δ′: `gamma/dpm[1,1]` is 24.1 to within 0.5%
+across the whole scan, so a converged Δ′ is both necessary and sufficient for a
+converged SLAYER growth rate on this surface.
+
+Third — and this is the blocker — the `implied` column (what
+`implied_knot_count` requests after the refined solve) grows monotonically
+relative to the grid actually used, from 1.7x to 2.2x. The two-pass scheme stops
+after pass 2, so tightening `psi_accuracy` moves it *further* from
+self-consistency rather than closer, and the warning's advice to "consider
+tightening psi_accuracy" is counterproductive in this regime.
+
+The same deck on a deterministic `ldp` grid, which skips the measure-and-re-form
+step entirely, converges:
+
+| mpsi | dpm[1,1] | dpm[2,2] | dpm[3,3] | gamma 2/1 (Hz) |
+|---|---|---|---|---|
+| 128  | 5.208 | -8.666 | -19.157 | 126.2 |
+| 256  | 7.097 | -5.187 | -15.729 | 171.2 |
+| 512  | 8.489 | -6.217 | -16.517 | 204.4 |
+| 1024 | 8.932 | -6.298 | -16.529 | 214.9 |
+| 2048 | 8.914 | -6.385 | -16.468 | 214.5 |
+
+Successive `dpm[1,1]` increments are +1.889, +1.392, +0.443, -0.018: the last
+doubling moves both Δ′ and the growth rate by 0.2%. So the q=2 Δ′ is **not**
+intrinsically ill-conditioned — it converges to ≈8.91, with gamma ≈214.5 Hz, and
+`dpm[3,3]` converges to ≈-16.5 on *both* grid families, which is what a physical
+value should do. The non-convergence under the auto grid is an artifact of the
+generator, not of the Δ′ extraction.
+
+Two consequences. A plateau criterion is implementable today against a fixed
+`ldp` grid, without waiting on the auto-grid work. And the auto grid's answers
+are biased in both directions relative to the converged value: at its default
+`psi_accuracy` it gave 6.39 (28% low), at its tightest 9.51 (7% high). Anything
+pinned on the auto grid should be read with that in mind.
+
+This is not implemented. Doing it properly needs:
+
+  - either a fixed `ldp` grid (which already converges, see above) or knot
+    refinement iterated to a fixed point (repeat the measure-and-re-form
+    step until `implied_knot_count` stops exceeding the grid in use), since
+    without it the scan target keeps moving;
+  - a scan driver over `psi_accuracy` × truncation (`psihigh`, `dmlim`, `qhigh`)
+    that records the Δ′ diagonal per configuration;
+  - a plateau/mode detector with an explicit stationarity tolerance, plus a
+    defined failure mode for surfaces where no plateau exists (the
+    near-separatrix surfaces are expected to fall in this class);
+  - a harness case that pins plateau values and their scan width, replacing the
+    single-point pins.
+
+Until then, treat the pinned diagonal as an order-of-magnitude and sign
+diagnostic only, and expect to re-pin it whenever the equilibrium grid changes.
+
 ### Regression Harness: Quick usage guide
 
 Set up an alias for convenience (optional):

docs/src/inner_layer.md

@@ -1,9 +1,18 @@
-# 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).
+
+## Layer Inputs
+
+Equilibrium/ForceFreeStates glue that assembles per-surface inner-layer inputs.
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Tearing]
+```
 
 ## InnerLayer
 
@@ -16,3 +25,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/utilities.md

@@ -14,6 +14,21 @@ The Utilities module currently provides:
 Modules = [GeneralizedPerturbedEquilibrium.Utilities, GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms]
 ```
 
+## Physical Constants
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Utilities.PhysicalConstants]
+```
+
+## Neoclassical Resistivity
+
+Parallel-resistivity closures (Spitzer, Spitzer-Härm, and the Sauter and Redl
+neoclassical models) used to set the Lundquist number in the tearing stack.
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Utilities.NeoclassicalResistivity]
+```
+
 ## IMAS Output
 
 ```@docs

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -0,0 +1,96 @@
+# 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
+
+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

examples/Solovev_ideal_example/gpec.toml

@@ -44,6 +44,11 @@ verbose = true                          # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
 reg_spot = 0.05                         # Regularization width for singular surfaces (0 = disabled)
 
+# Note: SLAYER tearing analysis is not run on the Solovev analytic equilibrium
+# — its Δ' / inner-layer dispersion does not yield meaningful tearing roots.
+# The SLAYER example and regression case use the DIII-D-like geqdsk instead
+# (see examples/DIIID-like_SLAYER_example).
+
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes