diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 3fda7d3ba..f7fb5ef39 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -333,6 +333,14 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect return identity.(species) # narrow to a concrete-eltype Vector{ResolvedNTVSpecies{...}} end +""" + has_ntv_fields(data::KineticProfileData) -> Bool + +Whether `data` carries the profiles the NTV view needs (`n_e`, `T_i`, `T_e` and `omega_E`), so that +`kinetic_profile_splines` and `resolve_ntv_species` can be built from it. +""" +has_ntv_fields(data::KineticProfileData) = all(!isnothing, (data.n_e, data.T_i, data.T_e, data.omega_E)) + """ kinetic_profile_splines(data::KineticProfileData; zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, diff --git a/src/ForceFreeStates/ForceBalance.jl b/src/ForceFreeStates/ForceBalance.jl index 29d70785d..f4db1a240 100644 --- a/src/ForceFreeStates/ForceBalance.jl +++ b/src/ForceFreeStates/ForceBalance.jl @@ -135,6 +135,8 @@ function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verb # A kinetic model (not supplied or synthetic matrices) is evaluated on the equilibrium's kinetic profiles. ctrl.kinetic === nothing || ctrl.kinetic isa Union{KineticMatrices,TestKineticMatrices} || equil.kinetic !== nothing || error("the KineticForces kinetic source needs kinetic profiles on the equilibrium: set kinetic_file or call attach_kinetic_profiles!") + ctrl.kinetic isa Union{Nothing,KineticMatrices,TestKineticMatrices} || Equilibrium.has_ntv_fields(equil.kinetic) || + error("the KineticForces kinetic source needs n_e, T_i, T_e and omega_E in the kinetic profiles; the kinetic file lacks at least one") ctrl.kinetic !== nothing && !(alg isa Forward) && error("$(nameof(typeof(alg)))() does not support kinetic runs; use Forward().") diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 8450df614..3be30b71b 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -166,6 +166,9 @@ function main_from_inputs( equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do Equilibrium.PlasmaEquilibrium(eq_config, additional_input) end + # The KineticForces stage needs the NTV view; stop before any solve rather than after the perturbed equilibrium. + deck.kinetic_forces !== nothing && equil.kinetic !== nothing && !Equilibrium.has_ntv_fields(equil.kinetic) && + error("[KineticForces] needs n_e, T_i, T_e and omega_E in the kinetic profiles; the kinetic file lacks at least one") if deck.force_balance !== nothing ffs = _timed_stage(runtimes, "force_free_states", "Force-Free States") do diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 010cc61ba..710aae592 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -170,11 +170,18 @@ function compute_calculated_kinetic_matrices( end # The force-balance solve's two-pass grid inputs: the kinetic profiles as the KineticForces model views them (its -# control when the source is one, else the default control) and the kinetic-resonance surfaces to pin. +# control when the source is one, else the default control) and the kinetic-resonance surfaces to pin. The model +# needs the NTV view, so a missing profile errors; otherwise profiles read for another consumer (e.g. SLAYER without +# omega_E) are left out and the grid follows the equilibrium alone. function ForceFreeStates.kinetic_grid_inputs(source, ctrl, equil, intr; rationals) kf_ctrl = source isa KineticForcesControl ? source : KineticForcesControl() - kinetic_profiles = ntv_profiles(kf_ctrl, equil) - pinned = resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, ntv_species(kf_ctrl, equil), equil, intr; rationals) + data = equil.kinetic + ntv = source isa KineticForcesControl || (data !== nothing && Equilibrium.has_ntv_fields(data)) + data !== nothing && !ntv && + @info "Two-pass psi grid: the kinetic profiles lack an NTV field (n_e, T_i, T_e, omega_E), so the grid follows the equilibrium alone" + kinetic_profiles = ntv ? ntv_profiles(kf_ctrl, equil) : nothing + species = ntv ? ntv_species(kf_ctrl, equil) : nothing + pinned = resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, species, equil, intr; rationals) return kinetic_profiles, pinned end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 7b56441b1..eaac9c7ff 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -164,6 +164,22 @@ using HDF5 ideal = GPEC.maybe_reform_equilibrium(pass1, sl_intr, sl_ctrl) pass1.kinetic = kinetic @test refined.rzphi_xs != ideal.rzphi_xs + # Profiles without omega_E (read for SLAYER, say) cannot form the NTV view: the grid follows the + # equilibrium alone, unless the run computes calculated kinetic matrices, which need it. + no_omega = GPEC.Equilibrium.KineticProfileData(; psi=kinetic.psi, n_i=kinetic.n_i, n_e=kinetic.n_e, + T_i=kinetic.T_i, T_e=kinetic.T_e, chi_e=kinetic.chi_e, chi_i=kinetic.chi_i, chi_phi=kinetic.chi_phi) + @test GPEC.Equilibrium.has_ntv_fields(kinetic) + @test !GPEC.Equilibrium.has_ntv_fields(no_omega) + pass1.kinetic = no_omega + no_omega_grid = GPEC.maybe_reform_equilibrium(pass1, sl_intr, sl_ctrl) + @test no_omega_grid.kinetic === no_omega + @test no_omega_grid.rzphi_xs == ideal.rzphi_xs + calc_deck = merge(sl_deck, Dict("ForceFreeStates" => merge(sl_deck["ForceFreeStates"], Dict("kinetic" => "KineticForces")))) + calc_ctrl = GPEC.force_balance_control(calc_deck, sl_dir) + @test_throws "missing required 'omega_E'" GPEC.maybe_reform_equilibrium(pass1, sl_intr, calc_ctrl) + # The KineticForces source stops up front, before any solve. + @test_throws "omega_E" solve(pass1, FFS.Forward(); nn=1, kinetic=GPEC.KineticForcesControl()) + pass1.kinetic = kinetic end @testset "galerkin solve returns a gal result with a unified delta_prime" begin