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

docs/src/equilibrium.md

@@ -55,10 +55,11 @@ 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
+## F and P from a g-file, i-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`).
+`PPRIME`); an IMAS equilibrium does the same (`f`, `pressure`, `f_df_dpsi`, `dpressure_dpsi`), and so does an
+i-file (`ldp_i`) whose optional trailing FF′ and p′ records are present (TokaMaker `save_ifile` writes them).
 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

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
+  - `ldp_i` / `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

@@ -106,6 +106,8 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin
         eq_input = read_chease_ascii(eq_config)
     elseif eq_type in ["chease", "chease_binary"]
         eq_input = read_chease_binary(eq_config)
+    elseif eq_type in ["ldp_i", "ifile"]
+        eq_input = read_ldp_i(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
@@ -487,6 +489,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, their sum (the residual), the residual normalized by the
+largest term, and the θ-integrated residual per surface.
 """
 function equilibrium_gse!(equil::PlasmaEquilibrium)
 
@@ -645,6 +651,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

@@ -42,7 +42,7 @@ specified in the input.
   - `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.
+  - `profile_source::String` - Which 1D arrays of a g-file, i-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.
@@ -380,8 +380,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_ldp_i`) — 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.

src/Equilibrium/ReadEquilibrium.jl

@@ -481,6 +481,69 @@ 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_ldp_i(config)
+
+Parses an inverse i-file (L. Don Pearlstein's format, also written by TokaMaker `save_ifile`) into an
+`InverseRunInput`. Sequential unformatted little-endian 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, then optionally `FF′(mx)` and
+`p′(mx)` per Wb/rad, which `profile_source = "derivatives"` integrates for F and P (see
+`file_profiles`). Port of Fortran `read_eq_ldp_i`.
+"""
+function read_ldp_i(config::EquilibriumConfig)
+    @info "Reading inverse i-file (ldp_i): $(config.eq_filename)"
+
+    psi, f, p, q, r, z, ffp, pp = 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)
+        derivs = eof(io) ? (zeros(mx), zeros(mx)) : (_read_fortran_reals(io, mx), _read_fortran_reals(io, mx))
+        return (profiles..., r, z, derivs...)
+    end
+    mx, my = length(psi), size(r, 1)
+    @info "Parsed from header: mx = $mx surfaces, my = $my poloidal points"
+
+    psio_signed = psi[end] - psi[1]
+    psio = abs(psio_signed)
+    xs = (psi .- psi[1]) ./ psio_signed
+
+    # FF′ and p′ are per Wb/rad; psio_signed converts them to ψ_norm (absent records read as zero)
+    f_nodes, p_nodes = file_profiles(config, xs, f, p, ffp .* psio_signed, pp .* psio_signed)
+    sq_fs = hcat(f_nodes, p_nodes .* 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,5 @@
+using Statistics
+
 
 @testset "Equilibrium Unit Tests" begin
 
@@ -197,6 +199,46 @@
         @test integrated[:, 3:4] == tabulated[:, 3:4]
     end
 
+    @testset "Load TokaMaker i-file (ldp_i)" 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="ldp", mpsi=64, psilow=0.01, psihigh=0.99, kw...)
+        eq_i = Eq.setup_equilibrium(cfg("ldp_i", "i33x65.ifile"))
+        eq_g = Eq.setup_equilibrium(cfg("efit", "g65.geqdsk"))
+        eq_s = Eq.setup_equilibrium(cfg("ifile", "i33x65_single.ifile"))
+        @test eq_i.ingest isa Eq.InverseIngest
+
+        # 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: small, and well below the g-file's
+        gse_median(eq) = (g = Eq.equilibrium_gse!(eq); loc = vec(maximum(g.error; dims=2)); median(loc[0.05 .< g.xs .< 0.7]))
+        @test gse_median(eq_i) < 5e-4
+        @test gse_median(eq_i) < gse_median(eq_g) / 5
+
+        # A real*4 file reads to single-precision agreement
+        @test isapprox(eq_s.ro, eq_i.ro; rtol=1e-6) && 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 FF′ and p′ records feed profile_source = "derivatives"; without them the tabulated F and P are used
+        table(file; kw...) = Eq.read_ldp_i(cfg("ldp_i", file; kw...)).ingest.sq_fs
+        integrated = @test_logs (:info, r"integrated from") match_mode = :any table("i33x65.ifile")
+        tabulated = table("i33x65.ifile"; profile_source="values")
+        @test integrated[end, 1:2] == tabulated[end, 1:2]
+        @test integrated[:, 1] != tabulated[:, 1]
+        @test maximum(abs.(integrated[:, 1] .- tabulated[:, 1]) ./ tabulated[:, 1]) < 1e-4
+        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 (@test_logs (:warn, r"absent or unusable") match_mode = :any Eq.read_ldp_i(cfg("ldp_i", stripped)).ingest.sq_fs) == tabulated
+    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 `ldp_i` reader test in `runtests_equil.jl`.