[matt-pharr]

This is the GPEC v2.0.0 consolidation. It reshapes the scripting API, the deck and the HDF5 output around one rule: a deck is a script written as TOML. Every object a script builds has one deck section named after it, holding exactly that constructor's keywords, and main runs the same public calls a script would. The physics is unchanged: across the 25 regression cases, every tracked quantity that both refs compute is identical (see Regression report).

It is a clean break. Nothing is kept for backward compatibility: no aliases, no key mapping, no accepted-and-ignored arguments. An old deck stops on its first moved or removed key, and the error says where the key lives now. The migration tables below list every change.

Closes #222
Closes #372
Closes #374
Closes #388

The plotting half of #326 moves to its own PR, stacked on this one. #613 (a run with a wall mixes a with-wall Λ and a no-wall L) is handled in #620, stacked on this one.

Contents

1. How to review this PR
2. Design rules
3. Changes by area
   • Driver and deck: main, Deck.jl, [Run], replay from gpec.h5
   • Output: save_to_hdf5, inputs recorded as deck sections
   • Equilibrium: analytic sections, kinetic data on the equilibrium
   • Vacuum: optional, typed walls
   • ForceFreeStates: ForceBalanceProblem, integrator objects, kinetic sources
   • LocalStability
   • InnerLayer and Tearing: one growth-rate convention, TearingProblem, GGJ tearing
   • PerturbedEquilibrium
   • KineticForces: kinetic_torque
   • ErrorFields: error_fields
   • Cleanup: dead code, duplicates, docstrings
   • Tests, examples, harness and docs
4. Migration: deck keys, API names, HDF5 paths
5. Release note
6. Regression report
7. Not in this PR

How to review this PR

The diff is large: about 280 files, +8.2k / −12.4k lines. It was built one change per commit, so reading it commit by commit is much easier than reading the full diff. Each area section below names its commits.

Every commit moves no numbers except 6c42f0364 (the growth-rate convention, a root shift of at most about 1e-5 relative, below every harness threshold).

The twelve <Module> - DOCS - Trim docstrings and comments ... commits (ddf311c6d to 9a2218338) change only docstrings and comments. A script that strips comments and docstrings confirmed every file is code-identical before and after, apart from one docstring moved onto the function it describes. These can be skimmed.

Suggested split by reviewer:

Reviewer	Areas
@logan-nc	Driver and deck, Output, ForceFreeStates, PerturbedEquilibrium, KineticForces, ErrorFields
@jhalpern30	Equilibrium, Vacuum
@d-burg	InnerLayer and Tearing, LocalStability

Design rules

These decide where every input lives, and the area sections apply them.

• solve(problem, algorithm) separates what from how. An input that can change the converged answer goes on the problem. An input that only changes accuracy or cost (method, tolerance, resolution, scan range) goes on the algorithm. solve is used only where there is a real choice of method:

  • the integrators: Forward(), Riccati(), Galerkin();
  • the inner-layer models: GGJModel(), SLAYERModel();
  • local stability: Ballooning().

  Everything else is a plain function that returns a result.

• The deck has the same shape as the API. Problems and stage functions get sections such as [ForceFreeStates], [Match], [Tearing] and [PerturbedEquilibrium]. Algorithms and models get sections named after their type, such as [Riccati], [Galerkin], [GGJ], [SLAYER] and [Ballooning]. A problem with a choice of method names it (integrator = "Galerkin", inner_model = "GGJ"). Including a section runs that stage; there are no enabled flags. Unknown keys stop the run.

• One save function. save_to_hdf5(target, result) has a method per result type. Computation never writes files. main saves after each stage, and a script saves what it wants.

• Every capability works from both a deck and the API, with a deck/API parity test (docs/src/developer_notes.md).

Changes by area

Driver and deck

main used to be a long monolithic driver. It is now a short script over the public API: for each section present in the deck, it builds the object, calls the same function a script would, and saves the result.

• src/Deck.jl (new) is the only place deck text becomes API objects. parse_deck turns every section into constructor keywords before anything runs. So a misspelled key, a moved key, or a section no stage reads stops the run at once with an error naming the fix. Moved keys say where they went (for example "gal_nx moved to [Galerkin] as nx").
• [Run] (new) holds the one output_file and the one verbose for the whole run. The per-section copies (HDF5_filename, output_filename, write_outputs_to_HDF5, verbose) are gone. A deck run always writes its file.
• Stopping early works by leaving out the later sections; force_termination is gone.
• Relative paths are resolved against the deck's directory once, in the parser. That covers the equilibrium, kinetic, forcing, coil and wall files. dir_path is gone from the problem and result types.
• The return value is one field per stage, nothing for stages that did not run: (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing). The old slayer, coil_sensitivities, monte_carlo, locking_risk and efc_couplings fields are now tearing and error_fields.
• Replay from gpec.h5: a parsed deck can be written back to the deck text that states it (deck_sections), the inverse of the parser. Rerun.jl replays any gpec.h5, including one written by a script, because inputs are recorded as deck sections (see Output).

Commits: cb523a196, fb53c7343, a0232cc79, 711a9945f, b92ec6b3e, d5b0ef1dc, and the review-sweep fix 9ecffadc6.

Output

• save_to_hdf5(target, result) is declared in Utilities (src/Utilities/Output.jl), and each module adds the method for its own result type. target is a path or an open file or group. The old driver writer write_outputs_to_HDF5 is gone.
• Inputs are recorded as deck sections. Each result keeps the problem and algorithm that produced it, and saving it writes them under Input/ in deck form. A script run and a deck run of the same workflow give the same gpec.h5, apart from Info/Runtimes, and either can be replayed.
• Per-stage writes: main writes each stage's results as soon as the stage finishes, so a crash later in the run keeps the earlier output.
• Dataset paths are unchanged except for the few listed under HDF5 changes.

Commits: c6536ef24, d5b0ef1dc, 36adc5bc1.

Equilibrium

• Analytic equilibria have their own sections: [Solovev], [LargeAspectRatio] and [TJAnalytic]. The section's presence selects the source, and the deck then sets no eq_type. [TJAnalytic] solver = "inverse" | "direct" replaces eq_type = "tj_analytic_direct". The analytic configs are immutable.
• The plasma's composition and transport defaults live on the kinetic data (equil.kinetic, a KineticProfileData), so every consumer reads one source. Previously SLAYER, KineticForces and matching each carried their own copies. These keys are now [Equilibrium] keys:
  • species: [[Equilibrium.ion_species]], impurity_charge, impurity_mass;
  • transport defaults and resistivity closure: zeff_default, chi_e_default, chi_i_default, chi_phi_default, tau_E, resistivity_model, lnLambda_form.
• Removed: the grid_type = "log_asymptotic" alias, equilibrium_solver_by_inversion's refine/β_r/β_z keywords, and the never-read EquilibriumParameters.b_norm and psi_axis_norm.

Commits: fae914299, daaea0236.

Vacuum

• Walls are optional and typed. Leaving out [Wall] means no wall. [Wall] shape = ... selects a WallSettings subtype (ConformalWall, EllipticalWall, DeeWall, ModDeeWall, FileWall), and each takes only its own parameters. WallShapeSettings and shape = "nowall" are gone.

Commit: 14a374549.

ForceFreeStates

• EulerLagrangeProblem is renamed ForceBalanceProblem. With kinetic terms the force operator is not self-adjoint, so the equation is no longer an Euler–Lagrange equation. Internal file and function names are unchanged.
• Integrator numerics live on the algorithm objects: Forward(; eulerlagrange_tolerance, renorm_threshold, numsteps_init, numunorms_init, save_interval). Riccati takes the same plus nchunks and extended_precision_bvp. Galerkin takes the old gal_* keys without the prefix. Physics and domain settings stay on the problem, and the deck selects the integrator with integrator = "Forward" | "Riccati" | "Galerkin".
• nn = 1 or nn = [1, 2, 3] (consecutive) replaces nn_low/nn_high, matching the API's Int or UnitRange.
• The kinetic source is chosen by type with the kinetic keyword: nothing (ideal), a KineticForcesControl (computed from equil.kinetic), a KineticMatrices (precomputed), or TestKineticMatrices(; scale) (synthetic, for tests). In a deck that is [ForceFreeStates] kinetic = "KineticForces". kinetic_factor is now only a scale (default 1), not an on/off switch, and kinetic_source is gone.
• kinetic_matrices(prob) evaluates the kinetic matrices once and returns a KineticMatrices. Passing it back as kinetic= reuses it across solves, for example in a kinetic_factor scan. The internal FKG type is renamed FKGMatrices.
• One frobenius_order on the problem (default 6) is used by all three integrators. gal_sing_order is gone, and [Galerkin] frobenius_order_ceiling keeps Galerkin's per-surface raise.
• Galerkin factorization is automatic (banded Cholesky, with LU fallback), so gal_solver is gone.
• MatchProblem defaults eta, rho and rotation from the equilibrium's kinetic data, through the same layer closures the tearing solve uses. With no kinetic data, the error names exactly which inputs are missing. This closes Add automatic computation of spitzer and neoclassical resistivities with profiles #222.
• psilim_from_layer_overlap works from scripts as well as decks. true uses a default SLAYERModel, and passing a SLAYERModel sets the closures. The scan runs, and ForceFreeStates/LayerOverlap is written, only when the cap is on.
• Removed: sing_start (its non-zero branches were never implemented), sing_get_dua, the allocating sing_matvec/sing_matmul wrappers, and never-read fields on OdeState, SingAsymptotics, FKGMatrices and ResistGeometry.

Commits: 4f9e4d6fa, de7d87d60, 924109e43, 464f5304f, b8f50de22, b63a64901.

LocalStability

• theta_k and alpha_boundary_scan move from Ballooning onto LocalStabilityProblem. Ballooning keeps only numerics.
• verbose is a solve keyword, defaulting to false.
• In a deck, [LocalStability] takes the problem keys and [Ballooning] the numerics.

Commit: 0671b86b0.

InnerLayer and Tearing

• One growth-rate convention. solve_inner(model, p, s) takes s = γ − iω in 1/s for every inner-layer model. layer_response(model, p, z) evaluates a model in its own scan variable (GGJ: s; SLAYER: Q).

  • SLAYER no longer conjugates its Δ. Its kernel works in ĝ = s·τ_k with ĝ − i·Q_e.
  • zero_winding is removed, and every *Winding count is +1 where it was −1.
  • Roots shift by at most about 1e-5 relative. The outer Δ′ diagonal has a ~1e-9 to 1e-7 imaginary part, so solving Δ = Δ′ rather than conj(Δ) = Δ′ moves deep roots slightly.

  6c42f0364.

• Tearing.Runner is folded into Tearing. The entry point is solve(TearingProblem(ffs; ...), model; config=TearingConfig(...)), and the result is a TearingResult. GGJ tearing growth rates now go through the same problem. 643a78df3 moves each GGJ layer into its surface's E×B frame, so GGJ roots are in the lab frame like SLAYER's.

• SLAYER's layer closures live on SLAYERModel: P_perp_model, D_c_type, D_R, D_geo, chi_perp_Dc, theta_sample and variant. Physics shared with other modules comes from equil.kinetic.

• The deck sections are reorganized:

  • [Tearing] holds the problem: inner_model, coupling_mode, msing_max, omega_E_kHz.
  • [Tearing.search] holds the root search (TearingConfig).
  • [Tearing.CriticalResonantField] holds the critical resonant field settings.
  • [SLAYER] and [GGJ] hold the models.

• Units in names are fixed. omega_Hz/gamma_Hz held rad/s and 1/s, so they are now omega/gamma. gap_kHz_threshold compared a rad/s value against kHz; it is now gamma_gap_threshold in 1/s. This closes SLAYER gap_kHz_threshold compares a rad/s quantity against a kHz threshold (2π error) #372.

• ggj_growth_rates(; omega_range) bounds ω, the same as TearingConfig and the result.

• Shared code: the GGJ ggj_rows are shared by _physical_uv and Ray.ode_matrix. Galerkin and GGJ share one Hermite FEM kernel (Utilities.HermiteFEM).

Commits: 345d1d107, 58cadc46e, 643a78df3, 6c42f0364, 8e13504e2.

PerturbedEquilibrium

• Response fields that were not reconstructed are left out of gpec.h5, rather than written as empty placeholders.
• Tearing/PerSurface/lu is renamed S (the Lundquist number).
• The metric labels now say g_ij/J, and b_psi_modes is labeled J·b^ψ. These were mislabeled before; the data are unchanged.
• regularization_width is a perturbed_equilibrium keyword. The PlasmaResponse field that was never set is gone.

Together with the HDF5 changes, this closes #388. Commit: 36adc5bc1.

KineticForces

• kinetic_torque(pe, ctrl) evaluates the NTV torque on a perturbed equilibrium and returns a KineticForcesResult, saved with save_to_hdf5.
• [KineticForces] nn is gone: n comes from the perturbed equilibrium. The species keys moved to [Equilibrium].

Commit: 6375f8874.

ErrorFields

• error_fields(resp, coil_sets, coil_config; ...) runs the whole assessment on a plasma response and returns one ErrorFieldsResult. That holds the sensitivities, dominant coupling, tolerances, Monte Carlo, locking risk, tolerance scan and correction-coil couplings, each nothing when not requested. One save_to_hdf5 writes them all at their existing paths.
• efc_couplings moves into ErrorFields, and run_error_fields is gone.

Commit: 7d19edf2d.

Cleanup

• Dead code removed: helpers and fields with no remaining reader (b63a64901, 54c2bbb6d).
• Duplicates removed: duplicated kernels left by the consolidation (8e13504e2). Each replacement was verified bit-identical. That covers the Hermite FEM basis and grid packing, the GGJ ODE rows, the Riccati renormalization, and the CHEASE reader.
• Docstrings and comments trimmed across src/ (ddf311c6d to 9a2218338, net −3.5k lines). Each docstring now says what the function does, gives the paper reference, and lists its fields, arguments and return value. Walkthroughs that restated the code, history notes and Fortran line references are gone. Docstrings found wrong along the way are fixed.
• Review-sweep fixes: a ten-agent review of the full diff found a few oversights, fixed in 9ecffadc6, 69325f463 and the test commit that follows them:
  • the deck now rejects sections no stage reads;
  • Forward lost an nchunks field that nothing read;
  • a broken benchmark and stale docs were fixed;
  • lost test coverage was restored.

Tests, examples, harness and docs

• Every example deck, test fixture, benchmark and harness case is rewritten in the new deck format. runtests_deck_roundtrip.jl parses every one of them.
• Tests and benchmarks build their setups through the same solve front half (dd1409f76). Finished one-off scripts are deleted.
• Two harness case changes: remove_sections lets a case drop deck sections, replacing the force_termination override, and the case TOMLs read Tearing/PerSurface/S.
• Deck and API parity is now a stated rule, with a parity testset in test/runtests_solve_api.jl (ec07adbb1).
• docs/src/api.md documents the script form of each stage.

Migration

The run errors on every old key, so a deck can be migrated by running it and following the messages. The tables give the full list.

Deck sections and keys

Old	New
[*] verbose, HDF5_filename, output_filename, write_outputs_to_HDF5	[Run] verbose, [Run] output_file
force_termination	leave out the later sections
[SOL_INPUT], [LAR_INPUT], [TJ_ANALYTIC_INPUT] + eq_type	[Solovev], [LargeAspectRatio], [TJAnalytic] (no eq_type)
eq_type = "tj_analytic_direct"	[TJAnalytic] solver = "direct"
grid_type = "log_asymptotic"	grid_type = "auto"
[Wall] shape = "nowall"	leave out [Wall]
[Wall] with unused parameters	only the chosen shape's parameters
[ForceFreeStates] nn_low, nn_high	nn = n or nn = [n1, n2, ...]
[ForceFreeStates] integrator = "forward" (lowercase)	"Forward", "Riccati", "Galerkin"
[ForceFreeStates] eulerlagrange_tolerance, renorm_threshold, numsteps_init, numunorms_init, save_interval	[Forward] or [Riccati]
[ForceFreeStates] nchunks, extended_precision_bvp	[Riccati]
[ForceFreeStates] kinetic_source	kinetic = "KineticForces"; kinetic_factor is only a scale
[ForceFreeStates] sing_start	removed
[ForceFreeStates] gal_nx, gal_nq, gal_pfac, gal_dx0, gal_dx1, gal_dx2, gal_cutoff, gal_tol, gal_gnstep, gal_dx1dx2_flag, gal_rpec_flag, gal_edge_onesided	[Galerkin] nx, nq, ... (no prefix)
[ForceFreeStates] gal_solver	removed (automatic)
[ForceFreeStates] gal_sing_order	[ForceFreeStates] frobenius_order
[ForceFreeStates] gal_sing_order_ceiling	[Galerkin] frobenius_order_ceiling
[ForceFreeStates] gal_match_flag	add a [Match] section
[ForceFreeStates] gal_ideal_flag	[Match] ideal
[ForceFreeStates] gal_eta, gal_rho, gal_rotation, gal_gamma	[Match] eta, rho, rotation, gamma (default from the kinetic data)
[ForceFreeStates] gal_inner_solver, gal_inner_xfac, gal_inner_nx, gal_inner_nq, gal_inner_cutoff, gal_inner_kmax	[GGJ] solver, xfac, nx, nq, cutoff, kmax
[DEBUG] gal_basis_output	[Galerkin] basis_output
[DEBUG] (the rest)	removed
[LocalStability] ode_tol, n_scan, max_alpha_scale, alpha_tol, psi_window	[Ballooning]
[KineticForces] ion_species, ion_charge, ion_mass	[[Equilibrium.ion_species]] with z, m
[KineticForces] impurity_charge, impurity_mass	[Equilibrium]
[KineticForces] kinetic_file, [SLAYER] profile_file, profile_group	[Equilibrium] kinetic_file, kinetic_group
[KineticForces] nn	removed (taken from the run)
[SLAYER] enabled	a [Tearing] section with inner_model = "SLAYER"
[SLAYER] inner_model, coupling_mode, msing_max, omega_E_kHz	[Tearing] (the SLAYER variant is [SLAYER] variant)
[SLAYER.scan_grid] and the scan keys (scan_mode, Q_re_range, nre, amr_passes, pole_threshold, polish_roots, ...)	[Tearing.search]
[SLAYER] gap_kHz_threshold	[Tearing.search] gamma_gap_threshold, in 1/s
[SLAYER.CriticalResonantField]	[Tearing.CriticalResonantField]
[SLAYER] zeff, chi_perp_e, chi_perp_i, chi_tor	[Equilibrium] zeff_default, chi_e_default, chi_i_default, chi_phi_default
[SLAYER] ion_mass, impurity_charge, tau_E, resistivity_model, lnLambda_form	[Equilibrium] (ion_mass via ion_species)
[SLAYER] bt	removed (B_T comes from the equilibrium's F)
[PerturbedEquilibrium] compute_response	removed

Relative paths in a deck are now resolved against the deck's directory, including the coil files (dat_dir, dat_file, h5_file).

API names

Old	New
EulerLagrangeProblem	ForceBalanceProblem
ForceBalanceProblem(...; kinetic_source, kinetic_factor > 0)	ForceBalanceProblem(...; kinetic=KineticForcesControl() / KineticMatrices / TestKineticMatrices(; scale))
ODE settings on the problem	Forward(; ...), Riccati(; ...)
Galerkin(; solver)	Galerkin() (automatic factorization)
Ballooning(; theta_k, alpha_boundary_scan, verbose)	LocalStabilityProblem(eq; theta_k, alpha_boundary_scan); solve(...; verbose)
WallShapeSettings	nothing or a WallSettings subtype
write_outputs_to_HDF5, write_to_hdf5!	save_to_hdf5(target, result)
Tearing.Runner, SLAYERResult, slayer_control_from_toml	solve(TearingProblem(ffs; ...), SLAYERModel(...) / GGJModel(...); config=TearingConfig(...)) -> TearingResult
result omega_Hz, gamma_Hz (and _secondary)	omega, gamma (rad/s, 1/s)
solve_inner(model, p, Q) with model-specific conventions	solve_inner(model, p, s) with s = γ − iω; layer_response(model, p, z)
zero_winding	removed; winding is +1
ggj_growth_rates(; omega_range) bounding −ω	bounds ω
ggj_inner_deltas	removed
ErrorFields stage driver (run_error_fields)	error_fields(resp, coil_sets, coil_config; ...) -> ErrorFieldsResult
KineticForces stage driver	kinetic_torque(pe, ctrl) -> KineticForcesResult
main(...) -> (; ffs, pe, slayer, local_stability, coil_sensitivities, monte_carlo, locking_risk, efc_couplings)	(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)
compute_torque_all_methods! species from the control	species from equil.kinetic (required)
Ballooning verbose = true by default	verbose = false by default
newly exported	kinetic_matrices, KineticMatrices, layer_overlap_scan, layer_eta, layer_rho, layer_omega_E, IonSpecies, CriticalResonantFieldControl, CriticalResonantFieldResult

HDF5 changes

Path	Change
Tearing/PerSurface/lu	renamed Tearing/PerSurface/S
PerturbedEquilibrium/Response/*	fields that were not reconstructed are absent, not empty
Tearing/**/*Winding	+1 where it was −1
ForceFreeStates/LayerOverlap	written only when psilim_from_layer_overlap caps the domain
KineticForces/Validity	written with the torque or the kinetic force-free-states solve
Info/Runtimes/match	new: wall-clock time of the [Match] stage
Input/	inputs recorded as deck sections, on the script path too
metric and b_psi_modes attributes	labels corrected to g_ij/J and J·b^ψ

Behavior changes to know about

• A kinetic run (a kinetic source set) with kinetic_factor = 0 still builds the kinetic matrices and pins kinetic-resonance knots in the ψ grid. To run ideal, leave out kinetic.
• With no kinetic source, the two-pass ψ grid refinement no longer reads the [KineticForces] profile factors.
• A deck stops on [LocalStability], [Ballooning], [ForcingTerms] or [KineticForces] without [ForceFreeStates]. It also stops on [ForcingTerms] without [PerturbedEquilibrium], and on a [KineticForces] that neither the torque nor the kinetic source uses.
• psilim_from_layer_overlap = true without [SLAYER] caps the domain with the default SLAYERModel (it used to warn and do nothing).

Release note

• Audience: users
• Numerical impact: none in any tracked quantity (harness @ b521d76); SLAYER roots shift ≤ ~1e-5 relative from the growth-rate convention, and *Winding counts are +1
• Migration: deck sections, keys and API names change throughout; see the migration tables in the PR body. An old deck stops on its first stale key with the fix in the message.

GPEC v2.0.0 makes the deck a script written as TOML: each object a script builds has one deck section, main runs the same public calls a script does, and every result saves with save_to_hdf5 and records its inputs so any gpec.h5 replays.

Regression report

Harness on ITER SDCC, develop @ 728a7e7d5 vs refactor/v2-consolidation @ b521d76be, the current head. Every value matches the earlier run at 54c2bbb6d.

Case                               Result
diiid_error_field                  10 unchanged
diiid_kinetic_multi_n              develop FAILED (deck uses V-only keys); V: et[1] = 0.5186 - 0.1543i
diiid_multi_n                      19 unchanged
diiid_n1                           53 unchanged
diiid_n1_riccati                   17 unchanged
diiid_n1_riccati_precompiled       17 unchanged
diiid_slayer_n1                    21 unchanged, 1 N/A (PerSurface/S, renamed from lu)
diiid_slayer_n1_chi_perp_e         21 unchanged, 1 N/A (same)
diiid_slayer_n1_coupled            5 changed  (override keys are V-only: develop runs uncoupled)
diiid_slayer_n1_toroidal           5 unchanged
diiid_slayer_n1_winding            21 unchanged, 1 N/A (same)
diiid_truncation_n1                12 unchanged
efit_fixedbdy_separatrix           5 unchanged
gal_resistive_diiid                11 unchanged
gal_resistive_diiid_galinner       1 changed  (override keys are V-only: develop ignores the [GGJ] solver)
gal_resistive_pe                   8 N/A on both (designed: stops before PE)
ggj_ray_q500i                      4 unchanged
ggj_reference                      4 unchanged
lar_reversed_shear_n1              22 unchanged
solovev_kinetic_calculated         20 unchanged
solovev_kinetic_multiion           6 unchanged
solovev_kinetic_ntv                7 unchanged
solovev_kinetic_nuzero             14 unchanged
solovev_multi_n                    15 unchanged
solovev_n1                         22 unchanged

The two "changed" cases are not like-for-like. Their overrides use keys that only this branch reads, so develop silently runs a different configuration: the uncoupled per-surface search, and the default GGJ solver. Every case where both refs run the same inputs is unchanged.

Not in this PR

• Plots take a result object or a gpec.h5, with a plot registry (Analysis UX: introduce a GPECResult object + precondition-aware plot registry #326). This is a separate PR stacked on this one.
• PerturbedEquilibrium: a run with a wall mixes a with-wall Λ and a no-wall L #613 (PerturbedEquilibrium - BUGFIX! - Reject a perturbed equilibrium solved with a wall #620, stacked on this one): a perturbed equilibrium with a wall errors, and the Solovev PE examples and fixture move to an equilibrium that is stable without a wall. It moves numbers.
• Renaming the internal gal_* functions into a GalerkinFEM submodule.
• Follow-ups the review sweep logged: the remaining deck/API parity test cases, IMAS output from main, and a few latent issues in code this PR does not change.

---

🤖 Generated with Claude Code

https://claude.ai/code/session_014wKUyAaxKdTJ5ch1boMFE1

[matt-pharr]

Sorry for this eltrich horror of a PR 💀 good for halloween I guess? @jhalpern30 @logan-nc @d-burg

This is wrapped up now and ready for your feedback. I will do a few cleanup passes to make sure that everything works, but no results were changed in this PR besides a few slayer roots because we got rid of the double conjugate convention with SLAYER that moves roots a tiny bit but is in theory 100% equivalent. This PR achieves a lot of things with open issues, including taking GPEC.jl down to 315 lines, decimating the docstrings, deduping utility functions, helpers, etc, ensuring parity between decks/scripting, fixing scripting syntax, making the hdf5 savable via scripting, and fixing a number of inaccuracies reflected in variable names (omega_Hz and friends were not actually in Hz, docstrings that claim false things, etc). Anything that was labeled as a co-author with claude is something that I did not read 100% of the lines of and mostly amounts to easy develop merges and test fixes, everything else got an attentive eye. This PR achieves 90% of the remainder of what I want to achieve before 2.0.0. Do not feel obligated to review the entire thing; if I requested you it is because you will likely be interested in what is changing/reorganized in this PR and you should look at it and offer any feedback you wish to include so I can implement it.

[d-burg]

@matt-pharr Claude and I reviewed the tearing side of this PR: Δ′, SLAYER, GGJ, η, χ and b_crit. The physics holds up. The kernels, layer inputs, critical Δ and b_crit algebra are unchanged, no tolerance moved, and #222 looks closed in the code. Dropping the conjugate also checks out: on the DIII-D-like SLAYER example at develop, Im/|·| is at most 7e-8 on the Δ′ diagonal and 6e-8 on the off-diagonal products.

Five things to look at, verified in the code:

1. GGJ and SLAYER use different ρ. GGJ uses m_i·n_i and SLAYER uses μ_i·m_p·n_e, despite "every layer model sees one plasma". For D with C⁶⁺ at Z_eff = 2, GGJ's ρ is 20% low because it drops the impurity mass, and a file without n_i stops GGJ but not SLAYER. Suggested fix: one helper for both models, from the same Z_eff that η uses, with quasi-neutrality and the impurity already on the data:

   ρ = m_p n_e [ μ_i (Z − Z_eff)/(Z − 1) + μ_Z (Z_eff − 1)/(Z (Z − 1)) ]

   This needs no n_i, and for D + C it equals SLAYER's value exactly, so the shipped decks should not move. It should error, not clamp, when Z_eff is outside [1, Z].

2. The overlap warning is gone. "Integration domain extends past the resistive-layer overlap point" no longer exists, and the scan runs only with the cap on. The body covers the scan and the output but not the warning. A SLAYER run past the overlap point now says nothing, though the matched-asymptotic treatment isn't defined there. examples/DIIID-like_truncation_example/gpec.toml:43 still says the scan is recorded either way.

3. Multi-ion runs now stop. _layer_ion_mass throws for more than one species, and IonSpecies.m is an Int, so an effective D-T mass of 2.5 can't be set. [SLAYER] ion_mass allowed it. This isn't in the migration notes.

4. GGJ and Match stop on χ they don't use. layer_eta and layer_omega_E go through SLAYERModel(), which validates chi_e, chi_i and chi_phi. A file with one non-positive chi_i point therefore kills a GGJ run that doesn't need it, and GGJ has no P_perp_model to avoid it.

5. The harness doesn't compare the new paths. The body notes this for diiid_slayer_n1_coupled, but diiid_slayer_n1_winding has the same problem: develop ignores Tearing.search.root_classification, so its "21 unchanged" is threshold against winding. No case runs GGJ tearing, so GGJ γ isn't tracked.

Traced but not run, for your cleanup pass:

• kinetic_factor without kinetic_source now silently gives an ideal Δ′.
• One Mercier-unstable surface aborts the GGJ stage after earlier scans finish, and no Tearing group is written.
• Coupled GGJ doesn't warn when the Doppler offsets leave omega_range, so a missed root reads as stable.
• test/runtests_slayer_riccati.jl still passes Q-like values to solve_inner, which now takes s in 1/s.

Landing order: this PR doesn't include #617 (the #614 fix, approved) or #619 (chi_perp_Dc follow-up), so ForceBalance.jl:167 still builds the NTV view unconditionally. A trial merge conflicts in 1 file for #617 and 4 for #619. I'm happy to port both onto this PR after it lands so you don't have to absorb them. Tell me if you'd rather take them first.

[logan-nc]

@d-burg most of these feel like mechanical things you and your claude can fix proactively. PleSe try to minimize dead time between human iteration by telling you claude to make proactive fixes where the path is clear. Even if a human decision is needed, at least post a clear proposed plan so the next human inthe loop can just say "yes, do as proposed" if they approve of the idea

[logan-nc]

Tiered review covering structure, physics, conventions and tests. Visual summary: https://claude.ai/artifact/EEN5xY6MXFEBAkftX4HrXD

Pushed 8dda8ca9c, no src/ changes, so the harness stamp holds:

• Deck and ForceBalance are now in the naming Area table. This was the only failing CI check.
• Three deck kinetic-source error checks now call parse_deck instead of GPEC.main.
• runtests_ffs_output.jl is folded into the deck round trip, saving one full pipeline run.
• Removed the struct-echo assertions and a replay testset that a merge had duplicated.
• Fixed three harness labels that still said Hz.
• Fixed the api.md text that still called the problem an Euler-Lagrange system.
• Tests run on alfven at 70fe92e3b; the later amends touch only docs and labels. Results: solve API 165/165, rerun 87/87, deck round trip 113/113.

Physics: the checks pass.

• 6c42f0364: the new SLAYER coefficients are the exact conjugates of the old ones, evaluated at ĝ instead of conj(ĝ). cole_delta maps to Δ(iQ) on both sides. The GGJ plane is s = −iQ/τ with ω relabeled. Winding +1 follows.
• b8f50de22: implements the plan posted on #222, with Sauter/Redl η from the kinetic data, ρ = m_i·n_i and f = Ω_E/2π.

Line comments each carry a fix. The two suggestion blocks and the [GGJ] validation change src/, so they should land together under one harness re-stamp. I can push that batch and re-run the harness on request.

Merging still needs approval from a third-party human reviewer.

[logan-nc]

nchunks is Riccati-only now that Forward dropped it, but this message still sends it to "[Forward] or [Riccati]". A [ForceFreeStates] nchunks would get pointed at [Forward], which then rejects it.

Suggested change

	const _ODE_DECK_KEYS = ("eulerlagrange_tolerance", "renorm_threshold", "numsteps_init", "numunorms_init", "save_interval", "nchunks")
	const _SEARCH_DECK_KEYS = ("scan_mode", "Q_re_range", "Q_im_range", "nre", "nim", "amr_passes", "amr_max_cells", "boxes",
	"multi_box_prescreen_n", "pole_threshold", "pole_threshold_adaptive", "filter_above_poles", "filter_outside_re",
	"gamma_gap_threshold", "polish_roots", "validity_rtol", "root_classification", "store_scan")
	const _KINETIC_FILE_HOME = "name the kinetic-profile file once, as kinetic_file (and kinetic_group) in [Equilibrium]"
	const _MOVED_KEYS = Dict{String,Dict{String,String}}(
	"ForceFreeStates" => Dict{String,String}(
	"nn_low" => "set nn = n, or nn = [n1, n2, ...] for consecutive toroidal modes",
	"nn_high" => "set nn = n, or nn = [n1, n2, ...] for consecutive toroidal modes",
	"extended_precision_bvp" => "moved to [Riccati]",
	const _ODE_DECK_KEYS = ("eulerlagrange_tolerance", "renorm_threshold", "numsteps_init", "numunorms_init", "save_interval")
	const _SEARCH_DECK_KEYS = ("scan_mode", "Q_re_range", "Q_im_range", "nre", "nim", "amr_passes", "amr_max_cells", "boxes",
	"multi_box_prescreen_n", "pole_threshold", "pole_threshold_adaptive", "filter_above_poles", "filter_outside_re",
	"gamma_gap_threshold", "polish_roots", "validity_rtol", "root_classification", "store_scan")
	const _KINETIC_FILE_HOME = "name the kinetic-profile file once, as kinetic_file (and kinetic_group) in [Equilibrium]"
	const _MOVED_KEYS = Dict{String,Dict{String,String}}(
	"ForceFreeStates" => Dict{String,String}(
	"nn_low" => "set nn = n, or nn = [n1, n2, ...] for consecutive toroidal modes",
	"nn_high" => "set nn = n, or nn = [n1, n2, ...] for consecutive toroidal modes",
	"extended_precision_bvp" => "moved to [Riccati]",
	"nchunks" => "moved to [Riccati]",

[logan-nc]

A [Match] the solve cannot serve (a Forward integrator, Riccati without vac_flag, Galerkin without rpec_flag, or several n) only fails inside MatchProblem, after the full force-free-states solve. The same check fits here with the other cross-section checks, so it stops at parse time like every other deck error. It mirrors the MatchProblem guards:

Suggested change

	match = match_from_deck(inputs)
	match = match_from_deck(inputs)
	if match !== nothing
	(; keywords, algorithm) = force_balance
	first(keywords[:nn]) == last(keywords[:nn]) || error("[Match] is single-n: set [ForceFreeStates] nn to one toroidal mode number")
	(algorithm isa Galerkin && algorithm.rpec_flag) || (algorithm isa Riccati && get(keywords, :vac_flag, false)) ||
	error("[Match] needs Δ′: use integrator = \"Galerkin\" with [Galerkin] rpec_flag = true, or \"Riccati\" with vac_flag = true")
	end

Then the [Match] ideal = true case in runtests_solve_api.jl can call parse_deck instead of GPEC.main, which drops one solve from the suite.

[logan-nc]

[GGJ] is the only section that skips section_keywords: its keys go straight into GGJModel options. A typo (for example ncx = 64) is not caught here. It surfaces as an "unsupported keyword argument" MethodError from the backend at the first layer solve, after the force-free-states solve. That breaks the promise in the PR body that unknown keys stop the run.

Proposal: give each backend a list of its option names (_ggj_options(::Val{:ray}) etc., from Base.kwarg_decl of the backend's solve_inner/solve_inner_profile/solve_ray methods, or a hand-written tuple) and validate here with section_keywords. This is small, but it changes src/, so it should go in with the [Match] check above under one harness re-stamp. I'm happy to push both if you'd rather not.

[logan-nc]

Structure question, not a blocker for this PR. ForceFreeStates' own solve lives in the top module because preparing the solve calls Tearing.layer_overlap_scan and the KineticForces model. So the scripting entry point of FFS is not in FFS, and a reader looking for it in src/ForceFreeStates/ won't find it. One way out for later: pass those two in as hooks (the pattern already used for calculated_source), so ForceFreeStates owns solve and ForceBalanceProblem and the top module only wires them. Worth a follow-up issue in a later milestone? If you agree I'll file it.

[logan-nc]

Remaining suggestions are minor non-blocking. @matt-pharr can quickly take them or leave them.

Very impressive!!