Repository navigation
ForceFreeStates/Tearing - API! - Pose inner-layer matching and tearing as problems solved on the published result - #507
Open
matt-pharr wants to merge 17 commits into
Open
matt-pharr wants to merge 17 commits into
matt-pharr wants to merge 17 commits into
Conversation
…nd read all kinetic stages from it
…red layer-parameter builder
…Problem on the published result
…ped inner-layer model
…st-solve matching
github-actions
Bot
added
api
Collaborator
Author
|
CC @jhalpern30 since you were party to the earlier steps of this refactoring |
matt-pharr
changed the title
github-actions
Bot
removed
the
changed-results
…d model config's way Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…e surface geometry in one place resist_eval re-integrated the same six flux-surface averages that resist_geometry already stores on every surface's restype, and build_ggj_inputs repeated layer_parameters' resistivity and mass-density closure. Both are replaced by one ggj_parameters(sing, equil; eta, rho, gamma) on top of restype, shared by the driven match and the GGJ tearing solve; layer_parameters takes the profile source as a keyword so both callers derive eta/rho through it. The two geometry ports agree to <= 2.4e-14 relative on every surface of the LAR, DIIID Galerkin-resistive and DIIID SLAYER decks (the residue tracks q - q_spline at the rational root). Matched outputs move at round-off only (<= 1.2e-10); ideal, Galerkin-ideal and tearing outputs are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…er argument of solve ForceFreeStates.GGJ/SLAYER duplicated InnerLayer's GGJModel/SLAYERModel as configuration wrappers, with forwarding methods and a translation step in the tearing solve. The InnerLayer types are now the inner-layer argument of solve(MatchProblem, ·) and solve(TearingProblem, ·) directly: GGJModel carries its backend keywords as options, forwarded to every solve under call-site keywords; SLAYER's chi fallbacks are read from the SLAYERControl that already held them. Removes Matching/Models.jl and the exported GGJ name that collided with the InnerLayer.GGJ submodule. Deck outputs are byte-identical to the previous commit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
matt-pharr
changed the title
github-actions
Bot
added
the
changed-results
…sistivity the inner-layer defaults everywhere One set of inner-layer defaults, held in one place: - The deck's gal_inner_* keys no longer carry their own grid values (1280/5/10, the Fortran deltac reference); unset keys keep the Galerkin backend's defaults (512/4/1), which agree with the certified ray Δ about 20x more closely on the DIIID matched surfaces at a fraction of the cost. - The ray backend is the default GGJ backend in tearing as well as matching: TearingProblem accepts GGJModel(), and the [SLAYER] deck gains inner_model = "ggj_ray". - layer_parameters and MatchProblem derive the Sauter neoclassical η by default, matching the tearing solve; SpitzerModel() stays selectable. No example deck sets these keys or derives η, so deck outputs are byte-identical to the previous commit; API calls relying on the old defaults move. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ttach kinetic profiles in one way _matched_result rebuilt ForceFreeStatesResult and GalerkinResult from every field positionally; a field-name copy replaces both and no longer breaks when a field is added. attach_kinetic_profiles! forwards its keywords to load_kinetic_profiles instead of restating them, and becomes the one way to attach profiles: the kinetic_file keyword on PlasmaEquilibrium, which restated the same eight keywords again and had no callers, is removed. Deck outputs are byte-identical to the previous commit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf
…seded by MatchProblem scan_resistivity_m2.jl and scan_rotation_m2.jl plotted per-point deck runs that each re-solved the outer region; MatchProblem reuses one outer solve, and the scan is a short loop (docs/api.md). scan_match_m2.jl, added earlier in this branch for that comparison, goes with them. REFACTOR_PLAN.md records the cut-down series. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf
… to what a reader needs User-facing docstrings (MatchProblem, TearingProblem, layer_parameters, MatchResult, attach_kinetic_profiles!) say what is computed, in what units, and how to call it; developer-only file headers, helper docstrings and inline comments shrink to a line or two. Design rationale and history move out of the source. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf
…deal δW from matched results The match product lived inside GalerkinResult, so a Riccati-fed match had nowhere to go and kept only bpen. It now sits on ForceFreeStatesResult beside delta_prime and bpen, for every formalism, and is written under SingularSurfaces/Match/ next to the Δ′ it was solved from. A :matched result drops wp and free_boundary: the ideal δW does not describe the matched plasma, so PerturbedEquilibrium warns and skips instead of mixing ideal and matched physics (as Fortran gpec deallocates its ideal solutions under gal_flag). The ideal reference keeps its δW. Matched-deck datasets are identical apart from the Match/ path. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf
…ofile errors, test that solves leave inputs untouched MatchProblem rejects multi-n results, since the forced eigenvalue and field scale use a single n. Missing-profile errors name the deck keys on the deck path and say to attach before solving (or to ffs.equil) on the API path. New tests check that a solve leaves the equilibrium untouched and that repeated MatchProblem solves leave the outer result untouched. api.md notes that Riccati matching needs vac_flag=true and that GGJ tearing γ are placeholders; stale run_slayer references fixed. REFACTOR_PLAN.md leaves the repository; §7D is to be re-planned. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf
Collaborator
Author
|
Also CC @StuartBenjamin for the new syntax which will be particularly useful for tearing suite |
matt-pharr
marked this pull request as draft
October 8, 2026 21:00
…netic and deck matching rules, fix stale docs
matt-pharr
marked this pull request as ready for review
October 9, 2026 15:16
Collaborator
Author
|
@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
requested changes
Oct 14, 2026
d-burg
left a comment
d-burg
left a comment
Collaborator
There was a problem hiding this comment.
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:
- 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.
- 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.
- 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.'
- “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.
- #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?)
Comment on lines
+15
to
+18
| Fields: `ffs`, `control::SLAYERControl`. | ||
| """ | ||
| # `ffs` is duck-typed (equil, surfaces, delta_prime, dir_path) so tests can use stand-ins. | ||
| struct TearingProblem{R} |
Collaborator
There was a problem hiding this comment.
The comment on line 17 separates the docstring from the struct, and Julia then drops it (checked in a scratch module on 1.11.6: no doc is registered). Folding the note into the docstring fixes it.
Suggested change
| Fields: `ffs`, `control::SLAYERControl`. | |
| """ | |
| # `ffs` is duck-typed (equil, surfaces, delta_prime, dir_path) so tests can use stand-ins. | |
| struct TearingProblem{R} | |
| Fields: `ffs` (duck-typed on `equil`, `surfaces`, `delta_prime`, `dir_path`, so tests can use stand-ins), `control::SLAYERControl`. | |
| """ | |
| struct TearingProblem{R} |
Comment on lines
+99
to
+100
| error("TearingProblem: no kinetic profiles — set profile_file, or attach profiles " * | ||
| "to the equilibrium with attach_kinetic_profiles!(equil, file)") |
Collaborator
There was a problem hiding this comment.
From claude: after an auto-grid solve ffs.equil is a re-formed copy, so attaching to the equilibrium I passed in has no effect and this error then tells me to do what I just did. Please name ffs.equil, as layer_parameters does.
Suggested change
| error("TearingProblem: no kinetic profiles — set profile_file, or attach profiles " * | |
| "to the equilibrium with attach_kinetic_profiles!(equil, file)") | |
| error("TearingProblem: no kinetic profiles on ffs.equil — set profile_file, or attach them with " * | |
| "attach_kinetic_profiles!(ffs.equil, file) (after an auto-grid solve ffs.equil is a re-formed copy, not the equilibrium passed to solve)") |
Comment on lines
+120
to
+121
| tear = solve(TearingProblem(ffs; coupling_mode=:coupled), SLAYERModel()) | ||
| tear.gamma_Hz, tear.rational_q # root-found rates per surface |
Collaborator
There was a problem hiding this comment.
In coupled mode the result holds one global root, so tear.gamma_Hz and tear.rational_q do not line up and the comment is wrong.
Suggested change
| tear = solve(TearingProblem(ffs; coupling_mode=:coupled), SLAYERModel()) | |
| tear.gamma_Hz, tear.rational_q # root-found rates per surface | |
| tear = solve(TearingProblem(ffs), SLAYERModel()) # coupling_mode=:coupled returns one global root instead | |
| tear.gamma_Hz, tear.rational_q # root-found rate at each rational surface |
| @test r_off.enabled == false | ||
| end | ||
|
|
||
| @testset "the default :ray backend drives the tearing solve" begin |
Collaborator
There was a problem hiding this comment.
From claude: this test returns at the empty-surface short-circuit, so it does not show that the :ray backend drives anything. Nothing in the suite runs _tearing_params or _profiles_from_equilibrium. Could you add one small case that solves via profile_file and via attached profiles and compares the per-surface parameters?
Suggested change
| @testset "the default :ray backend drives the tearing solve" begin | |
| @testset "TearingProblem accepts a GGJModel and short-circuits on an empty surface list" begin |
Comment on lines
63
to
64
| - `bpen::Matrix{ComplexF64}` - `(msing × numpert_total)` penetrated resonant field from the | ||
| inner layer, per surface and driving mode; all zeros under `:ideal` closure. |
Collaborator
There was a problem hiding this comment.
After the move, result.match and SingularSurfaces/Match/* are indexed by the matched subset, but nothing says which surfaces those are, and result.surfaces beside them has a different length whenever the Galerkin domain or m-band drops a rational. Please carry the matched surface list on MatchResult and write it to the group
Suggested change
| - `bpen::Matrix{ComplexF64}` - `(msing × numpert_total)` penetrated resonant field from the | |
| inner layer, per surface and driving mode; all zeros under `:ideal` closure. | |
| - `bpen::Matrix{ComplexF64}` - Penetrated resonant field from the inner layer, per surface and | |
| driving mode. `(msing × numpert_total)` zeros under `:ideal` closure; under `:matched` the rows | |
| follow the matched surfaces, which for Galerkin is the in-domain in-band subset of `surfaces`. |
| gal_inner_cutoff::Int = 5 # cells carrying the large solution as driving term | ||
| gal_inner_kmax::Int = 8 # large-x asymptotic series order (↔ order_pow) | ||
| # --- Inner-layer "galerkin" backend overrides (gal_inner_solver = "galerkin"; unset keeps the backend default) --- | ||
| gal_inner_xfac::Union{Nothing,Float64} = nothing # asymptotic-matching radius multiplier |
Collaborator
There was a problem hiding this comment.
From claude: for a deck with gal_inner_solver = "galerkin" and the knobs unset, this silently changes xfac 10→1, nx 1280→512 and nq 5→4. Please qualify the numerical-impact line and replace “Decks unchanged”. I’d also like the numbers behind “~20× closer to ray”, and a harness quantity on the Galerkin inner backend.
| rotation::Union{Nothing,AbstractVector{<:Real}}=nothing, | ||
| mu_i::Real=2.0, | ||
| zeff::Real=1.0, | ||
| resistivity_model::NeoResistivityModel=SauterNeoModel(), |
Collaborator
| end | ||
|
|
||
| # n_e in m⁻³, T_e in eV, and the E×B frequency at ψ from either profile container. | ||
| _layer_profile(kp::Equilibrium.KineticProfileSplines, ψ) = (kp.ne_spline(ψ), kp.Te_spline(ψ) / E_CHG, kp.omegaE_spline(ψ)) |
Collaborator
There was a problem hiding this comment.
The attached profiles carry zeff_spline, loglam_spline and ni_spline, but this reads only n_e, T_e and ω_E and then uses zeff=1.0 and a recomputed lnΛ. The default should take Zeff from equil.kinetic and treat the keywords as overrides
| ffs::ForceFreeStatesResult; | ||
| eta::Union{Nothing,AbstractVector{<:Real}}=nothing, | ||
| rho::Union{Nothing,AbstractVector{<:Real}}=nothing, | ||
| rotation::Union{Nothing,AbstractVector{<:Real}}=nothing, |
Collaborator
There was a problem hiding this comment.
From claude: for the advertised rotation scan the user must know the matched-surface count before posing the problem; today that takes a first solve or a dummy MatchProblem(ffs; ideal=true). I’d accept a scalar for each of eta/rho/rotation and export an accessor for the matched surfaces.
| return cof[1:2msing, :], cof[(2msing+1):4msing, :], residual | ||
| end | ||
|
|
||
| # Write `m` under SingularSurfaces/Match/ (debug-only: exempt from the metadata contract). |
Collaborator
There was a problem hiding this comment.
From claude: with δW dropped and PerturbedEquilibrium skipped, this group is where a matched run’s results live, so “debug-only, exempt from metadata” no longer fits. Please annotate the datasets, write reconnected_flux, and record the inner model, backend and options.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Release note
Matching and tearing are now problems solved on a finished force-free-states result.
examples/DIIID-like_ideal_example/files (~1.5 min).vac_flag=true; no matched ξ.TearingProblemwithGGJModel()runs; its γ are placeholders.What changed
MatchProblem+GGJModel(): matching after the solve.TearingProblem+SLAYERModel()/GGJModel(): growth rates.result.matchholds the match, any formalism.InnerLayermodel types are thesolveargument.Migration
ResistiveMatch(...)→MatchProblem(ffs; ...)+GGJModel().Runner.run_slayer→TearingProblem+SLAYERModel().build_ggj_inputs→layer_parameters+ggj_parameters.attach_kinetic_profiles!(eq, file).result.galerkin.match→result.match.GalerkinIntegration/Match/→SingularSurfaces/Match/.[SLAYER]gainsinner_model = "ggj_ray".gal_inner_solver = "galerkin": setgal_inner_xfac=10, nx=1280, nq=5for the old grid.New defaults
:rayeverywhere, tearing included.Regression report
6fc89d021vs289f6ea83.gal_resistive_diiid: 3 match quantities N/A on develop (path rename).diiid_slayer_n1: γ 0.02% (known root-finder scatter).gal_resistive_pe: N/A on both refs (Galerkin PE gap).diiid_error_field.Full report
Reviewing this PR
_compute_matchre-indents a block.gal_match_rpec→_compute_match:ctrl.gal_*→prob.*, the backendif→solve_inner_profile(model, …)._match_system, verbatim; a dead duplicate inresonant_match_rpecis deleted.GalMatchResult→MatchResult, plusreconnected_flux.write_galerkin!'sMatch/block →write_match!, verbatim.run_slayer→solve(::TearingProblem)+_tearing_params; only the profile source changed.build_ggj_inputsη/ρ →layer_parameters; τ_A/τ_R →ggj_parameters.resist_evaldeleted: a second port ofresist_geometry.kinetic_profilesargument →equil.kinetic, in 6 driver functions.MatchProblem,_matched_result/_with, the δW drop.TearingProblem,_profiles_from_equilibrium,_tearing_params(::GGJModel).GGJModeloptions and their forwarding (InnerLayer/GGJ/GGJ.jl).attach_kinetic_profiles!,equil.kinetic.run_force_free_states; new defaults.runtests_matching_models.jl, newruntests_solve_api.jltestsets.Notes for reviewers
src/is −185 lines net vs develop.bpendiffer near the edge.CoreTypes.jl).REFACTOR_PLAN.mdleaves the repo (−1117 lines).🤖 Generated with Claude Code
https://claude.ai/code/session_01VkufdJ1x3rjhWKMz9FzaMf