3D Vacuum Part 2: Passing structs instead of arrays into kernel, cleaning up spline interface, and matrix renaming

[jhalpern30]

This is part of a series of PRs for 3D vacuum that I am breaking up to be more manageable. This PR is entirely aesthetic - there are no modifications to the outputs, but mainly focuses on simplifications/renaming that will carry over to 3D.

TL;DR: kernel! now accepts either a PlasmaGeometry or WallGeometry struct and this led to a VERY nice cleanup, we now actually use the dx_dtheta and dz_dtheta stores in the structs instead of recomputing them, removed the isgn variable, and made some cosmetic renaming of variables. Note that Claude's concerns only relate to the untested pickup code so I am just going to leave it

In approximate order of impact:

kernel! now accepts PlasmaGeometry and WallGeometry structs!

Gone are the days of passing in j1 and j2 tokernel! and then going through some confusing logic to figure out which block of the gradient green's matrix you're in and how that affect the singularities. In Julia, we now pass in the struct,

function kernel!(
    grad_greenfunction::Matrix{Float64},
    greenfunction::Matrix{Float64},
    observer::Union{PlasmaGeometry,WallGeometry},
    source::Union{PlasmaGeometry,WallGeometry},
    n::Int
)

and can use type checking on the observer and source positional argument to obtain everything else out such as

• taking the view of the block of the grad_greenfunction we're in

col_index = (source isa PlasmaGeometry ? 1 : 2)
    row_index = (observer isa PlasmaGeometry ? 1 : 2)
    grad_greenfunction_block = view(
        grad_greenfunction,
        ((row_index-1)*mtheta+1):(row_index*mtheta),
        ((col_index-1)*mtheta+1):(col_index*mtheta)
    )

• whether or not the plasma is the source, which affects whether or not greenfunction even needs to be filled. This is a new feature for both clarify and (possibly) some computational improvement, but not critical
• adding the greenfunction singularity correction term if its the plasma-plasma block
• determining whether to flip the normal direction (more on isgn later)
• what the value of the grad greenfunction singularity residue is

    if is_closed_toroidal # Chance eq. 89
       residue = (observer isa WallGeometry) ? 0.0 : (source isa PlasmaGeometry ? 2.0 : -2.0)

Cleaning up spline interface

I simplified some of the existing logic in the interp_to_new_grid functions and periodic_deriv functions now that we have the Julia splines. The main improvement here comes from the use of the dx_dthetaanddz_dtheta` arrays in the geometry structs. Previously, these arrays were being filled but never touched, since we just rebuilt the splines and derivatives at the beginning of the kernel call

   # Used for Z'_θ and X'_θ in eq.(51)
    # Close the loop for periodic BC by appending first point at the end
    theta_closed = vcat(collect(theta_grid), theta_grid[end] + dtheta)
    x_closed = vcat(x_sourcepoints, x_sourcepoints[1])
    z_closed = vcat(z_sourcepoints, z_sourcepoints[1])
    spline_x = cubic_interp(theta_closed, x_closed; bc=PeriodicBC())
    spline_z = cubic_interp(theta_closed, z_closed; bc=PeriodicBC())
    # Create derivative views once, reuse for all evaluations (avoids allocation per call)
    d1_spline_x = deriv1(spline_x)
    d1_spline_z = deriv1(spline_z)
    theta_vec = collect(theta_grid)
    dx_dtheta = d1_spline_x.(theta_vec)
    dz_dtheta = d1_spline_z.(theta_vec)

Now, we still create the d1 splines since they're needed for the off-grid points during Gaussian interpolation, but since we doing everything on-grid during Simpson integration we can just use the precomputed derivatives here

G_n, gradG_n, gradG_0 = green(x_obs, z_obs, source.x[isrc], source.z[isrc], source.dx_dtheta[isrc], source.dz_dtheta[isrc], n)

Subnote for equal_arc_grid (previously mentioned in #141)

Note that this means that almost this entire section was unused since we were recomputing the derivatives

x_wall, z_wall, _, theta_grid, _ = distribute_to_equal_arc_grid(x_wall, z_wall, mtheta)
theta_grid .= theta_grid .* (2π)  # Scale to [0, 2π) - irregular spacing
# Close the loop for periodic BC by appending first point at the end
theta_closed = vcat(collect(theta_grid), 2π)
x_closed = vcat(x_wall, x_wall[1])
z_closed = vcat(z_wall, z_wall[1])
spline_x = cubic_interp(theta_closed, x_closed; bc=PeriodicBC())
spline_z = cubic_interp(theta_closed, z_closed; bc=PeriodicBC())
d1_spline_x = deriv1(spline_x)
d1_spline_z = deriv1(spline_z)
theta_vec = collect(theta_grid)
dx_dtheta = d1_spline_x.(theta_vec)
dz_dtheta = d1_spline_z.(theta_vec)

and essentially all that's happening here is we're redefining the theta angle to be the equal arc theta angle. So we can use the same derivative calculation for both

if wall_settings.equal_arc_wall && (wall_settings.shape != "nowall")
    @info "Re-distributing wall points to equal arc length spacing"
    !is_closed_toroidal && @error "Wall is not closed toroidally; equal arc length distribution assumes periodicity as cannot be safely used."
    x_wall, z_wall, _, _, _ = distribute_to_equal_arc_grid(x_wall, z_wall, mtheta)
end

# (Same for equal_arc true or false)
dx_dtheta = periodic_deriv(θ_grid, x_wall)
dz_dtheta = periodic_deriv(θ_grid, z_wall)

Removal of the isgn parameter

This closes #133 and adds a nice clarification. As discussed in that issue, this parameter is a sign flip on the direction of the normal to account for it being the outward unit normal. In the clockwise theta convention of VACUUM, this points toward increasing psi, which ends up pointing out of the vacuum region for the wall and into the vacuum region for the plasma surface. This is why there is a factor of -1 multiplying the grad_greenfunction matrix. I added a much more clear comment on where this comes from, and also condensed it into a single broadcasted matrix multiply, instead of multiplying every single term that gets added during integration.

Renaming/clarifications

These are truly cosmetic.

• Making the offsets within the fourier transforms more explicit via COS_COL_OFFSET and the like + better comments on this
• grri -> green_fourier
• greenfunction_temp -> green_temp
• grad_greenfunction_mat -> grad_greenfunction
• and so on

[jhalpern30]

@claude review this PR

[claude]

Claude finished @jhalpern30's task in 1m 50s —— View job

---

Code Review Complete ✓

Summary

This PR improves the Vacuum module code by:

1. Passing structs into kernel function instead of individual arrays - excellent improvement for clarity
2. Cleaning up spline interface - consolidated and renamed functions
3. Renaming matrices for clarity (grad_greenfunction_mat → grad_green, greenfunction_mat → green_temp, grri → green_fourier)
4. Fixing basis function naming (sin_ln_basis/cos_ln_basis → sin_mn_basis/cos_mn_basis) to match actual usage

Overall, these changes significantly improve code readability and maintainability. The refactoring is well-executed with only minor issues noted below.

---

Detailed Findings

✅ Strengths

1. Excellent struct-based refactoring in kernel! (src/Vacuum/Kernel2D.jl:75-226)

   • Passing PlasmaGeometry/WallGeometry structs instead of individual arrays is much cleaner
   • Type dispatch (observer isa PlasmaGeometry) is more idiomatic Julia than integer flags
   • Block indexing logic is now clearer with explicit calculations

2. Proper basis function naming (src/Vacuum/DataTypes.jl:44-45, 161-162)

   • Changed sin_ln_basis/cos_ln_basis → sin_mn_basis/cos_mn_basis
   • Correctly reflects that these are sin(mθ - nν) and cos(mθ - nν), not sin(lθ - nν)
   • Variables now use m (poloidal mode) consistently with physics notation

3. Vectorized basis computation (src/Vacuum/DataTypes.jl:161-162)

   sin_mn_basis = sin.((mlow .+ (0:(mpert-1))') .* θ_grid .- n .* ν)
   cos_mn_basis = cos.((mlow .+ (0:(mpert-1))') .* θ_grid .- n .* ν)

   This is much more concise and performant than the nested loops

4. Better matrix naming (src/Vacuum/Vacuum.jl:235-244)

   • grad_greenfunction_mat → grad_green (clearer, shorter)
   • grri → green_fourier (more descriptive of what it contains)
   • Added clear row/column offset constants for readability

5. Simplified initialize_wall (src/Vacuum/DataTypes.jl:213-222)

   • Early return for nowall case reduces nesting
   • Removed complex equal arc length spline logic, now uses consistent periodic_deriv

6. Function consolidation (src/Vacuum/MathUtils.jl:62-108)

   • interp_to_new_grid now has clearer signature and documentation
   • Renamed periodic_cubic_deriv → periodic_deriv (simpler, implementation detail hidden)
   • Removed unused offset parameters from old interp_to_new_grid

---

⚠️ Issues & Suggestions

1. Missing mhigh field in VacuumInput (CRITICAL)

Location: src/Vacuum/MathUtils.jl:483

The _calculate_potential_chi function references inputs.mhigh:

l_modes = (inputs.mlow:inputs.mhigh)

But VacuumInput struct (src/Vacuum/DataTypes.jl:18-28) doesn't have an mhigh field - only mlow and mpert. This will cause a runtime error.

Fix needed: Either add mhigh to VacuumInput or compute it as mlow + mpert - 1. The variable l_modes appears to be unused anyway (no subsequent references), so you might just delete line 483.

Fix this →

2. Inconsistent function naming (MINOR)

Location: src/Vacuum/MathUtils.jl:500

The function uses cos_ln_basis and sin_ln_basis:

g_real[i_theta, l_idx] = aval * plasma_surf.cos_ln_basis[i_theta, l_idx]
g_imag[i_theta, l_idx] = aval * plasma_surf.sin_ln_basis[i_theta, l_idx]

But these fields were renamed to cos_mn_basis and sin_mn_basis in the PlasmaGeometry struct. This will cause a field not found error.

Fix needed: Change cos_ln_basis → cos_mn_basis and sin_ln_basis → sin_mn_basis on lines 500-501.

Fix this →

3. Potential performance issue with nonsing_idx allocation (OPTIMIZATION)

Location: src/Vacuum/Kernel2D.jl:126

nonsing_idx = mod1.((j+2):(j+mtheta-2), mtheta)

This allocates a new vector on every iteration of the outer loop (mtheta times). For large mtheta (e.g., 256), this could be wasteful.

Suggestion: Pre-allocate this outside the loop and reuse with in-place operations, or use a view if possible. However, this is a micro-optimization and not critical.

4. Magic number for convergence check (CLARITY)

Location: src/Vacuum/MathUtils.jl:392

grid_info[i, j] = (fintjj > 0.1) ? 1 : 0

The threshold 0.1 appears arbitrary and lacks documentation. What does this check represent physically?

Suggestion: Extract to a named constant like INTERIOR_THRESHOLD = 0.1 with a comment explaining what it represents (presumably a winding number check).

5. Unused l_modes variable (CLEANUP)

Location: src/Vacuum/MathUtils.jl:483

As mentioned in issue #1, the l_modes variable is computed but never used. Remove it.

6. Comment style inconsistency (STYLE)

Throughout the changes, you use both Unicode symbols (θ, ν, 𝒢) and ASCII equivalents in comments. While Unicode is great for readability, ensure it's used consistently. For example:

• src/Vacuum/Kernel2D.jl:84 uses "dtheta" in code but "θ" in comments
• src/Vacuum/DataTypes.jl:151-152 consistently uses θ in both

Suggestion: Consider standardizing on Unicode in comments since Julia supports it well and it improves physics readability.

---

Testing Recommendations

Since this is a significant refactoring:

1. Run existing vacuum tests: Verify test/runtests_vacuum_julia.jl passes
2. Check against Fortran: Run test/runtests_vacuum_fortran.jl to ensure parity
3. Full integration test: Run test/runtests_fullruns.jl to verify no downstream breakage
4. Consider benchmarking: The vectorization changes should improve performance - worth measuring

---

Adherence to CLAUDE.md Standards

✅ Commit message format: "VACUUM - BUGFIX - no wall struct error" and refactoring commits follow the CODE - TAG - Message format

⚠️ No step numbering: Good - no "Step 1, Step 2" comments present

✅ Code comments: Generally clear and physics-focused

---

Action Items

Before merging, please address:

1. ❗ CRITICAL: Fix inputs.mhigh reference (line 483 in MathUtils.jl) - see issue Port over splines from fortran to Julia #1
2. ❗ CRITICAL: Fix cos_ln_basis/sin_ln_basis references (lines 500-501 in MathUtils.jl) - see issue Fix/ffi2 #2
3. Run test suite to catch any other field name inconsistencies

Optional improvements:

• Address magic number threshold (issue Refactor EQUIL #4)
• Consider pre-allocating nonsing_idx (issue internal eval and integral added. #3)
• Standardize Unicode usage in comments (issue Equil - NEW - the initial work's from Jaymyoung Lee #6)

---

Overall Assessment

Code Quality: ⭐⭐⭐⭐½ (4.5/5)

This is high-quality refactoring work that significantly improves code clarity. The struct-based approach is a big win for maintainability. The two critical bugs (#1 and #2) need fixing before merge, but these appear to be simple oversights during the rename process. Once fixed, this will be a solid improvement to the codebase.

Recommendation: Request changes to fix the two critical field reference errors, then approve after verification.

---

[jhalpern30]

Verification that I didn't break anything

Note that my changes are slightly more memory efficient (Solovev ideal example)
this branch: 2.695 s (12724201 allocations: 1.38 GiB)
develop: 2.624 s (12743413 allocations: 1.44 GiB)

Tests on develop

DIIID-like Ideal Example

Computing free boundary energies
[ Info: Using no wall
Least Stable Eigenmode Energies:
  Plasma = -3.471e-02 +1.912e-05i
  Vacuum = +1.742e+00 -1.172e-16i
  Total  = +1.707e+00 +1.912e-05i
All free-boundary modes stable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 1.637e+00 seconds
Normal termination.

Solovev Ideal Example

Computing free boundary energies
[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
[ Info: Re-distributing wall points to equal arc length spacing
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 +1.966e-17i
  Total  = -4.372e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 3.046e+00 seconds
Normal termination.

Solovev Ideal Example (equal_arc_wall = false)

Computing free boundary energies
[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 +3.297e-17i
  Total  = -4.373e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 3.203e+00 seconds
Normal termination.

Tests on this branch

DIIID-like Ideal Example

Computing free boundary energies
[ Info: Using no wall
Least Stable Eigenmode Energies:
  Plasma = -3.471e-02 +1.912e-05i
  Vacuum = +1.742e+00 +4.181e-17i
  Total  = +1.707e+00 +1.912e-05i
All free-boundary modes stable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 1.475e+00 seconds
Normal termination.

Solovev Ideal Example

Computing free boundary energies
[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
[ Info: Re-distributing wall points to equal arc length spacing
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 -1.033e-17i
  Total  = -4.372e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 2.693e+00 seconds
Normal termination.

Solovev Ideal Example (equal_arc_wall = false)

Computing free boundary energies
[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 -2.558e-17i
  Total  = -4.372e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 2.906e+00 seconds
Normal termination.

[jhalpern30]

@logan-nc might have a problem here, both my unit tests and docs test are failing and I think its due to FastInterpolations

[logan-nc]

How bad of a problem? like the tests weren't properly undated and try to use Fortran that isn't there? or like the tests have caught a serious and mysterious bug in the spline implementation somewhere in the source code?

[logan-nc]

How bad of a problem? like the tests weren't properly undated and try to use Fortran that isn't there? or like the tests have caught a serious and mysterious bug in the spline implementation somewhere in the source code?

[jhalpern30]

How bad of a problem? like the tests weren't properly undated and try to use Fortran that isn't there? or like the tests have caught a serious and mysterious bug in the spline implementation somewhere in the source code?

I think the former but I'm not 100% sure. It looks like it makes the unit tests fail due to some equilibrium issues, and the docs fail because there are spline docstrings that don't exist?

[jhalpern30]

@copilot fix the unit test and doc issues that are causing the checks on this PR to fail

[Copilot]

@jhalpern30 I've opened a new pull request, #150, to work on those changes. Once the pull request is ready, I'll request review from you.

[jhalpern30]

UPDATE

I made distribute_to_equal_arc_grid use FastInterpolations as well, so now everything in VACUUM except the Lagrange interpolation within the Gaussian quadrature uses FastInterpolations.jl in a consistent fashion. This let me remove the lagrange1D function, and now this PR will officially close #141 since all comments I had there are now addressed.

Tests stay the same with this update:

Computing free boundary energies
[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
[ Info: Re-distributing wall points to equal arc length spacing
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 -1.561e-19i
  Total  = -4.372e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to euler.h5
----------------------------------
Run time: 2.876e+00 seconds
Normal termination.

[jhalpern30]

@claude please summarize the changes in this PR made AFTER commit 376271c

[claude]

Claude finished @jhalpern30's task in 2m 6s —— View job

---

Summary of Changes After Commit 376271c

Todo List:

• Find commit 376271c in git history
• List all commits after 376271c
• Analyze the diff of changes
• Summarize key changes

---

After commit 376271c (the merge from develop), this PR includes 5 additional commits that make significant improvements to the Vacuum module. These changes removed 704 lines while adding only 277 lines, achieving substantial code simplification.

---

Commits Summary (in chronological order)

1. 162e2f8 - VACUUM - BUGFIX - fixing interp_to_new_grid

Fixed handling of periodic data that may or may not include endpoint duplication. Updated compute_greensfunctions_only interface.

Files changed: src/Vacuum/MathUtils.jl, src/Vacuum/VacuumFromEquilibrium.jl

---

2. 550f04f - VACUUM - IMPROVEMENT - removing compute_green_only functionality

Eliminated the separate compute_green_only interface and replaced it with a boolean flag compute_green_only in compute_vacuum_response(). This consolidates the two code paths into one.

Impact: Removed ~97 lines of duplicate logic

Files changed:

• src/ForceFreeStates/Free.jl (+3)
• src/PerturbedEquilibrium/SingularCoupling.jl (refactored calls)
• src/Vacuum/DataTypes.jl (removed flag)
• src/Vacuum/Vacuum.jl (unified interface)
• src/Vacuum/VacuumFromEquilibrium.jl (-124 lines, major cleanup)

---

3. bc61e08 - VACUUM - IMPROVEMENT - consolidating vacuum input generation

Unified the vacuum input generation interface between Free.jl and SingularCoupling.jl to use the same constructor pattern.

Key changes:

• Added VacuumInput(equil, ψ, mtheta, mpert, mlow, n; force_wv_symmetry) constructor
• Removed duplicate vacuum input setup code from ForceFreeStatesStructs.jl
• Simplified VacuumFromEquilibrium.jl initialization logic

Impact: Removed 176 lines of duplicate code

Files changed:

• src/ForceFreeStates/ForceFreeStatesStructs.jl (-57 lines, removed duplicate struct)
• src/ForceFreeStates/Free.jl (-81 lines)
• src/PerturbedEquilibrium/SingularCoupling.jl (updated to use new constructor)
• src/Vacuum/DataTypes.jl (added constructor)
• src/Vacuum/MathUtils.jl (simplified)
• src/Vacuum/VacuumFromEquilibrium.jl (major simplification)

---

4. 3ebfa96 - bugfix

Minor kernel bugfix (1 line change)

Files changed: src/Vacuum/Kernel2D.jl

---

5. d6971a7 - VACUUM - IMPROVEMENT - ruthless removal of unnecessary parameters and functions

This is the largest cleanup commit, refactoring geometry initialization functions as constructors and removing unused parameters.

Major changes:

1. Removed unused fields from VacuumInput:

   • Removed qa (safety factor) - not needed in vacuum calculations
   • Removed mtheta_eq - redundant with mtheta

2. Simplified PlasmaGeometry struct:

   • Removed delta field (was computed from ν but never actually needed)
   • Removed dx_dtheta and dz_dtheta storage (derivatives computed on-the-fly when needed)
   • Now only stores: x, z, ν

3. Refactored initialization functions as constructors:

   • initialize_plasma_surface() → PlasmaGeometry(inputs) constructor
   • initialize_wall() → WallGeometry(inputs, plasma_surf, wall_settings) constructor
   • Much cleaner, more Julian style

4. FourierTransforms simplification:

   • Changed from qa and delta parameters to direct ν parameter
   • Simplified from explicit loops to vectorized broadcasting
   • Renamed cslth/snlth → cos_mn_basis/sin_mn_basis (more descriptive)
   • Formula changed from cos(m*θ + n*qa*δ) to cos(m*θ - n*ν) (eliminated intermediate delta computation)

5. Removed helper functions:

   • Eliminated several intermediate helper functions in MathUtils.jl (-74 lines)
   • Streamlined initialization logic in VacuumFromEquilibrium.jl (-121 lines)
   • Simplified test code in runtests_vacuum_julia.jl (-47 lines)

Impact: Net removal of 159 lines, major simplification

Files changed (9 files):

• src/ForceFreeStates/Free.jl
• src/PerturbedEquilibrium/SingularCoupling.jl
• src/Utilities/FourierTransforms.jl (vectorized, removed qa/delta)
• src/Vacuum/DataTypes.jl (struct simplification + constructors)
• src/Vacuum/Kernel2D.jl
• src/Vacuum/MathUtils.jl (removed helpers)
• src/Vacuum/Vacuum.jl
• src/Vacuum/VacuumFromEquilibrium.jl (major cleanup)
• test/runtests_vacuum_julia.jl (simplified tests)

---

6. b0c07b8 - TEST - BUGFIX - fixing unit testing exception throw error

Fixed error handling in DataTypes.jl for unit tests.

Files changed: src/Vacuum/DataTypes.jl (+9, -4)

---

Overall Impact

Lines changed: +277 insertions, -704 deletions across 10 files

Key improvements:

1. Eliminated redundant code paths - unified compute_green_only functionality
2. Removed unnecessary parameters - qa, delta, mtheta_eq, stored derivatives
3. Simplified data structures - PlasmaGeometry now minimal (only x, z, ν)
4. Constructor pattern - More Julian style with PlasmaGeometry(inputs) and WallGeometry(...)
5. Vectorized Fourier transforms - Replaced loops with broadcasting
6. Unified interface - Same vacuum input generation for ForceFreeStates and PerturbedEquilibrium

These changes maintain numerical correctness (as verified by your tests showing identical results) while significantly improving code clarity and maintainability.

---

[jhalpern30]

@logan-nc ok this is now ready for review again. Its a bit bigger than originally, but I think its worth it due to the significant cleanups I found from the perturbed equilibrium PR. Claude's summary above gives a decent idea, I think the most notable are that I got rid of the interp_to_new_grid and periodic_deriv functions, incorporated @mgyoo86's new periodic exclusive endpoint BCs from #148 into the VACUUM code (which clean things up very well), removed a lot of unnecessary parameters from the VacuumInput struct, and renamed initialize_... functions to just be constructors in a more Julian style.

I think one of the bigger things that remain here is more file organization, reserved for a later commit. MathUtils.jl now only contains distribute_to_equal_arc_grid and the compute_vacuum_field subfunctions (which don't seem to appear in the initial perturbed equilibrium PR). Furthermore, the VacuumFromEquilibrium.jl file only contains extract_plasma_surface_at_psi. Honestly, I am thinking this function could go into the Equilibrium Module - thoughts on this? There's nothing in it that is vacuum specific, and I think would play nicer in an area where the user can more easily look up what the rzphi splines are and how that related to what we're extracting

[jhalpern30]

Results from this branch:
DIIID example (JPEC stuff is garbage but proving I didn't break its ability to run through)

Computing free boundary energies
[ Info: Using no wall
Least Stable Eigenmode Energies:
  Plasma = -3.471e-02 +1.912e-05i
  Vacuum = +1.742e+00 -1.248e-16i
  Total  = +1.707e+00 +1.912e-05i
All free-boundary modes stable for n = 1.
Writing saved data to jpec.h5
...
    Surface 1: q = 2.0, ψ = 0.5511327570855618
      Resonant mode: m = 2, n = 1
      Delta': -1.7364381747522152e16
[ Info: Using no wall
[ Info: Using no wall
[ Info: Using no wall
  Singular coupling calculation complete
    Max island half-width: 5.662e+06
    Max Delta': 1.736e+16
----------------------------------
Run time: 9.180e-01 seconds
PERTURBED EQUILIBRIUM COMPLETE
Writing perturbed equilibrium data to examples/DIIID-like_ideal_example/jpec.h5
  Perturbed equilibrium output complete

Solovev example

[ Info: Calculating conformal wall shape 8.52e-02 m from plasma surface.
[ Info: Re-distributing wall points to equal arc length spacing (assumes closed, toroidal wall).
Least Stable Eigenmode Energies:
  Plasma = -1.102e+01 +2.883e-06i
  Vacuum = +1.058e+01 +8.995e-18i
  Total  = -4.372e-01 +2.883e-06i
Free-boundary mode unstable for n = 1.
Writing saved data to jpec.h5

[logan-nc]

Beautiful

[claude]

Code Review

This PR successfully implements a useful architectural improvement by passing geometry structs to kernel! instead of raw arrays. The changes eliminate confusing index logic and make the code more maintainable. However, there are several critical issues that must be addressed before merge.

Critical Issues - Must Fix Before Merge

1. Missing dx_dtheta and dz_dtheta fields in geometry structs

The PR description states the code now uses precomputed derivatives stored in the structs. However, DataTypes.jl:103-136 shows PlasmaGeometry and WallGeometry only contain x, z, and ν fields. Yet Kernel2D.jl:128 references source.dx_dtheta[isrc] and source.dz_dtheta[isrc], which will cause a field-not-found error at runtime. The structs must be updated to include dx_dtheta and dz_dtheta arrays.

2. Missing is_closed_toroidal field in WallGeometry

Kernel2D.jl:210 uses is_closed_toroidal = true as a hardcoded variable, but Vacuum.jl:310 references wall.is_closed_toroidal. Either add this field to WallGeometry or use a consistent approach.

3. Broken pickup code references non-existent fields

MathUtils.jl:156 in _calculate_potential_chi references plasma_surf.dx_dtheta and plasma_surf.dz_dtheta, but these fields don't exist. If this code is unused (per PR description), either fix it or clearly mark it as deprecated/broken.

Design Issues

4. Inconsistent offset naming

Vacuum.jl:279-282 defines PLASMA_ROW_OFFSET = 0, WALL_ROW_OFFSET = mtheta, etc. The zero-valued offsets add no clarity. Either use all offsets consistently or eliminate redundant ones.

5. FourierTransforms documentation uses old terminology

The module docstring and several function docstrings still reference qa*delta phase shifts, but the implementation now uses ν (toroidal angle offset). Update all documentation to match the new API.

6. Removed compute_vacuum_inputs is a breaking change

Free.jl removes the compute_vacuum_inputs function. This is a breaking API change for any external callers. Either keep a deprecated wrapper or document the migration path.

Code Quality

7. Spline derivatives still recreated on each kernel! call

Kernel2D.jl:112-115 creates spline derivative functors on every kernel! invocation. If dx_dtheta and dz_dtheta are precomputed in the geometry structs (per PR goals), these spline setups should only be needed for off-grid Gaussian quadrature points. Consider whether this can be optimized further.

8. Test coverage reduced significantly

runtests_vacuum_julia.jl has 60 deletions. Verify that removed tests were truly redundant and that remaining coverage is adequate for the new struct-based interface.

9. Comment in Kernel2D.jl:127 could be clearer

The comment says derivatives are "needed for off-grid points during Gaussian interpolation" but doesn't specify which integration method uses off-grid vs on-grid points. Be more explicit.

Minor Issues

• SingularCoupling.jl:120: Use keyword argument syntax consistently - VacuumInput constructor call mixes styles
• DataTypes.jl:348: Verify all distribute_to_equal_arc_grid call sites updated for new 2-value return signature
• Kernel2D.jl:204-206: The normal direction sign correction is clearer now but could use a reference to which equation in Chance 1997 this corresponds to

Positive Aspects

• The struct-based kernel! interface is a clear improvement over index-based logic
• Removing isgn parameter and replacing with type-based dispatch is cleaner
• Consolidating VacuumInput constructor into DataTypes.jl is better organization
• FourierTransforms.jl cleanup using broadcasting is more idiomatic Julia
• Better separation of concerns between geometry setup and kernel computation

Recommendations

1. BLOCKER: Add dx_dtheta and dz_dtheta fields to PlasmaGeometry and WallGeometry, populate in constructors
2. BLOCKER: Fix is_closed_toroidal handling - either add to WallGeometry or use consistent approach
3. BLOCKER: Fix or clearly deprecate the pickup field calculation code in MathUtils.jl
4. Update FourierTransforms documentation to use ν consistently instead of qa*delta
5. Run full test suite including vacuum tests to catch any runtime field access errors
6. Consider adding a deprecated compute_vacuum_inputs wrapper for API compatibility

The architectural direction is sound and will benefit 3D vacuum development. However, the missing struct fields will cause runtime failures and must be fixed before merge.