Skip to content
Merged
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
123 changes: 52 additions & 71 deletions src/GeneralizedPerturbedEquilibrium.jl
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@ function main_from_inputs(
metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles;
species=kf_species,
overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path))
ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric; runtimes=runtimes)
ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric)

if ctrl.write_outputs_to_HDF5
write_outputs_to_HDF5(
Expand Down Expand Up @@ -286,56 +286,44 @@ function main_from_inputs(
push!(runtimes, "force_free_states" => ffs_dt)
@info "Force-Free States completed in $(@sprintf("%.3f", ffs_dt)) s"

# Early exit if user only requested force-free states (SLAYER still runs).
if ctrl.force_termination
slayer_result = run_slayer_stage(ffs_result, inputs, nothing; runtimes=runtimes)
# A non-nothing result means the Tearing/ group was appended.
slayer_result === nothing || (written_h5 = ctrl.HDF5_filename)
total_dt = time() - total_start
push!(runtimes, "total" => total_dt)
if written_h5 !== nothing
_write_runtimes!(joinpath(intr.dir_path, written_h5), runtimes)
end
@info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
return (; ffs=ffs_result, pe=nothing, slayer=slayer_result, coil_sensitivities=nothing, monte_carlo=nothing, locking_risk=nothing, efc_couplings=nothing)
end
# A force_termination run stops after ForceFreeStates; only SLAYER still runs off the result.
run_response = !ctrl.force_termination

# PE (and SLAYER, which appends after it) writes to PE's own output filename when set,
# falling back to the ForceFreeStates file.
pe_file = if "PerturbedEquilibrium" in keys(inputs)
pe_file = if run_response && "PerturbedEquilibrium" in keys(inputs)
pe_out = get(inputs["PerturbedEquilibrium"], "output_filename", "")
isempty(pe_out) ? ctrl.HDF5_filename : pe_out
else
ctrl.HDF5_filename
end

pe_start = time()
pe_state = run_perturbed_equilibrium(ffs_result, inputs, preloaded_coil_sets; runtimes=runtimes)
# Record only when the stage ran; an absent [PerturbedEquilibrium] section would otherwise
# stamp a ~0 s entry for work that never happened.
if "PerturbedEquilibrium" in keys(inputs)
push!(runtimes, "perturbed_equilibrium" => time() - pe_start)
pe_state = nothing
if run_response && "PerturbedEquilibrium" in keys(inputs)
pe_state = _timed_stage(runtimes, "perturbed_equilibrium", "Perturbed Equilibrium") do
run_perturbed_equilibrium(ffs_result, inputs, preloaded_coil_sets)
end
# Mirrors the write gate inside `perturbed_equilibrium` (write flag defaults to true).
if get(inputs["PerturbedEquilibrium"], "write_outputs_to_HDF5", true)
written_h5 = pe_file
end
end

kf_start = time()
run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species)
if "KineticForces" in keys(inputs)
push!(runtimes, "kinetic_forces" => time() - kf_start)
if run_response && "KineticForces" in keys(inputs)
_timed_stage(runtimes, "kinetic_forces", "KineticForces") do
run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species)
end
# Mirrors the write gate inside `run_kinetic_forces` (needs a PE state to contract against).
if pe_state !== nothing && kf_ctrl.write_outputs_to_HDF5
written_h5 = kf_ctrl.HDF5_filename
end
end

ef_start = time()
error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles)
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)
coil_sensitivities, monte_carlo, locking_risk, efc_couplings = nothing, nothing, nothing, nothing
if run_response && "ErrorFields" in keys(inputs)
coil_sensitivities, monte_carlo, locking_risk, efc_couplings = _timed_stage(runtimes, "error_fields", "ErrorFields") do
run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles)
end
# Mirrors the write gate inside `run_error_fields` (write flag defaults to true, filename
# falls back to the ForceFreeStates file since the result carries this same control).
ef_raw = inputs["ErrorFields"]
Expand All @@ -345,9 +333,15 @@ function main_from_inputs(
end
end

slayer_result = run_slayer_stage(ffs_result, inputs, pe_file; runtimes=runtimes)
tearing_start = time()
slayer_result = run_slayer_stage(ffs_result, inputs, pe_file)
tearing_dt = time() - tearing_start
# A non-nothing result means the Tearing/ group was appended.
slayer_result === nothing || (written_h5 = pe_file)
if slayer_result !== nothing
written_h5 = pe_file
push!(runtimes, "tearing" => tearing_dt)
@info "SLAYER completed in $(@sprintf("%.3f", tearing_dt)) s"
end

# ----------------------------------------------------------------
# Done
Expand All @@ -372,6 +366,22 @@ Compact label for the resolved toroidal mode range, `"1"` for a single `n` and `
"""
_mode_range_label(intr::ForceFreeStatesInternal) = intr.npert == 1 ? "$(intr.nlow)" : "$(intr.nlow):$(intr.nhigh)"

"""
_timed_stage(f, runtimes, key, label) -> f()

Run the pipeline stage `f` under a `label` section banner, log its wall-clock time and record
it in `runtimes` as `key => seconds` for the `Info/Runtimes` record.
"""
function _timed_stage(f, runtimes::Vector{Pair{String,Float64}}, key::String, label::String)
@info "\n $label\n$_SECTION"
start = time()
out = f()
dt = time() - start
push!(runtimes, key => dt)
@info "$label completed in $(@sprintf("%.3f", dt)) s"
return out
end

"""
resolve_mode_space!(intr, ctrl) -> intr

Expand Down Expand Up @@ -741,16 +751,14 @@ end

Run the formalism selected by `ctrl.integrator` — the standalone Galerkin solve, or the
Euler-Lagrange sweep with its free-boundary energies and Δ′ BVP — and publish its products as a
`ForceFreeStatesResult`. A `runtimes` collector, when given, receives the Galerkin solve's
wall-clock seconds as a `"galerkin" => dt` pair for the `Info/Runtimes` record.
`ForceFreeStatesResult`.
"""
function run_force_free_states(
ctrl::ForceFreeStatesControl,
equil::Equilibrium.PlasmaEquilibrium,
mats,
intr::ForceFreeStatesInternal,
metric;
runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing
metric
)
nstring = _mode_range_label(intr)

Expand All @@ -768,7 +776,6 @@ function run_force_free_states(
wv = ctrl.vac_flag ? first(ForceFreeStates.compute_scaled_wv(ctrl, equil, intr)) : nothing
gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv)
gal_dt = time() - gal_start
runtimes === nothing || push!(runtimes, "galerkin" => gal_dt)
@info "Galerkin solve completed in $(@sprintf("%.3f", gal_dt)) s"
else
# Integrate Euler-Lagrange Equation
Expand Down Expand Up @@ -976,15 +983,8 @@ outputs. Returns `nothing` when the deck carries no `[PerturbedEquilibrium]` sec
function run_perturbed_equilibrium(
result::ForceFreeStatesResult,
inputs::Dict{String,Any},
preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}};
runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing
preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}
)
# ----------------------------------------------------------------
# Perturbed Equilibrium
# ----------------------------------------------------------------
@info "\n Perturbed Equilibrium\n$_SECTION"
pe_start = time()

# Check for PerturbedEquilibrium section and run if present
pe_state = nothing
if "PerturbedEquilibrium" in keys(inputs)
Expand All @@ -993,17 +993,15 @@ function run_perturbed_equilibrium(
# The deck's forcing block is an unscaled RMPField, so the TOML path and the
# scripting API share one stage.
pe_state = perturbed_equilibrium(result, ForcingTerms.RMPField(ft_ctrl);
coil_sets=preloaded_coil_sets, runtimes=runtimes,
coil_sets=preloaded_coil_sets,
(Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])...)
end

@info "Perturbed Equilibrium completed in $(@sprintf("%.3f", time() - pe_start)) s"

return pe_state
end

"""
perturbed_equilibrium(ffs, rmp; coil_sets=nothing, runtimes=nothing, kwargs...) -> PerturbedEquilibriumState
perturbed_equilibrium(ffs, rmp; coil_sets=nothing, kwargs...) -> PerturbedEquilibriumState

Compute the plasma response to the external field `rmp` on top of a force-free-states solve
`ffs`, and write the perturbed-equilibrium outputs. `rmp` is an [`RMPField`](@ref); keyword
Expand All @@ -1015,7 +1013,7 @@ Galerkin and `:matched` results publish no free-boundary δW yet and are refused
response or the singular coupling is requested.

`coil_sets` injects already-built coil geometry (the gpec.h5 replay path), bypassing the
geometry read. A `runtimes` collector, when given, receives the forcing-mode materialization's wall-clock seconds as `"forcing_terms" => dt`.
geometry read.

```julia
pe = perturbed_equilibrium(ffs, RMPField("forcing.dat"))
Expand All @@ -1025,7 +1023,6 @@ function perturbed_equilibrium(
ffs::ForceFreeStatesResult,
rmp::ForcingTerms.RMPField;
coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing,
runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing,
kwargs...
)
ctrl = ffs.control
Expand All @@ -1041,7 +1038,7 @@ function perturbed_equilibrium(
pe_intr.coil_sets = copy(coil_sets)
end

pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(ffs, rmp, pe_ctrl, pe_intr; runtimes=runtimes)
pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(ffs, rmp, pe_ctrl, pe_intr)

# Write perturbed equilibrium outputs to same HDF5 file
if pe_ctrl.write_outputs_to_HDF5
Expand Down Expand Up @@ -1081,9 +1078,6 @@ function run_kinetic_forces(
# ----------------------------------------------------------------
("KineticForces" in keys(inputs)) || return nothing

@info "\n KineticForces\n$_SECTION"
kf_start = time()

# One boundary for the whole run, shared with the EL kinetic matrices and the Validity output.
axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, result.equil,
Float64[sng.psifac for sng in result.surfaces])
Expand Down Expand Up @@ -1136,8 +1130,6 @@ function run_kinetic_forces(
end
end

@info "KineticForces completed in $(@sprintf("%.3f", time() - kf_start)) s"

return nothing
end

Expand Down Expand Up @@ -1182,9 +1174,6 @@ function run_error_fields(
)
("ErrorFields" in keys(inputs)) || return nothing

@info "\n ErrorFields\n$_SECTION"
ef_start = time()

# [ErrorFields.MonteCarlo], [ErrorFields.Risk] and [ErrorFields.scenario] are nested tables,
# excluded from the control-struct splat.
ef_raw = inputs["ErrorFields"]
Expand Down Expand Up @@ -1277,8 +1266,6 @@ function run_error_fields(
@info "Results written to $output_file"
end

@info "ErrorFields completed in $(@sprintf("%.3f", time() - ef_start)) s"

return sens, monte_carlo, risk, couplings
end

Expand Down Expand Up @@ -1434,22 +1421,16 @@ end

Run the SLAYER tearing-mode analysis off the force-free-states `result`, appending its group to
`pe_file` (or the force-free-states output when PE did not run). Needs only the result, so it
runs in both the `force_termination = true` path and the full pipeline. A `runtimes` collector,
when given, receives the solver's wall-clock seconds as a `"tearing" => dt` pair for the
`Info/Runtimes` record. Errors propagate; every earlier stage has already written its output.
runs in both the `force_termination = true` path and the full pipeline. Errors propagate; every
earlier stage has already written its output.
"""
function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any}, pe_file::Union{String,Nothing};
runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing)
function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any}, pe_file::Union{String,Nothing})
("SLAYER" in keys(inputs)) || return nothing
slayer_ctrl = Runner.slayer_control_from_toml(inputs["SLAYER"])
slayer_ctrl.enabled || return nothing
@info "\n SLAYER\n$_SECTION"
slayer_start = time()
model = Runner._build_inner_model(slayer_ctrl.inner_model)
slayer_result = solve(Runner.TearingProblem(result, slayer_ctrl), model)
slayer_dt = time() - slayer_start
runtimes === nothing || push!(runtimes, "tearing" => slayer_dt)
@info "SLAYER completed in $(@sprintf("%.3f", slayer_dt)) s"
h5_filename = pe_file === nothing ? result.control.HDF5_filename : pe_file
h5_path = joinpath(result.dir_path, h5_filename)
# Append the Tearing/ group; create the file if no prior stage wrote
Expand Down
8 changes: 1 addition & 7 deletions src/HDF5Schema.jl
Original file line number Diff line number Diff line change
Expand Up @@ -281,18 +281,12 @@ const MAIN_H5_ANNOTATIONS = [
# Stages that did not run leave their path absent and are skipped by `annotate!`.
const RUNTIME_H5_ANNOTATIONS = [
"Info/Runtimes/equilibrium" => (; long_name="wall-clock time of the Equilibrium construction stage", units="s"),
"Info/Runtimes/galerkin" => (; long_name="wall-clock time of the Galerkin outer-region solve (measured inside force_free_states, not additional to it)", units="s"),
"Info/Runtimes/force_free_states" => (; long_name="wall-clock time of the ForceFreeStates stability stage", units="s"),
"Info/Runtimes/tearing" => (; long_name="wall-clock time of the Tearing stage (SLAYER inner-layer solve)", units="s"),
"Info/Runtimes/forcing_terms" => (;
long_name="wall-clock time of the ForcingTerms forcing-mode materialization, incl. coil Biot-Savart onto the plasma surface (measured inside perturbed_equilibrium, not additional to it)",
units="s"
),
"Info/Runtimes/perturbed_equilibrium" => (; long_name="wall-clock time of the PerturbedEquilibrium stage", units="s"),
"Info/Runtimes/kinetic_forces" => (; long_name="wall-clock time of the KineticForces (NTV) stage", units="s"),
"Info/Runtimes/error_fields" => (; long_name="wall-clock time of the ErrorFields stage (coil sensitivities, Monte Carlo, locking risk)", units="s"),
"Info/Runtimes/total" =>
(; long_name="wall-clock time of the full GPEC run (the per-stage values nest rather than partition it, so they do not sum to this)", units="s")
"Info/Runtimes/total" => (; long_name="wall-clock time of the full GPEC run, including the work between stages", units="s")
]

# Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants.
Expand Down
9 changes: 2 additions & 7 deletions src/PerturbedEquilibrium/PerturbedEquilibrium.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ export ResonantCoupling, DominantCoupling, dominant_coupling, rootarea_field, co
export TorqueResponse, torque_response_matrices, torque_profile

"""
compute_perturbed_equilibrium(ffs, forcing, ctrl, intr; runtimes=nothing)::PerturbedEquilibriumState
compute_perturbed_equilibrium(ffs, forcing, ctrl, intr)::PerturbedEquilibriumState

Main entry point for perturbed equilibrium calculations.

Expand All @@ -59,8 +59,6 @@ no free-boundary δW yet, so asking for the response or the coupling on one is a
- `forcing`: the external-field description — a `ForcingTermsControl` (TOML path) or any [`ForcingTerms.RMPField`](@ref)
- `ctrl`: Control parameters from [PerturbedEquilibrium] section
- `intr`: Internal state variables
- `runtimes`: optional collector; receives `"forcing_terms" => dt` for the forcing-mode
materialization (coil Biot-Savart onto the plasma surface), when that step runs

## Returns

Expand All @@ -70,8 +68,7 @@ function compute_perturbed_equilibrium(
ffs::ForceFreeStatesResult,
forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField},
ctrl::PerturbedEquilibriumControl,
intr::PerturbedEquilibriumInternal;
runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing
intr::PerturbedEquilibriumInternal
)::PerturbedEquilibriumState

state = PerturbedEquilibriumState()
Expand All @@ -97,10 +94,8 @@ function compute_perturbed_equilibrium(
# Unit-norm modes an internal caller already placed on `intr` skip materialization, so they
# are never re-converted or re-weighted.
if isempty(intr.forcing_modes)
ft_start = time()
modes, coil_sets = materialize_forcing_modes(ffs, forcing;
dir_path=intr.dir_path, preloaded_coil_sets=intr.coil_sets, verbose=ctrl.verbose)
runtimes === nothing || push!(runtimes, "forcing_terms" => time() - ft_start)
intr.forcing_modes = modes
isempty(coil_sets) || (intr.coil_sets = coil_sets)
end
Expand Down
4 changes: 3 additions & 1 deletion test/runtests_error_fields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -99,10 +99,12 @@ include("h5_metadata_check.jl")
h5open(h5path, "r") do h5
@test haskey(h5, "Info/Runtimes/error_fields")
@test read(h5["Info/Runtimes/error_fields"]) > 0
@test haskey(h5, "Info/Runtimes/forcing_terms")
@test haskey(h5, "Info/Runtimes/perturbed_equilibrium")
@test haskey(h5, "Info/Runtimes/kinetic_forces")
@test haskey(h5, "Info/Runtimes/total")
# Timed inside a stage, so no longer recorded.
@test !haskey(h5, "Info/Runtimes/forcing_terms")
@test !haskey(h5, "Info/Runtimes/galerkin")
end

@test sens.coil_names == ["hoop_tilted", "hoop_axi", "hoop_efc"]
Expand Down
Loading