EQUIL + ForceFreeStates - Extend equilibrium splines to separatrix for diverted plasmas by logan-nc · Pull Request #179 · OpenFUSIONToolkit/GPEC

logan-nc · 2026-03-01T03:03:40Z

Summary

Equilibrium: Edge inverse splines for diverted plasmas

Adds InverseCubicSpline{S} type wrapping an iota = 1/q inner spline; callable returns 1/inner(ψ), enabling well-behaved evaluation of q and dV/dψ as psin → 1 in diverted plasmas (where q → ∞)
Makes ProfileSplines mutable with a q_spline::Union{S,I} pointer field temporarily swapped to the edge inverse spline for ODE integration chunks that extend beyond psihigh
Adds direct_xpoint_find using Newton iteration for Br=0, Bz=0 with physically motivated initial guesses on the inboard side; correctly classifies diverted vs limited topology (critical fix: initial R guess at rs1 + 0.35*(ro-rs1) rather than midpoint)
For diverted plasmas, builds edge iota inverse splines over [0.9·psihigh, 1.0] with anchor iota(1.0) = 0 (CubicFit BC) enforcing the known separatrix physics
Sets qa = Inf for diverted plasmas; q_for_mode_range fallback prevents trunc(Int, Inf) in mode range computation
dV/dψ evaluation in Free.jl routes to dVdpsi_spline_inv when psilim > psihigh
New EquilibriumParameters fields: r_xpoint, z_xpoint, is_diverted
ProfileSplines fields: q_spline_direct, q_spline_iota_inverse, dVdpsi_spline_inv (all Nothing for limited plasmas → identical behavior to before)

ForceFreeStates: Integration beyond psihigh for diverted plasmas

Extends rational surface search in sing_find! above psihigh using the iota spline, stopping when surface spacing < edge_layer_width (new TOML parameter, default 1e-4)
Removed set_psilim_via_dmlim and dmlim — integration limit now determined entirely by psiedge < psilim logic; psilim for diverted plasmas is set from the last edge rational surface found
Moved sing_find! before sing_lim! in JPEC.jl so edge surfaces above psihigh are available when sing_lim! reads them
sing_lim! sets psilim = last_edge_surface + edge_layer_width/2 for diverted, giving a valid final ODE chunk beyond the last crossed surface
Fixed intr.mtheta → ctrl.mthvac bug in free_compute_wv_spline
Resizes dW_edge to match current step count in findmax_dW_edge! (integration can exceed initial numsteps_init allocation)
Handles SingularException in the edge dW scan — degenerate U₁ near the separatrix layer is skipped; those steps remain at -Inf and are ignored in the max search
transform_u! caps kfix = min(fixstep[ifix], odet.step) and breaks when jfix > odet.step to handle trimmed storage correctly after psiedge truncation
Fixed BoundsError when extracting psi_edge_scan: psi_store may be doubled by grow_storage!, so only the valid [1:odet.step] range is indexed before applying the edge mask

Diagnostic output: et[1] vs psi in the edge zone

Adds psi_edge_scan and et_edge_scan fields to OdeState; populated in findmax_dW_edge! from the full un-trimmed psi range before trim_storage! cuts to the peak-dW step
Saved to HDF5 as integration/psi_edge_scan and integration/et_edge_scan for post-run diagnostics

Benchmark script: `benchmarks/benchmark_edge_splines.jl`

Produces 7 diagnostic plots saved alongside jpec.toml:

iota (= 1/q) vs psin — inner spline descends smoothly to 0 at psin=1
q vs psin — direct spline (ExtendExtrap, diverges) vs edge inverse spline (well-behaved)
dV/dψ vs psin — same comparison showing asymptotic growth handled correctly
R(psin) and Z(psin) for fixed poloidal angles toward the edge
Flux surface cross-section with x-point marked at (R_x, Z_x), last few surfaces
Rational surface density in edge zone — vertical lines for surfaces found above psihigh
et[1] vs psi in the edge — reads integration/psi_edge_scan from jpec.h5, marks the peak dW truncation point

Backward compatibility

Limited plasmas: all edge fields are nothing, q_spline = q_spline_direct at startup, behavior unchanged
q_spline_iota_inverse and dVdpsi_spline_inv only activate when diverted x-point is detected
All test cases (Solovev, CHEASE, EFIT regression) remain unaffected

End-to-end verification (DIIID-like example)

X-point (lower null): R = 1.315 m, Z = -1.114 m
Edge inverse splines: psin ∈ [0.8937, 1.0] (28 nodes)
Extended edge search: 9 rational surfaces above psihigh = 0.993
ODE integration:     ψ = 0.993 → 0.999 (q = 5.3 → 14.5, crossing 9 edge surfaces)
Edge dW truncation:  peak at ψ = 0.995, q = 5.777

Test plan

runtests_equil.jl passes (includes new testsets for edge splines on DIIID-like EFIT and limited plasma regression)
runtests_solovev.jl passes
DIIID-like EFIT: pe.params.is_diverted == true, X-point found at R=1.315, Z=-1.114
Edge iota spline continuity: |q_direct(psihigh) - q_edge(psihigh)| < 1e-14, same for dV/dψ
iota(1.0) ≈ 0.0 (anchor condition)
Solovev (limited): q_spline_iota_inverse === nothing, behavior unchanged
Run benchmarks/benchmark_edge_splines.jl after a DIIID-like JPEC run to verify all 7 plots

🤖 Generated with Claude Code

… separatrix for diverted plasmas Adds infrastructure to integrate the ideal MHD ODE beyond psihigh toward the near-separatrix edge layer in diverted plasmas, where q→∞ and dV/dψ→∞. Key changes: - New InverseCubicSpline{S} type wrapping an iota=1/q inner spline with a self-contained hint; returns 1/inner(ψ) for well-behaved evaluation near psin=1 - ProfileSplines is now a mutable struct with a Union{S,I} q_spline pointer field that can be temporarily swapped to the edge inverse spline for ODE chunks beyond psihigh, with q_spline_direct always available for deriv2/deriv3 - X-point detection via Newton iteration (Br=0, Bz=0) with physically motivated initial guesses on the inboard side (rs1 + 0.35*(ro-rs1)); correctly identifies diverted vs limited plasma topology - For diverted plasmas: edge iota inverse splines for q and dV/dψ built over [0.9*psihigh, 1.0] with anchor iota(1.0)=0 (CubicFit BC) - qa=Inf for diverted; q_for_mode_range fallback in JPEC.jl avoids trunc(Int,Inf) - Extended rational surface search in sing_find! above psihigh using the iota spline, terminating when surface spacing < edge_layer_width (new TOML param) - sing_lim! updated for diverted: psilim from last edge surface, qlim=ctrl.qhigh - Pointer swap in integrate_el_region! before/after edge ODE chunks - eval_dVdpsi helper in check_for_zero_crossings! routes to dVdpsi_spline_inv when psi > psihigh for correct post-integration stability evaluation - Benchmark script: benchmarks/benchmark_edge_splines.jl with iota, q, dV/dψ comparison plots and rational surface density visualization Backward compatible: limited plasmas have edge fields=nothing and are unchanged. All existing tests pass. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…e spline feature - Fix dVdpsi_spline routing in Free.jl: dVdpsi_spline(psilim) in free_run! and free_compute_total now routes to dVdpsi_spline_inv when psilim exceeds psihigh for diverted plasmas, avoiding inaccurate ExtendExtrap at the edge - Add benchmark plots 4 (R,Z vs psin for several poloidal angles) and 5 (flux surface cross-section with x-point marked); rename rational surface density to plot 6; all plots use rzphi spline reconstruction R=ro+rfac*cos(η), Z=zo+rfac*sin(η) - Add edge_layer_width = 1e-4 to [ForceFreeStates] in DIIID-like example jpec.toml - Add formal test assertions for edge spline properties: limited plasma returns nothing, diverted plasma has correct anchor (iota(1.0)≈0), correct coverage, correct topology flags, qa=Inf, and continuity at psihigh Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…tend ODE integration past psihigh for diverted plasmas - Remove set_psilim_via_dmlim and dmlim control parameters entirely; integration limit is now determined solely by the psiedge < psilim logic in EulerLagrange.jl - Move sing_find! before sing_lim! in JPEC.jl so that edge rational surfaces above psihigh are available when sing_lim! sets psilim for diverted plasmas - sing_lim! now sets psilim = last_edge_surface + edge_layer_width/2 for diverted plasmas, giving a valid final ODE chunk beyond the last crossed surface - Fix intr.mtheta -> ctrl.mthvac bug in free_compute_wv_spline (Free.jl) - Resize dW_edge in findmax_dW_edge! to match current step count (integration may exceed initial numsteps_init allocation) - Handle SingularException in findmax_dW_edge! edge scan (degenerate U1 near the separatrix layer); those steps remain at -Inf and are ignored in max search - Fix transform_u! to cap kfix at odet.step and break early when jfix exceeds trimmed storage (prevents BoundsError after psiedge truncation) - Update DIIID-like example: remove dmlim, set psiedge = 0.995 to trigger edge dW scan above psihigh = 0.993 - Verified: DIIID-like end-to-end run integrates to psilim = 0.999 (q = 14.5), truncates at ψ = 0.995 (edge dW peak), completes FFS + PerturbedEquilibrium Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…d stability plot - Add psi_edge_scan / et_edge_scan fields to OdeState; populated in findmax_dW_edge! from the full psi_store before trim_storage! cuts it to the peak-dW step, so the complete stability profile from psiedge → psilim is retained - Save integration/psi_edge_scan and integration/et_edge_scan to jpec.h5 when non-empty - benchmark_edge_splines.jl: add plot 7 (et[1] vs psi in edge zone, reads from jpec.h5), add HDF5 import, update docstring to list all 7 plots Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ted psi_store psi_store may be doubled by grow_storage! so it is larger than odet.step. Index only the valid [1:odet.step] range before applying the edge_mask. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…e dW scan bugs; add F_iota_matched spline; GSE renames and plots Part A — Fix edge dW scan bugs in Free.jl: - Bug 1 (free_compute_total): dV/dψ was evaluated at intr.psilim instead of odet.psifac. Near the separatrix (psilim≈0.999), dV/dψ is ~100× larger than at the actual scan point, making the eigenvalue normalization wrong by ~100×. - Bug 2 (free_compute_wv_spline): VacuumInput used intr.psilim as the plasma boundary for all scan points. Now passes psi_array[i] so vacuum geometry corresponds to the surface actually being scanned. Part B — Benchmark Plot 7 now plots et[1] vs q (not psi): - Converts psi_es → q_es via iota inverse spline above psihigh - Adds hline at et=0 (stability boundary), vline at q(psihigh), and reference line at et=+1.707 (develop branch value at q≈5.2) Part C1 — GSE variable renames in equilibrium_gse!: - flux_fs→gs_div_term, flux_fsx→gs_div_dpsi, flux_fsy→gs_div_dtheta - s1→F_val, s1p→dF_dpsi, s2p→dP_dpsi - term→gse_terms_integrated, totali→gse_total_integrated - errori→gse_abs_error, errlogi→gse_log_error - Fixed pre-existing bug: gse.h5/gsei.h5 writes used nonexistent flux.xs/ys - Updated docstring to describe the GS equation and output groups Part C2 — F_iota_matched spline for self-consistent far-edge F(ψ): - ProfileSplines gains F_spline_direct, F_deriv_direct (stored originals), F_spline_iota_matched, F_deriv_iota_matched (Nothing for limited plasmas) - F_spline (and F_deriv) are now switchable pointers like q_spline - direct_build_F_iota_matched! computes J(ψ)=(1/2π)∫ J_coord/R² dθ from raw rzphi_nodes, scans n=1 rational surfaces for 10-pts-per-window grid density, then builds F_matched = q_target·2π/J; called inside if is_diverted block - integrate_el_region! switches both q_spline and F_spline to edge versions for ODE chunks that extend beyond psihigh, restores direct after Part C3 — Benchmark plots 8-10: - Plot 8: F(ψ) direct vs F_iota_matched with dF/dψ panel - Plots 9/10: GSE residual diagnostics via diagnose_src=true + HDF5 read Fixes stale dmlim/set_psilim_via_dmlim in two test TOML files. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ace sampling, vac_psi psihigh cap, always-save in edge zone - free_compute_wv_spline: switch from q-space to psi-space sampling (20 core + 20 edge points); avoids inverting the direct q spline beyond psihigh where it extrapolates linearly and returns psi >> 1.0, building the wvmat spline over the wrong domain - free_compute_wv_spline: cap vac_psi = min(psi_array[i], psihigh) when calling VacuumInput; the rzphi bicubic spline only covers [psilow, psihigh] and extrapolating beyond gives negative r^2 values → sqrt crash - integrator_callback!: always save every ODE step when integrator.t >= ctrl.psiedge (edge scan zone); with save_interval=10 only ~3 steps fall in [psiedge, psihigh] = [0.990, 0.993] by default, causing findmax_dW_edge! to miss the stability peak in that region Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…tion at psihigh, skip u_store above psihigh, fix F_spline chunk condition - eulerlagrange_integration: break out of chunk loop when edge scan is active and chunk.psi_start > psihigh; findmax_dW_edge! only scans [psiedge, psihigh] so chunks above psihigh contribute nothing but take exponentially more computation near the separatrix (q=9 alone requires 86k+ ODE steps → OOM) - integrator_callback!: return early (skip u_store save) when edge_scan_active and integrator.t > psihigh; the near-separatrix steps are never read by findmax_dW_edge! and storing them at save_interval=10 causes OOM - integrate_el_region!: change using_edge_spline condition from psi_END > psihigh to psi_START >= psihigh; the q=5→q=6 chunk straddles psihigh (0.986→0.996), so psi_end > psihigh was triggering F_spline_iota_matched for psi < 0.993 where it extrapolates a wrong constant F value → corrupt U₁ → SingularException in free_compute_total for all [psiedge, psihigh] steps Result: full DIIID-like example completes in 18.6 s (vs OOM/hours before); truncation at ψ=0.992, q=5.277; et[1] = +1.763 (stable, consistent with develop branch reference of +1.707 at q=5.2) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

logan-nc · 2026-03-02T13:22:33Z

Session update (2026-03-02): edge dW scan bugs fixed, full DIIID-like run verified

Five additional bugs found and fixed during verification of the edge dW scan:

Bug 3 (free_compute_wv_spline): Used q-space sampling with the direct q spline, which extrapolates linearly beyond psihigh and returns ψ>>1. Switched to ψ-space sampling (20 core + 20 edge points) with q_spline_iota_inverse for ψ > psihigh.

Bug 4 (free_compute_wv_spline): VacuumInput called with ψ values beyond psihigh where the rzphi bicubic spline fails (negative r² → sqrt crash). Added vac_psi = min(psi_array[i], psihigh) cap.

Bug 5 (integrate_el_region!): using_edge_spline checked psi_END > psihigh but should be psi_START >= psihigh. The q=5→q=6 chunk (ψ 0.986→0.996) straddles psihigh; using F_iota_matched (built only for ψ>0.993) for ψ<0.993 extrapolates a wrong constant F value → corrupt U₁ → SingularException in free_compute_total for all [psiedge, psihigh] steps.

Bug 6 (integrator_callback!): Steps above psihigh were saved to u_store (never needed). Added early return when edge_scan_active && integrator.t > psihigh.

Bug 7 (eulerlagrange_integration): ODE integrated all the way to psilim=1.0 even though findmax_dW_edge! only scans [psiedge, psihigh]. Near-separatrix ODE is extremely stiff (q=9 alone: 86k steps → OOM). Added break when edge_scan_active && chunk.psi_start > psihigh.

Verified result: Full DIIID-like example (examples/DIIID-like_ideal_example) now completes in 18.6 s with et[1] = +1.763 at ψ=0.992, q=5.277 — stable and consistent with the develop branch reference (+1.707 at q=5.2). Previously OOMed or ran for hours.

…and extend scan beyond psihigh Two bugs in free_compute_wv_spline / free_compute_total caused the edge stability scan to give et ≈ −14600 instead of the correct +1.763 for the DIIID-like diverted plasma. Bug 1 (dV/dψ at psilim instead of psifac, Free.jl ~line 211): free_compute_total was evaluating dV/dψ at intr.psilim (≈ 0.999 near the separatrix) instead of at odet.psifac (the current scan psi ≈ 0.990–0.993). With psilim ≈ 0.999, dV/dψ is ~100× larger than at psiedge, making the normalisation denominator huge and the eigenvalue tiny/wrong. Fix: use odet.psifac for the dV/dψ lookup in free_compute_total. Bug 2 (VacuumInput at psilim instead of psi_array[i], Free.jl ~line 172): free_compute_wv_spline was building the wvmat CubicSeriesInterpolant using intr.psilim ≈ 0.999 as the plasma boundary for all grid points in the spline. This pinned all vacuum responses to the near-separatrix geometry and produced norm_wv ≈ 26840 (wrong) instead of ≈ 5730 (correct). After the fix, the scan gives et ≈ +1.627–1.763 across [psiedge, psihigh], matching the develop-branch reference et = +1.707 (within 3.3%). Fix: use psi_array[i] (the actual scan point) as the plasma boundary. Additional improvements in this commit: - free_compute_wv_spline now builds a CORE-ONLY grid [psiedge, psihigh] (20 pts) with RAW vacuum response (no singfac). ExtendExtrap handles above-psihigh queries by returning the psihigh value. This avoids the derivative discontinuity at psihigh that previously shifted the scan peak. - singfac = (m - n·q(psifac)) is now applied analytically in free_compute_total, using the actual q at the scan psi. This keeps the spline well-conditioned. - Extend the edge scan beyond psihigh into [psihigh, psilim]: - EulerLagrange.jl: above-psihigh ODE chunks now run (up to 20k steps each) with sparse step saving (every 100×save_interval steps) to avoid OOM. - A per-chunk step budget cap (terminate!(integrator)) prevents exponential compute growth as q → ∞ near the separatrix. - integrator_callback! uses chunk_callback_count (not length(sol.t)) for near_start detection, since save_positions=(false,false) is now used. - OdeState gains chunk_callback_count, chunk_steps_start, chunk_is_above_psihigh fields for the above-psihigh tracking logic. - findmax_dW_edge! scans ALL stored steps (not just those ≤ psihigh). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…gnostic (Plot 11) Plot 7 (edge_spline_stability.png): The previous version plotted the full scan [psiedge, psilim], including above-psihigh steps where ExtendExtrap pins wv to the psihigh boundary geometry, producing et ≈ -1e8 that crushed the interesting core region (et ≈ +1.76 at q≈5.28) off-screen. Fixed: split into core mask (psi ≤ psihigh, wv interpolated) and above-psihigh (wv extrapolated, clamped to [-10, 10] for context). Core data now dominates the y-axis and the peak at et=+1.763, q=5.2774 is clearly visible alongside the develop-branch reference at et=+1.707. Plot 11 (edge_spline_gse_edge.png) — new: Three-panel edge GSE diagnostic showing the GS equation error from using the EFIT F(ψ) (extrapolated via ExtendExtrap as constant beyond psihigh) vs the iota-matched F: A. q_direct (EFIT) vs q_target = 1/ι(ψ): constant q vs diverging physical q → the EFIT F has no information about the near-separatrix q profile. B. Relative q deviation (q_direct − q_target)/q_target [%]: reaches -96% at ψ=0.9999, showing the EFIT F implies q~5.33 everywhere when the true q → ∞. C. GS source mismatch ΔFF′ = F_direct·F′_direct − F_matched·F′_matched: grows large with oscillatory structure near ψ→1 (from rational surface resolution in F_iota_matched). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ics to plots 9 and 10 For plots 9 (GS residual per theta) and 10 (theta-integrated GS error), overlay the far-edge (psi > psihigh) GS residual computed using F_iota_matched. Above psihigh, the rzphi splines use ExtendExtrap (geometry frozen at psihigh), so the GS residual reduces to the source term alone: source = J/(2pi*psi0*pi^2) * (F*F'/(2pi*R)^2 + P'). Key features: - Evaluation grid matches the F_iota_matched spline node construction (12 pts per rational surface window), limited to [psihigh, edge_surfaces[end]] where F_iota_matched is physically relevant (avoids the divergent q->inf region near psi=1) - Both NODES (spline knots) and inter-node MIDPOINTS are evaluated and plotted distinctly in Plot 10, providing a check that the cubic spline interpolates smoothly - Far-edge lines use same colors as core per-theta lines in Plot 9 (dashed) for direct visual comparison - Core per-theta lines now use explicit colors to match far-edge overlay colors - Added `using FastInterpolations` to imports Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…-edge rzphi extension Replace F_iota_matched approach with physically correct X-point asymptotic geometry. F'≈0 and P'≈0 near the separatrix means GS source ≈ 0 — no iota-matched F needed. Flux surfaces in the far edge (psin > psihigh) now scale from the X-point: R(ψ_n,θ) = R_X + s·(R_high−R_X), Z similarly, where s = √((1−ψ_n)/(1−psihigh)) Key changes: - DirectEquilibrium.jl: replace direct_build_F_iota_matched! with direct_build_xpoint_asymptotic!, which extends rzphi_nodes to psin=0.999999 using X-point scaling (10 pts per rational surface window) - direct_xpoint_find returns 6-tuple including ψ Hessian components at X-point - EquilibriumTypes.jl: add psi_hess_rr/zz/rz and psi_core_end to EquilibriumParameters; remove F_spline_iota_matched/F_deriv_iota_matched from ProfileSplines - Equilibrium.jl: equilibrium_gse!, equilibrium_separatrix_find!, equilibrium_global_parameters! use psi_core_end to index into core grid (pe.rzphi_xs is now the extended grid); fix gsei.h5 and gse.h5 to save psi_xs - EulerLagrange.jl: remove F_spline pointer swap (only q_spline swapped for edge chunks) - benchmark_edge_splines.jl: replace F comparison plot (Plot 8) with X-point geometry visualization; update far-edge GSE section; fix psi_gse_vals to use profiles.xs Result: core GS residual max = 3.71e-3 (below 1e-2); all 20 equilibrium tests pass. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…_xs instead of core grid equil.rzphi_xs is now the extended grid (250 nodes for diverted plasmas) after the X-point asymptotic geometry extension. make_metric was using length(equil.rzphi_xs) as mpsi, causing make_matrix to index profiles.xs at out-of-bounds positions (130+). Fix: use length(profiles.xs) for mpsi in make_metric, which gives the core psi grid size (129 nodes). The FGK matrices only need to cover the core equilibrium region; the far-edge extension is only needed for the ODE integration via the rzphi splines. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

… resonant spikes Previously the stability plot clipped at xlims = (4, q_psihigh * 1.02) ≈ (4, 5.44), hiding all above-psihigh data even though it was being plotted. The scan actually spans q = 5.16 → 7.99 (ψ = 0.990 → 0.998). Changes to Plot 7: - xlims expanded to full scan range: (4.0, q_max_scan * 1.02) - Above-psihigh data filtered to |et| ≤ 10 (removes resonant spikes where singfac → 0 near rational surfaces) and subsampled to ≤ 2000 points for rendering performance - ylims fixed to (-10, max(3, peak*1.2)) so spikes don't dominate the y scale - Added scan range and peak values to stdout summary - Title updated to show actual ψ scan bounds Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…g; add SingularException guard Key fixes: 1. EulerLagrange.jl: Move transform_u! to AFTER findmax_dW_edge!. The scan must use raw (pre-transform) u_store. U₁_raw is kept non-singular by the ODE fixups; the transformed U₁_core = U₁_raw × T₁ (where T₁ = G₁·…·Gₙ is the product of ALL Gaussian matrices) is ill-conditioned with many fixups, causing SingularException in free_compute_total for every core scan step. With transform before scan, only above-psihigh steps (where T_{n+1} = I) return values, forcing the peak to q≈8 instead of the correct q≈5.2. 2. FixedBoundaryStability.jl: Catch SingularException in check_for_zero_crossings!. U₂ can be degenerate near rational surfaces in the above-psihigh integration region. Propagate the previous crit_store value to avoid a spurious zero crossing. 3. Free.jl: Simplify free_compute_wv_spline; fix compute_vacuum_response to 4-value destructuring (wv_block, _, _, _); update docstring to clarify core-only grid design and rationale for ExtendExtrap above psihigh. Verified: et[1] = +1.763 at ψ=0.991, q=5.195 (develop reference: +1.707 at q=5.2). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…rease wv spline nodes Investigating why the edge dW scan gives only 6 valid points (|et|≤10) out of 11664 before jumping to et≈−343 at psi=0.99117. Changes in this commit: 1. free_compute_total: Replace vector @info with scalar diagnostic. Previous @info printed 68-element eigval vectors that were truncated as "68-element Vector{Float64}: …" in the log output, yielding no useful information. New diagnostic prints min eigenvalues of wt, wp, and wv separately, plus cond(U₁) raw and post-column-normalization. This will distinguish whether the jump at psi=0.99117 is caused by wp (U₁ conditioning) or wv (spline oscillation). 2. free_compute_wv_spline: Increase npsi_core from 20 to 50. With 20 nodes in [psiedge=0.990, psihigh=0.993] (spacing 0.000158), the scan step spacing (~0.000160) roughly equals the spline node spacing — any anomalous node causes oscillation across adjacent scan steps. 50 nodes (spacing 0.000061) gives ~3 scan steps per spline interval, providing better regularization. 3. ForceFreeStatesStructs + Fourfit: Add jmat_spline infrastructure. Builds a spline of Jacobian Fourier coefficients J(θ, ψ) vs ψ for potential use in per-step jmat normalization. Currently unused (jmat normalization deferred to free_run! where psilim is fixed) but retained for future investigation. Context: et[1]=+1.763 at the final truncation point is correct (matches develop branch +1.707). The scan values for steps 1-6 (+1.669→+1.800→+1.744) are smooth and physical. The jump to -343 at step 7 is the unresolved bug. cond(U₁)=1.49e90 at both the last valid step and the first invalid step — ruling out simple conditioning as the cause. Next step: run with new diagnostics to see whether evwv_min or evwp_min jumps at the transition. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

logan-nc · 2026-03-09T13:09:14Z

Edge dW scan investigation — current status ($(date '+%Y-%m-%d'))

What's working

Final et[1] = +1.763 at ψ=0.991, q=5.195 — correct, close to develop branch +1.707
Edge scan correctly truncates at the peak dW step (step 6 of the scan)
The 6 valid scan points (ψ ∈ [0.990, 0.991]) give smooth, physically reasonable values: +1.669 → +1.733 → +1.781 → +1.800 → +1.797 → +1.744

The bug

After psi=0.99101, et abruptly jumps to −343.85 and continues worsening to −6390 then slowly recovering. Only 6 of 11664 scan points are valid (|et| ≤ 10). This doesn't corrupt the final et[1] (since all invalid values are negative and are excluded by the max search), but the scan data is wrong for the benchmark plots.

Root cause investigation

What we know:

The transition psi=0.99117 corresponds to entering spline interval [node 7, node 8] of the 20-node (now 50-node) wv spline. Previously, spline node 7 was at psi=0.991105, one step before the transition.
cond(U₁) = 1.49e90 at the last valid step and 2.3e90 at the first invalid step — both astronomically large, yet step 6 gives +1.744. This rules out raw U₁ conditioning as the sole cause.
The column normalization in free_compute_total (normalizes U₁ and U₂ columns by the same U₁ column norms, preserving W = U₂·U₁⁻¹) may make the effective condition number much lower than the raw cond.
Pattern (−343 → −2803 → −5120 → −6390 peak → recovering) is smooth, suggesting a systematic error rather than pure noise — consistent with cubic spline oscillation around an anomalous node.

Two hypotheses:

wv spline: Node 7 (psi=0.991105) has a slightly wrong vacuum response value, causing cubic spline oscillation in the [7,8] interval. Fix: increase npsi_core (done: 20→50) or use Fortran-style singfac pre-scaling.
wp computation: The column-normalized U₁ still has poor conditioning in a way that causes wrong W=U₂·U₁⁻¹ at step 7 but not step 6. This could happen if a Gaussian reduction triggered between steps 6 and 7 creates a particular triangular structure that breaks the normalization.

What was committed today (`42beb70`)

Improved diagnostic in free_compute_total: now prints scalar min eigenvalues of wt, wp, and wv separately, plus raw and post-normalization cond(U₁). The previous diagnostic printed 68-element vectors truncated as "…" — useless.
wv spline density: 20 → 50 nodes in free_compute_wv_spline. Reduces the chance of a single bad node causing oscillation.
jmat_spline infrastructure: Builds a spline of J(θ, ψ) vs ψ in make_matrix. Currently unused but available for future per-step normalization experiments.

What to try next

Run JPEC with the new scalar diagnostics to get:

FCT_DIAG psi=0.991013: evwt_min=X  evwp_min=X  evwv_min=X  cU1_raw=X  cU1_norm=X  u1nrm=[min,max]
FCT_DIAG psi=0.991170: evwt_min=X  evwp_min=X  evwv_min=X  cU1_raw=X  cU1_norm=X  u1nrm=[min,max]

If evwv_min jumps (large change at psi=0.99117): wv spline is the cause → try Fortran-style singfac pre-application in the spline.

If evwp_min jumps (large change at psi=0.99117) while evwv_min stays normal: wp/U₁ is the cause → try SVD/pseudoinverse for W = U₂·pinv(U₁) instead of direct rdiv.

If cU1_norm is much smaller than cU1_raw: column normalization is helping, wp computation should be OK. Focus on wv.

If 50 nodes already fixes the scan (et values smooth throughout): it was a spline oscillation issue, no further work needed on the wv side.

…t normalization Three connected fixes for the edge dW scan and metric computation in diverted plasmas: 1. Fix findmax_dW_edge! scan ranking (Free.jl, EulerLagrange.jl): - free_compute_wv_spline now uses psi_array[i] as plasma boundary (not psihigh). Without this, wv(psihigh boundary) overestimates et for core steps relative to edge steps, producing a spurious peak at q≈3 instead of the physical peak near q≈5. - free_compute_total now normalizes by ffit.jmat * |ξ|² / dV_dψ(psifac), matching free_run! exactly. Since dV_dpsi varies ~44% across the scan range, unnormalized eigenvalues gave unfair advantage to core steps (smaller dV_dpsi → larger raw et). - Scan limited to psihigh + 0.5*(1-psihigh) for diverted plasmas: above this cutoff, singfac²=(m-nq)² grows unboundedly while wv is frozen at psihigh via ExtendExtrap. - Tracks ep_edge_scan, ev_edge_scan, evonly_edge_scan for diagnostic output to HDF5. - Verified: scan proxy (Total=+1.796 at ψ=0.943, q=4.198) matches actual free_run! output. 2. Limit Fourfit/make_matrix to psilim (Fourfit.jl, JPEC.jl, Sing.jl): - make_metric now accepts psilim kwarg; truncates rzphi_xs grid at psilim. Far-edge X-point geometry gives diverging FGK matrices; capping at psilim avoids ill-conditioned splines that corrupt the core edge integration. - make_matrix captures jmat at core_end_idx (last psi ≤ psihigh) instead of at the far-edge last node. The X-point metric is near-singular near the separatrix. - q and q' in make_matrix loop use iota inverse chain rule for psi > psihigh. - sing_lim! sets psilim just before the last edge surface using singfac_min buffer (replacing the previous edge_layer_width/2 offset that could trigger EL assertions). 3. Extend dense rzphi window past psilim (DirectEquilibrium.jl): - direct_build_xpoint_asymptotic! extends the dense per-window grid one full rational window beyond psilim so the spline covers the complete window psilim sits in. Also: update DIIID-like example psihigh to 0.992 and add ep/ev/evonly diagnostic plots to benchmark_edge_splines.jl (plots 12–14 showing plasma vs vacuum energy components). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

psilim for diverted plasmas is ~1.0 (near the separatrix), so the previous change to use rzphi_xs truncated at psilim inadvertently included all 77 far-edge nodes in the FGK matrix computation. Far-edge geometry near the X-point has J→∞, making the ODE extremely stiff above psihigh and causing very long integration times. Fix: compute min(psilim, psihigh) inside make_metric to always cap at psihigh. For limited plasmas (psilim < psihigh) this still correctly caps below psihigh. The jmat_at_core_end snapshot and q-via-iota-chain-rule branches are retained (correctly become no-ops for the psihigh-capped grid). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

logan-nc · 2026-03-11T12:13:44Z

Session summary (2026-03-11)

What was completed this session

Three bugs fixed in the edge dW scan and metric computation, now pushed as commits 0d6fa349 and ffbc3d4b:

1. `findmax_dW_edge!` scan ranking fixed (`Free.jl`, `EulerLagrange.jl`)

Two root causes for the false scan peak at q≈3 instead of the physical peak near q≈5:

Wrong plasma boundary for wv: free_compute_wv_spline was computing wv with psihigh as the plasma boundary for every scan step. Fixed to use psi_array[i] (the actual psifac at each scan point), making wv physically consistent with the plasma bounded by psifac.
Missing normalization: free_compute_total was returning a raw unnormalized eigenvalue. Since dV/dψ varies ~44% across the scan range, core steps (smaller dV/dψ) appeared more unstable than edge steps. Fixed by normalizing by ffit.jmat * |ξ|² / dV_dψ(psifac), matching free_run! exactly. Since jmat[mid](psihigh) ≈ dV/dψ(psihigh), this normalizes each step by dV/dψ(psihigh) / dV/dψ(psifac).

Scan now correctly returns Total = +1.796 at ψ = 0.9425 (q ≈ 4.198), confirmed to match actual free_run! output exactly. Also added tracking of ep_edge_scan, ev_edge_scan, evonly_edge_scan to HDF5.

2. `make_metric` grid capped at `psihigh` (`Fourfit.jl`)

The previous commit inadvertently used rzphi_xs truncated at psilim ≈ 1.0 for the FGK matrix grid, including all 77 far-edge X-point nodes (J→∞ near separatrix). This made the ODE extremely stiff above psihigh, causing very long runs or SIGTERM from macOS memory pressure. Fixed by capping at min(psilim, psihigh) inside make_metric.

3. `jmat` snapshot at `psihigh` (`Fourfit.jl`)

make_matrix now snapshots jmat at core_end_idx (last psi ≤ psihigh) rather than at the last loop iteration. This ensures ffit.jmat (used in free_run! and free_compute_total normalization) is evaluated at the physical plasma-vacuum interface, not at a far-edge node.

Verified end-to-end run

Full DIIID-like example runs cleanly with ffbc3d4b:

ψ = 0.551, q = 2.000, steps = 703
ψ = 0.811, q = 3.000, steps = 867
ψ = 0.926, q = 4.000, steps = 986
ψ = 0.986, q = 5.000, steps = 1105
Peak edge dW: ψ = 0.9425, q ≈ 4.198, et[1] = +1.796 (ep=0.153, ev=1.643)
19 edge surfaces found above psihigh = 0.992

Remaining TODOs before merge

Benchmark plots: Run benchmarks/benchmark_edge_splines.jl to generate the 11+ diagnostic plots and verify all look physically reasonable (especially plots 5/6: et vs q scan, plots 9/10: GSE residuals).
Non-Hermitian warning: After trimming to peak-dW step, evaluate_stability_criterion! sees 242 non-Hermitian steps (W matrix built from ill-conditioned transform products). This affects the Newcomb zero-crossing count but not et[1] (which uses free_run! with odet.u loaded before transform_u!). Low priority.
psiedge tuning: With psihigh=0.992, the scan covers ψ ∈ [0.800, 0.9425]. If only the near-edge region (q≈5) is of interest, raise psiedge from 0.800 to ~0.96 in jpec.toml.
psilim = 1.000 display: The run log shows psilim = 1.000 because 19 edge surfaces reach near ψ=1. This is correct behavior — the above-psihigh ODE runs through those surfaces with the 20k-step-per-chunk cap. No action needed unless runtime is a concern.
PR review: The 4 main feature commits on this branch are:
- f439a330: X-point asymptotic geometry for far-edge rzphi extension
- 33af3db2: Fix make_metric mpsi using extended rzphi_xs
- af132b29: Fix transform_u!/findmax_dW_edge! ordering; SingularException guard
- 42beb70d: Improve edge dW scan diagnostics; increase wv spline nodes
- 0d6fa349: Fix findmax_dW_edge! scan ranking via jmat normalization
- ffbc3d4b: Cap make_metric grid at psihigh

…_spline_extension

…hmark scripts can start from the actual examples and make any modifications dynamically - the replicated tomls and gfiles were already getting out of sync!

- benchmarks/benchmark_edge_splines.jl: update module name (JPEC→GeneralizedPerturbedEquilibrium), default toml path, and h5 filename (jpec.h5→gpec.h5) - test/runtests_equil.jl: fix JPEC→GeneralizedPerturbedEquilibrium and wrong geqdsk filename (TKMKR_D3Dlike_default_Hmode→TkMkr_D3Dlike_Hmode) - Project.toml: remove PlotlyJS dep (only used in analyze_example.jl scripts, incompatible with JSON 1.x — users should install it separately) - examples/DIIID-like_ideal_example/gpec.toml: raise psiedge 0.800→0.970 (consistent with verified PR run) - src/Equilibrium/DirectEquilibrium.jl: update JPEC→GPEC in comment Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

logan-nc · 2026-03-11T21:52:54Z

Session summary (2026-03-11, afternoon)

What was completed this session

GPEC rename cleanup (commit 28a4cdd9)

After merging develop (which renamed the package from JPEC → GeneralizedPerturbedEquilibrium and output from jpec.h5 → gpec.h5), several files on this branch still had stale references:

benchmarks/benchmark_edge_splines.jl: updated using JPEC → using GeneralizedPerturbedEquilibrium, default toml path (jpec.toml→gpec.toml), h5 filename references (jpec.h5→gpec.h5), and all JPEC.Equilibrium.* calls
test/runtests_equil.jl: fixed JPEC.Equilibrium → GeneralizedPerturbedEquilibrium.Equilibrium in the edge spline testset; also fixed wrong geqdsk filename (TKMKR_D3Dlike_default_Hmode.geqdsk → TkMkr_D3Dlike_Hmode.geqdsk)
src/Equilibrium/DirectEquilibrium.jl: updated JPEC→GPEC in comment
examples/DIIID-like_ideal_example/gpec.toml: raised psiedge 0.800 → 0.970 (consistent with the verified March 11 run)
Project.toml: removed PlotlyJS from deps — it's only used in analyze_example.jl scripts, incompatible with JSON 1.x (all PlotlyJS ≤ 0.18.17 cap JSON at 0.21.x), and was preventing using GeneralizedPerturbedEquilibrium from loading. Users who want analyze_example.jl should install PlotlyJS in their own environment.

Benchmark example dirs removed (commit 5bd1c1c8)

benchmarks/DIIID_ideal_example/ and benchmarks/Solovev_ideal_example/ deleted — orphaned dirs not referenced by any benchmark script (which already default to examples/*). Their toml configs were already out of sync (e.g. stale psiedge=1.0 vs examples psiedge=0.800).

Benchmark script ran successfully

benchmark_edge_splines.jl ran against the existing gpec.h5 and generated 9 diagnostic plots (geometry, q/iota splines, dV/dψ, R(psin)/Z(psin), cross-section, rational surfaces, GSE residuals, far-edge profiles). The et vs q scan plot (plot 7) is missing because the existing gpec.h5 predates the psi_edge_scan HDF5 fields — a fresh GPEC run with psiedge < psilim is needed to generate it.

Remaining TODOs before merge

Fresh GPEC run to get psi_edge_scan data: Background Julia runs are being SIGKILL'd (exit 137) by macOS memory pressure — this happens in the background runner but not interactively. Try running directly: julia --project=. examples/DIIID-like_ideal_example/run_example.jl and verify it completes, then re-run benchmark_edge_splines.jl to get plot 7 (et vs q).
Benchmark plots: Once the fresh h5 is available, verify plot 7 shows the peak near ψ=0.9425, q≈4.2 as expected from the verified March 11 run.
Non-Hermitian warning: ~242 non-Hermitian W steps after trim (does not affect et[1]).
Test suite: runtests_equil.jl geqdsk filename is now corrected — run it to confirm the edge spline testset passes.

…mat_spline at local psi Removes the erroneous infrastructure that froze the Jacobian Fourier coefficients at psihigh for normalization in free_run! and free_compute_total. This caused SingularException for all ODE steps above psihigh because v†·J(psihigh)·v was indefinite for eigenvectors characteristic of the far-edge region. Changes: - ForceFreeStatesStructs.jl: remove ffit.jmat field from FourFitVars - Fourfit.jl make_matrix: remove core_end_idx, jmat_at_core_end snapshot, and ffit.jmat assignment; jmat_flat still built over full grid for spline - Free.jl free_run!: evaluate jmat_at_psilim = ffit.jmat_spline(psilim) at the actual plasma boundary instead of frozen snapshot - Free.jl free_compute_total: evaluate jmat_local = ffit.jmat_spline(psifac) at each scan step's psi; remove wrong comment block Also includes prior session changes: - make_metric: cap at psilim (not psihigh) so FGK splines use real above-psihigh geometry rather than ExtendExtrap frozen at psihigh - direct_build_xpoint_asymptotic!: grid ends at psilim (~0.994), removing sparse tail to 0.999999 (reduces far-edge nodes from 77 to 66) - EulerLagrange.jl: fix profiles -> equil.profiles in integrator_callback!; edge scan stores all steps psiedge->psilim with NaN for failures - GeneralizedPerturbedEquilibrium.jl: new top-level edge_scan/ HDF5 group with intuitive field names (total_energy, plasma_energy, etc.) - benchmarks: updated HDF5 keys; NaN-aware peak finding and subsampling - Project.toml: re-add PlotlyJS dependency Result: 346/6831 valid edge scan steps (0 before fix above psihigh); peak et=+3.056 at q=5.487 (psi=0.994). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…real vacuum computations The free_compute_wv_spline previously limited the vacuum response grid to [psiedge, psihigh], using ExtendExtrap to freeze wv at psihigh for all above-psihigh scan steps. This was incorrect for diverted plasmas where the plasma boundary can meaningfully exceed psihigh. Changes: - free_compute_wv_spline: extend grid to min(psi_scan_end, rzphi_xs[end]) for diverted plasmas. Above-psihigh grid points invert iota(psi)=1/q via Newton's method on the iota inner spline (same as make_matrix). Real VacuumInput computations at each psi use the X-point asymptotic geometry, giving physically correct wv for the above-psihigh scan range. - findmax_dW_edge!: retain scan_psi_max = psihigh + 0.5*(1-psihigh) to exclude ODE fixup artifacts at edge rational surfaces (q=6,7,...) where u grows exponentially before each Gaussian elimination, producing spuriously large eigenvalues. The rationale is now fixup noise exclusion, not ExtendExtrap wv limitation. Updated comment accordingly. Result: peak et = +3.064 at q=5.487 (was +3.056 with frozen wv), a 0.3% change reflecting the actual vacuum geometry at psi=0.994 vs psi=0.992. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…uncation Replace global-maximum search in findmax_dW_edge! with a first-local-maximum algorithm: find the peak of et between the first ascending and first descending zero crossings. The physical instability peak (et≈3.06 at q≈5.5) lies in this window; artifacts from degraded ODE solution in the far above-psihigh region (et≈25 at q≈11) occur after et returns to negative, so the algorithm naturally excludes them without any hardcoded psi cap. Also add ucrit/10 above psihigh in compute_solution_norms! (via new is_above_psihigh parameter) for more frequent Gaussian fixups near the separatrix, reducing per-fixup artifact magnitude. Fix benchmark stability plot (Plot 7) x-axis: cap at qlim_bench (q at psilim) rather than q_max_scan which extended to q~30 with no physical content. Verified result: Peak et=+3.243 at ψ=0.994, q=5.488 (Plasma=+1.778, Vacuum=+1.465). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…cation; fix evonly plot Remove the first-local-maximum logic from findmax_dW_edge! — it was based on an incorrect assumption about the scan structure. The physical picture shows repeated oscillations (stable → unstable → stable → ...) at each rational surface, so a simple zero-crossing criterion cannot reliably identify the "first" physical peak. Users control the scan range via edge_layer_width. Also fixes the evonly (vacuum_eigenvalue) diagnostic panel in Plot 7: - Clip values with |vacuum_eigenvalue| > 5 before plotting so the O(1) physically-relevant region is visible (previously the axis was dominated by spurious values of -2988 at psi > rzphi_xs[end]) - Title reports how many steps were clipped - Rename all "evonly" labels to "vacuum_eigenvalue" to match the h5 key Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ifacts Negative vacuum_eigenvalue values are unphysical (evonly >= 0 by construction) and come from wv computed via ExtendExtrap beyond rzphi_xs[end]. Filter to [0, clip_ev] instead of |val|<=clip_ev, and pin ylims from 0 to max+15% so the plot stays in the physically meaningful range. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Panel B (log-scale |ev|) was a debugging artifact — the vacuum energy contribution is already visible in the top panel (et/ep/ev components) and the wv matrix quality is shown in the bottom panel (vacuum_eigenvalue). Remove it, leaving a clean two-panel figure. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

… sing_find! The edge rational surface scan was duplicated in direct_build_xpoint_asymptotic! (with its own stopping criterion, push-then-check) and sing_find! (the canonical source, check-before-push). These had gotten out of sync, causing rzphi_xs[end] to differ from psilim by ~0.0006 in psi, leaving a gap where wv was computed via ExtendExtrap with unphysical results. Changes: - direct_build_xpoint_asymptotic!: remove internal surface scan; accept edge_surface_psis::AbstractVector{Float64} as required parameter (empty → no-op) - setup_equilibrium: no longer calls direct_build_xpoint_asymptotic! (deferred) - New equilibrium_extend_rzphi!(equil, edge_surface_psis): recovers node values from existing bicubic splines via nodal_derivs.partials, calls direct_build_xpoint_asymptotic! with sing_find! surfaces, rebuilds the four rzphi splines in-place, updates equil.rzphi_xs - GeneralizedPerturbedEquilibrium.jl: call equilibrium_extend_rzphi! after sing_find!/sing_lim!, so the rzphi grid covers exactly [psilow, psilim] with the same edge surfaces used for psilim — single canonical source, no gap - free_compute_wv_spline: extrap=NoExtrap() — wv energy oscillates in each rational window so no smooth extrapolation is valid beyond psilim Confirmed psilim = 0.999785 for DIIID-like example (last sing_find! surface before gap < edge_layer_width=1e-4; formats as 1.000 at 3dp). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ation via eval_q - Remove `InverseCubicSpline` wrapper type entirely. Store raw iota spline as `profiles.iota_spline` (with `_iota_hint`) and raw 1/(dV/dψ) spline as `profiles.dVdpsi_inv_spline` (with `_dVdpsi_inv_hint`) on `ProfileSplines`. - Remove `q_spline` pointer field (was a Union for pointer-swap, now dead since the chunk-level spline switch was removed in prior commit). - Remove `ProfileSplines` type parameter `I` (no longer needed without the wrapper). - Add `eval_q(profiles, psi; deriv, output_lower_derivs, hint)` as the single unified interface for q evaluation at any psi. Dispatches to `iota_spline` above psihigh (with chain-rule derivatives via `deriv` keyword) or `q_spline_direct` in the core. `output_lower_derivs=true` returns all derivatives through `deriv` in a single pass, sharing iota evaluations. - Fix bug in `sing_find!`: finite-difference q' approximation (1e-6 step) replaced by analytic chain rule via `eval_q(profiles, psi; deriv=1)`. - Collapse all manual chain-rule implementations in Sing.jl, Fourfit.jl, Free.jl, EulerLagrange.jl into single `eval_q` calls. - Collapse `compute_sing_mmat!` pre-computation of `_d1/_d2/_d3_iota` functors into `eval_q(...; deriv=3, output_lower_derivs=true)`. - Update all `q_spline`/`q_spline_iota_inverse`/`dVdpsi_spline_inv` references across 13 files to use the new field names and `eval_q`. - Remove PlotlyJS from Project.toml (incompatible with JSON 1.x). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Now that the pointer-swap pattern is gone and eval_q owns edge dispatch, q_spline unambiguously refers to the direct cubic spline over [psilow, psihigh]. The _direct suffix was only needed to distinguish it from the Union pointer; removing it restores the natural name across all call sites. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

logan-nc and others added 2 commits February 28, 2026 22:16

logan-nc force-pushed the edge_spline_extension branch from 59e6c27 to 746034f Compare March 1, 2026 03:17

logan-nc and others added 6 commits March 1, 2026 09:33

logan-nc and others added 8 commits March 3, 2026 13:01

logan-nc and others added 2 commits March 11, 2026 02:34

logan-nc and others added 3 commits March 11, 2026 17:08

Merge branch 'develop' of github.com:OpenFUSIONToolkit/JPEC into edge…

51e5ad3

…_spline_extension

BENCHMARKS - CLEANUP - Removes example inputs in benchmarks. All benc…

5bd1c1c

…hmark scripts can start from the actual examples and make any modifications dynamically - the replicated tomls and gfiles were already getting out of sync!

logan-nc and others added 4 commits March 11, 2026 23:30

logan-nc and others added 5 commits March 12, 2026 12:51

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

EQUIL + ForceFreeStates - Extend equilibrium splines to separatrix for diverted plasmas#179

EQUIL + ForceFreeStates - Extend equilibrium splines to separatrix for diverted plasmas#179
logan-nc wants to merge 30 commits intodevelopfrom
edge_spline_extension

logan-nc commented Mar 1, 2026 •

edited

Loading

Uh oh!

logan-nc commented Mar 2, 2026

Uh oh!

logan-nc commented Mar 9, 2026

Uh oh!

logan-nc commented Mar 11, 2026

Uh oh!

logan-nc commented Mar 11, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

logan-nc commented Mar 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Equilibrium: Edge inverse splines for diverted plasmas

ForceFreeStates: Integration beyond psihigh for diverted plasmas

Diagnostic output: et[1] vs psi in the edge zone

Benchmark script: benchmarks/benchmark_edge_splines.jl

Backward compatibility

End-to-end verification (DIIID-like example)

Test plan

Uh oh!

logan-nc commented Mar 2, 2026

Session update (2026-03-02): edge dW scan bugs fixed, full DIIID-like run verified

Uh oh!

logan-nc commented Mar 9, 2026

Edge dW scan investigation — current status ($(date '+%Y-%m-%d'))

What's working

The bug

Root cause investigation

What was committed today (42beb70)

What to try next

Uh oh!

logan-nc commented Mar 11, 2026

Session summary (2026-03-11)

What was completed this session

1. findmax_dW_edge! scan ranking fixed (Free.jl, EulerLagrange.jl)

2. make_metric grid capped at psihigh (Fourfit.jl)

3. jmat snapshot at psihigh (Fourfit.jl)

Verified end-to-end run

Remaining TODOs before merge

Uh oh!

logan-nc commented Mar 11, 2026

Session summary (2026-03-11, afternoon)

What was completed this session

Remaining TODOs before merge

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

logan-nc commented Mar 1, 2026 •

edited

Loading

Benchmark script: `benchmarks/benchmark_edge_splines.jl`

What was committed today (`42beb70`)

1. `findmax_dW_edge!` scan ranking fixed (`Free.jl`, `EulerLagrange.jl`)

2. `make_metric` grid capped at `psihigh` (`Fourfit.jl`)

3. `jmat` snapshot at `psihigh` (`Fourfit.jl`)