Skip to content

VACUUM - Migrate surface current calculation into the VACUUM module - #345

Open
jhalpern30 wants to merge 16 commits into
developfrom
feature/vacuum-surface-inductance
Open

VACUUM - Migrate surface current calculation into the VACUUM module#345
jhalpern30 wants to merge 16 commits into
developfrom
feature/vacuum-surface-inductance

Conversation

@jhalpern30

@jhalpern30 jhalpern30 commented Jul 30, 2026

Copy link
Copy Markdown
Collaborator

This PR moves the surface inductance calculation (or at least the majority of it) into the VACUUM module. The surface inductance calculation is fundamentally a vacuum quantity. It was hacked into GPEC to avoid cross-code changes, but can more simply be integrated into the VACUUM module and passed into the PerturbedEquilibrium module.

This also makes it a boolean if it is computed or not. For the standard stability calculation where only Wv is needed, computing the surface current calculation costs an extra matrix inversion + storage since it requires forming and solving the interior kernel, which can add up especially when making the Wv spline in the edge. I didn't benchmark the speedup in 2D since its not massive, but this will be very helpful for 3D.

I have looked at the Fortran and checked that no not-yet-ported code will require the old way of doing this. And if anything new does, its more suitable to be coded up in the VACUUM module instead anyway.

Some notes on specific changes:

  • the grri/grre matrices are now internal to the VACUUM module only. The only quantity the rest of the code sees is the surface current matrix I_v
  • The machinery for compute_surface_current_from_greens is now drastically simplified
  • I removed the regularization logic surrounding the formation of the surface inductance matrix in SingularCoupling.jl. I looked at the Fortran and this was unfounded, and screamed "Claude did this because it had one failing test case while porting over the code and added try/catch logic". This same type of code appeared when we were doing the edge matrix splining for psiedge < 1 logic.
  • I've added more documentation on where code that handles the difference in theta orientation between GPEC and VACUUM appears and what it does
  • I've moved the Hermitization of the surface inductance into the VACUUM module

@jhalpern30
jhalpern30 marked this pull request as ready for review August 3, 2026 15:11
@jhalpern30

Copy link
Copy Markdown
Collaborator Author

Regression tests - this only really affects PerturbedEquilibrium so I ran DIIID

Regression Report: diiid_n1
================================================================================================================
Ref 1: local  @ local (2026-08-03)
Ref 2: develop  @ 167ada76 (2026-07-30)
----------------------------------------------------------------------------------------------------------------
Quantity                                      local            develop          Diff               Status       
----------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        7.923042e-01     7.923042e-01     0.0e+00            OK           
total energy Im(et[1])                        1.261466e-04     1.261466e-04     0.0e+00            OK           
plasma energy Re(ep[1])                       -1.375573e+00    -1.375573e+00    0.0e+00            OK           
vacuum energy Re(ev[1])                       2.167877e+00     2.167877e+00     0.0e+00            OK           
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-01     0.0e+00            OK           
plasma energy (all)                           [35 elem]        [35 elem]        0.0e+00            OK           
vacuum energy (all)                           [35 elem]        [35 elem]        0.0e+00            OK           
total energy (all)                            [35 elem]        [35 elem]        0.0e+00            OK           
ODE steps (saved)                             1353             1353             0.0e+00            OK           
ODE steps (total)                             1990             1990             0.0e+00            OK           
q0                                            1.204202e+00     1.204202e+00     0.0e+00            OK           
q95                                           4.781724e+00     4.781724e+00     0.0e+00            OK           
beta_t                                        1.327082e-02     1.327082e-02     0.0e+00            OK           
beta_n                                        1.372570e+00     1.372570e+00     0.0e+00            OK           
internal inductance li1                       8.842261e-01     8.842261e-01     0.0e+00            OK           
internal inductance li2                       7.080767e-01     7.080767e-01     0.0e+00            OK           
internal inductance li3                       7.304351e-01     7.304351e-01     0.0e+00            OK           
poloidal beta betap1                          6.681028e-01     6.681028e-01     0.0e+00            OK           
poloidal beta betap2                          5.350080e-01     5.350080e-01     0.0e+00            OK           
poloidal beta betap3                          5.519015e-01     5.519015e-01     0.0e+00            OK           
# singular surfaces                           5                5                0.0e+00            OK           
singular psi locations                        [5 elem]         [5 elem]         0.0e+00            OK           
singular q values                             [5 elem]         [5 elem]         0.0e+00            OK           
current beta betaj                            4.236772e-01     4.236772e-01     0.0e+00            OK           
plasma volume                                 1.829476e+01     1.829476e+01     0.0e+00            OK           
plasma current                                1.152130e+00     1.152130e+00     0.0e+00            OK           
mpert                                         35               35               0.0e+00            OK           
npert                                         1                1                0.0e+00            OK           
toroidal field bt0                            2.006573e+00     2.006573e+00     0.0e+00            OK           
wall field bwall                              3.880145e-01     3.880145e-01     0.0e+00            OK           
aspect ratio                                  2.845746e+00     2.845746e+00     0.0e+00            OK           
elongation kappa                              1.708322e+00     1.708322e+00     0.0e+00            OK           
q profile (checksum)                          6e5ed3ba20f0...  6e5ed3ba20f0...  identical          OK           
pressure profile (checksum)                   ba7ac449bf30...  ba7ac449bf30...  identical          OK           
Mercier D_I profile (checksum)                414eb69dd7de...  414eb69dd7de...  identical          OK           
resistive interchange D_R profile (checksum)  1978c3923c51...  1978c3923c51...  identical          OK           
ballooning Delta' profile (checksum)          908bad88b3c5...  908bad88b3c5...  identical          OK           
delta prime (BVP diagonal)                    [5 elem]         [5 elem]         0.0e+00            OK           
island half-widths                            [5 elem]         [5 elem]         1.5e-13            OK           
Chirikov parameter                            [5 elem]         [5 elem]         1.3e-12            OK           
||resonant area-weighted field||              3.953877e-04     3.953877e-04     6.4e-16            OK           
PE plasma energy                              3.444911e+00     3.444911e+00     3.3e-11            OK           
PE vacuum energy                              3.174509e+00     3.174509e+00     1.3e-11            OK           
PE surface energy                             5.962637e+00     5.962637e+00     1.196e-10 (0.00%)  ** CHANGED **
PE toroidal torque                            -5.031207e-02    -5.031207e-02    2.9e-13            OK           
NTV torque FGAR [N·m]                         5.772805e-01     5.772805e-01     1.3e-11            OK           
NTV kinetic energy dW FGAR [J]                7.235890e-02     7.235890e-02     1.6e-12            OK           
Runtime (s)                                   180.2s           170.1s                              --           
resonant area-weighted field b^r              [5 elem]         [5 elem]         2.1e-15            OK           
================================================================================================================
Summary: 1 changed, 47 unchanged

The small changes are due to slight variations in the order that Hermitization is performed. The 1 changed value is not physical, since it is just slightly above the arbitrary limit of 1e-10 that we set

@jhalpern30

Copy link
Copy Markdown
Collaborator Author

@claude review

@claude

claude Bot commented Aug 3, 2026

Copy link
Copy Markdown

Review: #345 — Migrate surface current calculation into VACUUM module

Verified the core refactor algebraically: the interior-kernel solve now uses A - 2I (subtracting 2 from the un-negated grad_green diagonal) instead of the old -A + 2I, and the downstream combination correspondingly switched from grri + grre to grre - grri. Working through both, grre_new - grri_new == grre_old + grri_old, so the new I_v-based path is algebraically equivalent to the removed compute_surface_inductance_from_greens — this is a real refactor, not a silent behavior change. The compute_Iv boolean is a sensible perf win: it avoids the interior solve + Green's-function storage on the common Wv-only path.

Issues

  1. VacuumData struct is missing the I_v field it's documented to need (src/ForceFreeStates/ForceFreeStatesStructs.jl). The PR removes grri/grre from VacuumData but doesn't add I_v, even though compute_vacuum_response!'s docstring (src/Vacuum/Vacuum.jl) lists I_v::AbstractMatrix{ComplexF64} as part of the required duck-typed contract "when compute_Iv=true". Today nothing breaks because the only call site against VacuumData (Free.jl:56) never passes compute_Iv=true, but the moment the in-place path is wired up for 3D (which the PR description flags as the actual motivation) it will throw has no field I_v. Worth adding the field now, or narrowing the docstring to say VacuumData doesn't support the I_v path yet.

  2. Test gap that hides Port over splines from fortran to Julia #1: the "in-place compute_vacuum_response! matches wrapper" test in test/runtests_vacuum.jl uses a hand-built NamedTuple (with I_v already present) rather than an actual VacuumData, so it wouldn't catch the struct/docstring mismatch above. Also, calc_surface_inductance (new public function in ResponseMatrices.jl) has no direct unit test — only the DIIID regression case (reported in PR comments) exercises it end-to-end. A small targeted test (e.g. checking the μ₀(2π)² normalization / inv() relation on a synthetic I_v) would pin this down independent of the full pipeline.

  3. Stale docs: docs/development/architecture.md (~lines 49, 133-134, 175) still describes grri/grre as the Vacuum module's outputs and as per-singular-surface data. This file wasn't touched by the PR and is now inconsistent with the I_v/compute_Iv design — docs/src/workflow.md and docs/src/vacuum.md were correctly updated, this one was missed.

  4. Minor: broken comment in SingularCoupling.jl: "Extract the (m_res, m_res) diagonal element is needed for singflx" reads as a leftover half-edit of "only the (m_res,m_res) diagonal element is needed" — should be reworded (e.g. "Only the (m_res, m_res) diagonal element is needed for singflx").

Regression harness

Good to see this was already run and posted (diiid_n1, local vs develop) per the CLAUDE.md mandate — 47/48 quantities identical, one (PE surface energy) differs by 1.2e-10 absolute, attributed to Hermitization-ordering changes and explained as noise. That's a reasonable result. Since the SingularCoupling diagonal-coupling path (island half-widths, Chirikov parameter) is one of the more direct consumers of the new I_vcalc_surface_inductance route, it'd be worth a second look at whether any existing/new case stresses that path more than diiid_n1 does, but nothing in the report looks alarming.

Removed try/catch regularization in surface inductance

Agree with removing the ad hoc regularization/try-catch fallback in the old compute_surface_inductance_from_greens. Per the PR description this had no Fortran basis and was a one-off patch for a specific failing case. Fine to let inv(I_v) fail loudly if it's ever singular rather than silently returning a placeholder.


No pull request should be merged without human review. Per this repo's CLAUDE.md, this is non-negotiable — please make sure a third-party human reviewer signs off before merging into develop, regardless of any automated review feedback (including this one).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant