From b3f412f352e73424af36f6dc93007a65263352f4 Mon Sep 17 00:00:00 2001 From: d-burg Date: Tue, 6 Oct 2026 00:23:57 -0400 Subject: [PATCH] =?UTF-8?q?Equilibrium=20-=20BUGFIX!=20-=20Build=20F=20and?= =?UTF-8?q?=20P=20of=20g-file=20and=20IMAS=20equilibria=20from=20the=20tab?= =?UTF-8?q?ulated=20FF=E2=80=B2=20and=20p=E2=80=B2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The g-file and IMAS readers took F and P from the tabulated profiles and left their derivatives to the cubic spline. F varies by only a few percent across the plasma, so rounding in the table (single precision in TokaMaker and classic EFIT output) leaves F′ nearly right in value but wrong in slope at the node scale, and Δ′ responds to that slope. `profile_source = "derivatives"` (the new default) integrates the file's FF′ and p′ inward from the boundary values of F and P, always for both profiles together. `profile_source = "values"` keeps the previous behaviour, and the reader falls back to it with a warning when the derivative arrays are absent or unusable. Co-Authored-By: Claude Opus 5.5 --- docs/src/conventions.md | 8 +-- docs/src/equilibrium.md | 25 +++++++ src/Equilibrium/EquilibriumTypes.jl | 13 +++- src/Equilibrium/ReadEquilibrium.jl | 106 ++++++++++++++++++++++++++-- test/runtests_equil.jl | 50 +++++++++++++ test/runtests_imas.jl | 29 ++++++++ 6 files changed, 216 insertions(+), 15 deletions(-) diff --git a/docs/src/conventions.md b/docs/src/conventions.md index cc1afd539..7deac2b42 100644 --- a/docs/src/conventions.md +++ b/docs/src/conventions.md @@ -106,11 +106,9 @@ helicity = bt_sign * Int(sign(crnt)) # src/ForcingTerms/CoilFourier.jl ## ``F = R B_\phi`` The poloidal-current function ``F = R B_\phi`` from the Grad-Shafranov equation is **forced -positive** via `abs` (`read_eq_efit` in `src/Equilibrium/ReadEquilibrium.jl`): - -```julia -abs.(fpol_data) -``` +positive**: the g-file and IMAS readers take its magnitude (`file_profiles` in +`src/Equilibrium/ReadEquilibrium.jl`), whether it is integrated from the file's ``FF'`` or read +from the tabulated values. The code always works with ``|F|``; the sign of ``B_t`` is carried separately (`fpol_sign`, and the `bt_sign`/`crnt` helicity inputs). diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index ecff97671..aa6c6a757 100644 --- a/docs/src/equilibrium.md +++ b/docs/src/equilibrium.md @@ -55,6 +55,31 @@ The Cartesian evaluation grid is clipped to the separatrix bounding box and its resolution is set adaptively from a bilinear interpolation error bound, so no manual tuning is needed. +## F and P from a g-file or IMAS equilibrium + +A g-file tabulates both the profiles (`FPOL`, `PRES`) and their flux derivatives (`FFPRIM`, +`PPRIME`); an IMAS equilibrium does the same (`f`, `pressure`, `f_df_dpsi`, `dpressure_dpsi`). +The stability matrices use ``F'`` directly and respond to how ``F'`` varies across each rational +surface, i.e. to the current gradient. ``F`` itself changes by only a few percent across the +plasma, so differentiating a tabulated ``F`` amplifies its rounding error by roughly +``1/(h^2\delta)`` in the slope of ``F'`` (``h`` the node spacing, ``\delta`` the fractional +variation of ``F``). On a 257-node file written in single precision that is an error of order +10 % in the current gradient, and several units in ``\Delta'``. + +`profile_source` in `[Equilibrium]` selects where the two profiles come from: + +| Value | Behaviour | +|---|---| +| `"derivatives"` (default) | ``F`` and ``P`` are the integrals of the file's ``FF'`` and ``p'``, anchored at the boundary values of the tabulated ``F`` and ``P`` | +| `"values"` | The tabulated ``F`` and ``P`` are used as written | + +With `"derivatives"` both profiles are always rebuilt together, so the pair stays consistent in +the force balance, and everything downstream still differentiates one spline. The reader logs the +largest interior difference between the rebuilt and the tabulated profiles and warns when it +exceeds `1e-4` of ``F`` or 5 % of the peak pressure. It falls back to `"values"` with a warning +when the derivative arrays are absent, all zero, or opposite in sign for only one of the two +profiles. The two end nodes are used as written. + ## Radial grid packing With `grid_type = "auto"` and `mpsi = 0` (the defaults; `"log_asymptotic"` is a legacy alias), the radial grid is diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index ac5a212fc..49a6d1367 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -41,6 +41,11 @@ specified in the input. - `etol::Float64` - Error tolerance for equilibrium solver - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations) - `use_galgrid::Bool` - Use the same grid as galerkin method + - `imas_cocos::Int` - COCOS convention of the input IMAS `dd.equilibrium` (11 = IMAS standard, 2 = internal) + - `profile_source::String` - Which 1D arrays of a g-file or IMAS equilibrium define F and P. + `"derivatives"` (default) integrates the file's FF′ and p′ inward from the boundary values of + F and P; `"values"` uses the tabulated F and P directly. The reader falls back to `"values"` + with a warning when the derivatives are absent or unusable. """ @kwdef struct EquilibriumConfig eq_type::String = "efit" @@ -72,8 +77,8 @@ specified in the input. force_termination::Bool = false use_galgrid::Bool = true - # IMAS-specific: expected COCOS convention of the input dd.equilibrium (11=IMAS standard, 2=GPEC internal) imas_cocos::Int = 11 + profile_source::String = "derivatives" """ Modified internal constructor that enforces self consistency within the inputs @@ -84,7 +89,7 @@ specified in the input. function EquilibriumConfig(eq_type, eq_filename, r0exp, b0exp, jac_type, _, _, _, _, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol, - force_termination, use_galgrid, imas_cocos) + force_termination, use_galgrid, imas_cocos, profile_source) if jac_type == "hamada" @info "Forcing hamada coordinate jacobian exponents: power_*" power_b = 0 @@ -147,10 +152,12 @@ specified in the input. @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0" end psihigh = min(psihigh, 1.0) + profile_source in ("derivatives", "values") || + error("Cannot recognize profile_source = $(profile_source); use \"derivatives\" or \"values\"") return new(eq_type, eq_filename, r0exp, b0exp, jac_type, power_bp, power_b, power_r, power_rc, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, etol, - force_termination, use_galgrid, imas_cocos) + force_termination, use_galgrid, imas_cocos, profile_source) end end diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index 755a91aca..f114bfba7 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -40,6 +40,90 @@ function _read_1d_gfile_format(lines_block::Vector{String}, num_values::Int) return parsed_values end +# Largest interior disagreement between the tabulated and the integrated profiles the reader accepts +# without a warning: relative for F, as a fraction of the peak pressure for P. +const PROFILE_F_MISMATCH_WARN = 1e-4 +const PROFILE_P_MISMATCH_WARN = 5e-2 + +""" + integrate_profile_derivatives(xs, f, p, ffprime, pprime) -> Union{Nothing,NamedTuple} + +Rebuild F and P on the nodes `xs` of a normalized flux coordinate (0 at the axis, 1 at the +boundary) by integrating the tabulated `ffprime` = d(F²/2)/dx and `pprime` = dP/dx inward from the +boundary values of `f` and `p`: + + F²(x) = F(1)² − 2∫ₓ¹ FF′ dx, P(x) = P(1) − ∫ₓ¹ P′ dx. + +F varies by only a few percent across the plasma, so rounding in a tabulated F is amplified by +roughly 1/(h²δ) in the slope of F′ (h the node spacing, δ the fractional variation of F), which is +the current gradient the stability matrices respond to. The tabulated derivatives carry that +information directly, and they are the source terms the Grad-Shafranov solution ψ(R,Z) was computed +from. F and P are rebuilt together so that the pair stays consistent in the force balance. + +Returns `(; f, p, f_mismatch, p_mismatch)` with `f` positive, where the mismatches are the largest +interior differences from the tabulated values (relative for F, as a fraction of the peak pressure +for P; the two end nodes are excluded). Returns `nothing` when the derivatives cannot be used: a +derivative array that is all zero or non-finite while its profile varies, a non-increasing `xs`, a +non-positive F², or derivative signs that match the tabulated profiles for only one of the two. +""" +function integrate_profile_derivatives(xs::AbstractVector{<:Real}, f::AbstractVector{<:Real}, p::AbstractVector{<:Real}, + ffprime::AbstractVector{<:Real}, pprime::AbstractVector{<:Real}) + n = length(xs) + (n >= 4 && length(f) == n && length(p) == n && length(ffprime) == n && length(pprime) == n) || return nothing + (all(isfinite, ffprime) && all(isfinite, pprime) && all(>(0), diff(xs))) || return nothing + f_abs = abs.(f) + f_has, p_has = any(!iszero, ffprime), any(!iszero, pprime) + # A zero derivative array is only believable for a profile that is itself flat + (f_has || maximum(f_abs) - minimum(f_abs) <= eps(maximum(f_abs))) || return nothing + (p_has || maximum(p) == minimum(p)) || return nothing + (f_has || p_has) || return nothing + + nodes = collect(Float64, xs) + tail(y) = (c = FastInterpolations.cumulative_integrate(cubic_interp(nodes, collect(Float64, y))); c[end] .- c) + f_tail, p_tail = tail(ffprime), tail(pprime) + interior = 2:(n-1) + p_scale = maximum(abs, p) + f_error(sgn) = maximum(i -> abs(sqrt(max(f_abs[end]^2 - 2 * sgn * f_tail[i], 0.0)) - f_abs[i]) / f_abs[i], interior) + p_error(sgn) = p_scale == 0 ? 0.0 : maximum(i -> abs(p[end] - sgn * p_tail[i] - p[i]), interior) / p_scale + + # Writers differ in the sign convention of the flux derivative, so take the sign that reproduces + # the tabulated profiles and require F and P to agree on it. + f_sign = f_error(1) <= f_error(-1) ? 1 : -1 + p_sign = p_error(1) <= p_error(-1) ? 1 : -1 + f_has && p_has && f_sign != p_sign && return nothing + sgn = f_has ? f_sign : p_sign + + f2 = f_abs[end]^2 .- 2 .* sgn .* f_tail + all(>(0), f2) || return nothing + return (; f=sqrt.(f2), p=p[end] .- sgn .* p_tail, f_mismatch=f_error(sgn), p_mismatch=p_error(sgn)) +end + +""" + file_profiles(config, xs, f, p, ffprime, pprime) -> (f_nodes, p_nodes) + +F (as a magnitude) and P on the 1D nodes of a file-based direct equilibrium, from the arrays that +`config.profile_source` selects. With `"derivatives"` the profiles come from +`integrate_profile_derivatives`, and the tabulated values are used with a warning when +that is not possible. +""" +function file_profiles(config::EquilibriumConfig, xs::AbstractVector{<:Real}, f::AbstractVector{<:Real}, p::AbstractVector{<:Real}, + ffprime::AbstractVector{<:Real}, pprime::AbstractVector{<:Real}) + config.profile_source == "values" && return abs.(f), collect(Float64, p) + built = integrate_profile_derivatives(xs, f, p, ffprime, pprime) + if built === nothing + @warn "profile_source = \"derivatives\", but the file's FF′ and p′ are absent or unusable; using its tabulated F and P" + return abs.(f), collect(Float64, p) + end + msg = "F and P integrated from the file's FF′ and p′; largest interior difference from the tabulated values: " * + "$(@sprintf("%.1e", built.f_mismatch)) of F, $(@sprintf("%.1e", built.p_mismatch)) of the peak pressure" + if built.f_mismatch > PROFILE_F_MISMATCH_WARN || built.p_mismatch > PROFILE_P_MISMATCH_WARN + @warn msg * ". The file's profiles and derivatives disagree; profile_source = \"values\" selects the tabulated profiles." + else + @info msg + end + return built.f, built.p +end + """ _read_efit(equil_in) @@ -91,10 +175,13 @@ function read_efit(config::EquilibriumConfig) psi_rz = reshape(psi_flat_vec, nw, nh) # --- Create 1D Profile Spline (sq_in) --- + psio_signed = sibry - simag psi_norm_grid = range(0.0, 1.0; length=nw) + # FFPRIM and PPRIME are derivatives in the file's ψ; psio_signed converts them to ψ_norm + f_nodes, p_nodes = file_profiles(config, psi_norm_grid, fpol_data, pres_data, ffprime_data .* psio_signed, pprime_data .* psio_signed) sq_fs_nodes = hcat( - abs.(fpol_data), - max.(pres_data .* mu0, 0.0), + f_nodes, + max.(p_nodes .* mu0, 0.0), qprof_data, sqrt.(psi_norm_grid) ) @@ -102,7 +189,6 @@ function read_efit(config::EquilibriumConfig) sq_in = cubic_interp(sq_xs, Series(sq_fs_nodes); extrap=ExtendExtrap()) # --- Process and Normalize 2D Psi Data --- - psio_signed = sibry - simag psi_proc = (sibry .- psi_rz) psio = abs(psio_signed) # Ensure psi at the magnetic axis is positive relative to the boundary @@ -495,11 +581,17 @@ function read_imas(config::EquilibriumConfig, dd) psi_norm_grid = range(0.0, 1.0; length=nw) # Build equilibrium source terms for spline interpolation - # abs(f_1d): F(ψ) can be negative depending on toroidal field direction convention; - # take absolute value since GPEC uses magnitude F = R·|Bt| + # F(ψ) can be negative depending on toroidal field direction convention; GPEC uses the + # magnitude F = R·|Bt|. f_df_dpsi and dpressure_dpsi are derivatives in the stored ψ, so the + # stored ψ span converts them to its normalized coordinate (absent arrays read as zero). + psi_stored = eqt.profiles_1d.psi + psi_span = psi_stored[end] - psi_stored[1] + stored_or_zero(name) = (v = getproperty(eqt.profiles_1d, name, Float64[]); length(v) == nw ? v .* psi_span : zeros(nw)) + f_nodes, p_nodes = file_profiles(config, (psi_stored .- psi_stored[1]) ./ psi_span, f_1d, p_1d, + stored_or_zero(:f_df_dpsi), stored_or_zero(:dpressure_dpsi)) sq_fs_nodes = hcat( - abs.(f_1d), - max.(p_1d .* mu0, 0.0), + f_nodes, + max.(p_nodes .* mu0, 0.0), q_1d, sqrt.(psi_norm_grid) ) diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index b4a257956..95063799c 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -147,6 +147,56 @@ (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl end + @testset "F and P from file derivatives" begin + Eq = GeneralizedPerturbedEquilibrium.Equilibrium + # Smooth profiles with exact derivatives; F varies by ~2 % like a tokamak F + xs = range(0.0, 1.0; length=257) + f_exact = @. sqrt(9.0 + 0.4 * (1 - xs)^2 - 0.1 * (1 - xs)^3) + ffprime = @. -0.4 * (1 - xs) + 0.15 * (1 - xs)^2 + f2_exact = @. (ffprime^2 / f_exact^2 * -1 + 0.4 - 0.3 * (1 - xs)) / f_exact # d²F/dx² + p_exact = @. 5e4 * (1 - xs)^2 + 300.0 + pprime = @. -1e5 * (1 - xs) + f_single = Float64.(Float32.(f_exact)) # a table written in single precision + + built = Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, pprime) + @test maximum(abs.(built.f .- f_exact) ./ f_exact) < 1e-7 # boundary anchor carries the rounding + @test maximum(abs.(built.p .- p_exact)) < 1e-6 * 5e4 + @test built.f_mismatch < 1e-6 && built.p_mismatch < 1e-9 + + # The slope of F′ is what single-precision rounding ruins and the derivative route keeps + slope_error(f) = maximum(abs.(deriv2(cubic_interp(collect(xs), f)).(xs[20:end-20]) .- f2_exact[20:end-20])) + @test slope_error(built.f) < 1e-2 * slope_error(f_single) + + # Either sign convention of the flux derivative gives the same profiles + flipped = Eq.integrate_profile_derivatives(xs, f_single, p_exact, -ffprime, -pprime) + @test flipped.f ≈ built.f && flipped.p ≈ built.p + @test Eq.integrate_profile_derivatives(xs, -f_single, p_exact, ffprime, pprime).f ≈ built.f + + # Unusable derivatives: absent, non-finite, or signs that suit only one of the two profiles + @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, zero(ffprime), zero(pprime)) === nothing + @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, zero(pprime)) === nothing + @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, fill(NaN, 257), pprime) === nothing + @test Eq.integrate_profile_derivatives(xs, f_single, p_exact, ffprime, -pprime) === nothing + + derivs = Eq.EquilibriumConfig(; profile_source="derivatives") + values = Eq.EquilibriumConfig(; profile_source="values") + @test_throws ErrorException Eq.EquilibriumConfig(; profile_source="spline") + @test Eq.file_profiles(values, xs, -f_single, p_exact, ffprime, pprime) == (f_single, p_exact) + @test (@test_logs (:warn, r"absent or unusable") Eq.file_profiles(derivs, xs, f_single, p_exact, zero(ffprime), zero(pprime))) == (f_single, p_exact) + @test (@test_logs (:info, r"integrated from") Eq.file_profiles(derivs, xs, f_single, p_exact, ffprime, pprime))[1] ≈ built.f + # A pressure table 20 % off its own derivative is reported, not silently accepted + @test_logs (:warn, r"disagree") Eq.file_profiles(derivs, xs, f_single, p_exact, ffprime, 1.2 .* pprime) + + # The g-file reader keeps the boundary values and moves the interior only at the rounding level + gfile = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode.geqdsk") + table(source) = Eq.read_efit(Eq.EquilibriumConfig(; eq_type="efit", eq_filename=gfile, profile_source=source)).ingest.sq_fs + tabulated, integrated = table("values"), table("derivatives") + @test integrated[end, 1:2] == tabulated[end, 1:2] + @test integrated[:, 1] != tabulated[:, 1] + @test maximum(abs.(integrated[:, 1] .- tabulated[:, 1]) ./ tabulated[:, 1]) < 1e-5 + @test integrated[:, 3:4] == tabulated[:, 3:4] + end + @testset "EFIT Method Consistency" begin # All three methods solve the same equilibrium — q-profiles should broadly agree. # Tolerance is 10% to allow for method-specific discretisation differences. diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl index 2d82d5b0e..41c47a032 100644 --- a/test/runtests_imas.jl +++ b/test/runtests_imas.jl @@ -72,6 +72,35 @@ using GeneralizedPerturbedEquilibrium.Equilibrium @test isapprox(result.psio, psi_bnd_int; rtol=1e-6) end + # read_imas — F and P come from the stored derivatives when the dd carries them + @testset "read_imas: profiles from stored derivatives" begin + for cocos in (11, 2) + dd, _ = make_mock_dd(; cocos=cocos) + prof = dd.equilibrium.time_slice[1].profiles_1d + nw = length(prof.psi) + # Derivatives in the stored ψ that match the mock's flat F and linear pressure + prof.f_df_dpsi = zeros(nw) + prof.dpressure_dpsi = fill(-1e4 / (prof.psi[end] - prof.psi[1]), nw) + # A pressure table carrying an error the derivatives do not have + exact_pressure = copy(prof.pressure) + prof.pressure = exact_pressure .* (1 .+ 1e-3 .* sinpi.(range(0, 8; length=nw))) + + mu0 = 4π * 1e-7 + read_mu0p(source) = Equilibrium.read_imas( + EquilibriumConfig(; eq_type="imas", eq_filename="mock", imas_cocos=cocos, profile_source=source), dd).ingest.sq_fs[:, 2] + @test read_mu0p("derivatives") ≈ mu0 .* exact_pressure rtol = 1e-9 atol = 1e-12 + @test read_mu0p("values") == mu0 .* prof.pressure + end + end + + # read_imas — a dd without the derivative arrays falls back to the tabulated profiles + @testset "read_imas: fallback without stored derivatives" begin + dd, _ = make_mock_dd() + config = EquilibriumConfig(; eq_type="imas", eq_filename="mock") + result = @test_logs (:warn, r"absent or unusable") match_mode = :any Equilibrium.read_imas(config, dd) + @test result.ingest.sq_fs[:, 1] == fill(5.0, 64) + end + # Test 3: read_imas — splines are callable after construction @testset "read_imas: splines callable" begin dd, _ = make_mock_dd()