Equilibrium - FEATURE - Read inverse i-files such as TokaMaker save_ifile output

docs/src/workflow.md

@@ -32,6 +32,7 @@ The single `gpec.toml` file supplies user-selected options to every module. The
 - Equilibrium data file in one of the supported formats:
   - `efit` — EFIT g-file from experimental reconstruction
   - `chease` / `chease2` — CHEASE equilibrium code output
+  - `ifile` — inverse i-file (e.g. TokaMaker `save_ifile`)
   - `lar` — Large aspect ratio analytical model
   - `sol` — Solov'ev analytical equilibrium
 - Kinetic profiles file (planned) — temperature and density profiles for the kinetic analysis path

src/Equilibrium/Equilibrium.jl

@@ -103,12 +103,15 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin
         eq_input = additional_input
         eq_input.config = eq_config
         eq_input.psihigh_resolved = eq_config.psihigh
+        resolve_psihigh!(eq_input)
     elseif eq_type in EFIT_KINDS
         eq_input = resolve_psihigh!(read_efit(eq_config))
     elseif eq_type in ["chease2", "chease_ascii"]
         eq_input = read_chease_ascii(eq_config)
     elseif eq_type in ["chease", "chease_binary"]
         eq_input = read_chease_binary(eq_config)
+    elseif eq_type == "ifile"
+        eq_input = resolve_psihigh!(read_ifile(eq_config))
     elseif haskey(ANALYTIC_EQ, eq_type)
         # Analytic kinds (sol/lar/tj_analytic[_direct]) dispatch off the ANALYTIC_EQ registry.
         # Their parameters live in the embedded `[*_INPUT]` section and are passed in as the
@@ -490,6 +493,10 @@ end
 Diagnoses the Grad-Shafranov solution by computing the residual of the
 Grad-Shafranov equation across the grid and writing diagnostic data to HDF5 files.
 Performs the same function as equil_out_gse in the Fortran code.
+
+Returns `(; xs, ys, flux_x, flux_y, source, total, error, errori)` on the `rzphi` (ψ, θ) nodes: the two
+flux-divergence terms, the source, the residual `total = flux_x - flux_y + source`, `error = |total|` over
+the largest term on the grid, and `errori`, the magnitude of the θ-integrated residual on each surface.
 """
 function equilibrium_gse!(equil::PlasmaEquilibrium)
 
@@ -648,6 +655,7 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
             file["errlogi"] = Float32.(errlogi)
         end
     end
+    return (; xs=equil.rzphi_xs, ys=equil.rzphi_ys, flux_x=flux_fsx[:, :, 1], flux_y=flux_fsy[:, :, 2], source, total, error, errori=vec(errori))
 end
 
 end # module Equilibrium

src/Equilibrium/EquilibriumTypes.jl

@@ -374,8 +374,8 @@ end
     InverseIngest
 
 The serializable raw arrays and scalars captured by an inverse-equilibrium reader
-(`read_chease_ascii`, `read_chease_binary`) — everything needed to rebuild an
-`InverseRunInput`'s splines without re-reading the original CHEASE file. Stored on
+(`read_chease_ascii`, `read_chease_binary`, `read_ifile`) — everything needed to rebuild an
+`InverseRunInput`'s splines without re-reading the original CHEASE file or i-file. Stored on
 `InverseRunInput.ingest` / `PlasmaEquilibrium.ingest` and reconstructed by
 `build_inverse_from_ingest`. See [`DirectIngest`](@ref) for the role this plays in
 the `gpec.h5` rerun snapshot.
@@ -503,8 +503,8 @@ A container struct for inputs to the `inverse_run` function.
   - `ingest::EquilibriumIngest` - captured raw arrays for the `gpec.h5` rerun snapshot
     (an [`InverseIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria)
   - `psihigh_resolved::Float64` - outer flux limit the equilibrium is formed on; see
-    [`DirectRunInput`](@ref). Equals `config.psihigh` for every inverse reader (CHEASE,
-    analytic); `efit_by_inversion` forwards the clamped value from its `DirectRunInput`.
+    [`DirectRunInput`](@ref). `config.psihigh`, clamped to the outermost input surface for i-files and
+    reruns (`resolve_psihigh!`); `efit_by_inversion` forwards the clamped value from its `DirectRunInput`.
 """
 mutable struct InverseRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND}
     config::EquilibriumConfig

src/Equilibrium/InverseEquilibrium.jl

@@ -41,6 +41,21 @@ end
 
 
 
+"""
+    resolve_psihigh!(input::InverseRunInput) -> InverseRunInput
+
+Clamp `input.psihigh_resolved` to the outermost input surface (below ψ_N = 1 for an i-file that
+carries the LCFS bounds), warning when it moves. `config.psihigh` keeps the request.
+"""
+function resolve_psihigh!(input::InverseRunInput)
+    psi_edge = input.rz_in_xs[end]
+    if input.psihigh_resolved > psi_edge
+        @warn "psihigh = $(input.psihigh_resolved) exceeds the outermost input surface ψ_N = $psi_edge; clamped to it"
+        input.psihigh_resolved = psi_edge
+    end
+    return input
+end
+
 function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing)
     @info "Starting Inverse Equilibrium Processing"
     # Extract input parameters
@@ -211,6 +226,10 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
     spl_fs = zeros(Float64, mtheta+1, 5)
     spl_post_buf = Vector{Float64}(undef, 5)   # scratch for post-SFL spline evaluation
 
+    # q: i-files use the field-line integral, as Fortran `inverse_run`; a tabulated q splined in ψ is
+    # inaccurate where edge surfaces are coarse (q ~ log(1-ψ)). Other inputs keep their input q (the
+    # field-line q is inaccurate in cylindrical LAR geometry).
+    q_from_geometry = config.eq_type == "ifile"
     f_sq_in_buf = Vector{Float64}(undef, size(sq_in.y, 2))
     sq_in_hint = Ref(1)
     hint2d = (Ref(1), Ref(1))
@@ -282,11 +301,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No
         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
+        sq_fs[ipsi+1, 4] = q_from_geometry ? spl_fsi[mtheta+1, 4] * sq_fs[ipsi+1, 1] / (2 * twopi * psio) : f_sq_in_buf[3]
     end
 
     sq = cubic_interp(sq_xs, Series(sq_fs); extrap=ExtendExtrap())

src/Equilibrium/ReadEquilibrium.jl

@@ -395,6 +395,84 @@ function read_chease_ascii(config::EquilibriumConfig)
 end
 
 
+"""
+    _read_fortran_reals(io, n) -> Vector{Float64}
+
+Read one little-endian Fortran sequential unformatted record of `n` reals, stored as real*4 or real*8
+(told apart by the record length).
+"""
+function _read_fortran_reals(io::IO, n::Int)
+    nbytes = ltoh(read(io, Int32))
+    T = nbytes == 8n ? Float64 : nbytes == 4n ? Float32 : error("Fortran record holds $nbytes bytes, expected $n reals")
+    data = ltoh.(read!(io, Vector{T}(undef, n)))
+    ltoh(read(io, Int32)) == nbytes || error("Fortran record end marker does not match its start")
+    return Float64.(data)
+end
+
+"""
+    read_ifile(config) -> InverseRunInput
+
+Read an inverse i-file, such as TokaMaker `save_ifile` writes: little-endian Fortran unformatted records,
+reals as real*8 or real*4, axis first, θ fastest.
+
+  - `mx, my` (int32); `psi(mx)` [Wb/rad], `f(mx)` = R·Bt [T·m], `p(mx)` [Pa], `q(mx)`
+  - `r(my, mx)`, `z(my, mx)` [m], with the periodic θ point duplicated
+  - optional `FF′(mx)`, `p′(mx)` per Wb/rad, skipped
+  - optional `[ψ_axis, ψ_LCFS]`: ψ_N is then normalized to this LCFS, so the last file surface (written just
+    inside it) sits below ψ_N = 1
+
+Port of Fortran `read_eq_ldp_i`.
+"""
+function read_ifile(config::EquilibriumConfig)
+    @info "Reading inverse i-file: $(config.eq_filename)"
+
+    psi, f, p, q, r, z, bounds = open(config.eq_filename, "r") do io
+        ltoh(read(io, Int32)) == 8 || error("i-file header record is not two int32")
+        mx, my = Int.(ltoh.(read!(io, Vector{Int32}(undef, 2))))
+        read(io, Int32)
+        profiles = [_read_fortran_reals(io, mx) for _ in 1:4]
+        r, z = (reshape(_read_fortran_reals(io, my * mx), my, mx) for _ in 1:2)
+        # Optional records, told apart by length: FF′ and p′ (mx each, unused), then the ψ bounds (2)
+        bounds = nothing
+        while !eof(io)
+            nbytes = ltoh(read(io, Int32))
+            if nbytes in (4mx, 8mx)
+                skip(io, nbytes + 4)
+            else
+                skip(io, -4)
+                bounds = _read_fortran_reals(io, 2)
+            end
+        end
+        return (profiles..., r, z, bounds)
+    end
+    mx, my = length(psi), size(r, 1)
+    @info "Parsed from header: mx = $mx surfaces, my = $my poloidal points"
+
+    psi_axis, psi_lcfs = bounds === nothing ? (psi[1], psi[end]) : bounds
+    psio_signed = psi_lcfs - psi_axis
+    psio = abs(psio_signed)
+    xs = (psi .- psi_axis) ./ psio_signed
+    bounds === nothing ||
+        @info "i-file ψ bounds: ψ_axis = $psi_axis, ψ_LCFS = $psi_lcfs; file surfaces span ψ_N ∈ [$(xs[1]), $(xs[end])]"
+
+    sq_fs = hcat(abs.(f), p .* mu0, q, sqrt.(xs))
+    sq_in = cubic_interp(xs, Series(sq_fs); extrap=ExtendExtrap())
+
+    R_data = Matrix(transpose(r))
+    Z_data = Matrix(transpose(z))
+    ro, zo = R_data[1, 1], Z_data[1, 1]
+    rz_in_ys = collect(range(0, 1; length=my))
+
+    opts2d = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
+    rz_in_R = cubic_interp((xs, rz_in_ys), R_data; opts2d...)
+    rz_in_Z = cubic_interp((xs, rz_in_ys), Z_data; opts2d...)
+    ingest = InverseIngest(xs, sq_fs, xs, rz_in_ys, R_data, Z_data, ro, zo, psio)
+
+    @info "Finished reading i-file. Magnetic axis at (ro=$(@sprintf("%.3f", ro)), zo=$(@sprintf("%.3f", zo))), psio=$(@sprintf("%.3e", psio))"
+    return InverseRunInput(config, sq_in, xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, ingest)
+end
+
+
 """
     build_direct_from_ingest(config::EquilibriumConfig, ingest::DirectIngest) -> DirectRunInput
 

test/runtests_equil.jl

@@ -1,3 +1,4 @@
+using Statistics
 
 @testset "Equilibrium Unit Tests" begin
 
@@ -192,6 +193,88 @@
         @test finite(eq_r)
     end
 
+    @testset "Load TokaMaker i-file" begin
+        Eq = GeneralizedPerturbedEquilibrium.Equilibrium
+        # One TokaMaker solve written by save_eqdsk (g65) and save_ifile (i33x65, real*8 and real*4)
+        tk_dir = joinpath(@__DIR__, "test_data", "TokaMaker_ifile")
+        cfg(type, file; kw...) = Eq.EquilibriumConfig(; eq_type=type, eq_filename=joinpath(tk_dir, file), jac_type="hamada",
+            grid_type="core_edge_packed", mpsi=64, psilow=0.01, psihigh=0.99, kw...)
+        eq_i = Eq.setup_equilibrium(cfg("ifile", "i33x65.ifile"))
+        eq_g = Eq.setup_equilibrium(cfg("efit", "g65.geqdsk"))
+        eq_s = Eq.setup_equilibrium(cfg("ifile", "i33x65_single.ifile"))
+
+        # Same equilibrium as the g-file of the same solve
+        @test isapprox(eq_i.ro, eq_g.ro; atol=1e-4) && isapprox(eq_i.zo, eq_g.zo; atol=1e-4)
+        @test isapprox(eq_i.psio, eq_g.psio; rtol=1e-6)
+        @test all(isapprox(eq_i.profiles.q_spline(x), eq_g.profiles.q_spline(x); rtol=1e-3) for x in (0.2, 0.5, 0.8))
+
+        # Grad-Shafranov residual on interior surfaces
+        g = Eq.equilibrium_gse!(eq_i)
+        @test median(vec(maximum(g.error; dims=2))[0.05 .< g.xs .< 0.7]) < 5e-4
+
+        # A real*4 file reads to single-precision agreement
+        @test isapprox(eq_s.psio, eq_i.psio; rtol=1e-6)
+        @test all(isapprox(eq_s.profiles.q_spline(x), eq_i.profiles.q_spline(x); rtol=1e-6) for x in (0.2, 0.5, 0.8))
+        @test maximum(abs.(eq_s.ingest.R_nodes .- eq_i.ingest.R_nodes)) < 1e-6
+
+        # The trailing FF′ and p′ records are ignored: a file without them reads the same
+        table(file) = Eq.read_ifile(cfg("ifile", file)).ingest.sq_fs
+        tabulated = table("i33x65.ifile")
+        raw = read(joinpath(tk_dir, "i33x65.ifile"))
+        mx = Int(reinterpret(Int32, raw[5:8])[1])
+        stripped = joinpath(mktempdir(), "no_derivatives.ifile")
+        write(stripped, raw[1:end-2*(8mx+8)])
+        @test table(stripped) == tabulated
+
+        # q is the field-line q of the file's geometry (Fortran inverse_run), not its tabulated q record
+        qrec = 16 + 3 * (8mx + 8) + 4
+        bumped = copy(raw)
+        bumped[(qrec+1):(qrec+8mx)] = reinterpret(UInt8, 1.1 .* reinterpret(Float64, raw[(qrec+1):(qrec+8mx)]))
+        bumped_file = joinpath(mktempdir(), "bumped_q.ifile")
+        write(bumped_file, bumped)
+        @test Eq.read_ifile(cfg("ifile", bumped_file)).ingest.sq_fs[:, 3] ≈ 1.1 .* tabulated[:, 3]
+        @test Eq.setup_equilibrium(cfg("ifile", bumped_file)).profiles.q_spline.y == eq_i.profiles.q_spline.y
+    end
+
+    @testset "TokaMaker i-file ψ bounds record" begin
+        Eq = GeneralizedPerturbedEquilibrium.Equilibrium
+        tk_dir = joinpath(@__DIR__, "test_data", "TokaMaker_ifile")
+        cfg(file; kw...) = Eq.EquilibriumConfig(; eq_type="ifile", eq_filename=file, jac_type="hamada", grid_type="core_edge_packed", mpsi=32,
+            psilow=0.01, psihigh=0.99, kw...)
+        raw = read(joinpath(tk_dir, "i33x65.ifile"))  # carries the FF′ and p′ records
+        mx = Int(reinterpret(Int32, raw[5:8])[1])
+        psi = reinterpret(Float64, raw[21:(20+8mx)])
+        stripped = raw[1:(end-2*(8mx+8))]
+        record(v) = vcat(reinterpret(UInt8, [Int32(8length(v))]), reinterpret(UInt8, Float64.(v)), reinterpret(UInt8, [Int32(8length(v))]))
+        dir = mktempdir()
+        variant(name, bytes) = (path = joinpath(dir, name); write(path, bytes); path)
+        ingest(path) = Eq.read_ifile(cfg(path)).ingest
+
+        base_s = ingest(variant("s.ifile", stripped))
+
+        # The true LCFS beyond the padded last surface: ψ_N shrinks by (1 - pad), psio grows by 1/(1 - pad).
+        # The record is read both after r, z (sb) and after FF′, p′ (db).
+        pad = 0.01
+        bounds = [psi[1], psi[1] + (psi[end] - psi[1]) / (1 - pad)]
+        sb = variant("sb.ifile", vcat(stripped, record(bounds)))
+        db = variant("db.ifile", vcat(raw, record(bounds)))
+        ing_sb = @test_logs (:info, r"ψ bounds") match_mode = :any ingest(sb)
+        ing_db = ingest(db)
+        for ing in (ing_sb, ing_db)
+            @test ing.sq_xs ≈ base_s.sq_xs .* (1 - pad)
+            @test ing.psio ≈ base_s.psio / (1 - pad)
+        end
+        @test ing_sb.sq_fs[:, 1:3] == ing_db.sq_fs[:, 1:3] == base_s.sq_fs[:, 1:3]
+
+        # psihigh beyond the last surface is clamped to it; a rerun from the ingest keeps the normalisation
+        eq = @test_logs (:warn, r"clamped") match_mode = :any Eq.setup_equilibrium(cfg(db; psihigh=0.995))
+        @test eq.params.psihigh_resolved == ing_db.sq_xs[end] ≈ eq.rzphi_xs[end]
+        @test all(isfinite, eq.profiles.q_spline.y) && all(isfinite, eq.eqfun_B.nodal_derivs.partials)
+        rerun = Eq.setup_equilibrium(cfg(db; psihigh=0.995), Eq.build_inverse_from_ingest(cfg(db), eq.ingest))
+        @test rerun.rzphi_xs == eq.rzphi_xs && rerun.psio == eq.psio
+        @test rerun.profiles.q_spline.y == eq.profiles.q_spline.y
+    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.

test/test_data/README_test_data.md

@@ -1,3 +1,8 @@
 Put relevant data for various tests here, most likely outputs from the Fortran code used to validate Julia outputs
 
 # TODO: store as hdf5 files instead?
+
+## TokaMaker_ifile
+One TokaMaker (OpenFUSIONToolkit) solve of a DIII-D-like H-mode, written as `g65.geqdsk` (`save_eqdsk`, 65×65)
+and as `i33x65.ifile` / `i33x65_single.ifile` (`save_ifile`, 33 surfaces × 65 angles, real*8 / real*4, with the
+trailing FF′ and p′ records). Used by the i-file reader and ψ bounds record tests in `runtests_equil.jl`.