Implement Bernstein's breakdown for the perturbed potential energy

[parkjk]

Summary

This PR adds post-processing diagnostics to GPEC that break the perturbed plasma potential energy δW_p of a DCON free-boundary eigenmode into Bernstein's terms, and compares each reconstruction with the DCON eigenvalue. It also adds a MATCH option, sol_num, that selects which DCON eigenfunction MATCH reconstructs.

Scope (definition of done): the four input flags below, their output files, the MATCH sol_num option, the derivation note in docs/tex/recon/, a Sphinx docs entry and one regression example. Every new flag defaults off, so default runs are unchanged. Further SNU work goes on its own purpose-scoped branches off develop.

The Julia port of the same diagnostics is OpenFUSIONToolkit/GPEC#465 ([EnergyDecomposition], options effective_field_form / standard_form). The flag names here match those options.

What is added

GPEC (gpec_output namelist, input/gpec.in)

Flag	Default	What it computes	Main routines
effective_field_flag	f	Effective-field (Bernstein) form 2μ₀δW_p = ∫J[|C|² − μ₀K|ξ_n|²] with C = Q + ξ_n μ₀ j×n̂ and K = K₁ + K₂ + K₃ (shear-current, parallel-current, pressure-curvature). Reports |C|² split into ψ/θ/ζ, the K terms, and the reconstructed δW_p against DCON ep(mode)	gpout_recon, gpeq_c, gpeq_shear, gpeq_curvature, gpeq_K, gpeq_epf, gpeq_dst
standard_form_flag	f	Standard (first) form |Q|²/μ₀ − j·(Q×ξ*) + (∇·ξ)*(ξ·∇p), with γp|∇·ξ|² omitted (Eq. (3.1) of docs/tex/recon/main.tex)	gpout_recon2, gpeq_firstform
effective_field_terms_flag	f	|C|² = |Q|² + 2Re(V*·Q) + |V|² reported as A, B, C∥, I⊥, with the identity check epf = A + B + C∥ + I⊥	gpout_recon3, gpeq_epfterms
curl_check_flag	f	Surface-by-surface check that (∇×C)·∇ψ = 0, judged on the residual relative to its cancelling terms (tolerance 1e-8)	gpeq_cveri
recon_out	t	Write the diagnostic files and gpec.log entries; f keeps only the terminal summary
recon_int	"spline"	Radial quadrature on the psifac grid for all three routes: "spline" or "trapezoid"

Usage: mode_flag=t, mode=<eigenmode>, coil_flag=f, reg_flag=f (and data_flag=f, harmonic_flag=f). reg_flag defaults to true, so it must be set explicitly. GPEC stops at input with a clear message if any of the three energy flags is set otherwise, because the comparison with ep(mode) is only meaningful for a bare DCON eigenmode.

Terminal / gpec.log summary: each route ends with its δW_p total, the same total in DCON normalization (dW_p(recon-normalize) and analogues), DCON ep(mode), and the ψ range it integrated over next to the DCON range:

  recon psi range      =   1.00000000E-004  9.91539010E-001
  DCON psi range       =   1.00000000E-004  9.91539010E-001

Outputs (per selected mode, sol<mode>): gpec_recon_{terms,integration,k,k_sigma_check,shear_fun,curvature,cwfun,cwmn,cvfun,cvmn,cv2fun,cv2mn,qvmn,cveri}_sol*.out, gpec_recon2_{terms,integration}_sol*.out and gpec_recon3_{terms,integration}_sol*.out, with R–Z resolved term densities, Fourier tables and cumulative radial integrals.

MATCH (input/match.in)

sol_num (default 1) selects which column of Wt (1 = most unstable) ideal_build uses as the edge condition. sol_num=1 writes solutions.bin as before; sol_num=n>1 writes solutions<n>.bin. It is bounds-checked against 1..mpert, mripple is range-checked, and the ripple path is unchanged.

Documentation

• docs/outputs.rst: an "Energy Decomposition Outputs" section (flags, preconditions, output files, printed summary) and a "MATCH Outputs" entry for solutions<n>.bin; it links docs/tex/recon/main.tex.
• docs/namelists.rst: includes input/match.in.
• docs/tex/recon/main.tex (with references.bib), "Bernstein–Chance decomposition of the perturbed potential energy in GPEC", is the derivation behind every routine: coordinate and basis conventions; the transformation of δW_F into the Bernstein form in DCON/SI units; the derivation of K; both reconstruction schemes; a code-variable ↔ symbol table; the exact integrand behind every logged quantity, with cross-checks between the three routes.

Regression example

docs/examples/regression_solovev_recon_example (Soloviev, inputs copied from regression_solovev_ideal_example, delta_mlow=delta_mhigh=8) runs dcon and gpec with effective_field_flag=t and curl_check_flag=t. check_recon.py passes when the recon ψ range equals the DCON range, dW_p(recon-normalize) matches DCON ep(1) within 2%, and the curl check reports PASS. Both CI jobs (gfortran and ifort) run it after the existing Soloviev examples.

Validation

• Clean make (gfortran 15, macOS arm64) with no new warnings: the same 16 as the merge base (type mismatches in stride/zvode1.f and xdraw linker paths). The regression example also passes in a build with the CI debug flags (-fcheck=all -finit-real=snan …).

• docs/examples/DIIID_ideal_example with default inputs gives the same outputs as the merge base with develop (902faac), apart from the imaginary part of e_P for high modes in gpec_response_n1.out, which also differs between two runs of the unmodified code (pre-existing nondeterminism). develop has moved 155 commits since then (including DCON/EQUIL/VACUUM changes), so a comparison against current develop differs for reasons outside this PR. The branch merges into develop without conflicts.

• Recon case on DIII-D (coil_flag=f, mode_flag=t, mode=1, reg_flag=f, all four flags on): the guard, rename, local-scratch, NEWUNIT and psi-range commits leave all 18 gpec_recon*.out files bit-identical; only the curl-check verdict lines (new relative criterion) and the two new psi-range lines change.

• Agreement with DCON ep(1):

  Case	effective field	standard form	effective-field terms	curl check (rel)
  DIII-D (mpert 34)	−0.91%	−0.23%	−0.91%	1.1e-10 PASS
  Soloviev example (mpert 32)	+0.41%	+0.41%	+0.42%	2.8e-13 PASS

  With the default 8 harmonics the Soloviev effective-field form misses by 21%: C is built in real space and Fourier truncated, so the example widens the band.

• recon_int: spline and trapezoid differ by ≤0.1% in all cases above, with no consistent winner, so the default stays "spline" (as in the Julia port). The remaining 0.4–0.9% is an edge-localized offset of the form itself, not quadrature (see #465).

Response to review

✳ = written with Claude Code (see "AI assistance" below).

Request	Change	Commit
What is docs/nametag? (2025-10)	Removed	985c69e
Document new flags and sol_num (2025-10)	input/gpec.in and input/match.in document every new input; the abandoned bernstien_flag is gone	6b27832, 0101e98, ccf446c, 6e2c2c2 ✳
Make the notes available or cite them (2025-10)	docs/tex/recon/main.tex is in the repository	0e1f875 ✳ and earlier DOCS commits
Finite objective; drop leftover diffs (2026-06)	input/dcon.in, dcon/dcon.F, equil/equil_out.f and docs/nametag no longer differ from develop; "most unstable" wording made consistent; sol_num/mripple bounds-checked	a49c79a ✳, 2c03f04 ✳, 3f2d03f ✳
Guard the recon flags against mode_flag/coil_flag/reg_flag (2026-10-03, required)	Stop at input unless mode_flag=t and reg_flag, coil_flag, data_flag, harmonic_flag are off; warn if curl_check_flag is set without effective_field_flag	c27437c ✳, b57c4f1 ✳
Division guards (2026-10-03, required)	rfac, |∇ψ|² and B² divisors floored at a private recon_tiny in the recon kernels; MAX(.,0) inside every recon SQRT(rzphi%f(1))	c27437c ✳, (next commit) ✳
Rename the numbered recon_flag# flags (2026-10-03)	recon_flag → effective_field_flag, recon_flag2 → standard_form_flag, recon_flag3 → effective_field_terms_flag, cveri_flag → curl_check_flag; no aliases. Theory note retitled	c27437c ✳, 13b97d8 ✳
Dev-history comments, non-ASCII in .f, gpeq_shear step header	Removed or fixed	c27437c ✳
Globals in the recon kernels	jac, jac1, rfac, eta, w, v are local in every recon kernel, and chi1 is computed locally in gpeq_firstform/gpeq_shear	4a70b91 ✳
Hard-coded units 72–99	All recon files opened with NEWUNIT=	47be71b ✳
cveri_tol is an absolute threshold on a dimensional quantity	Verdict on max|res| / max((|∂θC_ζ| + |∂ζC_θ|)/|J|), tolerance 1e-8	fe40fed ✳
Integration range may not match DCON's; decide the default recon_int	Each route prints its ψ range beside the DCON psilow–psilim range; default stays "spline" after the comparison above	df990fc ✳
Sphinx docs entry	docs/outputs.rst, docs/namelists.rst	aa88e64 ✳
Regression example	docs/examples/regression_solovev_recon_example + CI steps	2616e3f ✳
Update the PR body	This text	✳

AI assistance

Parts of this PR were written with Claude Code (Anthropic) and reviewed and committed by @jaesun57:

• All review-response commits from 2026-10: c27437c, 4a70b91, 47be71b, fe40fed, 13b97d8, b57c4f1, df990fc, aa88e64, 2616e3f (and the follow-up commit with the remaining SQRT guards).
• Earlier commits that carry a Co-Authored-By: Claude trailer, among them 3b14c8a, 6e2c2c2, 0e1f875, a49c79a, 2c03f04, 94547e9, 58173ef and 3f2d03f, plus several recon output and DOCS commits; git log --format='%h %(trailers:key=Co-Authored-By,valueonly)' develop..snu_offshoot lists them all.
• The validation runs above and this PR description.

[logan-nc]

What is this file?

[logan-nc]

Please document things as you add features and flags

[logan-nc]

Please document

[logan-nc]

Are these notes available to everyone? Can you add them to the repo or cite a publication?

[logan-nc]

@parkjk This branch appears to have completed a clear and finite feature addition. I recommend it be cleaned up no and merged into the develop branch for everyone to enjoy. Not doing so risks sliding into a never ending state of adding features but never finishing them rigorously - which could lead to a permanent split of the SNU code from the rest of the GPEC source.

[logan-nc]

@jaesun57, is this branch ready to merge? If not, what is needed and what it the expected timeline? I'd like to not let it live too much longer, as the mission will creep. PR's are best if finite so they do not diverge overmuch from the develop branch.

[jaesun57]

@logan-nc I'm mainly working on the recon branch right now, and it will take about two more weeks to finish. In the meantime, I'll make sure to keep both recon and snu_offshoot branches updated with the develop branch so they don't diverge too much. As for snu_offshoot, since I didn't create it myself, I'll need to check with @parkjk to see what else needs to be done there.

[logan-nc]

Thanks for the response @jaesun57.

What is the recon branch? Please open a Pull Request for it and make a clear statement of the intent as the opening comment. You can mark the PR and "draft" and that tells people not to bother reviewing the code yet, but it is important to signal what you are doing as it can impact what others are doing (we don't want multiple people accidentally doing the same thing or conflicting things). It looks like you branched off of snu_offshot - is your intent to merge back into snu_offshoot or directly into develop?

Similarly, it would be great if you could coordinate with the SNU team to edit the opening comment of this branch with a description of the finite objective, which will clarify whether it is ready to review and merge or not. We don't want infinitely diverging parallel develop branches.

Thanks!

[logan-nc]

Here is AI's summary of the branch:

Branch objective (finite)

Add a selectable solution-component parameter sol_num to MATCH, so the
ideal-MHD reconstruction can target any DCON ideal eigenfunction (column of
Wt) rather than only the most-unstable one.

What it does: sol_num (namelist match.in, default 1) selects which
column of Wt is used as the edge boundary condition in ideal_build, and
writes the result to solutions.bin (for sol_num=1) or solutions{sol_num}.bin
otherwise. With the default it reproduces existing behavior exactly, so the
change is backward-compatible.

Definition of done: It is ready to merge once
the diff is reduced to just that feature (plus its input template) and the
1 = least stable vs. 1 = most unstable wording is made consistent (see Q2).
Further SNU work belongs on its own purpose-scoped branch/PR, not here.

---

Questions needing clarification (@parkjk / @jaesun57)

1. Is sol_num the intended objective? It is the only substantive change
   surviving in the net diff (develop...snu_offshoot). The heavily-developed
   bernstein.f module (shear / geometry / curvature / metric calculations) was
   built up across many commits and then entirely deleted in 5ce5f7bb.

2. sol_num index description The new input/match.in comment says
   "1 is least stable" (good), but the match/match.f header comment says
   column 1 is the "most unstable ideal solution" (should be changed to match "least stable" terminology).

3. Leftover cruft The net diff also carries a few changes that
   look like leftovers from the abandoned experiments, unrelated to sol_num.
   Confirm whether these should be dropped before merge:

   • input/dcon.in: kin_flag t → f
   • dcon/dcon.F: comment xi' → xi (looks like a typo regression)
   • equil/equil_out.f: whitespace-only changes
   • docs/nametag: new branch-marker file

[parkjk]

Your points are good and well taken, though this touches on a big picture conversation we should have separately (a bit too much to get into here). Anyway, @jaesun57 has already committed to opening a PR within the next two weeks, so I'd ask for a little patience on that front. Once it's up, she'll follow the practice closely.

[jaesun57]

@logan-nc ,All addressed and pushed to snu_offshoot(now at ce90394)!

Q1 – sol_num: Yes, this is the objective. sol_num (match.in, default 1)
selects which DCON ideal eigenfunction (column of Wt) MATCH reconstructs —
so we can target the 2nd / 3rd / … solution instead of only the most-unstable
one. Default 1 reproduces existing behavior exactly. The old dcon/bernstein.f
(shear/geometry experiment) was an abandoned line, deliberately removed in
5ce5f7b — not part of the objective.

Q2 – index wording (fixed in 311922c): Standardized on "most unstable" in
both match.in and match.f:54, since sol_num=1 is the lowest-δW / most
-unstable column. (Went with "most unstable" rather than "least stable" as it's
the unambiguous description of column 1; happy to flip if you prefer the other
term.) Also fixed a "calculates" → "calculated" typo.

Q3 – leftovers (reverted to develop in 08ba8bd):

• input/dcon.in kin_flag f → t — restored to develop.
• dcon/dcon.F xi → xi' — restored the dropped prime (typo regression).
• equil/equil_out.f — whitespace reverted to develop.
• docs/nametag — already removed in 7ffde74; nothing left to drop.

All four leftover diffs vs develop are now zero.

That said, I still have a few other in-progress changes to wrap up, so I'll do
the final consolidation pass once those are all done.

[logan-nc]

Thanks @jaesun57 !

@parkjk this looks like a clean feature add rn. Our options are,

1. Update the PR name and merge this branch to develop → SNU continues to make new feature branches off develop as desired (e.g. recon → develop).
2. Branch a match_sol_num branch off on snu_offshoot now and merge that into develop as the finite feature, leave the snu_offshoot as an open ended sandbox for now and mark this PR as draft to it doesn't pop up when I am checking for stale PRs that need to get merged/worked-on.

Either way, I'm happy to chat about the big picture separately via slack, zoom, etc. or in person Aug 10.

[parkjk]

Option 1 sounds good. @jaesun57: Would this be a representing PR name at the moment?

"Implement Bernstein's breakdown for the perturbed potential energy"

or would you reserve it? Let me know the PR name otherwise and tell me when. Then I will change the name, merge pull, and let you go from there.

[jaesun57]

Would be a perfect name @parkjk and would be great to just close the pr now. And I'll keep merging regularly from now on to prevent the branch from drifting too far from develop. Thanks for the feedback, @logan-nc!

[jaesun57]

Hi @logan-nc
This is now ready for review. Since the June thread the branch consolidated into the full scope the title describes — Bernstein's breakdown of δW_p as GPEC post-processing diagnostics, with the MATCH sol_num selector riding along as discussed.

What I have done:

• recon_flag: reconstructs the effective perturbed field C = Q + V and the destabilizing kernel K (shear, parallel-current and curvature terms), integrates J|C|²/μ0 − JK|ξ_n|² over ψ, and compares the reconstructed 2δW_p against the DCON eigenvalue for the selected mode.
• recon_flag2: first-form plasma energy |Q|²/μ0 − j·(Q×ξ*) + (∇·ξ)*(ξ·∇p), each term evaluated directly (γp|∇·ξ|² omitted).
• recon_flag3: Bernstein |C|² decomposition into A, B, Cpar, Iperp with the independent identity check epf = A + B + Cpar + Iperp.
• cveri_flag: optional surface-by-surface verification that (curl C)·∇ψ = 0.
• Outputs: R-Z resolved term densities, Fourier tables, and cumulative radial integrals (spline or trapezoid via recon_int).
• Documentation: docs/tex/recon/main.tex writes up the full derivation — covariant/contravariant basis and coordinate conventions, the transformation of δW_F into the Bernstein form in DCON/SI units, the K derivation with explicit μ0 and its DCON-coordinate expression, the two reconstruction schemes (real-space on the (ψ,θ) grid and the radial matrix form of δW), a code-variable ↔ symbol mapping table, and the exact integrand behind every GPEC_RECON/2/3 log quantity with cross-checks between the three routes.
• MATCH: sol_num (default 1) selects which DCON eigenfunction ideal_build reconstructs, writing solutions.bin; bounds-checked, ripple path unchanged.

All new flags default off, so default runs are unchanged. The self-review comment above lists the cleanup items and known minor points.

[logan-nc]

@claude review this pr

[claude]

Claude finished @logan-nc's task in 1m 55s —— View job

---

Review of PR #213: Bernstein breakdown of δW_p (recon_flag*) and MATCH sol_num

I read the diff for the MATCH, input and GPEC changes and the core recon routines: gpeq_prep_c, gpeq_c, gpeq_cveri, gpeq_firstform, gpeq_epf, gpeq_dst, gpeq_shear, gpeq_K and gpout_recon. I skimmed gpout_recon2/3 and gpeq_epfterms. I did not build or run anything (no NetCDF toolchain here), and I did not check the physics in main.tex against the code.

origin/develop and this branch share no merge base in the checkout, so I compared the trees directly. Any file develop has touched since will show extra drift.

Overall

The scope is much better than the June state. Everything new is behind flags that default off, and the MATCH change is backward compatible. The recon code is well commented, and the Jacobian-weighting conventions are explained where they are used. I found no correctness bugs in the parts I read. The points below are mostly documentation, robustness and process.

MATCH (match/*, input/match.in)

• ✅ sol_num is bounds-checked against 1..mpert. mripple is now range-checked too. The ripple path forces ideal_write(1), so the output name stays solutions.bin. The default reproduces the old behavior.
• ✅ The wording is consistent ("most unstable") in match.in and match.f.
• Minor: the sol_num comment in input/match.in says "calculated by MATCH". It could instead say the output goes to solutions<n>.bin for n>1, since that naming is the user-visible effect.

GPEC recon

1. Energy comparison only makes sense for a DCON eigenmode. gpout_recon sets ep_index = mode only IF (mode_flag). Otherwise it silently compares against ep(1).
   • With coil_flag=t (or reg_flag=t) that comparison is meaningless, but it prints as "DCON ep(...)".
   • The gpec.in comment mentions coil_flag=f, mode_flag=t, reg_flag=f, but the code doesn't enforce it.
   • Suggest warning or stopping in gpec.f when a recon flag is set without mode_flag or with coil_flag or reg_flag.
2. Integration range may not match DCON's.
   • The main loop exits at psi > psilim, but the energy integral in gpout_recon runs over psifac points restricted to [rzphi%xs(0), rzphi%xs(mpsi)] (gpout.f ~7205–7243).
   • If rzphi%xs(mpsi) is below psilim, or xs(0) is above psifac(0), the "recon vs DCON" difference includes a truncation error. The output doesn't say so.
   • I'd print the integration limits next to the comparison.
3. Spline integration across singular surfaces. recon_int="spline" (the default) fits a natural cubic spline through epf/dst values sampled on the psifac grid. The integrands carry the 1/(m−nq) behavior near rational surfaces, and the solution is piecewise across them. Splines can ring around those jumps.
   • Consider making "trapezoid" the default.
   • Or document the expected difference, e.g. one line in gpec.in.
4. Division guards (per CLAUDE.md).
   • rfac=SQRT(rzphi%f(1)) is divided by in gpeq_c, gpeq_shear, gpeq_K and gpeq_firstform (e.g. jac*rfac, 1/(2*rfac*jac)).
   • shear_temp%fs(:,5) = fs(2)/fs(3) divides by |∇ψ|², and gpeq_K divides by bsq_val.
   • These are fine away from the axis, but the loop starts at ipsi=0. A cheap MAX(...,tiny) or an explicit check would match the repo guidelines.
   • gpeq_curvature uses SQRT(rzphi%f(1)) with no MAX(.,0), while gpeq_shear does use it.
5. Module-level scratch variables. The new routines assign the shared globals jac, jac1, rfac, eta, w and v (from gpglobal). They are not locals, apart from chi1/q/p1, which are shadowed locally.
   • It works today, because idcon_metric and the other routines recompute these before use.
   • It is fragile, though, and gpeq_firstform relies on the global chi1 being set in idcon. Please declare local copies, or at least note the dependency in the header comment.
6. Repeated work. gpeq_K (which calls gpeq_shear, including a spline fit per surface) runs once in the main loop and again inside gpeq_dst when recon_out=t, then once per psifac point in the energy pass. This is fine for a diagnostic. Reusing the first call's K_term* would roughly halve the cost.
7. cveri tolerance. cveri_tol = 1e-10 is an absolute threshold on a dimensional quantity, so PASS/WARN depends on units and normalization. Normalizing by max|∇·C| or by |C|/L would make the verdict meaningful.
8. Hard-coded Fortran units 72–99. The comment says they are reserved, and the routines run sequentially, so there is no clash today. Using NEWUNIT= or a shared parameter list would avoid collisions if gpout_recon2/3 ever run in the same call.

Style and process

• ✅ Executable lines are within 72 columns. The over-length lines are comments only, so they are harmless. Some contain UTF-8 (κ·∇ψ, ∂r/∂ξ_j), which gfortran accepts but ifort and others may not. Plain ASCII is safer in .f.
• Several comments describe dev history ("(xno/bno and theta-space c2v* were dev-time checks; removed)" in gpeq_c, "From recon_metric:" in gpeq_shear). These should go before merge.
• gpeq_shear Step 4 header says "Transform to mode space via FFT" but the step only deallocates.
• The PR description is still the placeholder text ("…the develop branch."). Per the June thread, it should state the finite objective, the new flags and the files touched, so reviewers know what "done" means.
• Docs: input/gpec.in is documented, which addresses the earlier review comment. But docs/ has no user-facing entry, such as the output tables in docs/outputs.rst or docs/releases.rst. That covers the new files (gpec_recon_*_sol<mode>.out) and the sol_num option. docs/tex/recon/main.tex is a great derivation reference, but it is not linked from the Sphinx docs. A 1.1 MB PNG in docs/tex/recon/ is OK but worth confirming.
• No regression test or example. I'd add a short run, e.g. in docs/examples, that exercises recon_flag=t and checks that dW_p(recon-normalize) matches DCON ep within a stated tolerance. The PR currently relies on manual verification.
• Please confirm a clean build with make (the DEPENDENCIES.inc change adds gpglobal.o as a dependency of gpout.o) and that default runs are bit-identical on one existing example.

Suggested before merge

1. Update the PR body.
2. Enforce or warn on the mode_flag/coil_flag/reg_flag preconditions for recon.
3. Print the integration limits next to the DCON comparison, and decide on the default recon_int.
4. Add division guards.
5. Strip the dev-history comments and the UTF-8 in comments.
6. Add a doc entry and one regression example.
   · branch snu_offshoot

[logan-nc]

I don't love the recon_flag# convention. It is not very descriptive for new users and introduces a new numbering scheme inconsistent with other flags. It would be better to name each one something more descriptive

[logan-nc]

@jaesun57 please do at least the flag guards and division guards in the above. The others seem like reasonable suggestions as well.
Note the fragility of the globals is always a big issue in the Fortran code, which is one of the reasons the Julia verzion is so much nicer now!

[logan-nc]

@jaesun57 what do you think of the above requests? Can you let me know if you plan to implement any? Or do you want to advocate for things to be merged just as they are?

[jaesun57]

@logan-nc Sorry for late replies. I updated codes with all of the requests above and also agree that some flag guards and division guards are needed. Could you take another look?