From f3583feded81ce2b67654f8c942ab154b0fcbfbe Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Fri, 9 Oct 2026 15:49:05 -0700 Subject: [PATCH 1/3] Equilibrium - BUGFIX - Support original, pow2 and rho grids for inverse equilibria grid_type = "original" set mpsi but never the grid, so it errored. It now uses the input's own psi nodes (Fortran inverse.f), axis included and unfloored; the axis node's 0/0 SFL and GSE terms are extrapolated from the innermost surfaces. Add pow2 (Fortran powspace, upper) and rho grids. Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/Equilibrium.jl | 1 + src/Equilibrium/EquilibriumTypes.jl | 1 + src/Equilibrium/InverseEquilibrium.jl | 70 +++++++++++++++++++-------- test/runtests_equil.jl | 28 +++++++++++ 4 files changed, 80 insertions(+), 20 deletions(-) diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 9bb0bd3ae..226c0be87 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -536,6 +536,7 @@ function equilibrium_gse!(equil::PlasmaEquilibrium) flux_fs[ipsi, itheta, 2] *= 2π * psio / f4 end end + equil.rzphi_xs[1] == 0.0 && extrap_axis_row!(flux_fs, equil.rzphi_xs, 1:2) # Create flux interpolants for Grad-Shafranov diagnostics @views flux_fs[:, end, :] .= flux_fs[:, 1, :] flux_opts = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap())) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 53f2ab8f2..339324e36 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -25,6 +25,7 @@ specified in the input. - `b0exp::Float64` - On-axis toroidal field normalization for CHEASE/EQDSK [T] - `grid_type::String` - Grid type for flux surface discretization ("auto" — two-pass measured-curvature refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; + inverse equilibria also "pow2", "rho" and "original" (the input file's own ψ nodes); "log_asymptotic" is a legacy alias for "auto") - `psilow::Float64` - Lower limit of normalized flux coordinate - `psihigh::Float64` - Requested upper limit of normalized flux coordinate; the value the diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 4a8d48bb6..8753af98a 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -39,6 +39,20 @@ function inverse_extrap(xx::Matrix{Float64}, ff::Matrix{Float64}, x::Float64)::V return f end +""" + extrap_axis_row!(fs, xs, ks; me=3) -> fs + +Fill the axis row (`xs[1] == 0`) of `fs[:, :, k]` for each `k` in `ks`, extrapolating the next +`me` surfaces in √ψ (∝ r), where the θ-dependent SFL and GSE terms are 0/0. +""" +function extrap_axis_row!(fs::Array{Float64,3}, xs::AbstractVector{Float64}, ks; me::Int=3) + rho = repeat(sqrt.(xs[2:(me+1)]), 1, size(fs, 2)) + for k in ks + fs[1, :, k] .= inverse_extrap(rho, fs[2:(me+1), :, k], 0.0) + end + return fs +end + function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing) @@ -124,14 +138,6 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No # Create 2D interpolants for r² and dη rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; itp_opts2d...) rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta; itp_opts2d...) - # c----------------------------------------------------------------------- - # c prepare new spline type for surface quantities. - # c----------------------------------------------------------------------- - - if grid_type in ["original", "orig"] - mpsi = size(sq_in._fs, 1) - 1 - end - # c----------------------------------------------------------------------- # c set up radial grid # c----------------------------------------------------------------------- @@ -165,6 +171,22 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No mpsi = 128 end sq_xs = [psilow + (psihigh - psilow) * (3(i / mpsi) - (i / mpsi)^3) / 2 for i in 0:mpsi] + elseif grid_type == "pow2" + if mpsi == 0 + mpsi = 128 + end + sq_xs = Utilities.powspace(psilow, psihigh, 2, mpsi + 1, "upper")[1, :] + elseif grid_type == "rho" + if mpsi == 0 + mpsi = 128 + end + # Uniform in √ψ from psihigh/(mpsi+1)² to psihigh; psilow is ignored (Fortran inverse.f) + sq_xs = [psihigh * i^2 / (mpsi + 1)^2 for i in 1:(mpsi+1)] + elseif grid_type in ("original", "orig") + # The input's own ψ nodes, axis included; psilow/psihigh are ignored (Fortran inverse.f, which + # takes the 1D input's nodes; they differ from the 2D ones only for efit_by_inversion) + sq_xs = copy(rz_in_xs) + mpsi = length(sq_xs) - 1 elseif grid_type == "uniform" if mpsi == 0 mpsi = 128 @@ -173,10 +195,10 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No else error("Unsupported grid_type: $grid_type") end - # Floor node spacing on the fixed grids (ldp/pow1 pack the edge below the integration-noise - # scale at high mpsi); the auto grid floors its refined grid in refined_psi_grid, and an - # override grid arrives already floored. - if override_psi_nodes === nothing && !(grid_type in ("auto", "log_asymptotic")) + # Floor node spacing on the fixed grids (ldp/pow1/pow2 pack the edge below the integration-noise + # scale at high mpsi); the auto grid floors its refined grid in refined_psi_grid, an override + # grid arrives already floored, and the original grid keeps the input's nodes. + if override_psi_nodes === nothing && !(grid_type in ("auto", "log_asymptotic", "original", "orig")) sq_xs = enforce_min_spacing(sq_xs, MIN_KNOT_SPACING) mpsi = length(sq_xs) - 1 end @@ -215,6 +237,14 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No for ipsi in 0:mpsi psifac = rzphi_xs[ipsi+1] sq_in(f_sq_in_buf, psifac; hint=sq_in_hint) + sq_fs[ipsi+1, 1] = f_sq_in_buf[1] * twopi + sq_fs[ipsi+1, 2] = f_sq_in_buf[2] + # Use the input q profile directly (from LAR ODE or CHEASE), matching the + # Fortran `inverse_chease4_run` convention (sq%fs(ipsi,4) = sq_in%f(3)). + # The field-line-integration-based q formula (spl_fsi * F / (2*twopi*psio)) + # is inaccurate for cylindrical LAR geometry. + sq_fs[ipsi+1, 4] = f_sq_in_buf[3] # q from input profile + psifac == 0.0 && continue # axis node: filled by extrapolation below spl_xs .= rzphi_ys for itheta in 0:mtheta theta = rzphi_ys[itheta+1] @@ -276,15 +306,13 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No spl_post(spl_post_buf, rzphi_ys[itheta+1]; hint=hint_post) @views rzphi_fs[ipsi+1, itheta+1, :] .= spl_post_buf[1:4] end - - sq_fs[ipsi+1, 1] = f_sq_in_buf[1] * twopi - sq_fs[ipsi+1, 2] = f_sq_in_buf[2] sq_fs[ipsi+1, 3] = spl_fsi[mtheta+1, 3] * twopi * pi # dV/d(psi) - # Use the input q profile directly (from LAR ODE or CHEASE), matching the - # Fortran `inverse_chease4_run` convention (sq%fs(ipsi,4) = sq_in%f(3)). - # The field-line-integration-based q formula (spl_fsi * F / (2*twopi*psio)) - # is inaccurate for cylindrical LAR geometry. - sq_fs[ipsi+1, 4] = f_sq_in_buf[3] # q from input profile + end + # An axis node (original grid) has r² = 0, where the SFL weights are 0/0: extrapolate it from the + # innermost surfaces as deta is above, θ-dependent fields in √ψ (∝ r) and dV/dψ in ψ. + if sq_xs[1] == 0.0 + extrap_axis_row!(rzphi_fs, sq_xs, 2:4; me) + sq_fs[1, 3] = inverse_extrap(reshape(sq_xs[2:(me+1)], :, 1), sq_fs[2:(me+1), 3:3], 0.0)[1] end sq = cubic_interp(sq_xs, Series(sq_fs); extrap=ExtendExtrap()) @@ -327,6 +355,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No v = zeros(Float64, 3, 3) for ipsi in 0:mpsi + sq_xs[ipsi+1] == 0.0 && continue # axis node: extrapolated below f_sq_vec = @view sq.y[ipsi+1, :] for itheta in 0:mtheta # Evaluate rzphi interpolants at grid points using nodal_derivs @@ -371,6 +400,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No eqfun_fs[ipsi+1, itheta+1, 3] = (v[2, 3] * v[3, 3] + f_sq_vec[4] * v[3, 3]^2) / (jacfac * eqfun_fs[ipsi+1, itheta+1, 1]^2) end end + sq_xs[1] == 0.0 && extrap_axis_row!(eqfun_fs, sq_xs, 1:3; me) # Create 2D interpolants for physics quantities (eqfun) eqfun_grid2d = (eqfun_xs, theta_range) @views eqfun_fs[:, end, :] .= eqfun_fs[:, 1, :] diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 2d33380c6..bf6cd676c 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -174,6 +174,34 @@ (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl end + @testset "Inverse grid types (original, pow2, rho)" begin + Eq = GeneralizedPerturbedEquilibrium.Equilibrium + cfg(grid) = Eq.EquilibriumConfig(; eq_type="chease_binary", eq_filename=joinpath(data_dir, "INP1_binary"), jac_type="boozer", + grid_type=grid, mpsi=16, psilow=0.01, psihigh=0.99, r0exp=6.8, b0exp=7.4) + splines(eq) = (eq.rzphi_rsquared, eq.rzphi_offset, eq.rzphi_nu, eq.rzphi_jac, eq.eqfun_B, eq.eqfun_metric1, eq.eqfun_metric2) + finite(eq) = all(isfinite, eq.profiles.dVdpsi_spline.y) && all(all(isfinite, s.nodal_derivs.partials) for s in splines(eq)) + eqs = Dict(grid => Eq.setup_equilibrium(cfg(grid)) for grid in ("original", "pow1", "pow2", "rho")) + + # original: the file's own ψ nodes, axis included, whatever psilow/psihigh/mpsi + eq_o = eqs["original"] + @test eq_o.rzphi_xs == Eq.read_chease_binary(cfg("original")).rz_in_xs + @test eq_o.rzphi_xs[1] == 0.0 + # The extrapolated axis node: B close to F/R₀ + b_axis = abs(eq_o.profiles.F_spline.y[1]) / (2π * eq_o.ro) + @test maximum(abs.(eq_o.eqfun_B.nodal_derivs.partials[1, 1, :] ./ b_axis .- 1)) < 1e-2 + + # pow1/pow2: Fortran powspace(psilow, psihigh, p, mpsi+1, "upper"); rho: psihigh·(i/(mpsi+1))², i = 1..mpsi+1 + t = (0:16) ./ 16 + @test eqs["pow1"].rzphi_xs ≈ 0.01 .+ 0.98 .* (3t .- t .^ 3) ./ 2 + @test eqs["pow2"].rzphi_xs ≈ 0.01 .+ 0.98 .* (15 / 8) .* (t .- 2t .^ 3 ./ 3 .+ t .^ 5 ./ 5) + @test eqs["rho"].rzphi_xs ≈ 0.99 .* ((1:17) ./ 17) .^ 2 + + for (grid, eq) in eqs + @test finite(eq) + @test isapprox(eq.profiles.q_spline(0.5), eqs["pow1"].profiles.q_spline(0.5); rtol=1e-3) + end + 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. From 44d73b2aff128e7ad5675342877456e651cf0694 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Mon, 12 Oct 2026 23:57:15 -0700 Subject: [PATCH 2/3] Equilibrium - MINOR - Drop the pow2 grid; the rho grid runs from psilow to psihigh Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/EquilibriumTypes.jl | 2 +- src/Equilibrium/InverseEquilibrium.jl | 11 +++-------- test/runtests_equil.jl | 11 ++++------- 3 files changed, 8 insertions(+), 16 deletions(-) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 339324e36..7c768e82e 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -25,7 +25,7 @@ specified in the input. - `b0exp::Float64` - On-axis toroidal field normalization for CHEASE/EQDSK [T] - `grid_type::String` - Grid type for flux surface discretization ("auto" — two-pass measured-curvature refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; - inverse equilibria also "pow2", "rho" and "original" (the input file's own ψ nodes); + inverse equilibria also "rho" (uniform in √ψ) and "original" (the input file's own ψ nodes); "log_asymptotic" is a legacy alias for "auto") - `psilow::Float64` - Lower limit of normalized flux coordinate - `psihigh::Float64` - Requested upper limit of normalized flux coordinate; the value the diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 8753af98a..7ad5b9969 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -171,17 +171,12 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No mpsi = 128 end sq_xs = [psilow + (psihigh - psilow) * (3(i / mpsi) - (i / mpsi)^3) / 2 for i in 0:mpsi] - elseif grid_type == "pow2" - if mpsi == 0 - mpsi = 128 - end - sq_xs = Utilities.powspace(psilow, psihigh, 2, mpsi + 1, "upper")[1, :] elseif grid_type == "rho" if mpsi == 0 mpsi = 128 end - # Uniform in √ψ from psihigh/(mpsi+1)² to psihigh; psilow is ignored (Fortran inverse.f) - sq_xs = [psihigh * i^2 / (mpsi + 1)^2 for i in 1:(mpsi+1)] + # Uniform in √ψ from psilow to psihigh + sq_xs = [(sqrt(psilow) + (sqrt(psihigh) - sqrt(psilow)) * i / mpsi)^2 for i in 0:mpsi] elseif grid_type in ("original", "orig") # The input's own ψ nodes, axis included; psilow/psihigh are ignored (Fortran inverse.f, which # takes the 1D input's nodes; they differ from the 2D ones only for efit_by_inversion) @@ -195,7 +190,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No else error("Unsupported grid_type: $grid_type") end - # Floor node spacing on the fixed grids (ldp/pow1/pow2 pack the edge below the integration-noise + # Floor node spacing on the fixed grids (ldp/pow1 pack the edge below the integration-noise # scale at high mpsi); the auto grid floors its refined grid in refined_psi_grid, an override # grid arrives already floored, and the original grid keeps the input's nodes. if override_psi_nodes === nothing && !(grid_type in ("auto", "log_asymptotic", "original", "orig")) diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index bf6cd676c..18b236dba 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -174,13 +174,13 @@ (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl end - @testset "Inverse grid types (original, pow2, rho)" begin + @testset "Inverse grid types (original, rho)" begin Eq = GeneralizedPerturbedEquilibrium.Equilibrium cfg(grid) = Eq.EquilibriumConfig(; eq_type="chease_binary", eq_filename=joinpath(data_dir, "INP1_binary"), jac_type="boozer", grid_type=grid, mpsi=16, psilow=0.01, psihigh=0.99, r0exp=6.8, b0exp=7.4) splines(eq) = (eq.rzphi_rsquared, eq.rzphi_offset, eq.rzphi_nu, eq.rzphi_jac, eq.eqfun_B, eq.eqfun_metric1, eq.eqfun_metric2) finite(eq) = all(isfinite, eq.profiles.dVdpsi_spline.y) && all(all(isfinite, s.nodal_derivs.partials) for s in splines(eq)) - eqs = Dict(grid => Eq.setup_equilibrium(cfg(grid)) for grid in ("original", "pow1", "pow2", "rho")) + eqs = Dict(grid => Eq.setup_equilibrium(cfg(grid)) for grid in ("original", "pow1", "rho")) # original: the file's own ψ nodes, axis included, whatever psilow/psihigh/mpsi eq_o = eqs["original"] @@ -190,11 +190,8 @@ b_axis = abs(eq_o.profiles.F_spline.y[1]) / (2π * eq_o.ro) @test maximum(abs.(eq_o.eqfun_B.nodal_derivs.partials[1, 1, :] ./ b_axis .- 1)) < 1e-2 - # pow1/pow2: Fortran powspace(psilow, psihigh, p, mpsi+1, "upper"); rho: psihigh·(i/(mpsi+1))², i = 1..mpsi+1 - t = (0:16) ./ 16 - @test eqs["pow1"].rzphi_xs ≈ 0.01 .+ 0.98 .* (3t .- t .^ 3) ./ 2 - @test eqs["pow2"].rzphi_xs ≈ 0.01 .+ 0.98 .* (15 / 8) .* (t .- 2t .^ 3 ./ 3 .+ t .^ 5 ./ 5) - @test eqs["rho"].rzphi_xs ≈ 0.99 .* ((1:17) ./ 17) .^ 2 + # rho: uniform in √ψ between psilow and psihigh + @test eqs["rho"].rzphi_xs ≈ (0.1 .+ (sqrt(0.99) - 0.1) .* (0:16) ./ 16) .^ 2 for (grid, eq) in eqs @test finite(eq) From f0497d3f73d24c9be01f286681740cc95001e81b Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 13 Oct 2026 01:29:27 -0700 Subject: [PATCH 3/3] Equilibrium - MINOR - Tighten the inverse grid-type comments and tests Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/EquilibriumTypes.jl | 3 ++- src/Equilibrium/InverseEquilibrium.jl | 14 +++++++------- test/runtests_equil.jl | 17 +++++++---------- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 7c768e82e..43d182fbd 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -25,7 +25,8 @@ specified in the input. - `b0exp::Float64` - On-axis toroidal field normalization for CHEASE/EQDSK [T] - `grid_type::String` - Grid type for flux surface discretization ("auto" — two-pass measured-curvature refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; - inverse equilibria also "rho" (uniform in √ψ) and "original" (the input file's own ψ nodes); + inverse equilibria also "rho" (uniform in √ψ from psilow to psihigh) and "original" (the input + file's own ψ nodes; ignores psilow and mpsi); "log_asymptotic" is a legacy alias for "auto") - `psilow::Float64` - Lower limit of normalized flux coordinate - `psihigh::Float64` - Requested upper limit of normalized flux coordinate; the value the diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 7ad5b9969..a551aa5d8 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -42,8 +42,8 @@ end """ extrap_axis_row!(fs, xs, ks; me=3) -> fs -Fill the axis row (`xs[1] == 0`) of `fs[:, :, k]` for each `k` in `ks`, extrapolating the next -`me` surfaces in √ψ (∝ r), where the θ-dependent SFL and GSE terms are 0/0. +Fill the axis row (`xs[1] == 0`, where the θ-dependent SFL and GSE terms are 0/0) of `fs[:, :, k]` +for each `k` in `ks`, extrapolating from the next `me` surfaces in √ψ (∝ r). """ function extrap_axis_row!(fs::Array{Float64,3}, xs::AbstractVector{Float64}, ks; me::Int=3) rho = repeat(sqrt.(xs[2:(me+1)]), 1, size(fs, 2)) @@ -175,11 +175,11 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No if mpsi == 0 mpsi = 128 end - # Uniform in √ψ from psilow to psihigh + # Uniform in √ψ from psilow to psihigh (Fortran's rho ignores psilow) sq_xs = [(sqrt(psilow) + (sqrt(psihigh) - sqrt(psilow)) * i / mpsi)^2 for i in 0:mpsi] elseif grid_type in ("original", "orig") - # The input's own ψ nodes, axis included; psilow/psihigh are ignored (Fortran inverse.f, which - # takes the 1D input's nodes; they differ from the 2D ones only for efit_by_inversion) + # The input's own ψ nodes; psilow and mpsi are ignored, psihigh still sets psilim (as Fortran, + # which takes the 1D input's nodes; they differ from these only for efit_by_inversion) sq_xs = copy(rz_in_xs) mpsi = length(sq_xs) - 1 elseif grid_type == "uniform" @@ -303,8 +303,8 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No end sq_fs[ipsi+1, 3] = spl_fsi[mtheta+1, 3] * twopi * pi # dV/d(psi) end - # An axis node (original grid) has r² = 0, where the SFL weights are 0/0: extrapolate it from the - # innermost surfaces as deta is above, θ-dependent fields in √ψ (∝ r) and dV/dψ in ψ. + # An axis node (ψ = 0) has r² = 0, where the SFL weights are 0/0: extrapolate it from the innermost + # surfaces as deta is above, θ-dependent fields in √ψ (∝ r) and dV/dψ in ψ. if sq_xs[1] == 0.0 extrap_axis_row!(rzphi_fs, sq_xs, 2:4; me) sq_fs[1, 3] = inverse_extrap(reshape(sq_xs[2:(me+1)], :, 1), sq_fs[2:(me+1), 3:3], 0.0)[1] diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 18b236dba..b8ce3f74f 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -180,23 +180,20 @@ grid_type=grid, mpsi=16, psilow=0.01, psihigh=0.99, r0exp=6.8, b0exp=7.4) splines(eq) = (eq.rzphi_rsquared, eq.rzphi_offset, eq.rzphi_nu, eq.rzphi_jac, eq.eqfun_B, eq.eqfun_metric1, eq.eqfun_metric2) finite(eq) = all(isfinite, eq.profiles.dVdpsi_spline.y) && all(all(isfinite, s.nodal_derivs.partials) for s in splines(eq)) - eqs = Dict(grid => Eq.setup_equilibrium(cfg(grid)) for grid in ("original", "pow1", "rho")) - # original: the file's own ψ nodes, axis included, whatever psilow/psihigh/mpsi - eq_o = eqs["original"] + # original: the file's own ψ nodes (this file starts at the axis), whatever psilow/mpsi + eq_o = Eq.setup_equilibrium(cfg("original")) @test eq_o.rzphi_xs == Eq.read_chease_binary(cfg("original")).rz_in_xs @test eq_o.rzphi_xs[1] == 0.0 + @test finite(eq_o) # The extrapolated axis node: B close to F/R₀ b_axis = abs(eq_o.profiles.F_spline.y[1]) / (2π * eq_o.ro) @test maximum(abs.(eq_o.eqfun_B.nodal_derivs.partials[1, 1, :] ./ b_axis .- 1)) < 1e-2 - # rho: uniform in √ψ between psilow and psihigh - @test eqs["rho"].rzphi_xs ≈ (0.1 .+ (sqrt(0.99) - 0.1) .* (0:16) ./ 16) .^ 2 - - for (grid, eq) in eqs - @test finite(eq) - @test isapprox(eq.profiles.q_spline(0.5), eqs["pow1"].profiles.q_spline(0.5); rtol=1e-3) - end + # rho: starts at psilow (Fortran's rho starts at psihigh/(mpsi+1)²) + eq_r = Eq.setup_equilibrium(cfg("rho")) + @test eq_r.rzphi_xs[[1, end]] ≈ [0.01, 0.99] + @test finite(eq_r) end @testset "EFIT Method Consistency" begin