Add D3D-like kinetic profiles and update example equilibrium

.gitignore

@@ -8,6 +8,8 @@
 *.i
 *.out
 *.h5
+# kinetic-profile inputs are tracked despite the blanket *.h5 (output) ignore
+!examples/**/*kinetic*.h5
 .vscode
 *.DS_Store
 *.pdf

docs/src/kinetic_forces.md

@@ -4,6 +4,27 @@ Kinetic torque and energy calculations for perturbed equilibria.
 Implements neoclassical toroidal viscosity (NTV) from the PENTRC formulation
 ([Logan & Park, 2013](citations.md); [Logan, 2015](citations.md)).
 
+## Kinetic profile file formats
+
+Kinetic profiles are read by `read_kinetic_file`, which dispatches on the file
+extension:
+
+  - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets
+    at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`,
+    `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional
+    `omega_tor`, `chi_e` (perpendicular heat diffusivity ``\chi_\perp``), and
+    `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``). Each dataset
+    carries a `units` attribute; the root carries `schema_version` and
+    `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
+    diffusivities m²/s. Write files with `write_kinetic_h5`.
+  - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table
+    `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
+    compatibility. Header rows are skipped.
+
+The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`,
+when present, are carried for the resistive-layer (SLAYER) analysis and ignored
+here.
+
 ## Profile Scaling Knobs
 
 Seven scaling factors are available on `KineticForcesControl` to modify kinetic

examples/DIIID-like_ideal_example/gpec.toml

@@ -7,11 +7,11 @@
 eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
 eq_type = "efit"               # Type of the input 2D equilibrium file
 jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc, park, custom)
-grid_type = "log_asymptotic"   # Radial grid packing type
+grid_type = "ldp"              # Radial grid packing (sin²: packs core+edge). Robust δW/Δ' vs log_asymptotic.
 psilow = 1e-4                  # Lower limit of normalized poloidal flux
-psihigh = 0.9995               # Upper limit of normalized poloidal flux
-mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
-psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+psihigh = 0.995                # captures q=6 with sas_flag/dmlim at this Ip
+mpsi = 256                     # Fixed radial grid (avoid mpsi=0 auto: it edge-starves → spurious δW/Δ')
+psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi (unused when mpsi>0)
 mtheta = 256                   # Number of poloidal grid points
 newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
 etol = 1e-10                   # Error tolerance for equilibrium solver
@@ -76,3 +76,24 @@ compute_singular_coupling = true        # Compute singular layer coupling metric
 verbose = true                          # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
 reg_spot = 0.05                         # Regularization width for singular surfaces (0 = disabled)
+
+[KineticForces]
+kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5"  # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_phi). Legacy .gpeckf/.kin still readable.
+nn = 1                                  # Toroidal mode number (matches nn_low/nn_high)
+nl = 6                                  # Bounce harmonics retained
+zi = 1                                  # Main-ion charge (deuterium)
+mi = 2                                  # Main-ion mass [proton masses]
+zimp = 6                                # Impurity charge (carbon)
+mimp = 12                               # Impurity mass [proton masses]
+electron = false                        # Ion-only NTV
+nutype = "harmonic"                     # Collision operator
+f0type = "maxwellian"                   # Equilibrium distribution
+moment = "pressure"                     # Pressure-moment NTV torque
+atol_xlmda = 1e-9                       # Inner pitch/energy integration tolerances
+rtol_xlmda = 1e-5
+atol_psi = 1e-2                         # Outer psi quadrature tolerances (relaxed from 1e-3:
+rtol_psi = 1e-2                         # NTV torque integrand has sharp resonant peaks at rational
+                                        # surfaces; tight tol over-refines. Proper fix = panel the
+                                        # psi-integral at rationals (deferred follow-up).
+write_outputs_to_HDF5 = true
+verbose = true

regression-harness/cases/diiid_n1.toml

@@ -347,6 +347,19 @@ label = "PE toroidal torque"
 noise_threshold = 1e-12
 order = 94
 
+# NTV (neoclassical toroidal viscosity) torque from the [KineticForces] block —
+# FGAR method, pressure moment. Diagnostic computed on top of the ideal FFS + PE
+# result. Complex total torque stored as [Re, Im]; tracked as a 2-element vector.
+# Absent (N/A) for refs without a [KineticForces] section. Threshold is a noise
+# floor for the adaptive-quadrature result; tune if run-to-run jitter appears.
+[quantities.ntv_torque_fgar]
+h5path = "kinetic_forces/fgar/total_torque"
+type = "real_vector"
+extract = "all_real"
+label = "NTV torque FGAR (Re,Im) [N·m]"
+noise_threshold = 1e-2
+order = 95
+
 # Profiles (checksums)
 [quantities.q_profile]
 h5path = "splines/profiles/q"

src/Equilibrium/Equilibrium.jl

@@ -26,7 +26,8 @@ include("KineticProfiles.jl")
 # --- Expose types and functions to the user ---
 export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParameters,
     ProfileSplines, GeometryProfileSplines, compute_geometry_profiles,
-    KineticProfileSplines, load_kinetic_profiles
+    KineticProfileSplines, load_kinetic_profiles,
+    KineticProfileData, read_kinetic_file, write_kinetic_h5
 export flux_surface_metric, flux_surface_area
 export compute_sqrt_jac_delpsi, compute_sqrtamat, rootarea_to_area_weight, area_to_rootarea_weight
 

src/Equilibrium/KineticProfiles.jl

@@ -1,13 +1,132 @@
 """
     KineticProfiles
 
-Reading and processing kinetic profile data (density, temperature, rotation)
-from ASCII tables. The output is a `KineticProfileSplines` of independent named
-splines that downstream physics modules (KineticForces NTV, etc.) can read
-directly without re-implementing data shimming.
+Reading and processing kinetic profile data (density, temperature, rotation,
+and optionally transport diffusivities) from the GPEC HDF5 kinetic format or
+from legacy 6-column ASCII tables (`.gpeckf`/`.kin`). The processed output is a
+`KineticProfileSplines` of independent named splines that downstream physics
+modules (KineticForces NTV, etc.) read directly without re-implementing data
+shimming.
 """
 
 using DelimitedFiles
+using HDF5
+
+"""
+    KineticProfileData
+
+Raw kinetic-profile columns read from a kinetic file (HDF5 or ASCII), before
+resampling and spline construction. `psi` (normalized poloidal flux) is always
+present; every profile field is `nothing` when the source omits it, so each
+consumer validates only the fields it needs. Units: densities m⁻³,
+temperatures eV, frequencies rad/s, diffusivities m²/s.
+
+| field       | meaning                              |
+|-------------|--------------------------------------|
+| `n_i`       | main-ion density                     |
+| `n_e`       | electron density                     |
+| `T_i`/`T_e` | ion / electron temperature           |
+| `omega_E`   | ExB rotation                         |
+| `omega_tor` | toroidal rotation (optional)         |
+| `chi_e`     | perpendicular heat diffusivity χ⊥    |
+| `chi_phi`   | toroidal momentum diffusivity χ_φ    |
+"""
+struct KineticProfileData
+    psi::Vector{Float64}
+    n_i::Union{Nothing,Vector{Float64}}
+    n_e::Union{Nothing,Vector{Float64}}
+    T_i::Union{Nothing,Vector{Float64}}
+    T_e::Union{Nothing,Vector{Float64}}
+    omega_E::Union{Nothing,Vector{Float64}}
+    omega_tor::Union{Nothing,Vector{Float64}}
+    chi_e::Union{Nothing,Vector{Float64}}
+    chi_phi::Union{Nothing,Vector{Float64}}
+    provenance::String
+end
+
+function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing,
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, provenance="")
+    _v(x) = x === nothing ? nothing : Float64.(collect(x))
+    return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e),
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), String(provenance))
+end
+
+const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
+
+# Dataset units written into / expected from the HDF5 kinetic schema.
+const _KINETIC_H5_UNITS = Dict("psi" => "normalized poloidal flux", "n_i" => "m^-3",
+    "n_e" => "m^-3", "T_i" => "eV", "T_e" => "eV", "omega_E" => "rad/s",
+    "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_phi" => "m^2/s")
+
+"""
+    read_kinetic_file(path; group="/") -> KineticProfileData
+
+Read a kinetic-profile file, dispatching on extension: HDF5 (`.h5`/`.hdf5`) via
+the GPEC kinetic schema, otherwise a 6-column ASCII table (`.gpeckf`/`.kin`/`.dat`:
+`psi_n n_i n_e T_i[eV] T_e[eV] omega_E`). Returns the raw columns; downstream
+consumers validate the fields they require.
+"""
+function read_kinetic_file(path::AbstractString; group::AbstractString="/")
+    isfile(path) || error("Kinetic profile file not found: $path")
+    ext = lowercase(splitext(path)[2])
+    return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path)
+end
+
+"""Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`."""
+function _read_kinetic_ascii(path::AbstractString)
+    psi, ni, ne, Ti, Te, omegaE = _read_kinetic_table(path)
+    return KineticProfileData(; psi=psi, n_i=ni, n_e=ne, T_i=Ti, T_e=Te, omega_E=omegaE,
+        provenance="ASCII 6-column: $(basename(path))")
+end
+
+"""Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are
+optional and surfaced as-is."""
+function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
+    h5open(path, "r") do f
+        g = group == "/" ? f : f[group]
+        haskey(g, "psi") || error("kinetic HDF5 '$path' group '$group': missing required dataset 'psi'")
+        rd(name) = haskey(g, name) ? Float64.(vec(read(g[name]))) : nothing
+        prov = ""
+        ats = attributes(g)
+        haskey(ats, "provenance") && (prov = string(read(ats["provenance"])))
+        return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))),
+            n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"),
+            omega_E=rd("omega_E"), omega_tor=rd("omega_tor"),
+            chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), provenance=prov)
+    end
+end
+
+"""
+    write_kinetic_h5(path, data; group="/", schema_version="1.0", provenance=data.provenance)
+
+Write a `KineticProfileData` to the GPEC HDF5 kinetic schema. Each present field
+is written as a dataset with a `units` attribute; absent (`nothing`) fields are
+skipped. The group carries `schema_version` and `provenance` attributes.
+"""
+function write_kinetic_h5(path::AbstractString, data::KineticProfileData;
+    group::AbstractString="/", schema_version::AbstractString="1.0",
+    provenance::AbstractString=data.provenance)
+    h5open(path, "w") do f
+        g = group == "/" ? f : create_group(f, group)
+        function put(name, v)
+            v === nothing && return
+            g[name] = collect(Float64, v)
+            attributes(g[name])["units"] = _KINETIC_H5_UNITS[name]
+        end
+        put("psi", data.psi)
+        put("n_i", data.n_i)
+        put("n_e", data.n_e)
+        put("T_i", data.T_i)
+        put("T_e", data.T_e)
+        put("omega_E", data.omega_E)
+        put("omega_tor", data.omega_tor)
+        put("chi_e", data.chi_e)
+        put("chi_phi", data.chi_phi)
+        attributes(g)["schema_version"] = schema_version
+        attributes(g)["provenance"] = provenance
+    end
+    return path
+end
 
 """
     load_kinetic_profiles(kinetic_file::AbstractString;
@@ -17,20 +136,21 @@ using DelimitedFiles
                           chi1::Union{Nothing,Float64}=nothing)
         → KineticProfileSplines
 
-Parse an ASCII kinetic profile file, interpolate onto a regular 101-point ψ
-grid, optionally apply profile scaling knobs, derive collisional / Z_eff
-diagnostics, and return a `KineticProfileSplines` with independent named cubic
-splines.
+Parse a kinetic profile file (ASCII or HDF5, dispatched by `read_kinetic_file`),
+interpolate onto a regular 101-point ψ grid, optionally apply profile scaling
+knobs, derive collisional / Z_eff diagnostics, and return a
+`KineticProfileSplines` with independent named cubic splines.
 
 # Expected file format
 
-Six whitespace-separated columns (header rows are filtered out):
+An HDF5 file following the GPEC kinetic schema (fields read by name), or a legacy
+six-column whitespace-separated ASCII table (header rows are filtered out):
 
     psi_n  n_i[m^-3]  n_e[m^-3]  T_i[eV]  T_e[eV]  omega_E[rad/s]
 
 # Arguments
 
-  - `kinetic_file`: Path to the ASCII kinetic profile file
+  - `kinetic_file`: Path to the ASCII or HDF5 kinetic profile file (see `read_kinetic_file`)
   - `zi`, `zimp`: Main ion and impurity charge numbers
   - `mi`, `mimp`: Main ion and impurity mass numbers (in proton masses)
   - `density_factor`: Density scaling factor (applied to ni, ne)
@@ -60,12 +180,17 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
     ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0,
     chi1::Union{Nothing,Float64}=nothing)
 
-    if !isfile(kinetic_file)
-        error("Kinetic profile file not found: $kinetic_file")
-    end
-
-    psi_input, ni_input, ne_input, Ti_input_eV, Te_input_eV, omegaE_input =
-        _read_kinetic_table(kinetic_file)
+    data = read_kinetic_file(kinetic_file)
+
+    # NTV requires n_e, T_i, T_e, omega_E; n_i defaults to n_e (quasineutrality)
+    # when absent. omega_tor / chi_e / chi_phi (if present) are ignored here.
+    _need(field, name) = field === nothing ? error("kinetic file '$kinetic_file' missing required '$name'") : field
+    psi_input = data.psi
+    ne_input = _need(data.n_e, "n_e")
+    Ti_input_eV = _need(data.T_i, "T_i")
+    Te_input_eV = _need(data.T_e, "T_e")
+    omegaE_input = _need(data.omega_E, "omega_E")
+    ni_input = data.n_i === nothing ? copy(ne_input) : data.n_i
 
     eV_to_J = 1.602e-19
     mp = 1.672_614e-27

src/KineticForces/Torque.jl

@@ -37,12 +37,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
               imx_override::Union{Nothing,Float64}=nothing,
               atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6)
 
-    if intr.verbose
-        println("torque - tpsi function, psi = ", psi)
-        println("  electron ", electron)
-        println("  ell ", l)
-    end
-
     # Enforce bounds
     if psi > 1
         tpsi_var[] = 0.0
@@ -196,15 +190,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
     wdhat = q^3 * wtran^2 / (4 * epsr * wgyro) * wdfac
     nueff = nu_s / (2 * epsr)
 
-    if intr.verbose
-        @printf("   eq values = %.1e %.1e %.1e %.1e %.1e %.1e %.1f %d\n",
-               wdian, wdiat, welec, wdhat, wbhat, nueff, q, 0)
-    end
-
-    if intr.verbose
-        println("  method = ", method)
-    end
-
     # Method selection — route on the registry dispatch tag (errors on unknown method)
     kind = method_kind(method)
     if kind == :fcgl
@@ -248,10 +233,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
                                    rex_override=rex_override, imx_override=imx_override)
     end
 
-    if intr.verbose
-        println("torque - end function, psi = ", psi)
-    end
-
     return nothing
 end
 

test/runtests_kinetic.jl

@@ -371,4 +371,59 @@
         @test mr.total_energy == 0.0 + 0.0im
         @test isempty(mr.records)
     end
+
+    # =========================================================================
+    # Kinetic profile I/O: HDF5 schema, ASCII back-compat, writer round-trip
+    # =========================================================================
+    @testset "kinetic profile I/O" begin
+        E = GeneralizedPerturbedEquilibrium.Equilibrium
+        exdir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example")
+        gpeckf = joinpath(exdir, "TkMkr_D3Dlike_Hmode_kinetic.gpeckf")
+        h5 = joinpath(exdir, "TkMkr_D3Dlike_Hmode_kinetic.h5")
+
+        @testset "extension dispatch + raw read" begin
+            da = E.read_kinetic_file(gpeckf)
+            db = E.read_kinetic_file(h5)
+            # legacy ASCII carries no optional fields; the example .h5 carries chi
+            @test da.chi_e === nothing && da.chi_phi === nothing
+            @test db.chi_e !== nothing && db.chi_phi !== nothing
+            # core columns identical (the .h5 was generated from the .gpeckf)
+            @test da.psi == db.psi
+            @test da.n_e == db.n_e
+            @test da.T_i == db.T_i
+            @test da.omega_E == db.omega_E
+        end
+
+        @testset "NTV splines identical: ASCII vs HDF5" begin
+            ka = E.load_kinetic_profiles(gpeckf)
+            kb = E.load_kinetic_profiles(h5)
+            for p in 0.0:0.1:1.0
+                @test ka.ne_spline(p) == kb.ne_spline(p)
+                @test ka.Ti_spline(p) == kb.Ti_spline(p)
+                @test ka.omegaE_spline(p) == kb.omegaE_spline(p)
+                @test ka.zeff_spline(p) == kb.zeff_spline(p)
+            end
+        end
+
+        @testset "writer round-trip" begin
+            db = E.read_kinetic_file(h5)
+            tmp = tempname() * ".h5"
+            E.write_kinetic_h5(tmp, db; provenance="roundtrip")
+            dc = E.read_kinetic_file(tmp)
+            @test dc.chi_e ≈ db.chi_e
+            @test dc.chi_phi ≈ db.chi_phi
+            @test dc.n_e ≈ db.n_e
+            @test dc.provenance == "roundtrip"
+            rm(tmp; force=true)
+        end
+
+        @testset "missing required field errors" begin
+            # HDF5 with only psi + n_e: NTV load must error on the missing T_i
+            tmp = tempname() * ".h5"
+            d = E.KineticProfileData(; psi=[0.0, 0.5, 1.0], n_e=[1e19, 9e18, 8e18])
+            E.write_kinetic_h5(tmp, d)
+            @test_throws ErrorException E.load_kinetic_profiles(tmp)
+            rm(tmp; force=true)
+        end
+    end
 end