TokaMaker: Native jphi-linterp Ip hold (via #267) + bootstrap scaler simplification#1
Open
d-burg wants to merge 18 commits into
Open
TokaMaker: Native jphi-linterp Ip hold (via #267) + bootstrap scaler simplification#1d-burg wants to merge 18 commits into
d-burg wants to merge 18 commits into
Conversation
Consolidates post-PR#195 improvements to bootstrap.py: - Add prominence filter to edge spike detection - Fix degenerate curve_fit bounds in analyze_bootstrap_edge_spike - Tighten Ip scale factor bracket from ±10% to ±1% - Add verbose flag to solve_with_bootstrap and find_optimal_scale - Move log Lambda calc inside relevant logic - Various bug fixes and diagnostics Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
- Replace `Itor_target` with `Ip_target` throughout code
OpenFUSIONToolkit#248) and related improvements
Address the well-known ~0.5-1% Ip discretization mismatch between
jphi_update's flux-averaged normalization (gs_flux_int) and
gs_comp_globals' physical-Ip integration (R*P' + 0.5*FFP/R).
Disagreement is structural -- the two formulas integrate
different quantities, agreeing only when <R>*<1/R> = 1 pointwise
(circular plasma). For D-shaped DIII-D plasma, ~0.76% mismatch.
Fix: outer iteration in tokamaker_solve wrapper, gated on
jphi_flux_func profile type detection. After each inner Picard
solve, compute Ip_phys via gs_comp_globals, compare to original
Itor_target, scale Itor_target by a damped multiplicative
correction, re-Picard. Continues until rel_err < 5e-4 or
max_outer (5) iterations.
Key design choices:
* jphi_update remains stateless from caller perspective --
avoids breaking caller-side iterative wrappers (eg. bouquet's
find_optimal_scale, _corrective_jphi_iteration) that depend
on solve being a deterministic FFP-from-jphi mapping for a
given equilibrium
* Damped correction (^0.7) + hard clamp [0.7, 1.3] for safety
* Only activates for jphi_flux_func profiles -- PP'/FF' users
unaffected
* Skipped for vacuum solves (no Ip target)
Validation (DIII-D 204441 @ 4.4s recon):
* Ip error: -0.7570% -> -0.0396% (19x improvement)
* PIN_JPHI draw 1 bnd-diag post-Ip-align: 3.68 mm -> 0.00 mm
* Recon's [Ip match] secant converges identically (deterministic
callable behaviour preserved)
* No build regressions
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…loop User-reported regression: with the outer loop in place, recon flow showed ~2% Ip drift after many internal solves. Cause: outer loop scales gseq%Itor_target multiplicatively but never restores it. Each subsequent mygs.solve() call starts from the previous call's final (inflated) Itor_target -- cumulative inflation across the ~dozens of solves a typical recon flow performs. Fix: save Itor_target magnitude at entry, restore at exit (preserving whatever sign jphi_update last set). The FFP is calibrated for the inflated target during the final inner Picard, so the equilibrium correctly has Ip ~ user_target; only the bookkeeping target is reset. Future mygs.solve() calls start from the clean user-set target and re-run the outer loop if needed. Validation (PIN_JPHI=1, sigma=0, N=3 draws): * Ip-align all 3 draws: actual_Ip == target (no-op corrections) * bnd-diag post-Ip-align: 0.00 mm (all draws) * bnd-diag post-homotopy: 0.39 mm (just iso-update + homotopy noise) * F-drift 0.00%, VSC drift 0.02% on every draw * 3/3 IN_SPEC * No cumulative drift across draws Before fix (original pre-outer-loop OFT): PIN_JPHI sigma=0: 3.69 mm RMS, Ip err -0.76% After fix (outer loop + restore): PIN_JPHI sigma=0: 0.39 mm RMS, Ip err +0.005% -- ~10x and 150x Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…d ip_phys Guard the multiplicative Itor_target correction against a corrupted ip_phys returned by gs_comp_globals after a nominally-successful (ierr=0) inner solve. Three failure modes, all observed under tight coil-bound homotopy at low j_phi pinning: ip_phys <= 0 (degenerate), ip_phys NaN/Inf (singular flux state), or |ip_phys| outside [0.3, 3.0]x target (QP went singular, plasma collapsed to device center, integral is bogus). Applying a correction from that state inflates Itor_target to nonsense and corrupts downstream solves; instead bail without correction and let the caller's rollback recover. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…j0 only The GS solve now holds Ip to target natively (jphi-linterp cut-cell fix + Ip-correction outer loop), so the find_optimal_scale(find_j0=False) Ip-scale secant in solve_with_bootstrap is redundant -> final_scale_Ip = 1.0. Strip the now-unused find_j0=False branch, get_Ip_error, and the find_j0/scale_j0 parameters from find_optimal_scale; it is now a clean core-j0 matcher. Callers updated. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
This PR updates TokaMaker’s jphi-linterp workflow to (1) hold plasma current Ip to the requested target during the Grad–Shafranov solve via a Fortran-side outer correction loop, and (2) simplify the Python bootstrap workflow now that Ip alignment no longer requires a downstream secant rescaling.
Changes:
- Add a
jphi-linterp-specific outer iteration intokamaker_solveto iteratively correctItor_targetso that physicalIp(fromgs_comp_globals) matches the original target. - Simplify
find_optimal_scaleto a single “core-j0matching” secant path and remove the Python-sideIprescaling secant insolve_with_bootstrap. - Improve edge-spike peak selection in
analyze_bootstrap_edge_spikeby adding a prominence filter and preferring the tallest peak in the far edge region.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 2 comments.
| File | Description |
|---|---|
src/python/wrappers/tokamaker_f.F90 |
Adds a jphi-linterp-only outer correction loop in tokamaker_solve to converge realized Ip to the user’s target. |
src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py |
Removes redundant Python Ip rescaling in bootstrap solve flow and simplifies scale optimization to core-j0 matching. |
Comment on lines
+210
to
+213
| # Find peak in the edge region (prominence filter suppresses noise | ||
| # oscillations that can be misidentified for small bootstrap spikes) | ||
| min_prominence = 0.05 * numpy.max(j_edge) if numpy.max(j_edge) > 0 else 0.0 | ||
| peaks, properties = find_peaks(j_edge, height=0., prominence=min_prominence) |
Comment on lines
+724
to
+730
| IF(ip_phys_bad)THEN | ||
| IF(oft_debug_print(1))WRITE(*,'(A,I0,A,2ES16.8)') & | ||
| ' [JPHI_IP] iter ', outer_it, & | ||
| ' SAFETY BAIL -- ip_phys out of sane range, ip_phys/target =', & | ||
| ABS(ip_phys), ip_phys_target | ||
| EXIT ! caller will see equilibrium as solver produced it; restore handled below | ||
| END IF |
- Fortran (identical to OpenFUSIONToolkit#290 update c68d383): signal failure on the jphi-linterp Ip safety-bail via error_str so a corrupted equilibrium no longer returns as a nominal success; fix the bail debug label. - bootstrap.py `analyze_bootstrap_edge_spike`: guard an empty edge slice (psi_N may not reach psi_N>=0.7) so numpy.max no longer raises; compute the edge max once.
Owner
Author
|
Addressed in 21fcd91:
|
# Conflicts: # src/physics/grad_shaf.F90 # src/python/wrappers/tokamaker_f.F90
…rop our outer-loop, keep bootstrap.py simplification
d-burg
changed the title
Owner
Author
|
@hansec I have some small cleanups in the wake of your Ip fix, which this builds on. The jphi-linterp Ip fix ( Actual changes in this PR (
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
TokaMaker: Native jphi-linterp Ip hold (via OpenFUSIONToolkit#267) + bootstrap scaler simplification
Summary
When TokaMaker solves with a
jphi-linterpprofile (caller specifies⟨j_φ⟩,solver back-solves
FF'), the integrated plasma current drifts fromIp_targetby ~0.5–1 % due to a cut-cell quadrature mismatch near the separatrix. This branch
makes the GS solve hold
Ipto target natively (~0.0004 %) and, now that thenative hold is reliable, removes the redundant Python
Ip-rescaling secant fromsolve_with_bootstrapand simplifiesfind_optimal_scaleto a core-j_0-onlyscaler.
This branch now contains two pieces:
gs_itor_nldamped relative-norm update + pressure-update ordering)bootstrap.py)What's included
jphi-linterp Ip handling (Fortran) — from OpenFUSIONToolkit#267:
jphi_updatenormalizes against the physical current (gs_itor_nl) under adamped multiplicative update relative to the previous norm
(
jphi_norm = (1 + Ip_target/itor_nl)·norm_last/2), which converges where theearlier absolute-estimate relaxation oscillated.
P%updateis moved ahead of the current integral (andp_scaleset) so thecurrent is computed against an up-to-date pressure profile.
Itor_target→Ip_targetrename andIp_target_skipflag.Bootstrap scaler simplification (
bootstrap.py) — this branch:solve_with_bootstrap: drop thefind_optimal_scale(find_j0=False)Ip-scalesecant →
final_scale_Ip = 1.0(Ip is now held by the solve). The core-j_0scale (
find_j0=True) is unchanged.find_optimal_scale: remove the now-unusedfind_j0=Falsebranch,get_Ip_error, and thefind_j0/scale_j0parameters — now a clean core-j_0matcher.
analyze_bootstrap_edge_spike: guard an empty edge slice sonumpy.maxnolonger raises.
Why
Downstream (the
bouquetperturbed-equilibrium workflow) carried several PythonIp-alignment secants purely to compensate for the jphi-linterp drift. With the
native hold those are redundant; this removes the OFT-side one and unblocks
removing the rest downstream.
Validation
jphi-linterp Ip hold across operating points (DIII-D-like H-mode, same mesh /
coils / isoflux; realized Ip via
get_stats):Holds Ip to ~0.0004 % of target in 12–15 nonlinear iterations across the range.
(For reference, our earlier outer-loop held ~0.01–0.05 % in 14–39 iterations on
the same cases; OpenFUSIONToolkit#267 is ~100× tighter and faster.)
find_optimal_scalecore-j_0behavior unchanged; the Python package importsand the merged tree builds cleanly.
Notes
mainand a merge of TokaMaker: Improve Ip matching with Jphi profiles OpenFUSIONToolkit/OpenFUSIONToolkit#267; canrebase to a clean linear PR (or reduce to just
bootstrap.pyonce TokaMaker: Improve Ip matching with Jphi profiles OpenFUSIONToolkit/OpenFUSIONToolkit#267 isupstream) on request.