From a5ab5882e5be3d1cbf1902e59687a469cbd3514f Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 19:29:36 +0200 Subject: [PATCH 1/7] Equilibrium - API! - Read the kinetic-profile file once into the equilibrium --- .../benchmark_diiid_ideal_ntv_torque.jl | 4 +- docs/src/api.md | 12 ++ docs/src/kinetic_forces.md | 4 +- src/Equilibrium/Equilibrium.jl | 8 +- src/Equilibrium/EquilibriumTypes.jl | 68 ++++++++- src/Equilibrium/KineticProfiles.jl | 130 +++++++----------- src/GeneralizedPerturbedEquilibrium.jl | 10 +- .../CalculatedKineticMatrices.jl | 2 +- test/runtests_kinetic.jl | 34 ++++- test/runtests_multiion.jl | 30 ++-- 10 files changed, 183 insertions(+), 119 deletions(-) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 430c353cc..d9de8fa81 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -307,8 +307,8 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) _p("\n--- Torque integration ---") t2 = time() - kinetic_profiles = Eq.load_kinetic_profiles( - kf_ctrl.kinetic_file; + kinetic_profiles = Eq.kinetic_profile_splines( + Eq.read_kinetic_file(kf_ctrl.kinetic_file); zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass) kf_state = KF.KineticForcesState() diff --git a/docs/src/api.md b/docs/src/api.md index bf734f762..eb81bee19 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -38,6 +38,18 @@ The four objects map one-to-one onto the pipeline stages: - `perturbed_equilibrium` runs the plasma-response stage — the `[PerturbedEquilibrium]` section. +## Kinetic profiles + +The equilibrium carries the kinetic profiles. Name the file with the `kinetic_file` (and, +for HDF5, `kinetic_group`) keyword, the same keys as the `[Equilibrium]` deck section, and it +is read once into `equil.kinetic`; every kinetic stage derives its own view from it. +[`attach_kinetic_profiles!`](@ref) attaches a file or in-memory `KineticProfileData` to an +equilibrium built without one, such as an analytic equilibrium. + +```julia +eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada", kinetic_file="kinetic.h5") +``` + ## Combining forcing sources `RMPField`s form a vector space: `+`, `-` and scalar `*` build lazy linear combinations that diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index ad866cf06..b2e9fff1c 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -95,8 +95,8 @@ profiles and physics parameters for sensitivity studies: ### Profile-loader knobs (`density_factor`, `temperature_factor`, `ExB_rotation_factor`, `toroidal_rotation_factor`) -These four knobs are applied during kinetic profile loading in -`load_kinetic_profiles`, before any physics evaluation. The physical model is: +These four knobs are applied when the NTV view of the kinetic profiles is built by +`kinetic_profile_splines`, before any physics evaluation. The physical model is: ```math \omega_\phi = \omega_E + \omega_{*n,i} + \omega_{*T,i} diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 51c92eaf9..f36f02069 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -30,8 +30,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, shift_exb_rotation, - KineticProfileData, read_kinetic_file, write_kinetic_h5 + KineticProfileSplines, kinetic_profile_splines, shift_exb_rotation, + KineticProfileData, read_kinetic_file, write_kinetic_h5, attach_kinetic_profiles! export flux_surface_metric, flux_surface_area export wants_two_pass, refined_psi_grid, merge_mandatory_nodes, bracket_mandatory_nodes, enforce_min_spacing, implied_knot_count export compute_sqrt_jac_delpsi, compute_sqrtamat, rootarea_to_area_weight, area_to_rootarea_weight @@ -86,7 +86,8 @@ callers with `mpsi=0` get the coarse first-pass grid unless they refine and re-f function setup_equilibrium(path::String="equil.toml") return setup_equilibrium(EquilibriumConfig(path)) end -function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothing; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing) +function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothing; override_psi_nodes::Union{Nothing,Vector{Float64}}=nothing, + kinetic::Union{Nothing,KineticProfileData}=nothing) eq_type = eq_config.eq_type @@ -142,6 +143,7 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin # Forward the captured ingest so the gpec.h5 writer can snapshot it (nothing for analytic). plasma_equilibrium.ingest = eq_input.ingest + plasma_equilibrium.kinetic = kinetic === nothing ? configured_kinetic_profiles(eq_config) : kinetic equilibrium_global_parameters!(plasma_equilibrium) equilibrium_qfind!(plasma_equilibrium) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index c195b3b54..e75dca716 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -42,6 +42,10 @@ specified in the input. - `flux_surface_rtol::Float64` - Relative tolerance of the ODE solves that trace each flux surface; their absolute tolerance is `equil_abstol(flux_surface_rtol)`, i.e. `flux_surface_rtol` capped at `EQUIL_ABSTOL_MAX` - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations) + - `imas_cocos::Int` - Expected COCOS convention of an IMAS `dd.equilibrium` input (11 = IMAS standard, 2 = GPEC internal) + - `kinetic_file::String` - Kinetic-profile file read once into `PlasmaEquilibrium.kinetic` (see [`read_kinetic_file`](@ref)); + empty for none. Every kinetic consumer (NTV, tearing, the layer-overlap scan, grid refinement) derives its view from it. + - `kinetic_group::String` - HDF5 group of `kinetic_file` holding the profiles """ @kwdef struct EquilibriumConfig eq_type::String = "efit" @@ -75,6 +79,9 @@ specified in the input. # IMAS-specific: expected COCOS convention of the input dd.equilibrium (11=IMAS standard, 2=GPEC internal) imas_cocos::Int = 11 + kinetic_file::String = "" + kinetic_group::String = "/" + """ Modified internal constructor that enforces self consistency within the inputs """ @@ -84,7 +91,7 @@ specified in the input. function EquilibriumConfig(eq_type, eq_filename, r0_norm, b0_norm, jac_type, _, _, _, _, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, flux_surface_rtol, - force_termination, imas_cocos) + force_termination, imas_cocos, kinetic_file, kinetic_group) if jac_type == "hamada" @info "Forcing hamada coordinate jacobian exponents: power_*" power_b = 0 @@ -151,7 +158,7 @@ specified in the input. return new(eq_type, eq_filename, r0_norm, b0_norm, jac_type, power_bp, power_b, power_r, power_rc, jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc, grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, flux_surface_rtol, - force_termination, imas_cocos) + force_termination, imas_cocos, kinetic_file, kinetic_group) end end @@ -185,6 +192,9 @@ function EquilibriumConfig(equil_dict::Dict{String,Any}, base_path::String="./") if haskey(config_data, "eq_filename") && !isempty(config_data["eq_filename"]) && !isabspath(config_data["eq_filename"]) config_data["eq_filename"] = normpath(joinpath(base_path, config_data["eq_filename"])) end + if !isempty(get(config_data, "kinetic_file", "")) && !isabspath(config_data["kinetic_file"]) + config_data["kinetic_file"] = normpath(joinpath(base_path, config_data["kinetic_file"])) + end return EquilibriumConfig(; symbolize_keys(config_data)...) end @@ -768,6 +778,51 @@ function GeometryProfileSplines(xs::Vector{Float64}, ) end +""" + 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 χ⊥ (shipped data: electron power-balance χ_e; see the stability docs) | +| `chi_phi` | toroidal momentum diffusivity χ_φ | +| `zeff` | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure | +| `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | +| `provenance` | short string recording the source file/format | +""" +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}} + zeff::Union{Nothing,Vector{Float64}} + species_densities::Union{Nothing,Dict{String,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, zeff=nothing, + species_densities=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), _v(zeff), species_densities, String(provenance)) +end + """ KineticProfileSplines @@ -894,6 +949,8 @@ This object provides a complete representation of the processed plasma equilibri - `ingest::EquilibriumIngest`: raw arrays forwarded from the equilibrium input for the `gpec.h5` rerun snapshot — a [`DirectIngest`](@ref)/[`InverseIngest`](@ref) for file-based equilibria, or `nothing` for analytic ones (regenerated from their TOML section on replay) + - `kinetic::Union{Nothing,KineticProfileData}`: the raw kinetic profiles of `config.kinetic_file` (or of + [`attach_kinetic_profiles!`](@ref)), on the file's own ψ grid; `nothing` when none are attached """ mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D<:FastInterpolations.CubicInterpolantND} config::EquilibriumConfig @@ -921,13 +978,14 @@ mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D psio::Float64 ingest::EquilibriumIngest + kinetic::Union{Nothing,KineticProfileData} end -# Solvers build the equilibrium before setup_equilibrium forwards eq_input.ingest, so allow -# construction without it; ingest defaults to nothing and is assigned post-construction. +# Solvers build the equilibrium before setup_equilibrium forwards eq_input.ingest and attaches the +# kinetic profiles, so allow construction without them; both are assigned post-construction. PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys, rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac, eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio) = PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys, rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac, - eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing) + eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing, nothing) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 3d4132d72..f344fec8d 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -12,7 +12,7 @@ shimming. using DelimitedFiles using HDF5 -# Physical constants and collisionality normalizations shared by `load_kinetic_profiles` and +# Physical constants and collisionality normalizations shared by `kinetic_profile_splines` and # `resolve_ntv_species` (Krook collision frequency, Logan & Park 2013 Eq. 6). const _EV_J = 1.602e-19 # eV → J const _E_CHARGE = 1.602e-19 # elementary charge [C] (PENTRC-parity value, not CODATA) @@ -25,11 +25,14 @@ const _ME = 9.109_1e-31 # electron mass [kg] # NRL-formulary Coulomb logarithm (natural log); n_e in m⁻³, T_e in J. _coulomb_log(ne, Te) = 17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(Te / _KEV_J) +# Names kinetic profiles in error messages by their provenance, when they carry one. +_describe(data) = isempty(data.provenance) ? "kinetic profiles" : "kinetic profiles ($(data.provenance))" + # Hard stop on an unphysical negative density (bad input or a cubic-resample overshoot). -function _assert_nonneg_density(kinetic_file, label, arr, psi_reg) +function _assert_nonneg_density(source, label, arr, psi_reg) i = findfirst(<(0), arr) i === nothing || - error("kinetic file '$kinetic_file': negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical") + error("$source: negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical") end """ @@ -50,51 +53,6 @@ struct ResolvedNTVSpecies{P} profiles::P end -""" - 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 χ⊥ (shipped data: electron power-balance χ_e; see the stability docs) | -| `chi_phi` | toroidal momentum diffusivity χ_φ | -| `zeff` | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure | -| `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | -| `provenance` | short string recording the source file/format | -""" -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}} - zeff::Union{Nothing,Vector{Float64}} - species_densities::Union{Nothing,Dict{String,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, zeff=nothing, - species_densities=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), _v(zeff), species_densities, String(provenance)) -end - const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") # Dataset units written into / expected from the HDF5 kinetic schema. @@ -154,6 +112,27 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") end end +""" + attach_kinetic_profiles!(equil, path; group="/") -> equil + attach_kinetic_profiles!(equil, data::KineticProfileData) -> equil + +Attach raw kinetic profiles to `equil` as `equil.kinetic`. Every kinetic consumer (the NTV torque, the +calculated kinetic matrices, tearing, the layer-overlap scan and the two-pass grid refinement) derives its +own view from them, so the file is read once. Equilibria built from a config with `kinetic_file` set +attach it on construction; use this for analytic or IMAS equilibria built without it, or in-memory data. +""" +function attach_kinetic_profiles!(equil::PlasmaEquilibrium, data::KineticProfileData) + equil.kinetic = data + return equil +end + +attach_kinetic_profiles!(equil::PlasmaEquilibrium, path::AbstractString; group::AbstractString="/") = + attach_kinetic_profiles!(equil, read_kinetic_file(path; group=group)) + +# The kinetic profiles a config names, read once at construction; `nothing` when it names none. +configured_kinetic_profiles(config::EquilibriumConfig) = + isempty(config.kinetic_file) ? nothing : read_kinetic_file(config.kinetic_file; group=config.kinetic_group) + """ write_kinetic_h5(path, data; group="/", schema_version="1.0", provenance=data.provenance) @@ -233,7 +212,7 @@ function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; end """ - resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{ResolvedNTVSpecies} + resolve_ntv_species(data::KineticProfileData, ion_species; electron, zimp, mimp, ...) -> Vector{ResolvedNTVSpecies} Resolve the **full NTV species set** for a multi-ion run: the main ions, the neutrality-closing impurity, and (if `electron`) the electrons. Returns `ResolvedNTVSpecies` descriptors — the @@ -244,16 +223,16 @@ Each descriptor's `profiles` view carries **that species' resonant density** `n_ and its **full-composition collision frequency** `ν_s` (as `nui_spline`); `ne/Te/ωE/Zeff/lnΛ` (and the electron `nue`) are shared. Main-ion density is resolved from `IonSpecies`: `fraction` ⇒ `n_s = fraction · n_i` (the kinetic file's `n_i` = total main-ion density); `density` ⇒ an explicit -`n_*`-named dataset in the HDF5 kinetic file (e.g. `"n_D"`), cubic-resampled to the working grid. +`n_*`-named dataset of the HDF5 kinetic data (e.g. `"n_D"`), cubic-resampled to the working grid. In a mixed fraction/density list the impurity content is set by the file's n_i/n_e deficit, not by a fraction shortfall. Composition (Zeff, zpitch, `n_imp`) comes from `multi_ion_composition`; per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. The `zpitch·n_main` field term is exact for z=1 main ions (D-T, D-H) and approximate otherwise (a z>1 main ion's charge enters only through Zeff in the impurity closure). The impurity (`zimp`,`mimp`, density = the quasineutrality `n_imp`) is included as its own resonant species. Reduces to -`load_kinetic_profiles` for one z=1 main ion with `fraction=1`. +[`kinetic_profile_splines`](@ref) for one z=1 main ion with `fraction=1`. """ -function resolve_ntv_species(kinetic_file::AbstractString, ion_species::AbstractVector; +function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVector; electron::Bool=false, zimp::Integer=6, mimp::Integer=12, density_factor::Float64=1.0, temperature_factor::Float64=1.0, ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0) @@ -262,8 +241,8 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract toroidal_rotation_factor == 1.0) || error("profile scaling factors are not yet supported together with a multi-ion ion_species list") - data = read_kinetic_file(kinetic_file) - _need(f, n) = f === nothing ? error("kinetic file '$kinetic_file' missing required '$n'") : f + source = _describe(data) + _need(f, n) = f === nothing ? error("$source missing required '$n'") : f psi_in = data.psi ne_in = _need(data.n_e, "n_e") Ti_in = _need(data.T_i, "T_i") @@ -278,8 +257,8 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* _EV_J Te = _cubic_resample(psi_in, Te_in, psi_reg) .* _EV_J omegaE = _cubic_resample(psi_in, omE_in, psi_reg) - _assert_nonneg_density(kinetic_file, "n_e", ne, psi_reg) - _assert_nonneg_density(kinetic_file, "n_i (total)", ni_total, psi_reg) + _assert_nonneg_density(source, "n_e", ne, psi_reg) + _assert_nonneg_density(source, "n_i (total)", ni_total, psi_reg) # Fraction bookkeeping: fractions are shares of the total n_i, so they can never sum above 1, # and when every species is fraction-specified they must account for all of n_i. (In a mixed @@ -313,14 +292,14 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract if has_dens (data.species_densities !== nothing && haskey(data.species_densities, s.density)) || error( - "ion_species[$si]: density profile `$(s.density)` not found in kinetic file " * + "ion_species[$si]: density profile `$(s.density)` not found in the kinetic profiles " * "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))" ) ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg) else ns[si] = s.fraction .* ni_total end - _assert_nonneg_density(kinetic_file, "ion_species[$si] density", ns[si], psi_reg) + _assert_nonneg_density(source, "ion_species[$si] density", ns[si], psi_reg) end # Shared composition + collisionality (natural-log Coulomb log), per grid point. @@ -341,7 +320,7 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract nue[i] = Te[i] > 0 ? (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * (Te[i] / _KEV_J)^1.5) : 0.0 end # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical. - _assert_nonneg_density(kinetic_file, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg) + _assert_nonneg_density(source, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg) # Approaching the zpitch pole at Zeff = zimp signals a dilute main-ion mix taken up by impurity. maximum(zeff) > 0.9 * zimp && @warn "max(Zeff) = $(maximum(zeff)) is within 10% of zimp = $zimp — zpitch closure near its pole; check ion_species densities/fractions" @@ -378,28 +357,23 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract end """ - load_kinetic_profiles(kinetic_file::AbstractString; - zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, - density_factor::Float64=1.0, temperature_factor::Float64=1.0, - ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0, - chi1::Union{Nothing,Float64}=nothing) + kinetic_profile_splines(data::KineticProfileData; + zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, + density_factor::Float64=1.0, temperature_factor::Float64=1.0, + ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0, + chi1::Union{Nothing,Float64}=nothing) → KineticProfileSplines -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 +The NTV view of raw kinetic profiles (e.g. `equil.kinetic`): 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. -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] +Needs `n_e`, `T_i`, `T_e` and `omega_E`; `n_i` defaults to `n_e`. See +[`read_kinetic_file`](@ref) for the file formats. # Arguments - - `kinetic_file`: Path to the ASCII or HDF5 kinetic profile file (see `read_kinetic_file`) + - `data`: Raw kinetic profiles (`KineticProfileData`) - `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) @@ -424,17 +398,15 @@ Collisionality is recomputed from the (possibly scaled) profiles. This differs f Fortran PENTRC, which computes collisionality from unscaled profiles. Use `nufac` (in `KineticForcesControl`) for independent collisionality scaling. """ -function load_kinetic_profiles(kinetic_file::AbstractString; +function kinetic_profile_splines(data::KineticProfileData; zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, density_factor::Float64=1.0, temperature_factor::Float64=1.0, ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0, chi1::Union{Nothing,Float64}=nothing) - 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 + _need(field, name) = field === nothing ? error("$(_describe(data)) 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") diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 251e6da84..2a8d39621 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -434,9 +434,9 @@ function load_kinetic_context( needs_kinetic_profiles = haskey(inputs, "KineticForces") || (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") if needs_kinetic_profiles - kinetic_file = joinpath(intr.dir_path, kf_ctrl.kinetic_file) - kinetic_profiles = Equilibrium.load_kinetic_profiles( - kinetic_file; + kinetic_data = Equilibrium.read_kinetic_file(joinpath(intr.dir_path, kf_ctrl.kinetic_file)) + kinetic_profiles = Equilibrium.kinetic_profile_splines( + kinetic_data; zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass, density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, @@ -445,7 +445,7 @@ function load_kinetic_context( if !isempty(kf_ctrl.ion_species) || kf_ctrl.electron speclist = isempty(kf_ctrl.ion_species) ? [KineticForces.IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)] : kf_ctrl.ion_species - species = Equilibrium.resolve_ntv_species(kinetic_file, speclist; + species = Equilibrium.resolve_ntv_species(kinetic_data, speclist; electron=kf_ctrl.electron, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.impurity_mass, density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, @@ -501,7 +501,7 @@ function maybe_reform_equilibrium( else nothing # fall back to re-reading the input file end - equil = Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) + equil = Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes, kinetic=equil.kinetic) implied = Equilibrium.implied_knot_count(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles) if implied > 1.5 * (length(psi_nodes) - 1) @warn "Two-pass psi grid: refined equilibrium implies $implied knots vs $(length(psi_nodes) - 1) used — " * diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index eaa35bae7..fa7a8e64a 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -45,7 +45,7 @@ section. carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the KineticForces kernel needs but ForceFreeStatesControl does not expose. - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- - profile splines loaded via `Equilibrium.load_kinetic_profiles`. + profile splines built by `Equilibrium.kinetic_profile_splines`. # Returns diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 4d141f626..12b15a271 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -1,3 +1,4 @@ +using TOML @testset "KineticForces Unit Tests" begin @@ -689,8 +690,8 @@ end @testset "NTV splines identical: ASCII vs HDF5" begin - ka = E.load_kinetic_profiles(gpeckf) - kb = E.load_kinetic_profiles(h5) + ka = E.kinetic_profile_splines(E.read_kinetic_file(gpeckf)) + kb = E.kinetic_profile_splines(E.read_kinetic_file(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) @@ -712,12 +713,31 @@ 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" + # Only psi + n_e: the NTV view must error on the missing T_i 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) + @test_throws "missing required 'T_i'" E.kinetic_profile_splines(d) + end + + @testset "attached to the equilibrium" begin + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + sol = E.SolovevConfig(deck["SOL_INPUT"]) + # A relative kinetic_file resolves against the deck directory, like eq_filename. + config = E.EquilibriumConfig(merge(deck["Equilibrium"], Dict{String,Any}("kinetic_file" => relpath(h5, template))), template) + @test config.kinetic_file == normpath(h5) + raw = E.read_kinetic_file(h5) + equil = E.setup_equilibrium(config, sol) + @test equil.kinetic.psi == raw.psi + @test equil.kinetic.T_i == raw.T_i + @test equil.kinetic.chi_e == raw.chi_e + # A re-form passes the data through instead of reading the file again. + @test E.setup_equilibrium(config, sol; kinetic=equil.kinetic).kinetic === equil.kinetic + + bare = E.setup_equilibrium(E.EquilibriumConfig(deck["Equilibrium"], template), sol) + @test bare.kinetic === nothing + E.attach_kinetic_profiles!(bare, gpeckf) + @test bare.kinetic.n_e == E.read_kinetic_file(gpeckf).n_e + @test E.attach_kinetic_profiles!(bare, raw).kinetic === raw end end end diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index 773174c9e..9b4cfde52 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -49,9 +49,9 @@ using HDF5 end end - # Single species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles. - single = EQ.load_kinetic_profiles(f; zi=1, zimp=6, mi=2, mimp=12) - sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) + # Single species [z=1,m=2,fraction=1] reproduces kinetic_profile_splines. + single = EQ.kinetic_profile_splines(EQ.read_kinetic_file(f); zi=1, zimp=6, mi=2, mimp=12) + sp1 = EQ.resolve_ntv_species(EQ.read_kinetic_file(f), [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) main = sp1[1] for ψ in (0.3, 0.6) @test main.profiles.ni_spline(ψ) ≈ single.ni_spline(ψ) rtol = 1e-8 @@ -60,7 +60,7 @@ using HDF5 end # D-T-e set = {D, T, impurity, electron}; ν scales as 1/√m at fixed z, T. - sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)]; + sp = EQ.resolve_ntv_species(EQ.read_kinetic_file(f), [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)]; electron=true, zimp=6, mimp=12) @test length(sp) == 4 @test count(s -> !s.electron && s.z == 1, sp) == 2 # D, T @@ -71,11 +71,11 @@ using HDF5 @test sp[iD].profiles.nui_spline(0.5) / sp[iT].profiles.nui_spline(0.5) ≈ sqrt(3 / 2) rtol = 1e-6 # Exactly one of {fraction, density} required. - @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2)]; zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(f), [KF.IonSpecies(; z=1, m=2)]; zimp=6, mimp=12) # Unified electron semantics: a single-ion + electron run resolves as # {main ion, impurity, electron} — electrons in ADDITION to the ions. - spe = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=true, zimp=6, mimp=12) + spe = EQ.resolve_ntv_species(EQ.read_kinetic_file(f), [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=true, zimp=6, mimp=12) @test length(spe) == 3 @test !spe[1].electron && spe[1].z == 1 # main ion unchanged @test spe[1].profiles.ni_spline(0.5) ≈ main.profiles.ni_spline(0.5) @@ -101,14 +101,14 @@ using HDF5 file["n_D"] = 0.6 .* Ni file["n_T"] = 0.4 .* Ni # unequal explicit shapes end - sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; + sp = EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; electron=false, zimp=6, mimp=12) iD = findfirst(s -> s.m == 2, sp) iT = findfirst(s -> s.m == 3, sp) @test sp[iD].profiles.ni_spline(0.5) ≈ 0.6 * 1.0e20 * (1 - 0.25) rtol = 1e-6 @test sp[iT].profiles.ni_spline(0.5) ≈ 0.4 * 1.0e20 * (1 - 0.25) rtol = 1e-6 # A missing named profile errors clearly. - @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_He")]; zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, density="n_He")]; zimp=6, mimp=12) rm(h5) end @@ -127,7 +127,7 @@ using HDF5 @test back.species_densities["n_D"] ≈ data.species_densities["n_D"] @test back.species_densities["n_T"] ≈ data.species_densities["n_T"] # The written file feeds the explicit-density multi-ion input directly. - sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; + sp = EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; zimp=6, mimp=12) @test count(s -> !s.electron && s.z == 1, sp) == 2 # Non-`n_*` species names are rejected at write time (they would not round-trip). @@ -157,20 +157,20 @@ using HDF5 end DT(f1, f2) = [KF.IonSpecies(; z=1, m=2, fraction=f1), KF.IonSpecies(; z=1, m=3, fraction=f2)] # An all-fraction list must account for the full n_i (the impurity is set by n_i/n_e, not a shortfall). - @test_throws ErrorException EQ.resolve_ntv_species(h5, DT(0.5, 0.4); zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), DT(0.5, 0.4); zimp=6, mimp=12) # Fractions are shares of n_i, so summing above 1 is fatal in any list, mixed included. - @test_throws ErrorException EQ.resolve_ntv_species(h5, DT(0.7, 0.5); zimp=6, mimp=12) - @test_throws ErrorException EQ.resolve_ntv_species(h5, + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), DT(0.7, 0.5); zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, fraction=1.2), KF.IonSpecies(; z=1, m=3, density="n_T")]; zimp=6, mimp=12) # A partial fraction alongside an explicit density is legal — the n_i/n_e deficit is the impurity. - sp = EQ.resolve_ntv_species(h5, + sp = EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, density="n_T")]; zimp=6, mimp=12) @test count(s -> !s.electron && s.z == 1, sp) == 2 # Main-ion z must stay below the impurity charge (zpitch closure pole at Zeff = zimp). - @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=6, m=12, fraction=1.0)]; zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=6, m=12, fraction=1.0)]; zimp=6, mimp=12) # A dilute main-ion mix drives Zeff toward zimp and warns about the zpitch pole. @test_logs (:warn, r"zpitch closure near its pole") match_mode = :any EQ.resolve_ntv_species( - h5, [KF.IonSpecies(; z=1, m=2, density="n_dilute")]; zimp=6, mimp=12) + EQ.read_kinetic_file(h5), [KF.IonSpecies(; z=1, m=2, density="n_dilute")]; zimp=6, mimp=12) # Vacuum point: neutral composition, no impurity, finite zpitch. @test EQ.multi_ion_composition([1], [1.0e19], 0.0; zimp=6, mimp=12) == (1.0, 1.0, 1.0e19, 0.0) rm(h5) From bbde42250cea4233abd5b594349c041e22cd35c4 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:08:34 +0200 Subject: [PATCH 2/7] KineticForces - API! - Derive the NTV views from the equilibrium's kinetic profiles and delete [KineticForces] kinetic_file --- .../benchmark_diiid_ideal_ntv_torque.jl | 4 +- .../benchmark_diiid_kinetic_stability.jl | 2 +- docs/src/kinetic_forces.md | 4 +- .../DIIID-like_error_field_example/gpec.toml | 2 +- examples/DIIID-like_ideal_example/gpec.toml | 2 +- .../Solovev_kinetic_NTV_example/gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- examples/a10_kinetic_example/gpec.toml | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 148 +++++++----------- src/KineticForces/KineticForcesStructs.jl | 1 - src/KineticForces/Utils.jl | 32 ++++ test/runtests_error_fields.jl | 8 +- test/runtests_kinetic.jl | 26 +++ .../gpec.toml | 1 + .../gpec.toml | 1 + 16 files changed, 133 insertions(+), 106 deletions(-) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index d9de8fa81..f137fb916 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -280,7 +280,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) kf_ctrl = KF.KineticForcesControl(; fgar_flag=true, tgar_flag=tgar_enabled, nn=1, bounce_harmonic_max=4, ion_mass=2, ion_charge=1, nutype="harmonic", f0type="maxwellian", - collisionality_factor=1, psilims=[0.0, 1.0], kinetic_file=inputs.kin_file, + collisionality_factor=1, psilims=[0.0, 1.0], verbose=false) # Initialize KF internal state @@ -308,7 +308,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) t2 = time() kinetic_profiles = Eq.kinetic_profile_splines( - Eq.read_kinetic_file(kf_ctrl.kinetic_file); + Eq.read_kinetic_file(inputs.kin_file); zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass) kf_state = KF.KineticForcesState() diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index 1c87d9814..2e619152c 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/benchmarks/benchmark_diiid_kinetic_stability.jl @@ -122,6 +122,7 @@ function build_benchmark_tomldir(eq_file::String, kin_file::String) [Equilibrium] eq_filename = "$eq_name" +kinetic_file = "$kin_file" eq_type = "efit" jac_type = "hamada" grid_type = "pow1" # Matches Fortran DIIID_kinetic_example (equil/grid.f90 powspace) @@ -160,7 +161,6 @@ renorm_threshold = 1e4 write_outputs_to_HDF5 = true [KineticForces] -kinetic_file = "$kin_file" nn = 1 nl = 4 zi = 1 diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index b2e9fff1c..b4395bae9 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -42,8 +42,10 @@ paths: the post-PE ψ-quadrature diagnostic and the self-consistent `kinetic_source = "calculated"` matrices. ```toml -[KineticForces] +[Equilibrium] kinetic_file = "kinetic.h5" # n_i column/dataset = TOTAL main-ion density + +[KineticForces] electron = true # add electron NTV in addition to the ion species impurity_charge = 6 # impurity charge (closes quasineutrality) impurity_mass = 12 # impurity mass diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index a4f303562..a9cfaa10b 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/examples/DIIID-like_error_field_example/gpec.toml @@ -10,6 +10,7 @@ [Equilibrium] eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +kinetic_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "auto" # Radial grid packing type ("auto" = two-pass measured-curvature grid) @@ -240,7 +241,6 @@ rotation_span_factor = 1.0 # The scan covers ±span_factor × max(| rotation_offset_factor = 2.0 # Rough neoclassical offset as a multiple of the ion temperature-gradient diamagnetic frequency ω_*T [KineticForces] -kinetic_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # GPEC HDF5 kinetic schema (psi,n_i,n_e,T_i,T_e,omega_E,+chi_e,chi_phi) nn = 1 # Toroidal mode number (matches nn_low/nn_high) bounce_harmonic_max = 6 # Bounce harmonics retained ion_charge = 1 # Main-ion charge (deuterium) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 448ac5d07..7e6fa7d93 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -5,6 +5,7 @@ [Equilibrium] eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "auto" # Radial grid packing type ("auto" = two-pass measured-curvature grid) @@ -81,7 +82,6 @@ write_outputs_to_HDF5 = true # Write ErrorFields/CoilSensitivities/ t verbose = true # Log per-coil-set progress and linearity diagnostics [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) bounce_harmonic_max = 6 # Bounce harmonics retained ion_charge = 1 # Main-ion charge (deuterium) diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index cd02c559f..c8dbd63dc 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -7,6 +7,7 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) grid_type = "core_edge_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux @@ -86,5 +87,4 @@ b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixe f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) [KineticForces] -kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index a6bce31b9..ebf68e397 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -8,6 +8,7 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "core_edge_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized flux coordinate @@ -66,7 +67,6 @@ write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = false # Enable verbose logging [KineticForces] -kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns [[KineticForces.ion_species]] z = 1 # Species charge (fundamental units) diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index bee59220e..c039dd0b2 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -7,6 +7,7 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) grid_type = "core_edge_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux @@ -88,7 +89,6 @@ f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, [KineticForces] # Multi-ion NTV: the analytic Solovev main ion is modelled as a 50/50 D-T mix (n_i = total), # so the summed NTV runs over D + T + the neutrality-closing impurity + electrons. -kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns electron = true # Add electron NTV in addition to the ion species impurity_charge = 6 # Impurity charge (fundamental units); closes quasineutrality impurity_mass = 12 # Impurity mass (proton masses) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index c7a32b36c..8511b7639 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -8,6 +8,7 @@ [Equilibrium] eq_filename = "fix_a100_k10_q2_bn010_prof1" # Path to equilibrium file +kinetic_file = "a10_prof1.txt" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "core_edge_packed" # Radial grid packing type @@ -47,7 +48,6 @@ write_outputs_to_HDF5 = true # Write outputs to the HDF5 file force_termination = true # Stop after force-free states — the kinetic stability eigenvalues are the result of interest [KineticForces] -kinetic_file = "a10_prof1.txt" # Path to the kinetic profile file (bundled) nn = 1 # Toroidal mode number bounce_harmonic_max = 6 # Bounce-harmonic number ion_charge = 1 # Ion charge (fundamental units) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 2a8d39621..76b91e9a3 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -222,9 +222,8 @@ function main_from_inputs( equil = Equilibrium.setup_equilibrium(eq_config, additional_input) - kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) - equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; - kf_ctrl=kf_ctrl, kf_species=kf_species) + kf_ctrl = kinetic_forces_control(inputs, ctrl, equil) + equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl; kf_ctrl=kf_ctrl) equil_dt = time() - equil_start push!(runtimes, "equilibrium" => equil_dt) @@ -257,8 +256,7 @@ function main_from_inputs( ffs_start = time() locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; - species=kf_species, + metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl; overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path)) ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric; runtimes=runtimes) @@ -276,11 +274,11 @@ function main_from_inputs( @info "Results written to $(ctrl.HDF5_filename)" end - # Validity describes the kinetic model, not the torque contraction, so it is written wherever - # kinetic profiles were loaded — including FFS-only runs that never reach run_kinetic_forces. - if ctrl.write_outputs_to_HDF5 && kinetic_profiles !== nothing + # Validity describes the kinetic model, not the torque contraction, so it is written for every + # run that reads the kinetic profiles through KineticForces — including FFS-only runs. + if ctrl.write_outputs_to_HDF5 && _needs_kinetic_profiles(inputs, ctrl) h5open(joinpath(ffs_result.dir_path, ctrl.HDF5_filename), "cw") do h5file - KineticForces.write_validity!(h5file, kf_ctrl, kf_species, kinetic_profiles, equil, + KineticForces.write_validity!(h5file, kf_ctrl, KineticForces.ntv_species(kf_ctrl, equil), KineticForces.ntv_profiles(kf_ctrl, equil), equil, Float64[sng.psifac for sng in intr.sing]) end end @@ -325,7 +323,7 @@ function main_from_inputs( end kf_start = time() - run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species) + run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl) if "KineticForces" in keys(inputs) push!(runtimes, "kinetic_forces" => time() - kf_start) # Mirrors the write gate inside `run_kinetic_forces` (needs a PE state to contract against). @@ -335,7 +333,7 @@ function main_from_inputs( end ef_start = time() - error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) + error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl) coil_sensitivities, monte_carlo, locking_risk, efc_couplings = error_fields === nothing ? (nothing, nothing, nothing, nothing) : error_fields if "ErrorFields" in keys(inputs) push!(runtimes, "error_fields" => time() - ef_start) @@ -407,80 +405,48 @@ function resolve_mode_space!(intr::ForceFreeStatesInternal, ctrl::ForceFreeState end """ - load_kinetic_context(inputs, intr, ctrl, equil) -> (kf_ctrl, kinetic_profiles, species) + kinetic_forces_control(inputs, ctrl, equil) -> kf_ctrl -Build the KineticForces control and load the kinetic profiles once for the whole run. -`kinetic_profiles` is `nothing` when no stage asks for them. +The `[KineticForces]` control of the run, defaults when the deck has no such section. A run that +needs kinetic profiles (a `[KineticForces]` section, or the calculated kinetic source) errors +here when the equilibrium carries none; every consumer derives its view from `equil.kinetic`. """ -function load_kinetic_context( - inputs::Dict{String,Any}, - intr::ForceFreeStatesInternal, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium -) - # The profiles are reused by the grid refinement, the stability kinetic callback (via - # `calculated_cb`), and the post-PE torque diagnostics block. The `"fixed"` kinetic source - # path in stability does not need kinetic_profiles, but the post-PE block always does, so we - # load whenever a [KineticForces] section is present or the stability path requests the - # calculated source. psio is invariant across grid re-formation. - kf_ctrl = - haskey(inputs, "KineticForces") ? - KineticForces.KineticForcesControl(; - (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) : - KineticForces.KineticForcesControl() - - kinetic_profiles = nothing - species = nothing - needs_kinetic_profiles = haskey(inputs, "KineticForces") || - (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") - if needs_kinetic_profiles - kinetic_data = Equilibrium.read_kinetic_file(joinpath(intr.dir_path, kf_ctrl.kinetic_file)) - kinetic_profiles = Equilibrium.kinetic_profile_splines( - kinetic_data; - zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, - mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass, - density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, - ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor, - chi1=2π * equil.psio) - if !isempty(kf_ctrl.ion_species) || kf_ctrl.electron - speclist = isempty(kf_ctrl.ion_species) ? - [KineticForces.IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)] : kf_ctrl.ion_species - species = Equilibrium.resolve_ntv_species(kinetic_data, speclist; - electron=kf_ctrl.electron, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.impurity_mass, - density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, - ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, - toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor) - end - end - - return kf_ctrl, kinetic_profiles, species +function kinetic_forces_control(inputs::Dict{String,Any}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) + kf_ctrl = haskey(inputs, "KineticForces") ? + KineticForces.KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) : + KineticForces.KineticForcesControl() + _needs_kinetic_profiles(inputs, ctrl) && equil.kinetic === nothing && + error("this run needs kinetic profiles ([KineticForces], or kinetic_source = \"calculated\"): set kinetic_file in [Equilibrium]") + return kf_ctrl end +# Whether a run reads the kinetic profiles through KineticForces (its stage, or the calculated EL source). +_needs_kinetic_profiles(inputs, ctrl) = haskey(inputs, "KineticForces") || (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") + """ - maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; - kf_ctrl=nothing, kf_species=nothing) -> equil + maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl; kf_ctrl=KineticForcesControl()) -> equil -Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic -profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory -input — no file re-read. Returns `equil` untouched when the configuration wants a single pass. +Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, and the NTV +view of its attached kinetic profiles), pin knots on rational surfaces, and re-form on the refined +grid from the in-memory input — no file re-read. Returns `equil` untouched when the configuration +wants a single pass. """ function maybe_reform_equilibrium( equil::Equilibrium.PlasmaEquilibrium, eq_config::Equilibrium.EquilibriumConfig, additional_input, intr::ForceFreeStatesInternal, - ctrl::ForceFreeStatesControl, - kinetic_profiles; - kf_ctrl::Union{Nothing,KineticForces.KineticForcesControl}=nothing, - kf_species=nothing + ctrl::ForceFreeStatesControl; + kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl() ) Equilibrium.wants_two_pass(eq_config) || return equil + kinetic_profiles = KineticForces.ntv_profiles(kf_ctrl, equil) mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) # Kinetic-resonance nodes are located here rather than at the call site so the validity # boundary that filters them is cleared against this pass's rationals (intr.sing is empty yet). - pinned = kf_ctrl === nothing ? Float64[] : - KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr; rationals=mandatory) + pinned = KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, KineticForces.ntv_species(kf_ctrl, equil), equil, intr; + rationals=mandatory) # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), # so the rational-surface brackets clear a zone large enough for every mode. n_min = minimum(abs(n) for n in intr.nlow:intr.nhigh if n != 0) @@ -629,7 +595,7 @@ function _layer_overlap_scan(overlap_inputs, intr::ForceFreeStatesInternal, equi end """ - prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; overlap_inputs=nothing) -> (metric, mats, overlap) + prepare_force_free_states!(intr, ctrl, equil, kf_ctrl; overlap_inputs=nothing) -> (metric, mats, overlap) Set up the force-free-states solve on `intr`: integration limits, the surviving singular surfaces and their GGJ coefficients, the poloidal mode range, and the metric plus @@ -639,9 +605,7 @@ function prepare_force_free_states!( intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - kf_ctrl::KineticForces.KineticForcesControl, - kinetic_profiles; - species=nothing, + kf_ctrl::KineticForces.KineticForcesControl; overlap_inputs=nothing ) # Resistive-layer overlap: locate where adjacent rational surfaces' layers run into each @@ -734,6 +698,8 @@ function prepare_force_free_states!( # Inject the KineticForces callback so the "calculated" source can # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). + kinetic_profiles = KineticForces.ntv_profiles(kf_ctrl, equil) + species = KineticForces.ntv_species(kf_ctrl, equil) # One boundary for the whole run: the widest-orbit species, moved clear of any rational # surface whose window the transition band would otherwise cut through. axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, @@ -953,19 +919,18 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat resolve_mode_space!(intr, ctrl) - # The API path never reads kinetic profiles, so the KineticForces control is only the - # placeholder `prepare_force_free_states!` threads into its (unused) callback. + # The problem has no [KineticForces] section, so the views of the attached kinetic profiles use the default control. kf_ctrl = KineticForces.KineticForcesControl() if Equilibrium.wants_two_pass(equil.config) && equil.ingest === nothing @warn "Two-pass auto grid needs the equilibrium's raw ingest, which analytic and IMAS equilibria do not carry; " * "solving on the single-pass grid. Set mpsi explicitly to choose the grid." else - equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl, nothing) + equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl) end locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) + metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl) result = run_force_free_states(ctrl, equil, mats, intr, metric) if ctrl.write_outputs_to_HDF5 @@ -1090,24 +1055,24 @@ function perturbed_equilibrium( end """ - run_kinetic_forces(inputs, result, pe_state, kf_ctrl, kinetic_profiles) + run_kinetic_forces(inputs, result, pe_state, kf_ctrl) Compute and write the neoclassical toroidal viscosity torque diagnostics when the deck carries a -`[KineticForces]` section. No-op when the perturbed-equilibrium state the operators contract +`[KineticForces]` section, from the kinetic profiles attached to `result.equil`. No-op when the perturbed-equilibrium state the operators contract against is missing. """ function run_kinetic_forces( inputs::Dict{String,Any}, result::ForceFreeStatesResult, pe_state, - kf_ctrl::KineticForces.KineticForcesControl, - kinetic_profiles, - species=nothing + kf_ctrl::KineticForces.KineticForcesControl ) # ---------------------------------------------------------------- # KineticForces (Neoclassical Toroidal Viscosity) # ---------------------------------------------------------------- ("KineticForces" in keys(inputs)) || return nothing + kinetic_profiles = KineticForces.ntv_profiles(kf_ctrl, result.equil) + species = KineticForces.ntv_species(kf_ctrl, result.equil) @info "\n KineticForces\n$_SECTION" kf_start = time() @@ -1121,7 +1086,6 @@ function run_kinetic_forces( if pe_state === nothing @info "Skipping NTV torque diagnostics: no perturbed-equilibrium data (e.g. kinetic_source=\"calculated\")." else - # kf_ctrl and kinetic_profiles were loaded once before the equilibrium was re-formed. kf_intr = KineticForces.KineticForcesInternal(result.equil; verbose=kf_ctrl.verbose) KineticForces.set_perturbation_data!(kf_intr, pe_state, result, result.equil, result.metric) @@ -1185,7 +1149,7 @@ function forcing_terms_control(inputs::Dict{String,Any}) end """ - run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) -> (sensitivities, monte_carlo, risk, couplings) or nothing + run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl) -> (sensitivities, monte_carlo, risk, couplings) or nothing Linearize every coil set's resonant drive with respect to its rigid shifts and tilts and write `ErrorFields/CoilSensitivities/` when the deck carries an `[ErrorFields]` section. When the @@ -1194,7 +1158,7 @@ Carlo on the full-window dominant mode with the `[ErrorFields.MonteCarlo]` setti `ErrorFields/MonteCarlo/`; with an `[ErrorFields.scenario]` table as well, evaluate the locking risk (and the tolerance scan when `[ErrorFields.Risk]` names `scan_scales`) and write `ErrorFields/Risk/`; with `[ErrorFields.NTV]` naming correction arrays, evaluate their overlap and NTV -torque couplings per kilo-ampere-turn (needs the kinetic context) and write `ErrorFields/NTV/`. +torque couplings per kilo-ampere-turn (needs the equilibrium's kinetic profiles) and write `ErrorFields/NTV/`. Stages not requested return `nothing`. Needs the perturbed-equilibrium state's singular-coupling matrix and coil-format forcing, and errors otherwise: a deck asking for error-field sensitivities without them is a misconfiguration, not a @@ -1205,8 +1169,7 @@ function run_error_fields( result::ForceFreeStatesResult, pe_state, preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}, - kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl(), - kinetic_profiles=nothing + kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl() ) ("ErrorFields" in keys(inputs)) || return nothing @@ -1286,7 +1249,7 @@ function run_error_fields( missing = setdiff(ntv_ctrl.efc_coils, [cs.name for cs in efc_sets]) isempty(missing) || error("[ErrorFields.NTV] efc_coils not among the run's coil sets: $(join(missing, ", "))") ntv_start = time() - couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, + couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, rotation_scan=ntv_ctrl.rotation_scan, scan_points=ntv_ctrl.rotation_scan_points, scan_max_points=ntv_ctrl.rotation_scan_max_points, scan_tolerance=ntv_ctrl.rotation_scan_tolerance, span_factor=ntv_ctrl.rotation_span_factor, offset_factor=ntv_ctrl.rotation_offset_factor) @info "NTV couplings of $(length(couplings)) correction arrays in $(@sprintf("%.1f", time() - ntv_start)) s" @@ -1311,7 +1274,7 @@ function run_error_fields( end """ - efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} + efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl; mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} The couplings of each correction coil array in `coil_sets` per kilo-ampere-turn: its dominant-mode overlap and resonant fraction from its spectrum on the control surface, and its @@ -1320,8 +1283,8 @@ torque is a quadratic form of the applied spectrum, so each is one plasma-respon of the unit-current spectrum (injected as forcing modes, nothing written) followed by the kinetic torque of `method`; both scale exactly with the square of the current. `cfg` is the run's `[ForcingTerms]` coil configuration, so the boundary grid matches the one the run's -sensitivities were evaluated on; `rc` must be the coupling of `ffs`'s own solve and -`kinetic_profiles` its kinetic context. The spectrum is normalized by the magnitude of the +sensitivities were evaluated on; `rc` must be the coupling of `ffs`'s own solve; the +torque uses the NTV view of the kinetic profiles attached to `ffs.equil`. The spectrum is normalized by the magnitude of the array's ampere-turns, `|nw| × max|I|`, keeping its winding sense and current pattern as the phase reference; the torque is stored with its sign. """ @@ -1331,8 +1294,7 @@ function efc_couplings( rc::PerturbedEquilibrium.ResonantCoupling, dom::PerturbedEquilibrium.DominantCoupling, cfg::ForcingTerms.CoilConfig, - kf_ctrl::KineticForces.KineticForcesControl, - kinetic_profiles; + kf_ctrl::KineticForces.KineticForcesControl; mode::Int=1, method::AbstractString="fgar", rotation_scan::Bool=true, @@ -1343,7 +1305,7 @@ function efc_couplings( offset_factor::Real=2.0, verbose::Bool=false ) - kinetic_profiles === nothing && error("efc_couplings needs kinetic profiles: add a [KineticForces] section with a kinetic_file") + ffs.equil.kinetic === nothing && error("efc_couplings needs kinetic profiles: set kinetic_file in [Equilibrium]") getfield(kf_ctrl, Symbol(method * "_flag")) || error("efc_couplings: KineticForces method \"$method\" is not enabled in the control") equil = ffs.equil grids = [(n, ForcingTerms.CoilForcingGrid(equil, cfg, n; psi=ffs.psilim)) for n in sort(unique(rc.n_modes))] @@ -1353,7 +1315,7 @@ function efc_couplings( # Reference rotation and scan span from the kinetic profiles, with the torque kernel's own # diamagnetic frequencies (Logan & Park 2013 Eq. 7): ω_φ = ω_E + ω_*n + ω_*T. - kp = kinetic_profiles + kp = KineticForces.ntv_profiles(kf_ctrl, equil) chi1 = 2π * equil.psio chrg = kf_ctrl.ion_charge * Utilities.E_CHG ψk = kp.xs diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 9014dd9a4..2152fed5f 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -175,7 +175,6 @@ builds a second control for its differing tolerance). save_records::Bool = false # Save detailed integration trajectories # Input files (relative to dir_path) - kinetic_file::String = "kinetic.dat" gpec_file::String = "gpec.h5" # Diagnostic flags diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 0c9615839..6a2cf0057 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -292,6 +292,38 @@ function clear_rational_windows(psi_c::Float64, rationals::Vector{Float64})::Flo return psi_c end +""" + ntv_profiles(kf_ctrl, equil) → Union{Nothing,KineticProfileSplines} + +The NTV view of the kinetic profiles attached to `equil` (`equil.kinetic`), with the main-ion and +impurity charges and masses and the profile scale factors of `kf_ctrl`; `nothing` when the +equilibrium carries no kinetic profiles. +""" +function ntv_profiles(kf_ctrl::KineticForcesControl, equil) + equil.kinetic === nothing && return nothing + return Equilibrium.kinetic_profile_splines(equil.kinetic; + zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass, + density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, + ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor, + chi1=2π * equil.psio) +end + +""" + ntv_species(kf_ctrl, equil) → Union{Nothing,Vector{ResolvedNTVSpecies}} + +The resolved NTV species set of the kinetic profiles attached to `equil` (see +`Equilibrium.resolve_ntv_species`) when `kf_ctrl` asks for a species sum (`ion_species` or +`electron`); `nothing` for the single main-ion run or when the equilibrium carries no kinetic profiles. +""" +function ntv_species(kf_ctrl::KineticForcesControl, equil) + (equil.kinetic !== nothing && (!isempty(kf_ctrl.ion_species) || kf_ctrl.electron)) || return nothing + speclist = isempty(kf_ctrl.ion_species) ? [IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)] : kf_ctrl.ion_species + return Equilibrium.resolve_ntv_species(equil.kinetic, speclist; + electron=kf_ctrl.electron, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.impurity_mass, + density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, + ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor) +end + """ axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil) → Float64 diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index ccd9e098e..3a9440dbe 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -80,7 +80,8 @@ include("h5_metadata_check.jl") cp(joinpath(@__DIR__, "test_data", "ErrorFields", "tolerances_two_hoops.toml"), joinpath(dir, "tolerances.toml")) # Kinetic profiles consistent with this Solovev pressure, for the NTV torque couplings. cp(joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat"), joinpath(dir, "kinetic.dat")) - inputs["KineticForces"] = Dict{String,Any}("kinetic_file" => "kinetic.dat", "verbose" => false, "write_outputs_to_HDF5" => true) + inputs["Equilibrium"]["kinetic_file"] = "kinetic.dat" + inputs["KineticForces"] = Dict{String,Any}("verbose" => false, "write_outputs_to_HDF5" => true) inputs["ErrorFields"] = Dict{String,Any}("verbose" => false, "tolerance_file" => "tolerances.toml", "MonteCarlo" => Dict{String,Any}("nsample" => 20_000, "nbatch" => 2, "seed" => 5, "nbins" => 100), "scenario" => Dict{String,Any}("n_e" => 12.0), @@ -531,7 +532,10 @@ include("h5_metadata_check.jl") @test length(with_profiles.subplots) == 3 # and the deposition profile its own panel @test with_profiles[3].series_list[1][:y] == c.torque_full_profile[:, 2] @test with_profiles[3].series_list[2][:y] == c.torque_residual_profile[:, 2] - @test_throws ErrorException GPEC.efc_couplings(ffs, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl(), nothing) + attached = ffs.equil.kinetic + ffs.equil.kinetic = nothing + @test_throws "needs kinetic profiles" GPEC.efc_couplings(ffs, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl()) + ffs.equil.kinetic = attached # Central differences: doubling the step moves the derivatives at O(h²). coarse = EF.compute_coil_sensitivities(sets, rc, ffs.equil, cfg, diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 12b15a271..db7560016 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -739,5 +739,31 @@ using TOML @test bare.kinetic.n_e == E.read_kinetic_file(gpeckf).n_e @test E.attach_kinetic_profiles!(bare, raw).kinetic === raw end + + @testset "NTV views of the attached profiles" begin + GPEC = GeneralizedPerturbedEquilibrium + KF = GPEC.KineticForces + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = E.setup_equilibrium(E.EquilibriumConfig(deck["Equilibrium"], template), E.SolovevConfig(deck["SOL_INPUT"])) + kf = KF.KineticForcesControl() + @test KF.ntv_profiles(kf, equil) === nothing + @test_throws "set kinetic_file in [Equilibrium]" GPEC.kinetic_forces_control( + Dict{String,Any}("KineticForces" => Dict{String,Any}()), GPEC.ForceFreeStates.ForceFreeStatesControl(), equil) + + raw = E.read_kinetic_file(h5) + E.attach_kinetic_profiles!(equil, raw) + view = KF.ntv_profiles(kf, equil) + ref = E.kinetic_profile_splines(raw; zi=kf.ion_charge, zimp=kf.impurity_charge, mi=kf.ion_mass, mimp=kf.impurity_mass, + chi1=2π * equil.psio) + for ψ in 0.0:0.1:1.0 + @test view.Ti_spline(ψ) == ref.Ti_spline(ψ) + @test view.nui_spline(ψ) == ref.nui_spline(ψ) + @test view.zeff_spline(ψ) == ref.zeff_spline(ψ) + end + @test KF.ntv_species(kf, equil) === nothing + with_e = KF.ntv_species(KF.KineticForcesControl(; electron=true), equil) + @test any(sp -> sp.electron, with_e) + end end end diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index a05caa035..233024d6d 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -4,6 +4,7 @@ [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "core_edge_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized poloidal flux diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index 5df1debce..5ca86fb47 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -5,6 +5,7 @@ # tracked separately by the solovev_kinetic_nuzero regression case. [Equilibrium] eq_type = "sol" # Type of the input 2D equilibrium file +kinetic_file = "kinetic.dat" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "core_edge_packed" # Radial grid packing type psilow = 1e-4 # Lower limit of normalized flux coordinate From 4f89be7785d9da796885dbf050e292eb95723ba0 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:16:49 +0200 Subject: [PATCH 3/7] Equilibrium - REFACTOR - Drive the HDF5 kinetic-profile reader and writer from one dataset table --- src/Equilibrium/KineticProfiles.jl | 31 +++++++++++------------------- test/runtests_kinetic.jl | 6 ++++++ 2 files changed, 17 insertions(+), 20 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index f344fec8d..bd318e0c2 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -55,9 +55,9 @@ 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", +# The optional profile datasets of the HDF5 kinetic schema and their units, one row per `KineticProfileData` +# field of the same name; the reader and writer both iterate it, so a new profile is a row here plus its field. +const _KINETIC_H5_PROFILES = ("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", "zeff" => "dimensionless") """ @@ -97,17 +97,14 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density # profiles for explicit multi-ion input; other names are reserved for future schema fields. - standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi", "zeff") extras = Dict{String,Vector{Float64}}() for k in keys(g) - (k in standard || !startswith(k, "n_")) && continue + (k == "psi" || any(p -> first(p) == k, _KINETIC_H5_PROFILES) || !startswith(k, "n_")) && continue v = read(g[k]) v isa AbstractArray && (extras[k] = Float64.(vec(v))) end 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"), zeff=rd("zeff"), + (Symbol(name) => rd(name) for (name, _) in _KINETIC_H5_PROFILES)..., species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end end @@ -146,21 +143,15 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData; provenance::AbstractString=data.provenance) h5open(path, "w") do f g = group == "/" ? f : create_group(f, group) - function put(name, v) + function put(name, v, units) v === nothing && return g[name] = collect(Float64, v) - attributes(g[name])["units"] = _KINETIC_H5_UNITS[name] + attributes(g[name])["units"] = units + end + put("psi", data.psi, "normalized poloidal flux") + for (name, units) in _KINETIC_H5_PROFILES + put(name, getfield(data, Symbol(name)), units) 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) - put("zeff", data.zeff) # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip. if data.species_densities !== nothing for (name, v) in data.species_densities diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index db7560016..802d37d67 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -712,6 +712,12 @@ using TOML rm(tmp; force=true) end + @testset "schema table matches the data fields" begin + # Every profile dataset the reader and writer iterate is a KineticProfileData field, and vice versa. + table = Set(Symbol(first(p)) for p in E._KINETIC_H5_PROFILES) + @test table == setdiff(Set(fieldnames(E.KineticProfileData)), Set([:psi, :species_densities, :provenance])) + end + @testset "missing required field errors" begin # Only psi + n_e: the NTV view must error on the missing T_i d = E.KineticProfileData(; psi=[0.0, 0.5, 1.0], n_e=[1e19, 9e18, 8e18]) From eee0393b4842545cb9b751e78be1cd02b69b06a2 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:32:51 +0200 Subject: [PATCH 4/7] Tearing - API! - Build the SLAYER and layer-overlap inputs from the equilibrium's kinetic profiles --- docs/src/api.md | 4 +- docs/src/stability.md | 12 ++--- examples/DIIID-like_SLAYER_example/gpec.toml | 2 +- .../DIIID-like_truncation_example/gpec.toml | 4 +- regression-harness/src/runner.jl | 2 +- src/ForceFreeStates/CoreTypes.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 17 ++++--- src/Tearing/Runner/Control.jl | 20 +++------ src/Tearing/Runner/Runner.jl | 2 +- src/Tearing/Runner/TearingProblem.jl | 16 ++++--- src/Tearing/Runner/run_slayer.jl | 45 +++++++------------ src/Utilities/KineticProfiles.jl | 4 +- test/runtests_layer_overlap_fitzpatrick.jl | 16 ++++--- test/runtests_slayer_runner.jl | 27 ++++++++--- 14 files changed, 85 insertions(+), 88 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index eb81bee19..85de34710 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -126,13 +126,13 @@ extraction is not implemented yet, so `GGJModel()` errors here. ```julia ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # Δ′ matrix for the dispersion -tear = solve(TearingProblem(ffs; profile_file="kin.h5"), SLAYERModel()) # coupling_mode=:coupled returns one global root instead +tear = solve(TearingProblem(ffs), SLAYERModel()) # coupling_mode=:coupled returns one global root instead tear.gamma_Hz, tear.rational_q # root-found rate at each rational surface ``` Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedure knobs (scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters); -kinetic profiles come from `profile_file`, relative to the run directory. +kinetic profiles come from the equilibrium the solve was posed on (see Kinetic profiles above). Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and remain TOML-driven. diff --git a/docs/src/stability.md b/docs/src/stability.md index b056753ab..03fdbe2cd 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -316,12 +316,12 @@ TokaMaker's `save_eqdsk` writes such files when `lcfs_pad` (``p``) is set with t ``\psi = 1/(1 - p)``. Writing with `truncate_eq=False` keeps the whole plasma; EFIT-written files place their boundary on the X-point flux. -The scan takes its kinetic profiles from the `[SLAYER]` `profile_file`, through the same control -and loader the SLAYER analysis uses, so both see the same profiles, the same HDF5 group and the -same χ⊥(ψ) — the layer widths cannot disagree between them. A deck without a `[SLAYER]` -`profile_file` does not get the scan. Because the file is resolved from the run's TOML, the scan -only runs on the `gpec.toml` entry point: the programmatic `solve(prob, alg)` path has no inputs to -resolve it from, so the cap cannot be applied there and setting the flag warns instead. Its result is always written to +The scan takes its kinetic profiles from the equilibrium (`[Equilibrium] kinetic_file`) through the +same view and `[SLAYER]` control the SLAYER analysis uses, so both see the same profiles and the +same χ⊥(ψ) — the layer widths cannot disagree between them. A deck without a `[SLAYER]` section or +kinetic profiles does not get the scan. Because the layer settings come from the run's `[SLAYER]` +section, the scan only runs on the `gpec.toml` entry point: the programmatic `solve(prob, alg)` path +has no such section, so the cap cannot be applied there and setting the flag warns instead. Its result is always written to `ForceFreeStates/LayerOverlap/`, but it only constrains the domain when `psilim_from_layer_overlap = true`. It then caps `qlim` before the `rational_offset_fraction` step, so `rational_offset_fraction` still selects the final surface from inside the cap. A cap beyond `psihigh` has no effect. diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index 08df13c11..e0f8d4759 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -6,6 +6,7 @@ [Equilibrium] eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +kinetic_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc) grid_type = "log_asymptotic" # Radial grid packing type @@ -83,7 +84,6 @@ polish_roots = true # refine each root to the true zero; enable store_scan = false # per-surface scans omitted from HDF5 (uncoupled has one per surface) # Standardized kinetic-profile file (GPEC HDF5 kinetic schema), shared with the # sibling ideal example; carries n_e, T_e, T_i, omega_E, chi_e (χ⊥), chi_phi (χ_φ). -profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled [SLAYER.scan_grid] Q_re_range = [-2.0, 2.0] # Scan box in the normalized Q plane, Re(Q) axis diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml index eb2180e03..29f5ec798 100644 --- a/examples/DIIID-like_truncation_example/gpec.toml +++ b/examples/DIIID-like_truncation_example/gpec.toml @@ -8,6 +8,7 @@ [Equilibrium] eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +kinetic_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (HDF5 GPEC kinetic schema or 6-column ASCII), read once into the equilibrium eq_type = "efit" # Type of the input 2D equilibrium file jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) grid_type = "auto" # Radial grid packing type ("auto" = two-pass measured-curvature grid) @@ -45,6 +46,5 @@ psilim_from_layer_overlap = true # Cap psilim where adjacent resistive layers [SLAYER] # The growth-rate analysis is not run here; the section is present because the layer-overlap scan -# takes its kinetic profiles, and its chi_perp(psi), from the SLAYER profile_file. +# builds its layer widths, and its chi_perp(psi), from the SLAYER settings and the equilibrium's kinetic profiles. enabled = false # Run the SLAYER tearing-mode analysis -profile_file = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5" # Kinetic-profile file (relative to the run dir), required when SLAYER is enabled diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index 2e7ec5415..ab6b88a35 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -33,7 +33,7 @@ With no `overrides`, the example deck runs in place (returns it untouched). With the deck is copied to a uniquely named throwaway sibling of the example dir and the named gpec.toml keys are patched there, so one shared example can serve several cases (e.g. a collisionless variant via `"KineticForces.nutype" => "zero"`). The copy sits beside the original so relative file -references in the deck (e.g. `SLAYER.profile_file = "..//..."`) still resolve. +references in the deck (e.g. `Equilibrium.kinetic_file = "..//..."`) still resolve. Override keys are dotted paths into the TOML; missing intermediate tables are created. Returns `(rundir, is_temp)`; the caller removes `rundir` when `is_temp`. """ diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index ef384af7b..c88df1030 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -134,7 +134,7 @@ gpec.toml. - `numunorms_init::Int` - Initial array size for solution normalization data - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced - `truncate_at_rational_offset::Bool` - Truncate the integration domain at `(last_rational_q + rational_offset_fraction) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`rational_offset_fraction = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `truncate_at_rational_offset = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `truncate_at_rational_offset = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + rational_offset_fraction / n" depends on which `n`); when `truncate_at_rational_offset = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`. - - `psilim_from_layer_overlap::Bool` - Cap the integration domain where adjacent rational surfaces' resistive layers overlap, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9: past that point no surface retains a well-separated inner region. Applied as an upper bound on `qlim` at the top of `sing_lim!`, so `rational_offset_fraction` / `qhigh` still select the final surface from inside it. The scan runs whenever the run's `[SLAYER]` section names a readable `profile_file` — the same profiles and χ⊥(ψ) the SLAYER analysis uses — and is recorded under `ForceFreeStates/LayerOverlap/` either way; this flag only decides whether it constrains the domain. It resolves that file from the run's TOML, so the programmatic `solve(prob, alg)` entry point cannot run it and the flag is inert there (with a warning). Off on every shipped deck; with the shipped DIII-D-like profiles the overlap point (ψ ≈ 0.9985) lies outside the 0.995 decks' domain but inside the SLAYER deck's 0.9995. See the integration-domain truncation section of the stability documentation. Default `false`. + - `psilim_from_layer_overlap::Bool` - Cap the integration domain where adjacent rational surfaces' resistive layers overlap, after Fitzpatrick, Nucl. Fusion 2025 Sect. 5.9: past that point no surface retains a well-separated inner region. Applied as an upper bound on `qlim` at the top of `sing_lim!`, so `rational_offset_fraction` / `qhigh` still select the final surface from inside it. The scan runs whenever the run has a `[SLAYER]` section and the equilibrium carries kinetic profiles — the same profiles and χ⊥(ψ) the SLAYER analysis uses — and is recorded under `ForceFreeStates/LayerOverlap/` either way; this flag only decides whether it constrains the domain. It resolves that file from the run's TOML, so the programmatic `solve(prob, alg)` entry point cannot run it and the flag is inert there (with a warning). Off on every shipped deck; with the shipped DIII-D-like profiles the overlap point (ψ ≈ 0.9985) lies outside the 0.995 decks' domain but inside the SLAYER deck's 0.9995. See the integration-domain truncation section of the stability documentation. Default `false`. - `rational_offset_fraction::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `truncate_at_rational_offset` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here. - `frobenius_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed). - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 76b91e9a3..c682f9a74 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -257,7 +257,7 @@ function main_from_inputs( locstab, ballooning_boundary = run_local_stability(ctrl, equil) metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl; - overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path)) + overlap_inputs=_overlap_slayer_inputs(inputs, equil.kinetic)) ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric; runtimes=runtimes) if ctrl.write_outputs_to_HDF5 @@ -532,17 +532,16 @@ function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.Pl return locstab, ballooning_boundary end -# Kinetic profiles for the overlap scan come from the [SLAYER] section, through the same control -# and loader the SLAYER analysis itself uses, so the scan sees the same profiles, the same HDF5 -# group and the same chi(psi). Returns `nothing` when there is no usable [SLAYER] profile_file. -function _overlap_slayer_inputs(inputs, dir_path::AbstractString) - (inputs isa AbstractDict && haskey(inputs, "SLAYER")) || return nothing +# Inputs of the overlap scan: the SLAYER view of the equilibrium's kinetic profiles under the run's +# [SLAYER] control, the same view the SLAYER analysis itself uses, so the scan sees the same +# profiles and the same chi(psi). Returns `nothing` without a [SLAYER] section or kinetic profiles. +function _overlap_slayer_inputs(inputs, kinetic::Union{Nothing,Equilibrium.KineticProfileData}) + (inputs isa AbstractDict && haskey(inputs, "SLAYER") && kinetic !== nothing) || return nothing ctrl = Runner.slayer_control_from_toml(inputs["SLAYER"]) - isempty(ctrl.profile_file) && return nothing loaded = try - Runner._load_profiles(ctrl, dir_path) + Runner._slayer_profiles(kinetic, ctrl) catch err - @warn "Layer-overlap scan skipped: could not read the [SLAYER] profile_file." exception = err + @warn "Layer-overlap scan skipped: the kinetic profiles do not serve the SLAYER layer." exception = err return nothing end # chi(psi) and Z_eff(psi) from the file where it carries them, else the control's scalar diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 94aed97f4..4ac3826dc 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -84,17 +84,12 @@ constructor. # Kinetic-profile source -SLAYER reads kinetic profiles through the shared `Equilibrium.read_kinetic_file` -reader (the same standardized object used by the kinetic/NTV physics), so -there is one consistent interface for resistive and kinetic profiles. - - - `profile_file` -- path to a kinetic-profile file (relative to the run - dir), required when SLAYER is enabled. HDF5 (`.h5`) files use the GPEC - kinetic schema and may carry `chi_e` (χ⊥) and `chi_phi` (χ_φ); ASCII - tables are also accepted but carry no χ. When a χ dataset is absent or - all-zero, the scalar `chi_perp`/`chi_tor` fallbacks below are used (so a - file can keep the χ keys set to 0 to defer to the scalars). - - `profile_group` -- group within the HDF5 file (default `"/"`) +SLAYER takes its kinetic profiles from the equilibrium (`[Equilibrium] kinetic_file`, +read once into `equil.kinetic`), the same object the kinetic/NTV physics reads. +They need `n_e`, `T_e`, `T_i` and `omega_E`. HDF5 files may also carry `chi_e` (χ⊥) +and `chi_phi` (χ_φ); ASCII tables carry no χ. When a χ dataset is absent or all-zero, +the scalar `chi_perp`/`chi_tor` fallbacks below are used (so a file can keep the χ +keys set to 0 to defer to the scalars). # Output control @@ -158,9 +153,6 @@ there is one consistent interface for resistive and kinetic profiles. # / failed-Δ'-BVP surface, not a real root. Flagged `:spurious`. validity_rtol::Float64 = 1e-3 - profile_file::String = "" - profile_group::String = "/" - store_scan::Bool = false end diff --git a/src/Tearing/Runner/Runner.jl b/src/Tearing/Runner/Runner.jl index 152bce0a2..6c48a66ea 100644 --- a/src/Tearing/Runner/Runner.jl +++ b/src/Tearing/Runner/Runner.jl @@ -31,7 +31,7 @@ import CommonSolve using FastInterpolations: cubic_interp using ..Utilities using ..Utilities: KineticProfiles -using ...Equilibrium: read_kinetic_file, KineticProfileData +using ...Equilibrium: KineticProfileData using ...ForceFreeStates: ForceFreeStatesResult using ..InnerLayer using ..InnerLayer: InnerLayerParameters, InnerLayerResponse, solve_inner, diff --git a/src/Tearing/Runner/TearingProblem.jl b/src/Tearing/Runner/TearingProblem.jl index 76b730f51..4d3b71b64 100644 --- a/src/Tearing/Runner/TearingProblem.jl +++ b/src/Tearing/Runner/TearingProblem.jl @@ -8,11 +8,12 @@ Tearing-mode growth rates from a finished force-free-states solve: hold the outer Δ′ fixed and find the growth rate at which the inner layer matches it. Solve with `solve(prob, SLAYERModel())`; GGJ growth-rate extraction is not implemented yet. Keywords are -[`SLAYERControl`](@ref) fields, the `[SLAYER]` deck knobs; kinetic profiles come from -`profile_file`. Without a full Δ′ matrix the per-surface Δ′ stubs are used, with a warning. +[`SLAYERControl`](@ref) fields, the `[SLAYER]` deck knobs; kinetic profiles come from the +equilibrium (`ffs.equil.kinetic`, see `attach_kinetic_profiles!`). Without a full Δ′ matrix the +per-surface Δ′ stubs are used, with a warning. ## Fields - - `ffs` - The finished solve; duck-typed on `equil`, `surfaces`, `delta_prime`, `dir_path`, so tests can use stand-ins. + - `ffs` - The finished solve; duck-typed on `equil` (with its `kinetic` profiles), `surfaces` and `delta_prime`, so tests can use stand-ins. - `control::SLAYERControl` - The tearing knobs. """ struct TearingProblem{R} @@ -77,7 +78,9 @@ function CommonSolve.solve(prob::TearingProblem, model::InnerLayer.InnerLayerMod control.enabled || return empty_slayer_result(control) isempty(surfaces) && return empty_slayer_result(control) - loaded = _load_profiles(control, ffs.dir_path) + equil.kinetic === nothing && + error("TearingProblem: the equilibrium carries no kinetic profiles; set kinetic_file in [Equilibrium] or attach_kinetic_profiles!") + loaded = _slayer_profiles(equil.kinetic, control) params = _tearing_params(model, equil, surfaces, loaded, control) @@ -110,9 +113,8 @@ function CommonSolve.solve(prob::TearingProblem, model::InnerLayer.InnerLayerMod rational_psi = Float64[surfaces[p.ising].psifac for p in params] rational_q = Float64[surfaces[p.ising].q for p in params] - # E×B rotation per unit n at each surface, for the coupled determinant's Doppler shifts; - # empty when the kinetic file carries none. - omega_E = loaded.has_omega_E ? Float64[loaded.profiles.omega_E(ψ) for ψ in rational_psi] : Float64[] + # E×B rotation per unit n at each surface, for the coupled determinant's Doppler shifts. + omega_E = Float64[loaded.profiles.omega_E(ψ) for ψ in rational_psi] return run_slayer_from_inputs(model, params, dp, control; rational_psi=rational_psi, rational_q=rational_q, omega_E=omega_E) end diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2143ed011..321a9df1b 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -4,36 +4,24 @@ # `run_slayer_from_inputs` (scan + growth-rate extraction). `TearingProblem` builds its inputs. # --------------------------------------------------------------------- -# Profile loading +# Kinetic-profile view # --------------------------------------------------------------------- -# Read kinetic profiles through the shared `Equilibrium.read_kinetic_file` -# reader and adapt them to the SLAYER inner layer. Returns the spline-based -# `KineticProfiles` plus χ⊥(ψ)/χ_φ(ψ) callables built from the file's -# `chi_e`/`chi_phi` (or `nothing` when the file carries no usable χ, in which -# case the caller falls back to the scalar `control.chi_perp`/`chi_tor`), and -# a Z_eff(ψ) callable: the file's own `zeff` profile when it carries one, else -# derived from its `n_i/n_e` (or `nothing` when the file carries neither, in which -# case the caller falls back to the scalar `control.zeff`), and whether the file -# carries an `omega_E` profile. -function _load_profiles(control::SLAYERControl, dir_path::AbstractString) - isempty(control.profile_file) && - error("run_slayer: [SLAYER] profile_file is empty — point it at a " * - "kinetic-profile file (HDF5 GPEC kinetic schema or ASCII table).") - path = isabspath(control.profile_file) ? control.profile_file : - joinpath(dir_path, control.profile_file) - data = read_kinetic_file(path; group=control.profile_group) - - for (name, v) in (("n_e", data.n_e), ("T_e", data.T_e), ("T_i", data.T_i)) +# The SLAYER view of raw kinetic profiles (`equil.kinetic`), on their own ψ grid. Returns the +# spline-based `KineticProfiles` plus χ⊥(ψ)/χ_φ(ψ) callables built from `chi_e`/`chi_phi` (or +# `nothing` when the data carry no usable χ, in which case the caller falls back to the scalar +# `control.chi_perp`/`chi_tor`), and a Z_eff(ψ) callable: the data's own `zeff` profile when they +# carry one, else derived from `n_i/n_e` (or `nothing` when they carry neither, in which case the +# caller falls back to the scalar `control.zeff`). The E×B rotation `omega_E` is required: the +# layer's Doppler shift has no default. +function _slayer_profiles(data::KineticProfileData, control::SLAYERControl) + source = isempty(data.provenance) ? "kinetic profiles" : "kinetic profiles ($(data.provenance))" + for (name, v) in (("n_e", data.n_e), ("T_e", data.T_e), ("T_i", data.T_i), ("omega_E", data.omega_E)) v === nothing && - error("run_slayer: kinetic file '$path' is missing required " * - "dataset '$name' for the SLAYER inner layer.") + error("run_slayer: $source missing required dataset '$name' for the SLAYER inner layer.") end - # `omega_E` carries the E×B rotation Ω_E (per unit n) when the file has it, else zero (no shift); # ω_*e/ω_*i are derived per surface from the density and temperature splines by `build_slayer_inputs`. - npsi = length(data.psi) - omega_E = data.omega_E === nothing ? zeros(npsi) : data.omega_E profiles = KineticProfiles(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, - T_i=data.T_i, omega_E=omega_E) + T_i=data.T_i, omega_E=data.omega_E) # χ⊥(ψ)/χ_φ(ψ) splines from the file. A χ array that is absent OR all-zero # is treated as "not provided" — χ must be positive (χ=0 ⇒ τ_⊥→∞), and the @@ -52,19 +40,18 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString) _usable(v) = v !== nothing && !all(iszero, v) zeff = if _usable(data.zeff) minimum(data.zeff) >= 1.0 - 1e-8 || - error("run_slayer: kinetic file '$path' has zeff = $(minimum(data.zeff)) < 1 — unphysical") + error("run_slayer: $source has zeff = $(minimum(data.zeff)) < 1 — unphysical") cubic_interp(psi_xs, collect(Float64, data.zeff)) elseif !_usable(data.n_i) nothing else z = [ne > 0 ? control.impurity_charge - (ni / ne) * (control.impurity_charge - 1.0) : control.impurity_charge for (ni, ne) in zip(data.n_i, data.n_e)] minimum(z) >= 1.0 - 1e-8 || - error("run_slayer: kinetic file '$path' has n_i > n_e, giving Z_eff = $(minimum(z)) < 1 — unphysical") + error("run_slayer: $source has n_i > n_e, giving Z_eff = $(minimum(z)) < 1 — unphysical") cubic_interp(psi_xs, z) end zeff_source = _usable(data.zeff) ? :file : zeff === nothing ? :scalar : :closure - return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source, - has_omega_E=data.omega_E !== nothing) + return (profiles=profiles, chi_perp=chi_perp, chi_tor=chi_tor, zeff=zeff, zeff_source=zeff_source) end # --------------------------------------------------------------------- diff --git a/src/Utilities/KineticProfiles.jl b/src/Utilities/KineticProfiles.jl index 5b900b859..b13497431 100644 --- a/src/Utilities/KineticProfiles.jl +++ b/src/Utilities/KineticProfiles.jl @@ -22,8 +22,8 @@ the normalized poloidal flux ψ ∈ [0, 1]. | `omega_E` | E×B rotation Ω_E, per unit n | rad/s | Construct via the keyword constructor `KineticProfiles(; psi, n_e, T_e, T_i, omega_E)` with -matched-length vectors. The SLAYER runner builds this object from a standardized -kinetic-profile file via `Equilibrium.read_kinetic_file`. Diamagnetic frequencies are not +matched-length vectors. The SLAYER runner builds this object from the equilibrium's +kinetic profiles (`equil.kinetic`). Diamagnetic frequencies are not carried here: `build_slayer_inputs` derives ω_*e and ω_*i per surface from the density and temperature splines. diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl index d532f205c..ea2e2b295 100644 --- a/test/runtests_layer_overlap_fitzpatrick.jl +++ b/test/runtests_layer_overlap_fitzpatrick.jl @@ -228,15 +228,15 @@ end # The scan must see the run's own plasma: every [SLAYER] setting the layer width depends on # reaches it, and changing one changes the widths. trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example") - slayer_section = TOML.parsefile(joinpath(trunc_dir, "gpec.toml"))["SLAYER"] + trunc_deck = TOML.parsefile(joinpath(trunc_dir, "gpec.toml")) + slayer_section = trunc_deck["SLAYER"] + kin = Equilibrium.read_kinetic_file(joinpath(trunc_dir, trunc_deck["Equilibrium"]["kinetic_file"])) slayer_inputs(over) = Dict{String,Any}("SLAYER" => merge(slayer_section, over)) - overlap_inputs(over; dir=trunc_dir) = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(slayer_inputs(over), dir) + overlap_inputs(over; kinetic=kin) = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(slayer_inputs(over), kinetic) oi = overlap_inputs(Dict("ion_mass" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter")) @test (oi.mu_i, oi.bt, oi.theta, oi.lnLambda_form) == (3.0, 1.7, 0.3, :sauter) # Z_eff(psi) comes from the kinetic file's own n_i/n_e with the [SLAYER] impurity charge. - kin_path = joinpath(trunc_dir, slayer_section["profile_file"]) - kin = Equilibrium.read_kinetic_file(kin_path) j = length(kin.psi) ÷ 2 @test oi.zeff(kin.psi[j]) ≈ 6.0 - kin.n_i[j] / kin.n_e[j] * 5.0 @test overlap_inputs(Dict("impurity_charge" => 10.0)).zeff(kin.psi[j]) ≈ 10.0 - kin.n_i[j] / kin.n_e[j] * 9.0 @@ -248,10 +248,12 @@ end base = (; psi=kin.psi, n_e=kin.n_e, T_i=kin.T_i, T_e=kin.T_e, omega_E=kin.omega_E, chi_e=kin.chi_e, chi_phi=kin.chi_phi) Equilibrium.write_kinetic_h5(joinpath(tmp, "own_zeff.h5"), Equilibrium.KineticProfileData(; base..., n_i=kin.n_i, zeff=own)) @test Equilibrium.read_kinetic_file(joinpath(tmp, "own_zeff.h5")).zeff == own - @test overlap_inputs(Dict("profile_file" => "own_zeff.h5", "impurity_charge" => 10.0, "zeff" => 2.5); dir=tmp).zeff(kin.psi[j]) ≈ 3.2 + @test overlap_inputs(Dict("impurity_charge" => 10.0, "zeff" => 2.5); kinetic=Equilibrium.read_kinetic_file(joinpath(tmp, "own_zeff.h5"))).zeff(kin.psi[j]) ≈ 3.2 Equilibrium.write_kinetic_h5(joinpath(tmp, "no_ni.h5"), Equilibrium.KineticProfileData(; base...)) - @test overlap_inputs(Dict("profile_file" => "no_ni.h5", "zeff" => 2.5); dir=tmp).zeff == 2.5 + @test overlap_inputs(Dict("zeff" => 2.5); kinetic=Equilibrium.read_kinetic_file(joinpath(tmp, "no_ni.h5"))).zeff == 2.5 end + # No kinetic profiles: no scan. + @test overlap_inputs(Dict{String,Any}(); kinetic=nothing) === nothing widths(over) = GeneralizedPerturbedEquilibrium._layer_overlap_scan(overlap_inputs(over), _fresh(), equil).width_dr w_ref = widths(Dict{String,Any}()) @@ -327,7 +329,7 @@ end # The scan's overlap point must therefore not depend on where the grid stops, and no # rational may fall between the in-grid search and the law. - oi = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(trunc_inputs, trunc_dir) + oi = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(trunc_inputs, long.kinetic) one_n() = (i = ForceFreeStatesInternal(); i.nlow = 1; i.nhigh = 1; i) scan_long = GeneralizedPerturbedEquilibrium._layer_overlap_scan(oi, one_n(), long) scan_short = GeneralizedPerturbedEquilibrium._layer_overlap_scan(oi, one_n(), short) diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index deeed6109..83fc4d28c 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -68,8 +68,7 @@ "max_cells" => 50_000), "growth_rate_filter" => Dict{String,Any}( "pole_threshold" => 1e5, - "filter_above_poles" => false), - "profile_file" => "profiles.h5" + "filter_above_poles" => false) ) c = slayer_control_from_toml(section) @test c.enabled @@ -88,7 +87,6 @@ @test c.amr_max_cells == 50_000 @test c.pole_threshold == 1e5 @test c.filter_above_poles == false - @test c.profile_file == "profiles.h5" # Unknown keys should raise bad = merge(section, Dict{String,Any}("mistyped_key" => 42)) @@ -98,15 +96,32 @@ @testset "TearingProblem solve forwards surfaces and Δ'" begin # A result with no singular surfaces short-circuits before any equilibrium access, # so a stand-in result is enough to pin the forwarding of surfaces / delta_prime. - c = SLAYERControl(; enabled=true, profile_file="unused.h5") + c = SLAYERControl(; enabled=true) no_surfaces = (equil=nothing, surfaces=GeneralizedPerturbedEquilibrium.ForceFreeStates.SingType[], delta_prime=nothing) r = solve(TearingProblem(no_surfaces, c), SLAYERModel()) @test isempty(r.params) # A disabled control never looks at the result at all. - r_off = solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=false, profile_file="unused.h5")), SLAYERModel()) + r_off = solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=false)), SLAYERModel()) @test r_off.enabled == false + + # With surfaces to analyse, the profiles must come from the equilibrium. + bare = (equil=(; kinetic=nothing), surfaces=[nothing], delta_prime=nothing) + @test_throws "carries no kinetic profiles" solve(TearingProblem(bare, c), SLAYERModel()) + end + + @testset "SLAYER view of the kinetic profiles" begin + kin_path = joinpath(dirname(@__DIR__), "examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode_kinetic.h5") + data = GeneralizedPerturbedEquilibrium.Equilibrium.read_kinetic_file(kin_path) + view = Runner._slayer_profiles(data, SLAYERControl()) + j = length(data.psi) ÷ 2 + @test view.profiles.n_e(data.psi[j]) ≈ data.n_e[j] + @test view.profiles.omega_E(data.psi[j]) ≈ data.omega_E[j] + @test view.chi_perp(data.psi[j]) ≈ data.chi_e[j] + # A missing E×B rotation is an error, never a silent zero. + no_rotation = GeneralizedPerturbedEquilibrium.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, T_i=data.T_i) + @test_throws "missing required dataset 'omega_E'" Runner._slayer_profiles(no_rotation, SLAYERControl()) end @testset "GGJ tearing is rejected until its growth-rate extraction exists" begin @@ -116,7 +131,7 @@ end no_surfaces = (equil=nothing, surfaces=GeneralizedPerturbedEquilibrium.ForceFreeStates.SingType[], delta_prime=nothing, dir_path=".") - @test_throws ErrorException solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=true, profile_file="unused.h5")), GGJModel()) + @test_throws ErrorException solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=true)), GGJModel()) end # Δ′ is unified across formalisms, so a Galerkin run feeds SLAYER exactly as a Riccati one From 0c0005e31b4674745b4d1cb6d681670d601718c3 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:55:25 +0200 Subject: [PATCH 5/7] ForceFreeStates - FEATURE - Run kinetic solves from the scripting API on an equilibrium with kinetic profiles --- docs/src/api.md | 5 +-- src/GeneralizedPerturbedEquilibrium.jl | 15 ++++---- test/runtests_solve_api.jl | 48 +++++++++++++++++++++++++- 3 files changed, 58 insertions(+), 10 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index 85de34710..294ea61b3 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -134,8 +134,9 @@ Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedur (scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters); kinetic profiles come from the equilibrium the solve was posed on (see Kinetic profiles above). -Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and remain -TOML-driven. +Kinetic runs (`kinetic_factor > 0`) take `Forward()`. The `"calculated"` source reads the +equilibrium's kinetic profiles, with the default `KineticForcesControl` species and collision +settings. ## Entry points diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index c682f9a74..baee69790 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -836,8 +836,8 @@ end dir_path=".", debug=DebugSettings(), kwargs...) The perturbed-plasma Euler-Lagrange problem posed on an equilibrium: the extremization of -the perturbed potential energy whose solutions are the force-free (and, via the TOML path, -kinetic) perturbed states. This is the WHAT of a stability solve; the integrator passed to +the perturbed potential energy whose solutions are the force-free (or, with +`kinetic_factor > 0`, kinetic) perturbed states. This is the WHAT of a stability solve; the integrator passed to [`solve`](@ref) is the HOW. A `PlasmaEquilibrium` hosts many possible problems — this type names this one, so `solve` stays unambiguous as other problem classes appear. Inner-layer matching is a separate problem posed on the finished solve: see `MatchProblem`. @@ -890,9 +890,10 @@ Solve the perturbed-plasma [`EulerLagrangeProblem`](@ref) with the formalism `al a `gpec.toml` run of `main` does and produces the same result object. The second form is sugar building the problem from an equilibrium and the problem keywords in one call. -Knobs owned by `alg` are rejected as `ForceFreeStatesControl` keywords. Kinetic -runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces]` profiles -and errors here. +Knobs owned by `alg` are rejected as `ForceFreeStatesControl` keywords. Kinetic runs +(`kinetic_factor > 0`) take `Forward()`; the `"calculated"` source needs kinetic profiles on the +equilibrium (`kinetic_file`, or `attach_kinetic_profiles!`) and uses the default +`KineticForcesControl` species and collision settings. ```julia eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") @@ -909,8 +910,8 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat ForceFreeStates._apply_alg!(ctrl_kwargs, alg) ctrl = ForceFreeStatesControl(; ctrl_kwargs...) - ctrl.kinetic_factor > 0 && - error("kinetic runs (kinetic_factor > 0) need the [KineticForces] profiles and are TOML-driven; run them through `main`") + ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" && equil.kinetic === nothing && + error("kinetic_source = \"calculated\" needs kinetic profiles on the equilibrium: set kinetic_file or call attach_kinetic_profiles!") intr = ForceFreeStatesInternal(; dir_path=prob.dir_path) intr.wall_settings = prob.wall diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 2a238b3a3..3667b866a 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -76,6 +76,51 @@ using TOML end end + @testset "calculated kinetic solve matches the TOML-driven run" begin + # The equilibrium carries the deck's kinetic profiles, so the API solves the same kinetic problem. + kin_template = joinpath(@__DIR__, "test_data", "regression_solovev_kinetic_calculated") + kdeck = TOML.parsefile(joinpath(kin_template, "gpec.toml")) + kequil = GPEC.Equilibrium.setup_equilibrium( + GPEC.Equilibrium.EquilibriumConfig(kdeck["Equilibrium"], kin_template), + GPEC.Equilibrium.SolovevConfig(kdeck["SOL_INPUT"]) + ) + @test kequil.kinetic !== nothing + kwall = GPEC.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in kdeck["Wall"])...) + kkwargs = Dict(Symbol(k) => v for (k, v) in kdeck["ForceFreeStates"] if !(k in ("integrator", "nn_low", "nn_high"))) + mktempdir() do dir + for name in readdir(kin_template) + cp(joinpath(kin_template, name), joinpath(dir, name)) + end + reference = GPEC.main([dir]).ffs + api = solve(kequil, Forward(); nn=1, wall=kwall, dir_path=dir, kkwargs...) + @test api.mats.kinetic !== nothing + @test api.free_boundary.et[1] ≈ reference.free_boundary.et[1] rtol = 1e-10 + @test api.diagnostics.nzero == reference.diagnostics.nzero + @test api.kinetic.kmsing == reference.kinetic.kmsing + end + end + + @testset "attached kinetic profiles refine the two-pass grid" begin + # main and solve both re-form through maybe_reform_equilibrium, which reads equil.kinetic, so + # the attached profiles shape the grid the same way on either path. The DIII-D-like profiles + # add knots over the ideal ones. + sl_dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_SLAYER_example") + sl_deck = TOML.parsefile(joinpath(sl_dir, "gpec.toml")) + sl_config = GPEC.Equilibrium.EquilibriumConfig(sl_deck["Equilibrium"], sl_dir) + pass1 = GPEC.Equilibrium.setup_equilibrium(sl_config) + @test pass1.kinetic !== nothing + sl_ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in sl_deck["ForceFreeStates"])...) + sl_intr = FFS.ForceFreeStatesInternal() + GPEC.resolve_mode_space!(sl_intr, sl_ctrl) + refined = GPEC.maybe_reform_equilibrium(pass1, sl_config, nothing, sl_intr, sl_ctrl) + @test refined.kinetic === pass1.kinetic # carried through the re-form, not re-read + kinetic = pass1.kinetic + pass1.kinetic = nothing + ideal = GPEC.maybe_reform_equilibrium(pass1, sl_config, nothing, sl_intr, sl_ctrl) + pass1.kinetic = kinetic + @test refined.rzphi_xs != ideal.rzphi_xs + end + @testset "galerkin solve returns a gal result with a unified delta_prime" begin mktempdir() do dir api = solve(equil, Galerkin(; nx=32); nn=1, wall=wall, dir_path=dir, ffs_kwargs...) @@ -188,7 +233,8 @@ using TOML end @testset "rejected keyword combinations" begin - @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5) + @test_throws "needs kinetic profiles on the equilibrium" solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5, + kinetic_source="calculated") @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., integrator="riccati") @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., nchunks=8) @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., nn_low=2) From 559d96d20ed0b9aa53dec16336c35f39ca29263f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 29 Oct 2026 10:10:23 -0400 Subject: [PATCH 6/7] Rerun - BUGFIX - Attach the stored deck's kinetic profiles to the restored equilibrium The replay constructor still passed the pre-kinetic positional arguments, so every post-hoc entry point built on equilibrium_from_h5 failed with a MethodError. The stored deck's kinetic_file is re-read against the h5's directory; a file that has since moved warns and leaves the profiles off. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- src/Rerun.jl | 10 +++++++++- test/runtests_frozen_equilibrium.jl | 14 ++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/Rerun.jl b/src/Rerun.jl index 94fbfdfca..5ba420ffc 100644 --- a/src/Rerun.jl +++ b/src/Rerun.jl @@ -321,6 +321,14 @@ function rebuild_equilibrium_inputs(inputs::Dict{String,Any}, ingest, output_dir return eq_config, additional_input end +# The stored deck's kinetic profiles, re-read from its file; a file that has since moved leaves the equilibrium without them. +function _replay_kinetic_profiles(config::Equilibrium.EquilibriumConfig) + isempty(config.kinetic_file) && return nothing + isfile(config.kinetic_file) && return Equilibrium.configured_kinetic_profiles(config) + @warn "The stored deck's kinetic_file is not readable; the restored equilibrium carries no kinetic profiles" file = config.kinetic_file + return nothing +end + """ equilibrium_from_h5(h5path) -> (; equil, inputs, psilim) @@ -365,7 +373,7 @@ function equilibrium_from_h5(h5path::AbstractString) geometry = Equilibrium.compute_geometry_profiles(xs, ys, rzphi[1], rzphi[2], rzphi[4], params.ro) analytic = haskey(Equilibrium.ANALYTIC_EQ, config.eq_type) ? build_analytic_config(config.eq_type, inputs) : nothing equil = Equilibrium.PlasmaEquilibrium(config, params, profiles, geometry, xs, ys, rzphi..., eqfun..., params.ro, params.zo, params.psio, - read_equilibrium_ingest(f), analytic) + read_equilibrium_ingest(f), analytic, _replay_kinetic_profiles(config)) return (; equil, inputs, psilim=Float64(read(f, "Info/psilim"))) end end diff --git a/test/runtests_frozen_equilibrium.jl b/test/runtests_frozen_equilibrium.jl index 3829c3c4a..23a52beeb 100644 --- a/test/runtests_frozen_equilibrium.jl +++ b/test/runtests_frozen_equilibrium.jl @@ -115,6 +115,20 @@ using GeneralizedPerturbedEquilibrium @test analytic_snapshot.analytic isa Eq.SolovevConfig @test analytic_snapshot.analytic.elongation == 1.8 @test analytic_snapshot.profiles.q_spline.y == profile_values[4] + # The stored deck's kinetic_file resolves against the h5's directory and is re-read; a missing file warns and leaves none. + kinetic_path = joinpath(dir, "kinetic.dat") + cp(joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat"), kinetic_path) + analytic_inputs["Equilibrium"]["kinetic_file"] = "kinetic.dat" + h5open(path, "r+") do f + delete_object(f, "Input/gpec_toml_raw") + f["Input/gpec_toml_raw"] = sprint(io -> TOML.print(io, analytic_inputs)) + end + with_kinetic = GPEC.equilibrium_from_h5(path).equil + @test with_kinetic.kinetic.n_e == Eq.read_kinetic_file(kinetic_path).n_e + @test with_kinetic.config.kinetic_file == kinetic_path + rm(kinetic_path) + without = @test_logs (:warn, r"kinetic_file is not readable") match_mode = :any GPEC.equilibrium_from_h5(path).equil + @test without.kinetic === nothing h5open(path, "r+") do f delete_object(f, "Equilibrium/Profiles/q") end From 255586db517a69ab68b9465a02919bc7ebbeb13d Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 29 Oct 2026 10:10:23 -0400 Subject: [PATCH 7/7] Repo - BUGFIX - Stop decks that still name the kinetic-profile file under KineticForces or SLAYER and fix the docs cross-reference A deck carrying the removed KineticForces kinetic_file or SLAYER profile_file/profile_group keys now stops with the new home named, in main and on replay, and the example-deck hook rejects the SLAYER keys. The scripting-API page qualifies the attach_kinetic_profiles! reference, which only the Equilibrium module exports. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- .pre-commit-config.yaml | 2 +- docs/src/api.md | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 20 ++++++++++++++------ test/runtests_equil.jl | 14 ++++++++++++++ 4 files changed, 30 insertions(+), 8 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d82cfe6c1..d82d998a2 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -77,5 +77,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag|profile_file|profile_group)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/docs/src/api.md b/docs/src/api.md index f6ffb9e03..854a2e274 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -52,7 +52,7 @@ The objects map one-to-one onto the pipeline stages: The equilibrium carries the kinetic profiles. Name the file with the `kinetic_file` (and, for HDF5, `kinetic_group`) keyword, the same keys as the `[Equilibrium]` deck section, and it is read once into `equil.kinetic`; every kinetic stage derives its own view from it. -[`attach_kinetic_profiles!`](@ref) attaches a file or in-memory `KineticProfileData` to an +[`Equilibrium.attach_kinetic_profiles!`](@ref) attaches a file or in-memory `KineticProfileData` to an equilibrium built without one, such as an analytic equilibrium. ```julia diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 382dc9427..873f112e8 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -100,10 +100,14 @@ const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_ "gal_flag", "nstep", "diagnose_ca", "diagnose", "ksing", "local_stability_flag") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc", "use_galgrid") const _DEPRECATED_KF_KEYS = ("moment", "eq_out", "theta_out", "xlmda_out", "eqpsi_out", "save_records", "gpec_file", - "fnml_flag", "ellip_flag", "diag_flag", "force_xialpha") + "fnml_flag", "ellip_flag", "diag_flag", "force_xialpha", "kinetic_file") const _DEPRECATED_PE_KEYS = ("fixed_boundary", "output_eigenmodes", "filter_modes", "singular_point_method", "compute_response") +const _DEPRECATED_SLAYER_KEYS = ("profile_file", "profile_group") const _DEPRECATED_KEYS = ("Equilibrium" => _DEPRECATED_EQUIL_KEYS, "ForceFreeStates" => _DEPRECATED_FFS_KEYS, - "KineticForces" => _DEPRECATED_KF_KEYS, "PerturbedEquilibrium" => _DEPRECATED_PE_KEYS) + "KineticForces" => _DEPRECATED_KF_KEYS, "PerturbedEquilibrium" => _DEPRECATED_PE_KEYS, "SLAYER" => _DEPRECATED_SLAYER_KEYS) +# The kinetic-profile file is named once, in [Equilibrium]; these keys stop with that remedy rather than the generic one. +const _MOVED_KINETIC_KEYS = ("KineticForces" => ("kinetic_file",), "SLAYER" => _DEPRECATED_SLAYER_KEYS) +const _MOVED_KINETIC_REMEDY = "name the kinetic-profile file once, as kinetic_file (and kinetic_group) in [Equilibrium]" # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated. function _stop_on_deprecated_keys(table, deprecated_keys, section::String; remedy::String="remove from gpec.toml") @@ -236,6 +240,9 @@ function main_from_inputs( intr = ForceFreeStatesInternal(; dir_path=path) ffs_table = inputs["ForceFreeStates"] _stop_on_deprecated_keys(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates") + for (section, keys) in _MOVED_KINETIC_KEYS + haskey(inputs, section) && _stop_on_deprecated_keys(inputs[section], keys, section; remedy=_MOVED_KINETIC_REMEDY) + end haskey(inputs, "KineticForces") && _stop_on_deprecated_keys(inputs["KineticForces"], _DEPRECATED_KF_KEYS, "KineticForces") haskey(inputs, "PerturbedEquilibrium") && _stop_on_deprecated_keys(inputs["PerturbedEquilibrium"], _DEPRECATED_PE_KEYS, "PerturbedEquilibrium") ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...) @@ -460,9 +467,10 @@ needs kinetic profiles (a `[KineticForces]` section, or the calculated kinetic s here when the equilibrium carries none; every consumer derives its view from `equil.kinetic`. """ function kinetic_forces_control(inputs::Dict{String,Any}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) - kf_ctrl = haskey(inputs, "KineticForces") ? - KineticForces.KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) : - KineticForces.KineticForcesControl() + kf_ctrl = + haskey(inputs, "KineticForces") ? + KineticForces.KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) : + KineticForces.KineticForcesControl() _needs_kinetic_profiles(inputs, ctrl) && equil.kinetic === nothing && error("this run needs kinetic profiles ([KineticForces], or kinetic_source = \"calculated\"): set kinetic_file in [Equilibrium]") return kf_ctrl @@ -1300,7 +1308,7 @@ function efc_couplings( ω_E = [kp.omegaE_spline(ψ) for ψ in ψk] ω_star_T = [-2π * kp.Ti_deriv(ψ) / (chrg * chi1) for ψ in ψk] ω_star_n = [ - (n=kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk + (n = kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk ] ω_φ = ω_E .+ ω_star_n .+ ω_star_T # Momentum weight n_i·dV/dψ with ⟨R²⟩ taken as R₀², inside the plasma only. diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index f8afac63c..fc7b01980 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -194,6 +194,20 @@ using TOML # A listed deck key must not still be a control field, or listing it would stop runs that set a live key. @test !any(k -> hasfield(GPEC.ForceFreeStates.ForceFreeStatesControl, Symbol(k)), GPEC._DEPRECATED_FFS_KEYS) @test !any(k -> hasfield(GPEC.KineticForces.KineticForcesControl, Symbol(k)), GPEC._DEPRECATED_KF_KEYS) + @test !any(k -> hasfield(GPEC.Runner.SLAYERControl, Symbol(k)), GPEC._DEPRECATED_SLAYER_KEYS) + # The kinetic-profile file moved to [Equilibrium]; the stop names the new home instead of "remove". + mktempdir() do deck_dir + template_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + foreach(name -> cp(joinpath(template_dir, name), joinpath(deck_dir, name)), readdir(template_dir)) + inputs = TOML.parsefile(joinpath(deck_dir, "gpec.toml")) + inputs["SLAYER"] = Dict{String,Any}("enabled" => false, "profile_file" => "kinetic.h5") + open(io -> TOML.print(io, inputs), joinpath(deck_dir, "gpec.toml"), "w") + @test_throws r"\[SLAYER\] sets `profile_file`, which is no longer a key; name the kinetic-profile file once, as kinetic_file" GPEC.main([deck_dir]) + delete!(inputs, "SLAYER") + inputs["KineticForces"] = Dict{String,Any}("kinetic_file" => "kinetic.dat") + open(io -> TOML.print(io, inputs), joinpath(deck_dir, "gpec.toml"), "w") + @test_throws r"\[KineticForces\] sets `kinetic_file`, which is no longer a key; name the kinetic-profile file once" GPEC.main([deck_dir]) + end pe_keywords = union((Base.kwarg_decl(m) for f in (GPEC.plasma_response, GPEC.perturbed_equilibrium) for m in methods(f))...) @test !any(k -> Symbol(k) in pe_keywords, GPEC._DEPRECATED_PE_KEYS) end