Skip to content

Commit 8dab0ed

Browse files
d-burgclaude
andcommitted
COILS - BUGFIX - Use psilim on the external-forcing path too, not just the coil path
The coil branch of PerturbedEquilibrium already passed psi=ffs_intr.psilim (PerturbedEquilibrium.jl:98), but the other branch — external forcing data read from file — called convert_forcing_normalization! with no psi, so it fell back to the default surface. That path therefore normalized on the wrong surface silently, with no error, whenever dmlim/qhigh/psiedge moved psilim inward of psihigh (e.g. 0.998902 vs 0.9995 on a DIII-D case: a 6e-4 gap, same order as the psi_N=1 extrapolation this PR already fixed). Both production paths now use psilim, matching Fortran (gpec/gpec.f:431 field_bs_psi(psilim,...)). psilim is the ODE's last grid point to within eps=1e-10 (the final chunk ends at psilim*(1-eps), EulerLagrange.jl:504) and is snapped to psi_store[end] exactly in the dW-peak truncation path (Riccati.jl:1375,1870, EulerLagrange.jl:230). The psihigh default on the three signatures is retained and re-documented as the bare-PlasmaEquilibrium fallback — reachable only by direct callers with no ForceFreeStates solve, for whom psilim does not exist yet — not as the intended control surface. Full test suite passes locally. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
1 parent b2ec04c commit 8dab0ed

2 files changed

Lines changed: 16 additions & 8 deletions

File tree

src/ForcingTerms/CoilFourier.jl

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -81,13 +81,17 @@ extrapolated distance is a fixed 1 − psihigh while the final spline interval s
8181
packing, so the error grows with resolution — which is why this surfaced only once the pass-1
8282
auto grid was refined.
8383
84-
NOTE: strictly faithful is `psilim`, the *integration* limit, not `psihigh`, the *equilibrium
85-
spline* limit. They are equal unless `dmlim`/`qhigh`/`psiedge` truncation fires, in which case
86-
`psilim < psihigh`. PPPL shipped a fix for exactly this confusion (`docs/releases.rst:281`:
87-
"Fixes inappropriate uses of psihigh, which may not be the end of integration psilim if
88-
sas_flag, qhigh, or peak_flag are used"). `psihigh` is used as the default here because it is
89-
all a bare `PlasmaEquilibrium` knows; callers holding a ForceFreeStates solve should pass
90-
`psi = intr.psilim` explicitly, as the benchmarks do.
84+
NOTE ON THE DEFAULT: the physically correct control surface is `psilim`, the *integration*
85+
limit, not `psihigh`, the *equilibrium spline* limit. They are equal unless
86+
`dmlim`/`qhigh`/`psiedge` truncation fires, in which case `psilim < psihigh`. PPPL shipped a fix
87+
for exactly this confusion (`docs/releases.rst:281`: "Fixes inappropriate uses of psihigh, which
88+
may not be the end of integration psilim if sas_flag, qhigh, or peak_flag are used").
89+
90+
**Both production paths pass `psi = ffs_intr.psilim` explicitly** (PerturbedEquilibrium.jl, coil
91+
and forcing-file branches), as do both Fortran-comparison benchmarks. This default of `psihigh`
92+
therefore applies only to direct callers holding a bare `PlasmaEquilibrium` with no
93+
ForceFreeStates solve — for which it is the outermost surface that exists, and always a strict
94+
improvement on extrapolating to ψ_N = 1. Pass `psi` explicitly whenever `psilim` is known.
9195
9296
The coupling surface therefore moves inward by (1 − psihigh): 0.05 % of flux at the default
9397
psihigh = 0.9995, 0.5 % on the example, which lowers psihigh to 0.995 to capture q=6.

src/PerturbedEquilibrium/PerturbedEquilibrium.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -107,8 +107,12 @@ function compute_perturbed_equilibrium(
107107
modes_n = filter(m -> m.n == n, intr.forcing_modes)
108108
isempty(modes_n) && continue
109109
m_vals = [m.m for m in modes_n]
110+
# Same control surface as the coil branch above: psilim, the integration
111+
# limit (Fortran: gpec/gpec.f:431 `field_bs_psi(psilim, ...)`). Without it
112+
# the normalization was taken on the equilibrium-spline limit, which differs
113+
# whenever dmlim/qhigh/psiedge truncation moves psilim inward.
110114
convert_forcing_normalization!(modes_n, norm_tag, equil, n,
111-
minimum(m_vals), maximum(m_vals))
115+
minimum(m_vals), maximum(m_vals); psi=ffs_intr.psilim)
112116
end
113117
end
114118
end

0 commit comments

Comments
 (0)