Repository navigation
Conversation
There was a problem hiding this comment.
Copilot review overview
🟡 Changes recommended
The review found one critical integration issue and two moderate correctness/schema issues.
Get a fresh assessment by requesting another Copilot review.
Review effort: Lite
Findings: 1
Open (2)
What changed in this PR
Adds an optional energy-decomposition diagnostic for free-boundary eigenmodes, with effective-field and standard-form calculations, HDF5 output, tests, regression tracking, and documentation.
Changes:
- Adds decomposition APIs and pipeline/TOML integration.
- Adds regression quantities and example configuration.
- Adds tests, documentation, and HDF5 schema guidance.
- Refactors displacement regularization to a keyword argument.
Unresolved findings: critical integration gating issue (4 votes), moderate disabled-field HDF5 shape issue (3 votes), and moderate truncated-domain integration issue (1 vote).
| File | Description |
|---|---|
test/runtests.jl |
Registers energy-decomposition tests. |
test/runtests_energy_decomposition.jl |
Tests calculations, toggles, metadata, and output. |
src/PerturbedEquilibrium/PerturbedEquilibrium.jl |
Includes and exports the new APIs. |
src/PerturbedEquilibrium/FieldReconstruction.jl |
Refactors regularization argument handling. |
src/PerturbedEquilibrium/EnergyDecomposition.jl |
Implements decomposition, diagnostics, and HDF5 writing. |
src/GeneralizedPerturbedEquilibrium.jl |
Integrates the optional stage and scripting API. |
regression-harness/cases/solovev_n1.toml |
Tracks decomposition diagnostics. |
regression-harness/cases/diiid_n1.toml |
Tracks decomposition diagnostics. |
examples/Solovev_ideal_example/gpec.toml |
Enables the feature in the example. |
examples/DIIID-like_ideal_example/gpec.toml |
Enables the feature in the example. |
docs/src/energy_decomposition.md |
Documents theory, usage, outputs, and API. |
docs/src/citations.md |
Adds supporting citations. |
docs/src/api.md |
Documents the scripting entry point. |
docs/make.jl |
Adds the documentation page to navigation. |
docs/development/references.md |
Records implementation references. |
docs/development/hdf5-conventions.md |
Documents the new output schema. |
docs/development/architecture.md |
Updates the architecture overview. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| ForceFreeStates.require_solution(result, "energy decomposition") || return nothing | ||
| ForceFreeStates.require(result, :free_boundary, "energy decomposition") || return nothing |
| :effective_b_curl_residual_density, :b_squared_density, :current_coupling_density, :pressure_compression_density) | ||
|
|
||
| # Empty result arrays are written zero-extent (the file's not-computed sentinel). | ||
| _h5_value(x::AbstractArray) = isempty(x) ? similar(x, 0) : x |
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
Consolidated review from three independent passes (structural, performance, code clarity), each claim checked against the source. Ordered by what I'd want settled before merge. The physics itself was not reviewed here — see the end.
1. Blocking: the stage aborts the run instead of skipping
run_energy_decomposition (GeneralizedPerturbedEquilibrium.jl, new hunk) gates only require_solution and require(:free_boundary), which warn and return. But decompose_energy throws ArgumentError on nlow != nhigh and mats.kinetic !== nothing (EnergyDecomposition.jl:698–699). A legal multi-n or kinetic deck that carries [EnergyDecomposition] therefore kills the whole run after FFS — before PE, KineticForces or SLAYER — contradicting the stage docstring ("warned and skipped, like every optional stage") and docs/src/energy_decomposition.md. Copilot flagged the same gap. ~10 lines: warn-and-return for both conditions in the driver, keep the throws for scripted use.
2. Decide before merge: where the output lives
- The group
ForceFreeStates/EnergyDecomposition/bakes "free-boundary eigenmodes only" into the schema, while the producer is a PE-module stage. Free-boundary enters the code through justwt[:,k](:342) andep[k](:814); adecompose_energy(ffs, edge::Vector; reference=nothing)entry would serve fixed-boundary modes and forced responses at ~half a day, andPerturbedEquilibrium/…would then be the natural home. HDF5 renames are clean breaks here with no shim, so the placement is worth deciding now even if the generalization waits. eigenmode_index = copy(eigenmodes)(:719) echoes a control value into an output group (hdf5-conventions.md: never echo control-struct values outsideInput/). It is defensible as the coordinate of theeigenmodedimension — in which case attach it as a dimension scale; onlypsiis attached today (:926–927).
3. Structural: the geometry block shadows Equilibrium
surface_geometry!(:160–166) reproducesEquilibrium.flux_surface_metric(FluxSurfaceMetrics.jl:16–23) line for line — the function whose docstring declares it "the single source of truth for the |∇ψ| computation reused across the FFS and PE modules" — then extends it withw21/w22/w31/w32and the tangent basis. Extendflux_surface_metric(or add a sibling in Equilibrium) and call it, rather than carrying a second implementation in PE.- The same block re-interpolates geometry with 13 bicubic evaluations per point at the equilibrium knots, where
rzphi_*.nodal_derivs.partials[1:3, i, j]already holds f, ∂f/∂ψ, ∂f/∂θ andmake_metric(Fourfit.jl:193–197) andGeometryProfiles.jl:44read them directly. Only the appended edge node needs the interpolant. The feature is cheap overall (order 0.1–1 s), so this is about consistency with the existing pattern more than speed. metric_bsq_interpolant(:213–230) is a private, persistent B² that competes withequil.eqfun_Band is rebuilt on every call (:717). Whether metric-consistent B² is the right physics is for the physics review; either way its owner should be Equilibrium ormake_metric, so a future consumer (LocalStability shear, KineticForces) does not grow a third B².eigenmode_modes(:361–367) re-encodes −A⁻¹(Bξ′ + Cξ) by hand;compute_node_xi_s!(ForceFreeStates/EulerLagrange.jl:1397, used atResult.jl:174) already does exactly that.surface_geometry!andcurvature_kernel!depend only on(equil, psi)but sit inside the eigenmode loop (:755–770): all geometry is computed twice foreigenmodes = [1, 2]. Hoist per surface.
4. Robustness
surface_geometry!(…; bsq=nothing)leavesgeom.bsqzero, andcurvature_kernel!divides by it (:285). Makebsqrequired.similar(x, 0)(:921–922) collapses disabled 2-D/3-D outputs to rank-1 while the metadata table declares(psi, eigenmode)/(psi, theta, eigenmode)— the writer cannot attach the declared dims (Copilot). The zero-extent sentinel has precedent, but only if the rank is kept.- Single-n is load-bearing, not just declared:
N = ffs.numpert_total(:349) is written intoxi_s[i, :]of widthmpert(:367), aDimensionMismatchthe guard at:698is what prevents. - Constructors hard-code field counts (
for _ in 1:16,1:6,1:13,1:6,1:8at:131/248/400/451/539);fieldcount(T)stops them rotting when a field is added.
5. Cleanup
- Code comments are clean — no run exceeds two lines. The narration is in docstrings:
:208–211and:409–412carry session history (keep the one-sentence reason; drop the rest), andtest:176,:208record observed percentages that belong in the docs, not a comment. docs/src/energy_decomposition.md:64says "better than 1 %"; the tests record 1.4 %, 1.2 % and 3.3 % at the mpsi = 16 fixture (test:176,:208). "≲1–3 % at mpsi = 16, converging with resolution" matches the evidence.test/runtests_energy_decomposition.jl(292 lines) would be ~200 with the same coverage:_ed_metadata_violations(:216–231) duplicatestest/h5_metadata_check.jl— include it;:32–50tests thereg_spotkeyword, not this kernel;:70–73asserts kwdef defaults;:86–90re-derive K1/K2/K3 from:293–295(tautological);:99–128hand-plumbs the curl residual that:173already pins;decompose_energyruns eight times — hoist the shared results.:287–291writes into the shared fixture'sgpec.h5and only works because it is the last testset — write to a temp copy.- Five structs describe their fields in prose rather than a
## Fieldsblock (:104,:232,:376,:435,:516); the unnamed1e-10at:315; μ₀j^θ and μ₀j^ζ computed three times (:286,:475,:572); the 13 hand-named accumulators plus three mirrored copy blocks (:737–749,:796–837) could be driven from one name tuple, as_allocate_result!already half-does.
6. Before merge
- The branch is ~250 commits behind and develop has since rewritten
FieldReconstruction.jl(225 lines) and the driver (+455); your 6-linereg_spotchange and 47-line hook land on rewritten code.DIRTYhere is semantic, not textual (zero conflict markers) — worth re-reading both hunks after the merge-forward. - The regression report was a laptop run (arm64, 4 threads) against a develop that has moved; please re-run on the shared machine after merging forward and refresh the stamp.
- This review deliberately did not cover the physics (metric-B² vs
eqfun_B, the pointwise covariant lowering your docstring at:408–412justifies, the ξ_s recomputation at knots). Those deserve a physics pass and I'd ask for one before approval.
Credit where due: Copilot's two inline comments — the gating gap and the similar(x, 0) rank collapse — were both real and are folded in above.
🤖 Generated with Claude Code
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
Physics pass, to close the gap my earlier review left open. Checked term by term against the Fortran counterpart and the theory note on origin/recon of the Fortran repository (gpout_recon/gpout_recon2 in gpout.f, gpeq_c/gpeq_K/gpeq_shear/gpeq_curvature/gpeq_cveri/gpeq_firstform in gpeq.f, docs/tex/recon/main.tex) — not on Fortran develop, which is worth stating in the docs since "fidelity to Fortran" means fidelity to that branch. Every claim marked verified was read or derived; the Fortran line references were checked by me as well.
Verdict: correct with caveats. No missing terms, no placeholders; every kernel either matches origin/recon line for line or is a defensible, more accurate deviation. One thing should be settled before merge — see §6 — because it is the untested explanation for the 1.4–3.3 % mismatch.
1. Both δW forms — verified
- Standard form (
:565–585): derived from −½∫ξ*·F(ξ) with the full ξ (the ξ∥ pieces of the J×Q and ∇·ξ terms cancel, so using the full Jξ^θ, Jξ^ζ is right).det_termat:572–575is J·det[j, b, ξ*] with j^ψ = 0; ∇·ξ at:577is (1/J)[∂ψ(Jξ^ψ) + ∂θ(Jξ^θ) − 2πin Jξ^ζ] in the exp(2πi(mθ−nζ)) convention. Term for term identical togpeq_firstform. One divergence: ∂θ(Jξ^θ) is spectral here (:415–418) and a periodic spline there — fine, but annotate it as a Fortran deviation. - Effective field b_eff = b + ξ_n μ₀ j×n̂ (
:455–480): all six component formulas verified by hand and identical togpeq_c(cwt/cwz/cvp/cvt/cvz). μ₀j^θ = −F′/J and μ₀j^ζ = −(μ₀p)′/χ′ − qF′/J follow from ∇×B and (j×B)·∇ψ = p′, givenF_deriv= d(2πF)/dψ_N (inferred from the Fortransq%f1(1)usage; please confirm). - K = K1 + K2 + K3 (
:293–295): identical togpeq_Kandmain.tex; units consistent (K|ξ_n|² is Pa); K3 = 2p′κ_ψ checks against the interchange term. Caveat: K1, K2 and the completing-the-square identity were verified againstmain.texandgpeq_Konly — Bernstein 1983, Chance 1992 and Boozer 2004 are not indocs/resources/, so the equation numbers this PR adds toreferences.md/citations.mdare unverified. Please either add the PDFs or cite the theory note. - ½, Re, conjugation, normalization: Re[(Jc^i)* c_i]/(μ₀J) is Hermitian-correct; the ½ matches Fortran
dw_raw = 0.5*total; the 2μ₀/ψ₀² scale (:602) with unit-power-normalizedwt[:,k]as edge BC (:342) is the same convention as Fortran'sdcon_facafteridcon.fundoes its rescaling;ep = diag(wt†·wp·wt)withwp = (u2/u1)/ψ₀²matchesfree.f. No vacuum or surface term belongs in the comparison —epis the plasma boundary term Ξ†(FΞ′+KΞ), i.e. the volume δW_p;evis separate. Verified in concept; the DCON ½/μ₀ scale was not re-derived from Glasser 2016, but the 1e-3 agreement rules out a missing O(1) factor.
2. Metric-B² vs eqfun_B — your choice is right; say so at the call site
Fortran gpeq_K and gpeq_curvature take B² from eqfun (bsq_val = eqfun%f(1)**2, with fx/fy for its derivatives). make_metric builds only g_ij/J, J, J_ψ — no eqfun — so the Euler–Lagrange operator and the b reconstruction carry the metric B², and the C∥ = K2|ξ_n|² cancellation holds only with that B². Using it is operator-consistent and an improvement over Fortran. The ~1 % is the mismatch between the equilibrium's direct |B| and the inverse-map metric at mpsi = 16 — coordinate-transformation accuracy converging with equilibrium resolution, not hidden physics. It needs the project's one-line deviation note at :285, :291, :294, not only in the docstring.
3. Pointwise covariant lowering — the claim holds, with two corrections
gpeq_c back-transforms band-limited modes before squaring, so ∮|C|² and ∮K2|ξ_n|² each lose out-of-band energy that does not cancel between them; pointwise products on the mtheta grid are Parseval-exact (aliasing aside). So the Julia effective-field numbers will differ from Fortran recon reference values at the few-% level, and Julia is the more accurate. Two corrections to the docstring at :408–412: "resolution-independent" should read "ψ-resolution-independent, mpert-dependent" (the deficit shrinks with the band); and this is a hybrid — Jξ^ψ is formed pointwise (:416, :426) but Jξ^θ, Jξ^ζ still come from the truncated jmat convolution in compute_contra_displacements — the docstring implies all-pointwise.
4. ξ_s at the knots — verified, no physical content in knot-vs-node
−A⁻¹(Bξ′ + Cξ) (:361–367) matches Fortran gpeq_sol and Julia's own compute_node_xi_s! (EulerLagrange.jl:1397, same raw A; Cholesky vs LDLᵀ immaterial). It is the algebraic Euler–Lagrange relation for Ξ_s, so evaluating at knots rather than ODE nodes changes only the interpolation accuracy of (ξ, ξ′), nothing physical. That argues for calling the existing helper rather than re-encoding it (as my earlier review noted).
5. The curl residual — right check, state its scope
(:495–515) is identical to gpeq_cveri and is the Ξ_s equation itself (∇ψ·∇×C = 0 ⇔ δj·∇ψ = 0 closure), so at the knots it holds by construction. It validates the ξ_s algebra, the g22/g23/g33 lowering and the b_θ, b_ζ assembly against A, B, C — not the ψ-ODE and not b_ψ or g11/g12/g13/|∇ψ|². A roundoff-level value says nothing about the eigenmode's correctness. Worth one sentence in the docs so it isn't read as a global check.
6. Convergence — the one thing to settle before merge
Expected: ~h⁴ for the cubic-spline ψ-quadrature on a smooth integrand, spectral in θ. Observed from the PR's own numbers (mpsi 16 → 128): effective form ~h^0.9, standard form ~h^1.8. That is not quadrature-order behaviour.
The leading explanation is a mis-documented Fortran deviation. gpout.f says verbatim, just above the energy check: "Perform the energy check on the DCON solution grid. Do not use rzphi%xs for this comparison; it is a coarse geometry grid." — and integrates on psifac. decomposition_grid (:309–316) integrates on the rzphi_xs equilibrium knots, and the docstring at :305 says gpout_recon "loops over the same knots" — true of its diagnostics loop, not of its energy integral. With the 1/(m − nq) structure in ξ_s and Jξ^θ, the coarse-knot integral is the natural suspect for both the size of the mismatch and the sub-quadrature rate; the slower effective-form rate is consistent with K depending on spline derivatives (bsq_psi, dratio) that lose an order. A confound: ep itself moves with mpsi.
Discriminating test, cheap to run: fix mpsi, integrate on ODE nodes ∪ knots instead of knots alone, and report how the relative error changes. If it collapses, fix the grid and the docstring; a separate mpert scan would then bound the band-truncation floor at mpsi = 128.
7. Fortran correspondence and the annotations owed
gpout_recon/gpeq_c/gpeq_K/gpeq_shear/gpeq_curvature/gpeq_cveri/gpeq_epf/gpeq_dst/gpeq_firstform ↔ decompose_energy/effective_field!/curvature_kernel!/curl_residual!/surface_densities!. Term-level divergences, each needing the project's one-line deviation comment at the call site: (a) metric B² vs eqfun B²; (b) pointwise vs band-truncated C_i, ξ_n; (c) spectral vs spline ∂θ; (d) radial grid knots vs psifac; (e) reg_spot = 0 vs Fortran's regularized xmp1/xmt path in gpeq_firstform (inferred from variable names, not traced). Also: the long_name on current_coupling and pressure_compression (:940–941) says "imaginary part is roundoff" — true of the volume totals, not of the per-surface profiles; reword.
Not checked: the Bernstein/Chance/Boozer equation numbers; the DCON ½/μ₀ convention from Glasser 2016 directly; the edge amplitude in sum_eigenmode_contributions; gpeq_firstform's regularization path; gpout_recon2's radial grid; aliasing error of the mtheta Riemann sum.
🤖 Generated with Claude Code
| - Published: Physics of Plasmas **24**, 032505 (2017) | ||
| - Describes: Self-consistent coupling with neoclassical effects | ||
|
|
||
| The energy decomposition of free-boundary eigenmodes (`EnergyDecomposition.jl`) follows: |
There was a problem hiding this comment.
I am not sure if these should be here.
- This does not conform to the standard module markdown section formatting (minor)
- This
references.mdis intended to be references for people to cite our code - the papers by us, the developers (not just every paper about relevant physics).
I am torn since this new feature does not have a paper yet, but I think we should wait for @jaesun57 to publish something using it and then add that paper here. In the meantime, these papers can be cited in source code annotations to help developers looking into the details without confusing the general users into thinking Bernstein and Boozer are GPEC developers 😆
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
I flagged a few general organizational things for now. A more in depth AI assisted review can come in a bit
| xwp_modes, xwt_modes, xwz_modes, xmt_modes, xmz_modes = compute_contra_displacements( | ||
| xi_psi_modes, clebsch_psi1, clebsch_alpha, | ||
| psi_grid, equil, ffs, metric, ctrl | ||
| psi_grid, equil, ffs, metric; reg_spot=ctrl.reg_spot |
There was a problem hiding this comment.
This has me a bit worried. We just had a number of issues/PRs about regularization from @jhalpern30 and @matt-pharr and some deeper thought than I can give skimming this during the ITPEA meeting is needed to ensure any use here is consistent with what we are doing elsewhere.
| push!(runtimes, "force_free_states" => ffs_dt) | ||
| @info "Force-Free States completed in $(@sprintf("%.3f", ffs_dt)) s" | ||
|
|
||
| run_energy_decomposition(ffs_result, inputs) |
There was a problem hiding this comment.
@matt-pharr this is already quite clean - but do you need it tweaked into a "solve" type format to conform to the new standard API you are enforcing? Is there an open issue describing it that we can point claude to and ask it to conform to for things like this?
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
Review of #465: energy decomposition of free-boundary eigenmodes
@jaesun57, first: this is a lot of careful work, and much of it is good.
You ported a hard diagnostic from the Fortran.
You checked the result against ep[k] in two independent forms.
You also chased the convergence question further than most people would.
The finding that splining ξ through the saved ODE steps hid a cancellation, so the old 1.7e-5 agreement was luck and not convergence, is a real result.
solution_state_at is a small, well-argued, reusable piece.
You also answered every item of the last round of reviews.
Thank you.
Even so, I'm requesting changes, and they are not small.
The physics can stay.
The way it is built into the code cannot.
The rest of this review explains why, then gives a concrete plan for the rework.
Important
This PR (and its rework) must not be merged without approval from a third-party human reviewer. That is non-negotiable.
#465 has no requested reviewer, and the author is the only assignee.
Please request a human reviewer, or mark the PR as Draft while the rework is under way.
Big picture
1. It does not belong in the main run chain
The decomposition is a verification diagnostic.
It reconstructs δW_p of eigenmodes that ForceFreeStates (FFS) has already computed, and checks the result against ep[k].
Nothing downstream consumes its output.
Even so, it adds:
- a new
[EnergyDecomposition]TOML section with four keys - a new stage in
main_from_inputs - a new
gpec.h5group with about 38 annotated datasets - the section in two example decks
- five rows in each of two regression-harness cases
All of that is upkeep the whole team carries, for a calculation a user runs on purpose, a few times, while investigating a mode.
That is what a scripted diagnostic is for: call one function on a ForceFreeStatesResult in a script or notebook and plot what comes back.
The PR already contains that entry point (energy_decomposition(ffs; eigenmodes=...)), so the TOML section and the stage add nothing a user can't already do.
The placement is also inconsistent:
- It is configured as its own stage.
- The code lives in
PerturbedEquilibrium. - It writes under
PerturbedEquilibrium/even when PerturbedEquilibrium never runs. - It is a diagnostic of FFS eigenmodes.
The rework below moves it into Analysis, which already loads after PerturbedEquilibrium and is where our post-processing lives.
2. It recomputes geometry the code already has
Equilibrium.flux_surface_basis hand-expands the Fortran gpeq_c v and w matrices from a 12-scalar positional argument list.
-
The
vmatrix already exists four times. The same tangent-basis expressions, built from the samerzphi_*.nodal_derivs, are in:ForceFreeStates.make_metric(src/ForceFreeStates/Fourfit.jl)build_kinetic_metric_matrices(same file)src/Equilibrium/Equilibrium.jl(two places)src/ForceFreeStates/Surfaces/Resist.jl
This PR adds a fifth copy.
-
ffs.metricalready holds the result on exactly the grid you integrate on.make_metricstores e_i·e_j/J, J and ∂J/∂ψ at every equilibrium knot. The decomposition already readsffs.metric(ineigenmode_modes). -
The
wmatrix carries no new information. It is exactly the cofactor matrix ofvdivided by J (I checked all six entries). The contravariant pieces therefore follow from the covariant metric g_ij:- |∇ψ|² = (g₂₂g₃₃ − g₂₃²)/J²
- ∇ψ·∇θ = (g₁₃g₂₃ − g₁₂g₃₃)/J²
- ∇ψ·∇ζ = (g₁₂g₂₃ − g₁₃g₂₂)/J²
-
B² comes from the same metric. The "metric B²" is (χ'/J)²(g₂₂ + 2q g₂₃ + q² g₃₃), so it doesn't need a new exported Equilibrium interpolant (
metric_bsq_interpolant) either.
This is the pattern we most need to avoid in the Julia code base: porting a Fortran routine's private scaffolding verbatim next to an equivalent that already exists.
3. Too much machinery for the size of the calculation
By the PR's own measurement the whole decomposition takes about 0.15 s, but it is built like a performance-critical kernel:
- Five all-
Vectorscratch structs.SurfaceGeometry,CurvatureKernel,SurfaceFields,EffectiveFieldandSurfaceDensitiesare built by reflection (_zeros_struct) and filled field by field. For example,surface_geometry!copies 14 fields out of a NamedTuple one at a time. - Hand-managed threading. Per-
threadid()buffer arrays under@threads :static, and per-thread hint tuples. - Global state changed inside a library function.
BLAS.set_num_threads(1)is set and then restored in afinallyblock. - Fields addressed by
Symbol.getfield/setfield!loops over_ED_*_PROFILESand names built withSymbol(name, :_density). A typo becomes a runtime error, and the loops are type-unstable inside the threaded loop. - A 35-field result struct. It is filled with zero-extent sentinels so two
Booltoggles can switch a form off. Several fields can be derived from others:dW_cumulative*,dW_plasma_reference(a copy ofep[k]), both*_relative_error, and three statistics for one curl check. - Dead code. The repeated-node bookkeeping in
_radial_integralnever runs, becausedecomposition_gridalready returns strictly increasing nodes. - Integrals done twice.
print_energy_decomposition_summaryre-integrates profiles thatdecompose_energyintegrated a few lines earlier. - Validation done twice. Both
run_energy_decompositionanddecompose_energycheck the integrator, n range and kinetic flag.
4. Most of the tests can't fail
These assertions only restate the implementation:
- Sums restated. The component sum equals the total, and
surface_integralequals the same sum written out again. - One integral checked against itself. The last cumulative value equals the total, which comes from the same
cumulative_integratecall. - Densities checked against themselves. Densities re-integrated equal the profiles they were integrated into.
- Arithmetic, not physics.
rms ≤ max, andreference ≈ ep. - Bookkeeping. The name tuples cover the struct fields, and the return value of
mainhas no "energy" key.
The 1e-8 bound on the curl residual is a fair unit test of the b_eff formulas.
But the PR itself notes that (∇×b_eff)·∇ψ = 0 is an algebraic identity on the knots, so it is not a physics diagnostic.
It should not go to gpec.h5, and it should not be in the harness: a 1e-9 noise threshold on a roundoff-level number will mostly raise false alarms.
The tight test is the one described in your own PR body: on the knots, the standard-form integrand equals the DCON Lagrangian density to about 1e-11.
That assertion is independent of the implementation and would catch real bugs.
It isn't in the suite.
What should survive
These parts are good and should carry over:
ForceFreeStates.solution_state_at- widening
compute_node_xi_s!toAbstractVecOrMat - the
reg_spotkeyword refactor ofcompute_contra_displacements - the physics: both forms, K₁/K₂/K₃, the pointwise (not band-truncated) products, and ξ at the knots
- the citations
- your convergence results, which belong in the docs page
Rework plan
The plan below is written the way I would write it for myself (or for Claude Code) before starting.
Treat it as the spec for the rework.
If you disagree with a step, say so in this thread before writing code.
Context
PR #465 ports the Fortran gpout_recon/gpout_recon2 energy reconstruction.
The physics is verified: both forms reproduce ep[k] to 1e-4–1e-2, with the error converging as mpsi and etol are refined.
The integration is wrong:
- a main-pipeline stage and its TOML section
- duplicated geometry (
flux_surface_basis,metric_bsq_interpolant) - about 950 lines of scratch structs, symbol-driven result plumbing and an HDF5 writer
- tests that mostly restate the implementation
Goal: the same physics as a scripted diagnostic in Analysis.
It reads all geometry from ffs.metric, returns plain arrays, and is checked by a few independent tests.
Non-goals:
- no TOML section and no main stage
- no
gpec.h5output and no harness rows - no new Equilibrium exports
- no threading
- no fixed-boundary or forced-response generalization
Target API
using GeneralizedPerturbedEquilibrium
res = GeneralizedPerturbedEquilibrium.main(["path/to/run"]) # forward integrator, vac_flag = true
ed = Analysis.EnergyDecomposition.decompose_plasma_energy(res.ffs; eigenmodes=[1, 2])
Analysis.EnergyDecomposition.plot_energy_decomposition(ed; eigenmode=1)decompose_plasma_energy returns one NamedTuple:
| field | shape | meaning |
|---|---|---|
psi |
(npsi,) |
decomposition nodes (equilibrium knots inside the ξ domain, plus ψ_edge) |
theta |
(mtheta,) |
θ grid, in [0, 1) |
R, Z |
(npsi, mtheta) |
node positions, for 2-D maps |
terms |
NTuple{7,Symbol} |
(:b_eff_sq, :K1, :K2, :K3, :b_sq, :j_b_xi, :divxi_gradp) |
density |
(npsi, mtheta, 7, nmodes) |
pointwise integrands, already scaled by 2μ₀/ψ₀² |
profile |
(npsi, 7, nmodes) |
∮J·density dθ |
dW_effective, dW_standard |
(nmodes,) |
½∫ of the two combinations |
There are no toggles and no sentinels: both forms are always computed, because it costs a fraction of a second.
The densities are always returned.
Store them as ComplexF64, which is about 4 MB per mode at mpsi = 128 and mtheta = 256.
The caller derives everything else: cumulative profiles come from cumulative_integrate, and the relative error is (dW − real(ffs.free_boundary.ep[k])) / abs(ep[k]).
Files
New
src/Analysis/EnergyDecomposition.jl, the submoduleAnalysis.EnergyDecomposition. Target at most about 300 lines including docstrings.test/runtests_energy_decomposition.jl, rewritten. Target at most about 120 lines.
Edit
src/Analysis/Analysis.jl: add the include and the module docstring entry.docs/src/analysis.md: add an@autodocsblock for the new submodule.docs/src/energy_decomposition.md: trim it and keep the physics and convergence tables. Alternatively, fold it intoanalysis.md.
Revert to develop
src/Equilibrium/FluxSurfaceMetrics.jland theEquilibrium.jlexport line._grad_psi_poloidalonly existed to share code withflux_surface_basis.src/GeneralizedPerturbedEquilibrium.jl: removerun_energy_decomposition,energy_decomposition, the stage call and the export.src/PerturbedEquilibrium/PerturbedEquilibrium.jl: remove the include and exports.examples/*/gpec.toml: remove the[EnergyDecomposition]sections.regression-harness/cases/{solovev_n1,diiid_n1}.toml: remove the five new rows.docs/development/hdf5-conventions.mdanddocs/development/architecture.md: revert the PR's one-line additions.
Keep as they are on the branch
ForceFreeStates.solution_state_at(src/ForceFreeStates/Result.jl)compute_node_xi_s!widened toAbstractVecOrMat(src/ForceFreeStates/EulerLagrange.jl)- the
reg_spotkeyword change (src/PerturbedEquilibrium/FieldReconstruction.jl) - the citation additions in
docs/src/citations.mdanddocs/development/references.md
The cleanest history is a fresh branch from develop.
Cherry-pick 66c60af (reg_spot), then bring over solution_state_at and the compute_node_xi_s! change.
Then either force-push feature/energy-decomposition or open a new PR linked from here.
Algorithm
All geometry comes from ffs.metric (make_metric, src/ForceFreeStates/Fourfit.jl).
metric.fs[:, :, 1:6] hold e_i·e_j/J in the order (ψψ, θθ, ζζ, θζ, ζψ, ψθ).
Index 7 is J and index 8 is ∂J/∂ψ.
The docstring labels them g^ij·J; the code computes the covariant quantities.
Confirm this numerically in milestone M0.
-
Nodes. Use the equilibrium knots in
[psi_store[1], psi_store[step]], plus ψ_edge when the last knot falls short of it (keepdecomposition_grid's logic). -
Geometry per node. Read the eight columns at each node.
- At a knot, index the row directly.
- At ψ_edge, the one node off the knots, evaluate a single
CubicSeriesInterpolantalong ψ over the 8 × mtheta columns. That is a native op; write no new geometry function. - Build one series interpolant for all nodes, so there is a single code path (it reproduces the knot rows exactly), and take ∂ψ of the columns from the same object.
- Derive |∇ψ|², ∇ψ·∇θ and ∇ψ·∇ζ from the cofactor identities above.
- Derive B² = (χ'/J)²(g₂₂ + 2q g₂₃ + q² g₃₃). Take ∂ψB² from the series interpolant and ∂θB² from a periodic spline or the spectral derivative on the surface.
- Take R and Z from existing helpers (
grepAnalysis and Equilibrium first) or from therzphinodal values. Add no exported function.
-
ξ and b per eigenmode. Keep the PR's chain unchanged:
solution_state_at, thenel_derivatives!, then thewt[:, k]edge condition, thencompute_node_xi_s!, thenPerturbedEquilibrium.compute_perturbed_field_modesandcompute_contra_displacements(...; reg_spot=0.0). Then inverse-transform to θ withUtilities.FourierTransforms. -
Densities. Write one plain loop over (ψ, θ) that computes the seven terms with local scalars:
- b_eff from the effective-field relations
- K₁, K₂, K₃
- |b|², j·(b×ξ*) and (∇·ξ)*(ξ·∇p)
Use no scratch structs. Allocate per surface where needed. Use no threading and leave BLAS settings alone.
-
Integrals.
profile = sum(density .* J; dims=θ) / mtheta, written withsum/dot, not a helper.
For the totals, useFastInterpolations.integrate(cubic_interp(psi, profile)).
For complex terms, takerealbefore the ψ-integral and say why in one line: the imaginary part integrates to roundoff. -
Plot. Write one function,
plot_energy_decomposition(ed; eigenmode). It shows radial profiles of the seven terms and the two cumulative curves, followingdocs/development/plotting.md.
Comments are one line, with no step numbers.
A deviation from the Fortran gets a single short comment stating the choice made.
Docstrings state what a function computes, not how we got there.
Milestones
Each milestone has an acceptance check.
Commit after each one, using Area - TAG - Imperative summary subjects.
- M0: Pin down the geometry identities (scratch script, not committed).
On the Solovev fixture, compare the PR'sflux_surface_basis(at branch tip) withmetric.fs× J and the cofactor identities, at every knot.
Accept when the maximum relative difference is at most 1e-12 for g_ij, J, ∂J/∂ψ, |∇ψ|², ∇ψ·∇θ, ∇ψ·∇ζ and B².
If any of them fails, stop and report: the plan's premise would then be wrong. - M1: Baseline numbers (scratch, not committed).
RecorddW_plasmaanddW_plasma_standard_formfrom the current branch on the Solovev fixture and onexamples/Solovev_ideal_example, for eigenmodes [1, 2].
These are the reference values for M3. - M2: Strip the branch.
Do the reverts listed under Files.
Accept whenmainon both example decks produces agpec.h5identical todevelop,Pkg.testpasses, and Documenter builds. - M3: The new diagnostic.
AddAnalysis.EnergyDecomposition.
Accept when both totals match the M1 baselines to at most 1e-10 relative difference on both decks, the source file is at most about 300 lines, and there are nogetfield/setfield!loops, nothreadid, and noBLAS.set_num_threadscalls. - M4: Tests (below).
Accept when every assertion could fail for a physics or algebra bug and none of them restates the implementation.
Each testset should run in a few seconds after the fixture is built. - M5: Docs and PR.
Trim the docs page, add the@autodocsblock and update the PR body.
The body should name the rework, and its release note should say "Numerical impact: none", with the harness run againstdevelop.
Every existing row should be unchanged and there should be no new rows.
Request a human reviewer.
Tests (replaces the current file)
Use the existing Solovev fixture (mpsi = 16, forward integrator, vac_flag = true).
Run one main and one decompose_plasma_energy(ffs; eigenmodes=[1, 2]), and share them.
- Standard form equals the DCON Lagrangian, pointwise.
At every knot,profile[:, :b_sq] − profile[:, :j_b_xi] + profile[:, :divxi_gradp]must equal the Euler-Lagrange quadratic form in ξ and ξ′, built from theffs.matsF, K and G splines with the same scaling.
Usertol = 1e-9.
This is the main correctness test, and it is independent of the θ-space code. - Effective form equals standard form, integrated.
dW_effective ≈ dW_standardat a tolerance set from the measured fixture offset with a stated margin (not an arbitrary 0.03), plus a one-line comment citing where the offset comes from: equilibriumetol. - Totals equal
ep[k].
Checkabs(dW_standard[k] − real(ep[k])) / abs(ep[k])against the measured fixture value plus a stated margin. - Curl identity.
Check that (∇×b_eff)·∇ψ / max|∂θb_eff,ζ/J| < 1e-10 at the knots.
It is an algebraic identity, so it is a unit test of the b_eff construction only.
It is not returned or stored. - Shear surface average.
Check ⟨J S⟩ = χ'² q′ (keep the PR's analytic check). solution_state_at.
A saved step must come back unchanged, and integrating from the left and from the right must meet (keep the PR's two checks).- The off-knot geometry node.
At ψ_edge, the series-interpolatedmetric.fsmust match a directrzphi-spline evaluation to the interpolation tolerance.
Drop everything else: shape checks, haskey checks, toggle tests, sum-equals-components checks, and the checks on what main returns.
Risks and open questions
- ψ_edge accuracy. If the series interpolation of
metric.fsat ψ_edge is much less accurate than therzphispline evaluation, test 7 will show it. The fallback is to end the radial integral at the last knot and report how much energy that drops. Decide that from data, not up front. - The
make_metricdocstring mislabels its arrays as g^ij·J. Fixing the docstring is a separate one-lineDOCScommit and should not be folded into this PR. - Size. If M3 grows well past about 300 lines, stop and post the reason here before adding structure.
I know this throws away a lot of code you worked hard on.
The physics and the verification are the valuable part, and they carry over almost entirely.
What changes is the packaging.
I'm happy to go through the plan on a call before you start.
|
|
|
Merged develop, retained the new geometry helpers in |
… ψ from the nearest saved step Add ForceFreeStates.solution_state_at(ffs, psi), which returns the stored basis state of a forward run at an arbitrary ψ. A saved step comes back as stored; otherwise the linear Euler-Lagrange system is integrated to psi from the nearest saved step that no singular surface separates from it, so the state does not depend on how densely the sweep saved. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…enmodes on the (ψ,θ) grid Add Analysis.EnergyDecomposition, a scripted diagnostic that evaluates δW_p of the free-boundary DCON eigenmodes of a forward, ideal, single-n solve in the effective-field form ½∫[|b_eff|²/μ₀ − (K₁+K₂+K₃)|ξ_n|²] and the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)], scaled to compare with ep[k]. decompose_plasma_energy returns a NamedTuple of the seven pointwise integrands, their surface integrals and both totals; plot_energy_decomposition draws the radial profiles and cumulative δW_p. All geometry comes from ffs.metric: the contravariant products follow from the covariant metric by cofactors and B² is the metric field strength. ξ at each knot comes from ForceFreeStates.solution_state_at. It refuses kinetic, multi-n, non-forward and fixed-boundary results. On the same solve it reproduces the former PerturbedEquilibrium stage to 1e-13 in both totals. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…rangian and the plasma energy On the Solovev fixture: the standard-form surface integral equals 2/ψ₀² times the Euler-Lagrange quadratic form on every knot; both totals match each other and ep[k] within the measured quadrature error; the (∇×b_eff)·∇ψ identity, the shear surface average, the off-knot metric and solution_state_at are checked; the plot is drawn. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…rinciple references Add the EnergyDecomposition section to the Analysis page (both forms, the K split, usage, grid and accuracy, Fortran deviations, the mpsi convergence table) with its @autodocs block, and cite Bernstein et al. 1958, Bernstein 1983, Chance et al. 1992 and Boozer 2004 in citations.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…docstring, not the citations page The citations page lists the papers to cite for GPEC itself. Move Bernstein et al. 1958, Bernstein 1983, Chance et al. 1992 and Boozer 2004 out of citations.md and the Analysis page text into a References section of the EnergyDecomposition module docstring, next to the code that uses them. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ismatch of each decomposed eigenmode decompose_plasma_energy now returns and checks two success measures per eigenmode: relative_error_effective/relative_error_standard, the mismatch of each δW_p total against ep[k] (warned above 5 %; it is the ψ quadrature of the knot grid and flags an mpsi that is too coarse), and curl_residual, the largest relative (∇×b_eff)·∇ψ over the knots (warned above 1e-8; an algebraic identity there that checks the b_eff assembly). The tests bound the returned values instead of recomputing the curl, and the Analysis page lists the checks. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
f6f5739 to
fb2ab23
Compare
…tart and normalize its curl check as Fortran does The ψ integrals started at the first saved ODE step. The forward sweep stores that step one ODE step after its start, so the integrals dropped [start, first knot]: 1.25 % of ep for a qlow = 2.5 start (Solovev fixture, mode 2, mpsi 16 to 64). They now start at max(psilow, first profile knot), where initialize_el_at_axis! starts. curl_residual was normalized surface by surface by max|∂θ b_eff,ζ/J|, which inflates it where the mode is weak. It is now the largest residual on the knots over the largest (|∂θ b_eff,ζ| + |∂ζ b_eff,θ|)/J, as in Fortran gpeq_cveri. On the DIII-D-like example with mpsi = 512 that is 2.2e-5 instead of 3.4e-3. The tests could not see the F′ terms or the g_θζ and g_ψζ terms: the Solovev fixture has F′ = 0 and uses PEST coordinates, and 5 of 33 injected errors survived. Changes to the tests: - The fixture now runs in Boozer coordinates. - A DIII-D-like case checks the standard form against the DCON Lagrangian and the curl residual. - A qlow case checks where the integration starts. 4 of the 5 are now caught. The fifth scales (∇·ξ)*(ξ·∇p), which vanishes for DCON displacements. The docs now say what each success check can and cannot detect. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… ODE steps as well as the equilibrium knots The ψ integrals ran over the equilibrium knots only. The knots follow the equilibrium profiles, not the eigenmode, so near the rational surfaces a spline through the knot densities misses the structure of ξ. Against ep, the standard form was off by: - 0.3–0.5 % on the DIII-D-like example, at every knot grid; - 35 % for eigenmode 1 of the LAR reversed-shear example at mpsi = 128; - 5 % on the a10 example at mpsi = 64. On DIII-D the whole error was knot quadrature: Re(ξ†u₂) at the edge equals ep, the jumps at the five rational surfaces sum to 5e-6, and ∫L over the saved ODE steps equals ep to 1e-5. The nodes are now the ODE start and every saved step, as in Fortran gpout_recon, plus the knots, where the geometry is exact and the curl check runs. The standard form now matches ep to: - 1e-5 on DIII-D; - 6e-5 on a10 (mpsi 64); - 1e-4 on TJ (mpsi 256); - 7e-4 on LAR (mpsi 128). A saved step needs no re-integration, so the extra nodes are cheap. The density arrays grow with them; DIII-D goes from 281 to 1444 nodes. Test changes: - The fixture runs at mpsi = 32, since between its own 16 knots the splined metric moves the totals by 3 %. - The success-check tolerances tighten to the new measured values. - A DIII-D check guards the fix. The effective-field total still differs from the standard-form one by a grid-independent amount: −0.4 % on DIII-D and −9.5 % for LAR eigenmode 1. The docs now say so. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Release note
Analysis.EnergyDecomposition.decompose_plasma_energy(ffs; eigenmodes)decomposes the plasma potential energy δW_p of free-boundary eigenmodes on the (ψ,θ) grid, in the effective-field form ½∫[|b_eff|²/μ₀ − (K₁+K₂+K₃)|ξ_n|²] and the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)], scaled to compare withffs.free_boundary.ep[k].plot_energy_decompositiondraws the radial profiles and the cumulative δW_p.Regression report
Run on a laptop (arm64, Julia 1.11.6, pinned manifest 44d8c447, 4 threads), develop @ 728a7e7: diiid_n1 53 unchanged, solovev_n1 22 unchanged, no new rows.
Notes for reviewers
This is the rework @logan-nc asked for in the 2026-10-19 review. The branch was rebuilt from develop; the physics is unchanged, the packaging is new.
mainstage, nogpec.h5group, no harness rows, no example-deck sections. A scripted diagnostic inAnalysis.ffs.metric: the contravariant products by cofactors of the covariant metric, B² = (χ'/J)²(g₂₂ + 2q g₂₃ + q² g₃₃).flux_surface_basisandmetric_bsq_interpolantare gone; Equilibrium is identical to develop. Verified on both Solovev decks thatffs.metricreproduces the formerflux_surface_basis(g_ij, J, ∂J/∂ψ, |∇ψ|², ∇ψ·∇θ, ∇ψ·∇ζ, B²) to ≤1.1e-15 at every knot.BLAS.set_num_threadsorgetfield/setfield!loops. Returns one NamedTuple:psi,theta,R,Z,terms,eigenmodes,density (npsi, mtheta, 7, neig),profile (npsi, 7, neig),dW_effective,dW_standard,relative_error_effective,relative_error_standard,curl_residual.qlowrun;solution_state_at; plot. Of 33 single-term errors injected into the module, the earlier version of these tests missed 5; the current tests catch 4 of them, and the fifth scales (∇·ξ)*(ξ·∇p), which vanishes for DCON displacements. The sum-restating, self-integral, arithmetic and bookkeeping tests are gone.decompose_plasma_energyitself throws on kinetic, multi-n, non-forward and fixed-boundary results.references.mdreferences.mdorcitations.md; Bernstein 1958/1983, Chance 1992 and Boozer 2004 are cited in a References section of the module docstring, as suggested inline.docs/src/analysis.md(physics, usage, Fortran deviations, convergence table); no separate page.gpec.h5or the harness: the δW_p mismatch against ep[k] (warned above 5 %) flags anmpsitoo coarse for the ψ quadrature, which does happen on real grids, and the curl residual (normalized and warned above 1e-8 as in Fortrangpeq_cveri; measured 7e-14 on the test fixture and 4e-12 on the DIII-D-like example) flags an inconsistent b_eff assembly.Kept from the old branch:
ForceFreeStates.solution_state_at(ξ at each knot integrated from the nearest saved step). Not kept: thecompute_node_xi_s!widening — ξ_s is linear in (ξ', ξ), so the module calls the existing matrix signature on the basis and applies the edge combination afterwards.Port check: on the same
ffsthe new module reproduces the formerPerturbedEquilibrium.decompose_energyto ≤1e-13 in both totals and ≤1e-10 in every density (fixture andexamples/Solovev_ideal_example).The ψ integrals run over every saved ODE step as well as the equilibrium knots, and start where the forward sweep starts,
max(psilow, first knot).qlow = 2.5integrating from that step instead had biased mode 2 by −1.25 %.Relative error against ep on the test fixture in its PEST coordinates (
mpert = 8):mpsiAt mpsi = 16 the metric splined between the knots is too coarse. From mpsi = 32 the two forms agree to 4e-5 on this fixture.
Open: the effective form differs from the standard form by a grid-independent amount on some equilibria, while the standard form matches ep:
The gap is unchanged when the integral runs over the saved ODE steps, so it is not quadrature. On DIII-D it also does not change with
mtheta,mpert, the coordinates or the B² source.Open, in ForceFreeStates and not in this PR: when the forward sweep gets very close to a rational surface, the stored states left of the surface carry a spurious component. On
examples/Solovev_ideal_example, Re(ξ†u₂) of eigenmode 1 jumps across the q = 3 crossing by:singfac_mindown to 1e-7eulerlagrange_tolerance = 1e-11singfac_min = 1e-8The decomposition flags both cases and the warning fires. Integrated over the knots only, the same corruption had shown up as a shift of about 1e-3.
🤖 Generated with Claude Code