ForceFreeStates/Tearing - API! - Pose inner-layer matching and tearing as problems solved on the published result

[matt-pharr]

Release note

• Audience: users
• Numerical impact: default (ray) matching moves ≤ 1.1e-11 vs develop; gal_inner_solver = "galerkin" decks move Δ 6.9e-4, toward ray (harness @ 20887ed)
• Migration: see below

Matching and tearing are now problems solved on a finished force-free-states result.

using GeneralizedPerturbedEquilibrium

eq = PlasmaEquilibrium("TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada", grid_type="ldp",
    psilow=1e-4, psihigh=0.993, mpsi=128, mtheta=256)

# Force-free states; matching needs the coil-response block (rpec_flag)
prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=4, qlow=1.02, set_psilim_via_dmlim=false)
ffs = solve(prob, Galerkin(; pfac=0.03, rpec_flag=true, cut_solution=true))

# Driven matching with the GGJ inner layer (ray backend); scalars apply to every surface
matched = solve(MatchProblem(ffs; eta=8e-8, rho=3.3e-7, rotation=1.0), GGJModel())
matched.bpen                      # penetrated resonant field, surface × coil

# A rotation scan reuses the one outer solve
bpens = [solve(MatchProblem(ffs; eta=8e-8, rho=3.3e-7, rotation=f), GGJModel()).bpen for f in (1.0, 4.0, 16.0)]

# Tearing growth rates with the SLAYER slab layer
tear = solve(TearingProblem(ffs; profile_file="TkMkr_D3Dlike_Hmode_kinetic.h5"), SLAYERModel())
tear.gamma_Hz

• Runs as written in examples/DIIID-like_ideal_example/ (~70 s).
• Per-surface inputs follow ffs.surfaces, core to edge.
• SLAYER does tearing only: slab, tearing parity.
• Riccati also matches: add vac_flag=true; no matched ξ.
• GGJ tearing errors until its γ extraction lands.

What changed

• MatchProblem + GGJModel(): matching after the solve.
• TearingProblem + SLAYERModel(): growth rates.
• Scans reuse one outer solve.
• result.match holds the match and its model, any formalism.
• Match/ is annotated; adds reconnected_flux, inner_model.
• Tearing/ records the model that ran.
• Matched results drop ideal δW; PE skips them.
• InnerLayer model types are the solve argument.
• One GGJ-parameter builder, from stored surface geometry.

Migration

• ResistiveMatch(...) → MatchProblem(ffs; eta, rho, rotation) + GGJModel().
• Runner.run_slayer → TearingProblem(ffs; profile_file) + SLAYERModel().
• build_ggj_inputs → ggj_parameters, with explicit η and ρ.
• result.galerkin.match → result.match.
• gpec.h5: GalerkinIntegration/Match/ → SingularSurfaces/Match/.
• [SLAYER] inner_model = "ggj_*" now errors.
• Decks with gal_inner_solver = "galerkin": set gal_inner_xfac=10, nx=1280, nq=5 for the old grid.
• Matching rejects multi-n, already-matched, and Galerkin solves that dropped a surface.

New defaults

• GGJ matching backend: :ray.
• Galerkin inner grid: backend's 512/4/1 (was 1280/5/10).
• vs ray, 512/4/1 is 26–29× closer in Δ and cout, ~9× in bpen.

Regression report

• 8 cases, develop 278057d96 vs 20887ed8c.
• gal_resistive_diiid (+ new _galinner): match quantities N/A on develop (path rename).
• One-off develop-vs-head, same deck and packages:
  • ray inner: Δ 1.1e-11, cout 1.0e-13, bpen 1.1e-12.
  • Galerkin inner, old grid: Δ 1.5e-7 (inputs agree to 1e-15; backend amplifies).
• diiid_slayer_n1: layer parameters identical; roots within noise.
• gal_resistive_pe: N/A on both refs (Galerkin PE gap).
• Everything else unchanged.

Full report

================================================================
Case: gal_resistive_diiid — DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)
================================================================

Regression Report: gal_resistive_diiid
====================================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
--------------------------------------------------------------------
Quantity                 278057d96     20887ed8c     Diff     Status
--------------------------------------------------------------------
gal # singular surfaces  4             4             0.0e+00  OK
gal singular q values    [4 elem]      [4 elem]      0.0e+00  OK
gal PEST3 Δ diagonal     [4 elem]      [4 elem]      0.0e+00  OK
||gal Δ′ matrix||        1.242969e+04  1.242969e+04  0.0e+00  OK
||gal Δ_coil block||     2.415908e+02  2.415908e+02  0.0e+00  OK
gal D_I per surface      [4 elem]      [4 elem]      0.0e+00  OK
gal α per surface        [4 elem]      [4 elem]      0.0e+00  OK
||gal inner-layer Δ||    N/A           1.474569e+06  N/A      N/A
||gal match cout||       N/A           3.062664e-02  N/A      N/A
||gal match bpen||       N/A           3.748382e+00  N/A      N/A
gal match residual       N/A           1.806146e-16  N/A      N/A
Runtime (s)              833.7s        137.0s                 --
====================================================================
Summary: 7 unchanged, 4 missing/N/A


================================================================
Case: gal_resistive_diiid_galinner — DIII-D-like, n=1, Galerkin Δ′ matched with the Galerkin inner-layer backend
================================================================

Regression Report: gal_resistive_diiid_galinner
====================================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
--------------------------------------------------------------------
Quantity                 278057d96     20887ed8c     Diff     Status
--------------------------------------------------------------------
gal # singular surfaces  4             4             0.0e+00  OK
gal singular q values    [4 elem]      [4 elem]      0.0e+00  OK
gal PEST3 Δ diagonal     [4 elem]      [4 elem]      0.0e+00  OK
||gal Δ′ matrix||        1.242969e+04  1.242969e+04  0.0e+00  OK
||gal Δ_coil block||     2.415908e+02  2.415908e+02  0.0e+00  OK
gal D_I per surface      [4 elem]      [4 elem]      0.0e+00  OK
gal α per surface        [4 elem]      [4 elem]      0.0e+00  OK
||gal inner-layer Δ||    N/A           1.474558e+06  N/A      N/A
||gal match cout||       N/A           3.062664e-02  N/A      N/A
||gal match bpen||       N/A           3.748382e+00  N/A      N/A
gal match residual       N/A           1.688777e-16  N/A      N/A
Runtime (s)              156.2s        129.0s                 --
====================================================================
Summary: 7 unchanged, 4 missing/N/A


================================================================
Case: gal_resistive_pe — DIII-D-like, n=1, DRIVEN/RPEC: gal-matched resistive ξ → PerturbedEquilibrium (coil-driven singular coupling)
================================================================

Regression Report: gal_resistive_pe
==========================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
----------------------------------------------------------
Quantity                278057d96  20887ed8c  Diff  Status
----------------------------------------------------------
rational q values       N/A        N/A        N/A   N/A
island half-widths      N/A        N/A        N/A   N/A
Chirikov parameter      N/A        N/A        N/A   N/A
||resonant b field||    N/A        N/A        N/A   N/A
resonant b field (all)  N/A        N/A        N/A   N/A
penetrated b field      N/A        N/A        N/A   N/A
PE Δ' per surface       N/A        N/A        N/A   N/A
||C resonant b field||  N/A        N/A        N/A   N/A
Runtime (s)             154.3s     149.8s           --
==========================================================
Summary: 8 missing/N/A


================================================================
Case: ggj_reference — GGJ inner-layer Galerkin solver, Glasser & Wang 2020 Eq. 55, γ = 1+i
================================================================

Regression Report: ggj_reference
==========================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
----------------------------------------------------------
Quantity     278057d96      20887ed8c      Diff     Status
----------------------------------------------------------
Re(Δ_odd)    -1.901769e-01  -1.901769e-01  0.0e+00  OK
Im(Δ_odd)    -4.159478e-01  -4.159478e-01  0.0e+00  OK
Re(Δ_even)   -3.892003e-01  -3.892003e-01  0.0e+00  OK
Im(Δ_even)   -3.606845e-01  -3.606845e-01  0.0e+00  OK
Runtime (s)  5.6s           5.6s                    --
==========================================================
Summary: 4 unchanged


================================================================
Case: ggj_ray_q500i — GGJ inner-layer :ray backend, q=4 surface, Q = 500i (beyond :galerkin)
================================================================

Regression Report: ggj_ray_q500i
========================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
--------------------------------------------------------
Quantity     278057d96     20887ed8c     Diff     Status
--------------------------------------------------------
Re(Δ_odd)    2.472008e+00  2.472008e+00  0.0e+00  OK
Im(Δ_odd)    1.335402e+01  1.335402e+01  0.0e+00  OK
Re(Δ_even)   1.374966e-01  1.374966e-01  0.0e+00  OK
Im(Δ_even)   7.427542e-01  7.427542e-01  0.0e+00  OK
Runtime (s)  13.5s         14.1s                  --
========================================================
Summary: 4 unchanged


================================================================
Case: solovev_kinetic_ntv — Solovev analytical equilibrium, n=1, ideal + PE + NTV torque quadrature
================================================================

Regression Report: solovev_kinetic_ntv
======================================================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
--------------------------------------------------------------------------------------
Quantity                                   278057d96     20887ed8c     Diff     Status
--------------------------------------------------------------------------------------
NTV torque fgar [Re, Im]                   [1 elem]      [1 elem]      0.0e+00  OK
NTV ψ quadrature evaluations               840           840           0.0e+00  OK
root-area-weighted total energy Re(et[1])  4.419680e-01  4.419680e-01  0.0e+00  OK
# singular surfaces                        2             2             0.0e+00  OK
singular psi locations                     [2 elem]      [2 elem]      0.0e+00  OK
q0                                         1.900006e+00  1.900006e+00  0.0e+00  OK
Runtime (s)                                190.2s        191.2s                 --
======================================================================================
Summary: 6 unchanged


================================================================
Case: diiid_error_field — DIII-D-like error-field assessment, n=1: ErrorFields stage with the NTV-limited correction couplings
================================================================

Regression Report: diiid_error_field
============================================================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
--------------------------------------------------------------------------------------------
Quantity                                         278057d96     20887ed8c     Diff     Status
--------------------------------------------------------------------------------------------
EFC overlap per kAt                              [1 elem]      [1 elem]      0.0e+00  OK
EFC whole-field NTV torque per kAt² [N·m]        [1 elem]      [1 elem]      0.0e+00  OK
EFC residual-field NTV torque per kAt² [N·m]     [1 elem]      [1 elem]      0.0e+00  OK
EFC reference rotation ω_ref [rad/s]             [1 elem]      [1 elem]      0.0e+00  OK
EFC rotation-scan shifts [rad/s]                 [9 elem]      [9 elem]      0.0e+00  OK
EFC residual torque against rotation [N·m/kAt²]  [9 elem]      [9 elem]      0.0e+00  OK
dominant-coupling singular values                [3 elem]      [3 elem]      0.0e+00  OK
|delta_nominal| of coil set 1                    1.407189e-06  1.407189e-06  0.0e+00  OK
|Σ δ_as_designed| of the error-field coils       1.502909e-18  1.502909e-18  0.0e+00  OK
Monte Carlo mean |δ| of the error-field coils    4.701483e-05  4.701483e-05  0.0e+00  OK
Runtime (s)                                      830.5s        627.9s                 --
============================================================================================
Summary: 10 unchanged


================================================================
Case: diiid_slayer_n1 — DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis (uncoupled per-surface, AMR, validity-gated)
================================================================

Regression Report: diiid_slayer_n1
=========================================================================
Ref 1: 278057d96  @ 278057d96 (2026-10-15)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
Ref 2: 20887ed8c  @ 20887ed8c (2026-10-16)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest 3ec0d96d (pinned), 6 threads/6 BLAS, jit build
-------------------------------------------------------------------------
Quantity                            278057d96  20887ed8c  Diff     Status
-------------------------------------------------------------------------
SLAYER surface indices              [6 elem]   [6 elem]   0.0e+00  OK
SLAYER poloidal m                   [6 elem]   [6 elem]   0.0e+00  OK
SLAYER toroidal n                   [6 elem]   [6 elem]   0.0e+00  OK
SLAYER minor radius rs              [6 elem]   [6 elem]   0.0e+00  OK
SLAYER r-based shear                [6 elem]   [6 elem]   0.0e+00  OK
SLAYER Lundquist S                  [6 elem]   [6 elem]   0.0e+00  OK
SLAYER D_norm                       [6 elem]   [6 elem]   0.0e+00  OK
SLAYER P_perp                       [6 elem]   [6 elem]   0.0e+00  OK
SLAYER tauk                         [6 elem]   [6 elem]   0.0e+00  OK
SLAYER iota_e                       [6 elem]   [6 elem]   0.0e+00  OK
SLAYER Q_e                          [6 elem]   [6 elem]   0.0e+00  OK
SLAYER Q_i                          [6 elem]   [6 elem]   0.0e+00  OK
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   1.2e-05  OK
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   3.4e-01  OK
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   8.4e-02  OK
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00  OK
SLAYER enabled flag                 1          1          0.0e+00  OK
Runtime (s)                         227.8s     186.9s              --
=========================================================================
Summary: 17 unchanged

Reviewing this PR

• Review the full diff, not commit by commit.
• Turn on "Hide whitespace": _compute_match re-indents a block.
• Moved (check only the listed edits):
  • gal_match_rpec → _compute_match: ctrl.gal_* → prob.*, the backend if → solve_inner_profile(model, …).
  • Its 4·msing system → _match_system, verbatim; a dead duplicate in resonant_match_rpec is deleted.
  • GalMatchResult → MatchResult, plus model.
  • write_galerkin!'s Match/ block → write_match!, plus annotations.
  • run_slayer → solve(::TearingProblem) + _tearing_params, verbatim (develop's Z_eff included).
  • build_ggj_inputs τ_A/τ_R → ggj_parameters.
  • resist_eval deleted: a second port of resist_geometry.
• New (where attention belongs):
  • MatchProblem, _matched_result/_with, the δW drop.
  • TearingProblem, and the GGJ error.
  • GGJModel options and their forwarding (InnerLayer/GGJ/GGJ.jl).
  • Deck match block in run_force_free_states; new defaults.
  • Tests: runtests_matching_models.jl, new runtests_solve_api.jl testsets.

Notes for reviewers

• src/ is −260 lines net vs develop.
• Deck outputs byte-identical commit to commit, except:
  • default matching moves at round-off (above);
  • SLAYER roots scatter run to run (known).
• Kinetic profiles on the equilibrium: moved to the next PR.
• GGJ tearing γ extraction is future work.
• Matched PE needs resistive δW: future physics work.
• Fortran gpec does the same: drops ideal solutions.
• Riccati vs Galerkin bpen differ near the edge.
• Galerkin Δ′/Δ_coil edge instability: validate vs Fortran.
• |bpen| rising at low rotation is known GGJ behavior.
• Develop merged three times; conflicts in CoreTypes.jl and Tearing (Z_eff).
• REFACTOR_PLAN.md leaves the repo (−1117 lines).

🤖 Generated with Claude Code

https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf

[matt-pharr]

CC @jhalpern30 since you were party to the earlier steps of this refactoring

[matt-pharr]

Also CC @StuartBenjamin for the new syntax which will be particularly useful for tearing suite

[matt-pharr]

@d-burg please prioritize this as I want the shape of this to be final before I go through and refactor GeneralizedPerturbedEquilibrium.jl and this is desired to be in for v2.0.0

[d-burg]

Thanks for this, the direction is right overall and the moved matching logic is equivalent to the deleted code (signs, transposes, parity and the bpen formula all match by my eye). Before the shape is frozen for v2 it'd be great to have these settled:

1. The inner layer run must be included in the result. The model is now the solve argument, but SLAYERControl still carries inner_model and the result records the control, so the two can disagree. MatchResult also records no model, backend or options.
2. Regarding GGJ tearing: either the GGJ γ are marked unmistakably as unvalidated on the object and in the file, or solve(TearingProblem, GGJModel()) should error until extraction is implemented. I have local scripts for this that work robustly, but need to be cleaned up and probably need their own PR.
3. From claude: '"Profiles live on the equilibrium” currently means the NTV container (SI units, 101 points, and no χ). The tearing layer needs eV, the file’s own grid and χ. The bridge is where the unit and χ differences below come from.'
4. “Round-off” probably needs qualifying -- the harness shows the three match quantities as N/A on develop, and decks on the Galerkin inner backend change grid.
5. #463 and #501 change code in the two run_slayer methods this PR deletes, so git cannot carry their hunks. I’d prefer those two land first and this rebases onto them.

Claude also raised some lower priority questions:

• Attaching profiles before an auto-grid solve now changes the ψ grid on the API path.
• A deck with [KineticForces] and [SLAYER] but no profile_file used to fail the SLAYER stage and now runs.
• gamma_Hz/omega_Hz versus the HDF5 units (1/s and rad/s) for the same arrays.
• TearingProblem accepts multi-n results; MatchProblem rejects them.
• The default rotation is ω_E from one profile container and toroidal rotation from the other.
• Riccati matching is effectively new in this PR and is known to differ from Galerkin near the edge (is this correct?)

[matt-pharr]

Thank you for the thoughtful review @d-burg -- are #463 and #501 things that we can get wrapped up in the next day or two if I can help with the reviews, or are those longer outlook? If the latter, are they at least going to make it in with a week+ to spare before APS? I have more refactoring to do that will depend on this branch's changes.

[matt-pharr]

Thanks Daniel — addressed in d1754db28, c4682c696, 76caa9ee2, 24b03813c, plus a develop merge (20887ed8c). PR body restamped.

1. Model on the result. MatchResult.model and SLAYERResult.model hold the model actually solved (with its options). Match/ and Tearing/ record it as an inner_model attribute (Match/ also inner_model_options). SLAYERControl.inner_model stays only as the deck's request.

2. GGJ tearing. solve(TearingProblem(ffs), GGJModel()) now errors; your extraction PR can replace the error.

3. Profiles on the equilibrium. Pulled out of this PR entirely. Tearing reads profile_file through _load_profiles again, as on develop (now with develop's Z_eff). Matching takes explicit eta/rho/rotation. The next PR stores the file's KineticProfileData on the equilibrium and has each consumer derive its own view. It will also cover the unit constant, χ, Z_eff, the η default (Spitzer until Utilities - BUGFIX! - Compute the trapped-particle fraction from the pitch-angle integral #497), and the profile_file-vs-attached equivalence test.

4. Round-off. Develop-vs-head one-off on gal_resistive_diiid (old path read on develop, same packages):

   • ray inner: Δ 1.1e-11, cout 1.0e-13, bpen 1.1e-12;
   • Galerkin inner, old grid: Δ 1.5e-7. Inputs agree to 1e-15; the backend amplifies them.

   The Galerkin inner default grid change moves those decks' Δ by 6.9e-4, toward ray: 26–29× closer in Δ and cout, ~9× in bpen. bpen is now tracked, and a gal_resistive_diiid_galinner case covers the Galerkin inner backend. The release note says this.

Lower-priority questions:

• Auto-grid with attached profiles: gone with item 3; the next PR documents it.
• [KineticForces] + [SLAYER] without profile_file: fails the SLAYER stage again, as on develop.
• gamma_Hz/omega_Hz naming: inherited from develop's SLAYERResult; the HDF5 units are right. Happy to rename in your extraction PR.
• Multi-n: TearingProblem accepts it, as run_slayer did; matching is single-n because the Galerkin Δ′ is.
• Rotation containers: gone. Matching takes rotation explicitly, and tearing uses the file's ω_E as on develop. Unifying it (including ExB_rotation_factor) is in the next PR.
• Riccati matching: correct. It's new here and differs from Galerkin near the edge. Galerkin is the one with the edge instability; both need validating against Fortran rmatch (whose Riccati path has no Δ_coil).

[d-burg]

Thanks Matt!