From c05b2babc43df7e5c1a87a098058051a8411ca57 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 29 Oct 2026 23:52:54 -0400 Subject: [PATCH 1/2] =?UTF-8?q?Equilibrium=20-=20BUGFIX=20-=20Build=20the?= =?UTF-8?q?=20two-pass=20=CF=88=20grid=20from=20kinetic=20profiles=20that?= =?UTF-8?q?=20have=20no=20omega=5FE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The grid refinement and the EL kinetic stage built the NTV view of the attached profiles whenever any were attached, and that view requires omega_E, so a deck whose kinetic file has no omega_E (read for SLAYER) stopped in the equilibrium stage. The view is now built only when a calculated kinetic source needs it, where a missing field still errors, or when the profiles carry the NTV fields; otherwise the grid follows the equilibrium alone, as before the profiles moved onto the equilibrium. Closes #614 Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/KineticProfiles.jl | 8 ++++++++ src/GeneralizedPerturbedEquilibrium.jl | 19 +++++++++++++++---- test/runtests_solve_api.jl | 13 +++++++++++++ 3 files changed, 36 insertions(+), 4 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index e30b3af83..17d044bcc 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -347,6 +347,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/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index c791fd610..f1a0f2fdb 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -479,6 +479,16 @@ 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") +# The NTV view and species of the attached kinetic profiles for the EL stage and the grid. A calculated +# kinetic source needs them, so a missing profile errors; otherwise they are built only when the profiles +# carry the NTV fields, and profiles read for another consumer (e.g. SLAYER without omega_E) are left alone. +function _ntv_views(kf_ctrl::KineticForces.KineticForcesControl, equil::Equilibrium.PlasmaEquilibrium, ctrl::ForceFreeStatesControl) + needed = ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" + data = equil.kinetic + (data === nothing || !(needed || Equilibrium.has_ntv_fields(data))) && return (nothing, nothing) + return KineticForces.ntv_profiles(kf_ctrl, equil), KineticForces.ntv_species(kf_ctrl, equil) +end + """ maybe_reform_equilibrium(equil, intr, ctrl; kf_ctrl=KineticForcesControl()) -> equil @@ -495,12 +505,14 @@ function maybe_reform_equilibrium( ) eq_config = equil.config Equilibrium.wants_two_pass(eq_config) || return equil - kinetic_profiles = KineticForces.ntv_profiles(kf_ctrl, equil) + kinetic_profiles, species = _ntv_views(kf_ctrl, equil, ctrl) + (equil.kinetic === nothing || kinetic_profiles !== nothing) || + @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" 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 = KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, KineticForces.ntv_species(kf_ctrl, equil), equil, intr; + pinned = KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, species, 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. @@ -705,8 +717,7 @@ 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) + kinetic_profiles, species = _ntv_views(kf_ctrl, equil, ctrl) # 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, diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 0124fe2e2..7b6fff630 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -149,6 +149,19 @@ 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(no_omega) && GPEC.Equilibrium.has_ntv_fields(kinetic) + 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["ForceFreeStates"], Dict("kinetic_factor" => 1.0, "kinetic_source" => "calculated")) + calc_ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in calc_deck)...) + @test_throws "missing required 'omega_E'" GPEC.maybe_reform_equilibrium(pass1, sl_intr, calc_ctrl) + pass1.kinetic = kinetic end @testset "galerkin solve returns a gal result with a unified delta_prime" begin From fc9d9bb74c8123775b6ad043c82ae03ff7994dcf Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 30 Oct 2026 12:51:30 -0400 Subject: [PATCH 2/2] Equilibrium - BUGFIX - Stop NTV runs up front on profiles without omega_E and report the suppression they lose A [KineticForces] deck with omega_E-less profiles now errors in kinetic_forces_control rather than after the solve, a fixed kinetic source logs that the near-axis suppression is off, and the calculated-source test lives in one helper. Review changes. Co-Authored-By: Claude Opus 5.5 --- src/GeneralizedPerturbedEquilibrium.jl | 29 +++++++++++++++++--------- test/runtests_solve_api.jl | 7 ++++++- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index f1a0f2fdb..dd8c721da 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -464,28 +464,35 @@ end 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`. +here when the equilibrium carries none or they lack an NTV field; 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() - _needs_kinetic_profiles(inputs, ctrl) && equil.kinetic === nothing && - error("this run needs kinetic profiles ([KineticForces], or kinetic_source = \"calculated\"): set kinetic_file in [Equilibrium]") + if _needs_kinetic_profiles(inputs, ctrl) + equil.kinetic === nothing && + error("this run needs kinetic profiles ([KineticForces], or kinetic_source = \"calculated\"): set kinetic_file in [Equilibrium]") + Equilibrium.has_ntv_fields(equil.kinetic) || + error("this run needs n_e, T_i, T_e and omega_E in the kinetic profiles ([KineticForces], or kinetic_source = \"calculated\"); the kinetic file lacks at least one") + end return kf_ctrl end +# Whether the EL stage computes its kinetic matrices from the profiles. +_calculated_kinetic(ctrl) = ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" + # 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") +_needs_kinetic_profiles(inputs, ctrl) = haskey(inputs, "KineticForces") || _calculated_kinetic(ctrl) # The NTV view and species of the attached kinetic profiles for the EL stage and the grid. A calculated # kinetic source needs them, so a missing profile errors; otherwise they are built only when the profiles # carry the NTV fields, and profiles read for another consumer (e.g. SLAYER without omega_E) are left alone. function _ntv_views(kf_ctrl::KineticForces.KineticForcesControl, equil::Equilibrium.PlasmaEquilibrium, ctrl::ForceFreeStatesControl) - needed = ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" data = equil.kinetic - (data === nothing || !(needed || Equilibrium.has_ntv_fields(data))) && return (nothing, nothing) + (data === nothing || !(_calculated_kinetic(ctrl) || Equilibrium.has_ntv_fields(data))) && return (nothing, nothing) return KineticForces.ntv_profiles(kf_ctrl, equil), KineticForces.ntv_species(kf_ctrl, equil) end @@ -493,8 +500,8 @@ end maybe_reform_equilibrium(equil, intr, ctrl; kf_ctrl=KineticForcesControl()) -> equil 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 with `Equilibrium.reform` from the source the equilibrium remembers. Returns `equil` untouched +view of its attached kinetic profiles when they carry the NTV fields), pin knots on rational surfaces, +and re-form on the refined grid with `Equilibrium.reform` from the source the equilibrium remembers. Returns `equil` untouched when its configuration wants a single pass. """ function maybe_reform_equilibrium( @@ -506,7 +513,7 @@ function maybe_reform_equilibrium( eq_config = equil.config Equilibrium.wants_two_pass(eq_config) || return equil kinetic_profiles, species = _ntv_views(kf_ctrl, equil, ctrl) - (equil.kinetic === nothing || kinetic_profiles !== nothing) || + equil.kinetic !== nothing && kinetic_profiles === nothing && @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" mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) @@ -718,6 +725,8 @@ function prepare_force_free_states!( # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). kinetic_profiles, species = _ntv_views(kf_ctrl, equil, ctrl) + equil.kinetic !== nothing && kinetic_profiles === nothing && kf_ctrl.axis_validity_suppression && + @info "Kinetic matrices: the kinetic profiles lack an NTV field (n_e, T_i, T_e, omega_E), so the near-axis validity suppression is off" # 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, @@ -932,7 +941,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat ForceFreeStates._apply_alg!(ctrl_kwargs, alg) ctrl = ForceFreeStatesControl(; ctrl_kwargs...) - ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" && equil.kinetic === nothing && + _calculated_kinetic(ctrl) && equil.kinetic === nothing && error("kinetic_source = \"calculated\" needs kinetic profiles on the equilibrium: set kinetic_file or call attach_kinetic_profiles!") ctrl.kinetic_factor > 0 && !(alg isa ForceFreeStates.Forward) && error("$(nameof(typeof(alg)))() does not support kinetic runs (kinetic_factor > 0); use Forward().") diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 7b6fff630..8cff6149a 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -153,7 +153,8 @@ using HDF5 # 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(no_omega) && GPEC.Equilibrium.has_ntv_fields(kinetic) + @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 @@ -161,6 +162,10 @@ using HDF5 calc_deck = merge(sl_deck["ForceFreeStates"], Dict("kinetic_factor" => 1.0, "kinetic_source" => "calculated")) calc_ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in calc_deck)...) @test_throws "missing required 'omega_E'" GPEC.maybe_reform_equilibrium(pass1, sl_intr, calc_ctrl) + # An NTV consumer stops up front, before any solve; a deck without one runs. + @test_throws "omega_E" GPEC.kinetic_forces_control(Dict{String,Any}("KineticForces" => Dict{String,Any}()), sl_ctrl, pass1) + @test_throws "omega_E" GPEC.kinetic_forces_control(Dict{String,Any}(), calc_ctrl, pass1) + @test GPEC.kinetic_forces_control(Dict{String,Any}(), sl_ctrl, pass1) isa GPEC.KineticForces.KineticForcesControl pass1.kinetic = kinetic end