Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions src/Equilibrium/KineticProfiles.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
42 changes: 31 additions & 11 deletions src/GeneralizedPerturbedEquilibrium.jl
Original file line number Diff line number Diff line change
Expand Up @@ -464,27 +464,44 @@ 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)
data = equil.kinetic
(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

"""
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(
Expand All @@ -495,12 +512,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"

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: the (a || b) || @info guard is a double negative; equil.kinetic !== nothing && kinetic_profiles === nothing && @info ... reads directly.


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.
Expand Down Expand Up @@ -705,8 +724,9 @@ 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)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With kinetic_source = "fixed" and profiles lacking omega_E, this now gives axis_psi_c = 0, so near-axis suppression is silently off, while the same deck with omega_E gets a nonzero boundary. Not a regression (it errored before), but worth an @info here or a line in the PR description so it isn't silent.

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,
Expand Down Expand Up @@ -921,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().")
Expand Down
18 changes: 18 additions & 0 deletions test/runtests_solve_api.jl
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,24 @@ 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["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

@testset "galerkin solve returns a gal result with a unified delta_prime" begin
Expand Down
Loading