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()