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
@@ -70,8 +71,8 @@ variation of ``F``). On a 257-node file written in single precision that is an e
 
 | 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 |
+| `"values"` (default) | The tabulated ``F`` and ``P`` are used as written |
+| `"derivatives"` | ``F`` and ``P`` are the integrals of the file's ``FF'`` and ``p'``, anchored at the boundary values of the tabulated ``F`` and ``P`` |
 
 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

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

@@ -38,6 +38,8 @@ const mu0 = 4π * 1e-7
 
 # efit-family equilibrium kinds that share the separatrix-clamp / inversion handling.
 const EFIT_KINDS = ("efit", "efit_arclength", "efit_by_inversion")
+# eq_type names of the inverse i-file reader (`read_ldp_i`).
+const LDP_I_KINDS = ("ldp_i", "ifile")
 
 """
     AnalyticEqSpec
@@ -100,12 +102,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 in LDP_I_KINDS
+        eq_input = resolve_psihigh!(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 +492,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)
 
@@ -645,6 +654,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,10 +42,10 @@ 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.
-    `"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.
+  - `profile_source::String` - Which 1D arrays of a g-file, i-file or IMAS equilibrium define F and P.
+    `"values"` (default) uses the tabulated F and P directly; `"derivatives"` integrates the file's
+    FF′ and p′ inward from the boundary values of F and P, falling back to `"values"` with a
+    warning when the derivatives are absent or unusable.
 """
 @kwdef struct EquilibriumConfig
     eq_type::String = "efit"
@@ -78,7 +78,7 @@ specified in the input.
     use_galgrid::Bool = true
 
     imas_cocos::Int = 11
-    profile_source::String = "derivatives"
+    profile_source::String = "values"
 
     """
     Modified internal constructor that enforces self consistency within the inputs
@@ -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.
@@ -509,8 +509,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
@@ -209,6 +224,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 in LDP_I_KINDS
     f_sq_in_buf = Vector{Float64}(undef, size(sq_in.y, 2))
     sq_in_hint = Ref(1)
     hint2d = (Ref(1), Ref(1))
@@ -280,11 +299,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

@@ -48,9 +48,9 @@ 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`:
+Rebuild F and P on the nodes `xs` of a normalized flux coordinate (0 at the axis) by integrating the
+tabulated `ffprime` = d(F²/2)/dx and `pprime` = dP/dx inward from the values of `f` and `p` at the
+last node (written x = 1 below):
 
     F²(x) = F(1)² − 2∫ₓ¹ FF′ dx,    P(x) = P(1) − ∫ₓ¹ P′ dx.
 
@@ -101,7 +101,7 @@ 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
+F (as a magnitude) and P on the 1D nodes of a file-based 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.
@@ -481,6 +481,86 @@ 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) -> InverseRunInput
+
+Read an inverse i-file (L. Don Pearlstein's format, also written by TokaMaker `save_ifile`): 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, used by `profile_source = "derivatives"` (see `file_profiles`)
+  - 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_ldp_i(config::EquilibriumConfig)
+    @info "Reading inverse i-file (ldp_i): $(config.eq_filename)"
+
+    psi, f, p, q, r, z, ffp, pp, 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), then the ψ bounds (2)
+        ffp, pp, bounds = zeros(mx), zeros(mx), nothing
+        while !eof(io)
+            nbytes = ltoh(read(io, Int32))
+            skip(io, -4)
+            if nbytes in (4mx, 8mx)
+                ffp, pp = _read_fortran_reals(io, mx), _read_fortran_reals(io, mx)
+            else
+                bounds = _read_fortran_reals(io, 2)
+            end
+        end
+        return (profiles..., r, z, ffp, pp, 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])]"
+
+    # 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