diff --git a/docs/src/api.md b/docs/src/api.md index 4cb8938ddc..5510b21407 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -10,7 +10,7 @@ eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada", kinetic_file="kineti ffs = solve(eq, Forward(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) # [ForceFreeStates]: Forward, no wall, for the response below ls = solve(LocalStabilityProblem(eq; alpha_boundary_scan=true), Ballooning()) # [LocalStability] resp = plasma_response(ffs; compute_singular_coupling=true) # [PerturbedEquilibrium] -pe = perturbed_equilibrium(resp, RMPField("coils.dat")) # with the [ForcingTerms] drive +pe = perturbed_equilibrium(resp, ForcingField("coils.dat")) # with the [ForcingTerms] drive kf = kinetic_torque(pe, KineticForcesControl()) # [KineticForces] ef = error_fields(resp, pe.coil_sets, ForcingTerms.CoilConfig(ft_ctrl); # [ErrorFields]; ft_ctrl: the [ForcingTerms] control tolerances=ErrorFields.read_tolerance_toml("tolerances.toml")) @@ -29,8 +29,8 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear)) | ForceFreeStates | `[ForceFreeStates]` + `[Forward]`, `[Riccati]` or `[Galerkin]` | `solve(ForceBalanceProblem(eq; ...), Forward() \| Riccati() \| Galerkin())` | `ForceFreeStatesResult` | | LocalStability | `[LocalStability]` + `[Ballooning]` | `solve(LocalStabilityProblem(eq; ...), Ballooning(; ...))` | `LocalStabilityResult` | | Match | `[Match]` + `[GGJ]` | `solve(MatchProblem(ffs; ...), GGJModel())` | `ForceFreeStatesResult` with `closure == :matched` | -| Drive | `[ForcingTerms]` | `RMPField(ft_ctrl)` (a `ForcingTermsControl`), `RMPField(file)`, `RMPField(modes)` | `RMPField` | -| PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, rmp; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` | +| Drive | `[ForcingTerms]` | `ForcingField(ft_ctrl)` (a `ForcingTermsControl`), `ForcingField(file)`, `ForcingField(modes)` | `ForcingField` | +| PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, field; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` | | KineticForces | `[KineticForces]` | `kinetic_torque(pe, KineticForcesControl(; ...))` | `KineticForcesResult` | | ErrorFields | `[ErrorFields]` | `error_fields(resp, pe.coil_sets, CoilConfig(...); ...)` | `ErrorFieldsResult` | | Tearing | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `solve(TearingProblem(ffs; ...), SLAYERModel() \| GGJModel(); config=TearingConfig(...))` | `TearingResult` | @@ -107,16 +107,16 @@ The `kinetic_file` (and, for HDF5, `kinetic_group`) keyword reads the kinetic pr ## Combining forcing sources -`RMPField`s form a vector space: `+`, `-` and scalar `*` build lazy linear combinations that compute nothing until `perturbed_equilibrium` materializes them, and because the response is linear in the forcing, driving with a combination equals combining the materialized fields; a complex scalar phase-rotates a source. +`ForcingField`s form a vector space: `+`, `-` and scalar `*` build lazy linear combinations that compute nothing until `perturbed_equilibrium` materializes them, and because the response is linear in the forcing, driving with a combination equals combining the materialized fields; a complex scalar phase-rotates a source. ```julia -nominal = RMPField("nominal_efc.dat") -weld_field = RMPField("weld_fields.h5") +nominal = ForcingField("nominal_efc.dat") +weld_field = ForcingField("weld_fields.h5") pe = perturbed_equilibrium(resp, 2.0 * nominal + 0.5 * weld_field) ``` This is the substrate for error-field workflows: build per-unit sources (a coil set per ampere, a displacement field per millimetre), then assemble physical cases by weighting and summing, with one `plasma_response` serving every combination. -`RMPField(modes)` wraps a control-surface spectrum computed elsewhere, in the unit-norm ``\Phi_x`` convention, as a leaf. +`ForcingField(modes)` wraps a control-surface spectrum computed elsewhere, in the unit-norm ``\Phi_x`` convention, as a leaf. The weights are linear-combination coefficients, not physical amplitudes: a coil set with a failed conductor is `nominal - failed_coil`, not `0.9 * nominal`. ## Choosing an integrator diff --git a/docs/src/forcing_terms.md b/docs/src/forcing_terms.md index 8ad9ae5a0f..2cd2ea5915 100644 --- a/docs/src/forcing_terms.md +++ b/docs/src/forcing_terms.md @@ -1,7 +1,7 @@ # Forcing Terms `ForcingTerms` specifies the external magnetic perturbation that drives the perturbed equilibrium: the spectrum of the normal field on the control surface at `psilim`, given either as a file of (m, n) amplitudes or computed by Biot–Savart from coil geometry. -The `[ForcingTerms]` section is read when `[PerturbedEquilibrium]` runs; its result is the `RMPField` that `perturbed_equilibrium` applies to the plasma response. +The `[ForcingTerms]` section is read when `[PerturbedEquilibrium]` runs; its result is the `ForcingField` that `perturbed_equilibrium` applies to the plasma response. ## Normalization diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index b5b9d27a93..2c7ddfddf0 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/docs/src/perturbed_equilibrium.md @@ -7,7 +7,7 @@ It is the stage that turns a stability calculation into the quantities of error- The force-free states of [Stability](stability.md) give the displacement inside the plasma for any boundary value, and the [Vacuum](vacuum.md) response gives the field outside it. On the control surface at `psilim` the two are joined by requiring continuity of the normal field, which fixes the plasma response to an applied spectrum uniquely ([Park et al. (2007a)](citations.md#PerturbedEquilibrium), [Park et al. (2009)](citations.md#PerturbedEquilibrium)). -`plasma_response(ffs)` builds the operators of that response once, with no reference to a drive; `perturbed_equilibrium(resp, rmp)` applies them to one external field from [Forcing Terms](forcing_terms.md), so a scan over fields reuses one response. +`plasma_response(ffs)` builds the operators of that response once, with no reference to a drive; `perturbed_equilibrium(resp, field)` applies them to one external field from [Forcing Terms](forcing_terms.md), so a scan over fields reuses one response. `compute_perturbed_equilibrium` is the lower-level form of the same call. The resonant part of the response is reported at every rational surface as the singular coupling: the pitch-resonant field and current, the island half-width and the Chirikov overlap parameter. @@ -58,7 +58,7 @@ A kinetic forward solve (`integrator = "Forward"`, `kinetic = "KineticForces"`) With `compute_torque_response = true`, `plasma_response` refers that quadratic form back to the applied root-area-weighted field ``\tilde b`` and stores it at every node as the torque response matrix ``T_{xe}(\psi)``, cumulative in ψ, under `PerturbedEquilibrium/TorqueResponse/`; it is the `T_xe` of the Fortran GPEC, in the basis of the stored response matrices, in N·m/T². ``T_{xe}`` is built on the ``+n`` harmonic alone and is therefore twice the physical value, so the cumulative complex torque of an applied spectrum is ``T(\psi) = \tilde b^\dagger T_{xe}(\psi)\, \tilde b / 2`` (`torque_profile` applies the 1/2), with ``\operatorname{Re} T`` the toroidal torque inside ψ, equal to `Energies/toroidal_torque` at the control surface, and ``\operatorname{Im} T = 2n\,\delta W``. Only the Hermitian part ``(T_{xe} + T_{xe}^\dagger)/2`` contributes to the torque; its eigenvectors are the applied spectra of extremal torque, and the anti-Hermitian part carries the energy. -Each `perturbed_equilibrium(resp, rmp)` call contracts ``T_{xe}`` with its field into a `TorqueResponse`: the applied torque ``T_\mathrm{applied}`` and the source form ``T_\mathrm{source} = M^\dagger T M`` with one column of ``M`` per forcing source, so that driving the sources at scale factors ``s`` gives ``s^\dagger T_\mathrm{source}\, s / 2``. +Each `perturbed_equilibrium(resp, field)` call contracts ``T_{xe}`` with its field into a `TorqueResponse`: the applied torque ``T_\mathrm{applied}`` and the source form ``T_\mathrm{source} = M^\dagger T M`` with one column of ``M`` per forcing source, so that driving the sources at scale factors ``s`` gives ``s^\dagger T_\mathrm{source}\, s / 2``. ## Governing equations @@ -114,11 +114,11 @@ In a script the same two steps are using GeneralizedPerturbedEquilibrium resp = plasma_response(ffs; compute_singular_coupling=true) -pe = perturbed_equilibrium(resp, rmp; regularization_width=0.05) # rmp: an RMPField, see the Scripting API +pe = perturbed_equilibrium(resp, field; regularization_width=0.05) # field: a ForcingField, see the Scripting API save_to_hdf5("gpec.h5", pe) ``` -with the `RMPField` built from coils or a file as on the [Scripting API](api.md) page; computing the response does not write a file. +with the `ForcingField` built from coils or a file as on the [Scripting API](api.md) page; computing the response does not write a file. ## Outputs diff --git a/src/Analysis/CoilForcing.jl b/src/Analysis/CoilForcing.jl index 353a562782..11bcbd4934 100644 --- a/src/Analysis/CoilForcing.jl +++ b/src/Analysis/CoilForcing.jl @@ -232,18 +232,18 @@ function plot_coil_geometry_3d(pe::PerturbedEquilibriumState; save_path=nothing) return plot_coil_geometry_3d(pe.coil_sets; equil=pe.response.ffs.equil, save_path) end -# The coil sets of every coil term of `rmp`, each set's currents scaled by its term's weight (by the weight's magnitude +# The coil sets of every coil term of `field`, each set's currents scaled by its term's weight (by the weight's magnitude # when it is complex, which shifts the toroidal phase), with the number of terms that have no geometry. -function _rmp_coil_sets(rmp::FT.RMPField; n=1) +function _field_coil_sets(field::FT.ForcingField; n=1) sets, skipped, complex_weight = FT.CoilSet[], 0, false - for t in FT._rmp_terms(rmp) - if !(t isa FT.RMPSource && t.ctrl.forcing_data_format == "coil") + for t in FT._field_terms(field) + if !(t isa FT.ForcingCoils) skipped += 1 continue end complex_weight |= !isreal(t.scale) w = isreal(t.scale) ? real(t.scale) : abs(t.scale) - for cs in FT.load_coil_sets(FT.CoilConfig(t.ctrl), n) + for cs in FT.load_coil_sets(t.config, n) push!(sets, FT.CoilSet(cs.name, cs.ncoil, cs.s, cs.nw, cs.nsec, cs.x, cs.y, cs.z, w .* cs.currents)) end end @@ -251,15 +251,15 @@ function _rmp_coil_sets(rmp::FT.RMPField; n=1) end """ - plot_sources(rmp; n=1, color=:current, save_path=nothing) + plot_sources(field; n=1, color=:current, save_path=nothing) -Plot the coils of an `RMPField` in 3D at equal data aspect, with no plasma: every coil term of a sum together, each +Plot the coils of a `ForcingField` in 3D at equal data aspect, with no plasma: every coil term of a sum together, each term's currents scaled by its weight. Mode-table and spectrum terms have no geometry and are counted in the title. Use GLMakie to rotate it. ### Arguments - - `rmp`: An `RMPField` with at least one coil term + - `field`: A `ForcingField` with at least one coil term ### Keyword arguments @@ -271,10 +271,10 @@ GLMakie to rotate it. A Makie `Figure`. """ -function plot_sources(rmp::FT.RMPField; n::Int=1, color::Symbol=:current, save_path=nothing) - _require(plot_sources, rmp) +function plot_sources(field::FT.ForcingField; n::Int=1, color::Symbol=:current, save_path=nothing) + _require(plot_sources, field) color in (:current, :set) || throw(ArgumentError("plot_sources: color must be :current or :set")) - d = _rmp_coil_sets(rmp; n) + d = _field_coil_sets(field; n) notes = String[] d.skipped > 0 && push!(notes, "$(d.skipped) mode-table term$(d.skipped == 1 ? "" : "s") not drawn") d.complex_weight && push!(notes, "complex weights drawn at their magnitude") @@ -463,7 +463,7 @@ register!(plot_bn_contour, (PerturbedEquilibriumState,), (:forcing_modes,), "app register!(plot_coil_geometry_3d, (PerturbedEquilibriumState,), (:coil_sets,), "3D coils and plasma boundary") register!(plot_coil_geometry_rz, (PerturbedEquilibriumState,), (:coil_sets,), "coils over the plasma boundary in (R, Z)") register!(plot_mode_spectrum, (PerturbedEquilibriumState,), (:forcing_modes,), "|Φ_x(m, n)| forcing spectrum") -register!(plot_sources, (FT.RMPField,), (:coil_source,), "the coil sources in 3D, coloured by current or by set") +register!(plot_sources, (FT.ForcingField,), (:coil_source,), "the coil sources in 3D, coloured by current or by set") export plot_coil_geometry_3d, plot_coil_geometry_rz export plot_bn_contour, plot_mode_spectrum, plot_control_surface, plot_sources diff --git a/src/Analysis/Registry.jl b/src/Analysis/Registry.jl index 03940daa32..38d83d046e 100644 --- a/src/Analysis/Registry.jl +++ b/src/Analysis/Registry.jl @@ -61,7 +61,7 @@ const REASONS = Dict{Symbol,String}( :forcing_modes => "needs forcing modes: give ForcingTerms a forcing", :ef_sensitivities => "needs the coil sensitivities: add an [ErrorFields] section", :dominant_mode => "needs the dominant resonant-coupling mode: run PerturbedEquilibrium with compute_singular_coupling = true", - :coil_source => "needs a coil source: give RMPField coil sets (a mode table or spectrum has no geometry to draw)", + :coil_source => "needs a coil source: give ForcingField coil sets (a mode table or spectrum has no geometry to draw)", :coupling_matrices => "needs the resonant coupling matrices: run plasma_response with compute_singular_coupling = true", :ef_monte_carlo => "needs the tolerance Monte Carlo: give ErrorFields a tolerance file", :ef_risk => "needs the locking risk: give ErrorFields a tolerance file and a scenario", @@ -97,8 +97,7 @@ end capabilities(r::LocalStabilityResult) = isempty(r.alpha) ? Set{Symbol}() : Set([:ballooning_scan]) -capabilities(r::FT.RMPField) = - any(t -> t isa FT.RMPSource && t.ctrl.forcing_data_format == "coil", FT._rmp_terms(r)) ? Set([:coil_source]) : Set{Symbol}() +capabilities(r::FT.ForcingField) = any(t -> t isa FT.ForcingCoils, FT._field_terms(r)) ? Set([:coil_source]) : Set{Symbol}() function capabilities(r::PlasmaResponse) caps = Set{Symbol}() diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17ce..d449ba9e58 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -603,6 +603,20 @@ function _forcing_table(ctrl::ForcingTerms.ForcingTermsControl) return table end +# [ForcingTerms] of a single forcing source: a file by path, coils by their settings and the non-default keys of each set. +_forcing_table(f::ForcingTerms.ForcingFile) = Dict{String,Any}("forcing_data_file" => f.path, "forcing_data_format" => f.format) + +function _forcing_table(f::ForcingTerms.ForcingCoils) + cfg, default = f.config, ForcingTerms.CoilSetConfig() + table = Dict{String,Any}("forcing_data_format" => "coil", "machine" => cfg.machine, "mtheta_coil" => cfg.mtheta_coil, "nzeta_coil" => cfg.nzeta_coil) + cfg.dat_dir == ForcingTerms.BUNDLED_COIL_DIR || (table["dat_dir"] = cfg.dat_dir) + table["coil_set"] = Any[ + Dict{String,Any}(String(k) => _deck_out(getfield(sc, k)) for k in fieldnames(ForcingTerms.CoilSetConfig) if !isequal(getfield(sc, k), getfield(default, k))) + for sc in cfg.coil_sets + ] + return table +end + # [ErrorFields] with its sub-tables, from the settings an error-field assessment takes. function _error_fields_sections(settings, tolerance_file::AbstractString) table = _deck_table(settings.ctrl) @@ -646,7 +660,7 @@ function deck_sections(pe::PerturbedEquilibrium.PerturbedEquilibriumState) sections = Dict{String,Any}("PerturbedEquilibrium" => table) forcing === nothing && return sections sections["ForcingTerms"] = - forcing isa ForcingTerms.RMPSource && forcing.scale == 1 ? _forcing_table(forcing.ctrl) : + forcing isa Union{ForcingTerms.ForcingFile,ForcingTerms.ForcingCoils} && forcing.scale == 1 ? _forcing_table(forcing) : Dict{String,Any}("forcing_data_format" => "hdf5", "forcing_data_file" => _RECORDED_FORCING) return sections end diff --git a/src/Display.jl b/src/Display.jl index 662cbbd476..59ad5b3c9e 100644 --- a/src/Display.jl +++ b/src/Display.jl @@ -11,7 +11,7 @@ using .PerturbedEquilibrium: PerturbedEquilibriumState using .KineticForces: MethodResult, KineticForcesState using .ErrorFields: MonteCarloResult, RiskResult, EFCCoupling using .Tearing: GrowthRateResult -using .ForcingTerms: RMPSource, RMPSpectrum, RMPFieldSum +using .ForcingTerms: ForcingFile, ForcingCoils, ForcingSpectrum, ForcingLeaf, ForcingFieldSum _fmt(x::Real) = @sprintf("%.4g", x) _fmt(::Nothing) = "—" @@ -285,21 +285,21 @@ end _weight(w) = string(iszero(real(w)) && !iszero(imag(w)) ? _fmt(imag(w)) * "im" : _fmt(w), " × ") -function _rmp_label(f::RMPSource) - c = f.ctrl - c.forcing_data_format == "coil" || return string(basename(c.forcing_data_file), " (", c.forcing_data_format, ")") - sets = [isempty(get(d, "name", "")) ? get(d, "source", "coil set") : d["name"] for d in c.coil_sets_raw] +_field_label(f::ForcingFile) = string(basename(f.path), " (", f.format, ")") +function _field_label(f::ForcingCoils) + c = f.config + sets = [!isempty(sc.name) ? sc.name : !isempty(sc.source) ? sc.source : "coil set" for sc in c.coil_sets] return string(isempty(c.machine) ? "" : c.machine * " ", "coils ", _list(sets)) end -_rmp_label(f::RMPSpectrum) = string("spectrum of ", _plural(length(f.modes), "mode"), ", n = ", _irange([m.n for m in f.modes])) +_field_label(f::ForcingSpectrum) = string("spectrum of ", _plural(length(f.modes), "mode"), ", n = ", _irange([m.n for m in f.modes])) -Base.show(io::IO, f::Union{RMPSource,RMPSpectrum}) = print(io, "RMPField(", f.scale == 1 ? "" : _weight(f.scale), _rmp_label(f), ")") -Base.show(io::IO, f::RMPFieldSum) = print(io, "RMPField(sum of ", _plural(length(f.terms), "term"), ")") +Base.show(io::IO, f::ForcingLeaf) = print(io, "ForcingField(", f.scale == 1 ? "" : _weight(f.scale), _field_label(f), ")") +Base.show(io::IO, f::ForcingFieldSum) = print(io, "ForcingField(sum of ", _plural(length(f.terms), "term"), ")") -function Base.show(io::IO, ::MIME"text/plain", f::RMPFieldSum) - print(io, "RMPField sum of ", _plural(length(f.terms), "term")) +function Base.show(io::IO, ::MIME"text/plain", f::ForcingFieldSum) + print(io, "ForcingField sum of ", _plural(length(f.terms), "term")) for t in f.terms - print(io, "\n ", _weight(t.scale), _rmp_label(t)) + print(io, "\n ", _weight(t.scale), _field_label(t)) end return nothing end diff --git a/src/ErrorFields/NTVLimits.jl b/src/ErrorFields/NTVLimits.jl index 4504a8b457..598826d889 100644 --- a/src/ErrorFields/NTVLimits.jl +++ b/src/ErrorFields/NTVLimits.jl @@ -525,7 +525,7 @@ function efc_couplings( "scan ±$(@sprintf("%.3e", span)) rad/s on $(length(grid0)) initial points (max $(scan_max_points))" function response_of(modes::Vector{ForcingMode}) - pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(resp, ForcingTerms.RMPField(modes); regularization_width, verbose) + pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(resp, ForcingTerms.ForcingField(modes); regularization_width, verbose) kf_intr = KineticForces.KineticForcesInternal(equil; verbose=verbose) KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric) return kf_intr diff --git a/src/ForcingTerms/CoilGeometry.jl b/src/ForcingTerms/CoilGeometry.jl index e13277bdee..6d27f68a54 100644 --- a/src/ForcingTerms/CoilGeometry.jl +++ b/src/ForcingTerms/CoilGeometry.jl @@ -156,15 +156,16 @@ Base.@kwdef struct CoilConfig coil_sets::Vector{CoilSetConfig} = CoilSetConfig[] end +# The coil geometry files shipped with the package, the default `dat_dir`. +const BUNDLED_COIL_DIR = joinpath(@__DIR__, "coil_geometries") + """ CoilConfig(ft_ctrl::ForcingTermsControl) Construct a `CoilConfig` from a `ForcingTermsControl` populated from TOML. """ function CoilConfig(ft_ctrl::ForcingTermsControl) - dat_dir = isempty(ft_ctrl.dat_dir) ? - joinpath(@__DIR__, "coil_geometries") : - ft_ctrl.dat_dir + dat_dir = isempty(ft_ctrl.dat_dir) ? BUNDLED_COIL_DIR : ft_ctrl.dat_dir coil_sets = CoilSetConfig[_parse_coil_set_config(d) for d in ft_ctrl.coil_sets_raw] diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index c8ee6eab05..01da33e910 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -253,7 +253,7 @@ function load_forcing_from_h5_group!(forcing_modes::Vector{ForcingMode}, group): end """ - RMPField + ForcingField Abstract supertype of every external-forcing description the scripting API accepts: a resonant-magnetic-perturbation field, described by where it comes from and how strongly it @@ -261,30 +261,32 @@ drives (a forcing-mode file, a coil set with currents, or a precomputed spectrum Construction checks the description and throws an `ArgumentError` naming the problem: a forcing file must exist and parse as a mode table with a known normalization, and every coil set must resolve to geometry and carry currents. The -modes are materialized against an equilibrium only when the perturbed-equilibrium stage runs, so one `RMPField` can +modes are materialized against an equilibrium only when the perturbed-equilibrium stage runs, so one `ForcingField` can drive several solves. ## Construction - RMPField(path; format=..., scale=1.0, kwargs...) - RMPField(coil_sets; scale=1.0, kwargs...) - RMPField(ctrl::ForcingTermsControl; scale=1.0) - RMPField(modes::Vector{ForcingMode}; scale=1.0) + ForcingField(path; format=..., scale=1.0) + ForcingField(coil_sets; scale=1.0, kwargs...) + ForcingField(cfg::CoilConfig; scale=1.0) + ForcingField(ctrl::ForcingTermsControl; scale=1.0) + ForcingField(modes::Vector{ForcingMode}; scale=1.0) The first form points at a forcing-mode file, `format` defaulting to `"hdf5"` for an `.h5` or `.hdf5` extension and `"ascii"` otherwise. The second form takes one coil set or a vector of them, each a NamedTuple or Dict with the `[[ForcingTerms.coil_set]]` keys, and selects the coil format: ```julia -RMPField((name="c", currents=[1e3, 500, -500, -1e3, -500, 500]); machine="d3d") # the DIII-D C-coil, in amperes -RMPField([(name="iu", currents=iu), (name="il", currents=il)]; machine="d3d") -RMPField((source="pf_hoop", radius=12.0, height=1.2, currents=[2e3])) +ForcingField((name="c", currents=[1e3, 500, -500, -1e3, -500, 500]); machine="d3d") # the DIII-D C-coil, in amperes +ForcingField([(name="iu", currents=iu), (name="il", currents=il)]; machine="d3d") +ForcingField((source="pf_hoop", radius=12.0, height=1.2, currents=[2e3])) ``` -Remaining `kwargs` are -[`ForcingTermsControl`](@ref) fields. These three return a single-source leaf -([`RMPSource`](@ref)). The last form takes a control-surface spectrum computed elsewhere, already in -the unit-norm `Φ_x` convention, and returns a spectrum-literal leaf ([`RMPSpectrum`](@ref)). +Its remaining `kwargs` are the [`ForcingTermsControl`](@ref) coil settings (`machine`, `dat_dir`, `mtheta_coil`, +`nzeta_coil`). The third form takes the same coil description already parsed, and the fourth a deck's `[ForcingTerms]` +section. The file forms return a [`ForcingFile`](@ref) leaf and the coil forms a [`ForcingCoils`](@ref) leaf. The last form +takes a control-surface spectrum computed elsewhere, already in the unit-norm `Φ_x` convention, and returns a +spectrum-literal leaf ([`ForcingSpectrum`](@ref)). ## Algebra @@ -293,8 +295,8 @@ perturbation) build lazy linear combinations. The response is linear in the forc is materialized on the control surface and the mode amplitudes are summed. ```julia -nominal = RMPField("nominal_efc.dat") -weld_field = RMPField("weld_fields.h5"; scale=0.5) +nominal = ForcingField("nominal_efc.dat") +weld_field = ForcingField("weld_fields.h5"; scale=0.5) total = 2.0 * nominal + weld_field # lazy: records terms and weights, computes nothing pe = perturbed_equilibrium(plasma_response(ffs), total) # materializes each term, sums, drives PE ``` @@ -305,55 +307,76 @@ in a linear combination, applied uniformly to the materialized spectrum of every A source that is not a scalar multiple of another needs its own description: a coil set with a failed conductor is `nominal - failed_coil`, not `0.9 * nominal`. """ -abstract type RMPField end +abstract type ForcingField end """ - RMPSource + ForcingFile -A single-source [`RMPField`](@ref) leaf: one forcing description plus a complex weight. -Built by the `RMPField` constructors; scalar multiplication rescales the weight. +A forcing-mode file [`ForcingField`](@ref) leaf: a table of (n, m) amplitudes plus a complex weight. +Built by `ForcingField(path)`; scalar multiplication rescales the weight. ## Fields - - `ctrl::ForcingTermsControl` - The forcing source: format, file path or machine, and the raw coil-set blocks. + - `path::String` - Absolute path of the forcing-mode table. + - `format::String` - `"ascii"` or `"hdf5"`; see [`load_forcing_data!`](@ref). - `scale::ComplexF64` - The source's weight in a linear combination, applied to the materialized control-surface spectrum (complex = phase rotation). Not a physical amplitude. """ -struct RMPSource <: RMPField - ctrl::ForcingTermsControl +struct ForcingFile <: ForcingField + path::String + format::String scale::ComplexF64 end """ - RMPSpectrum + ForcingCoils -A spectrum-literal [`RMPField`](@ref) leaf: control-surface forcing modes computed elsewhere, in the +A coil [`ForcingField`](@ref) leaf: coil sets with their currents, whose Biot-Savart field is sampled on the control +surface, plus a complex weight. Built by `ForcingField(coil_sets)`; scalar multiplication rescales the weight. + +## Fields + + - `config::CoilConfig` - The coil sets and the boundary grid their field is sampled on. + - `scale::ComplexF64` - The source's weight in a linear combination, as for [`ForcingFile`](@ref). +""" +struct ForcingCoils <: ForcingField + config::CoilConfig + scale::ComplexF64 +end + +""" + ForcingSpectrum + +A spectrum-literal [`ForcingField`](@ref) leaf: control-surface forcing modes computed elsewhere, in the unit-norm `Φ_x` convention the response consumes, plus a complex weight. Built by -`RMPField(modes::Vector{ForcingMode})`; scalar multiplication rescales the weight. +`ForcingField(modes::Vector{ForcingMode})`; scalar multiplication rescales the weight. ## Fields - `modes::Vector{ForcingMode}` - The control-surface spectrum, copied at construction. - - `scale::ComplexF64` - The leaf's weight in a linear combination, as for [`RMPSource`](@ref). + - `scale::ComplexF64` - The leaf's weight in a linear combination, as for [`ForcingFile`](@ref). """ -struct RMPSpectrum <: RMPField +struct ForcingSpectrum <: ForcingField modes::Vector{ForcingMode} scale::ComplexF64 end +# A single-source ForcingField: one forcing description plus its weight. +const ForcingLeaf = Union{ForcingFile,ForcingCoils,ForcingSpectrum} + """ - RMPFieldSum + ForcingFieldSum -A lazy linear combination of [`RMPSource`](@ref) and [`RMPSpectrum`](@ref) leaves, built by `+`/`-` on -[`RMPField`](@ref)s. Holds the flattened term list; scalar multiplication distributes onto +A lazy linear combination of [`ForcingFile`](@ref), [`ForcingCoils`](@ref) and [`ForcingSpectrum`](@ref) leaves, built by `+`/`-` +on [`ForcingField`](@ref)s. Holds the flattened term list; scalar multiplication distributes onto the leaves. Materialization evaluates each term against the equilibrium and sums the mode amplitudes — valid because the perturbed-equilibrium response is linear in the forcing. ## Fields - - `terms::Vector{Union{RMPSource,RMPSpectrum}}` - The flattened weighted leaves. + - `terms::Vector{Union{ForcingFile,ForcingCoils,ForcingSpectrum}}` - The flattened weighted leaves. """ -struct RMPFieldSum <: RMPField - terms::Vector{Union{RMPSource,RMPSpectrum}} +struct ForcingFieldSum <: ForcingField + terms::Vector{ForcingLeaf} end """ @@ -366,34 +389,32 @@ _infer_format(path::AbstractString) = lowercase(splitext(path)[2]) in (".h5", ". const _NORMALIZATIONS = ("normal_field_T", "sfl_flux_Wb", "unit_flux") # Read a forcing-mode file now, so a missing or malformed file fails when the source is built. -function _check_forcing_file(ctrl::ForcingTermsControl) - path = ctrl.forcing_data_file - isfile(path) || throw(ArgumentError("RMPField: no forcing file at $path")) +function _check_forcing_file(path::AbstractString, format::AbstractString) + isfile(path) || throw(ArgumentError("ForcingField: no forcing file at $path")) modes = ForcingMode[] tag = try - load_forcing_data!(modes, path, ctrl.forcing_data_format) + load_forcing_data!(modes, path, format) catch err throw( ArgumentError( - "RMPField: $path is not a forcing-mode table (format $(ctrl.forcing_data_format)): $(sprint(showerror, err)). " * + "ForcingField: $path is not a forcing-mode table (format $format): $(sprint(showerror, err)). " * "A mode table holds rows of n, m, Re, Im (see `load_forcing_data!`); a coil geometry file goes in a coil set, " * - "as in RMPField((name=\"set\", dat_file=path, currents=[...]))." + "as in ForcingField((name=\"set\", dat_file=path, currents=[...]))." ) ) end - isempty(modes) && throw(ArgumentError("RMPField: $path holds no forcing modes")) - tag in _NORMALIZATIONS || throw(ArgumentError("RMPField: $path declares normalization \"$tag\"; use one of $(join(_NORMALIZATIONS, ", "))")) + isempty(modes) && throw(ArgumentError("ForcingField: $path holds no forcing modes")) + tag in _NORMALIZATIONS || throw(ArgumentError("ForcingField: $path declares normalization \"$tag\"; use one of $(join(_NORMALIZATIONS, ", "))")) return nothing end # Resolve every coil set now, so an unknown key, a missing geometry file or a missing current fails when the source is # built; window-pane sets placed by standoff need the plasma surface and are resolved when materialized. -function _check_coil_sets(ctrl::ForcingTermsControl) - isempty(ctrl.coil_sets_raw) && throw(ArgumentError("RMPField: a coil source needs at least one coil set")) - cfg = CoilConfig(ctrl) +function _check_coil_sets(cfg::CoilConfig) + isempty(cfg.coil_sets) && throw(ArgumentError("ForcingField: a coil source needs at least one coil set")) for sc in cfg.coil_sets label = isempty(sc.name) ? "a coil set" : "coil set \"$(sc.name)\"" - all(iszero, sc.currents) && throw(ArgumentError("RMPField: $label has no currents; give `currents` in amperes, one per conductor")) + all(iszero, sc.currents) && throw(ArgumentError("ForcingField: $label has no currents; give `currents` in amperes, one per conductor")) isnan(sc.standoff) || continue cs = try _build_raw_coil_set(cfg, sc) @@ -401,70 +422,74 @@ function _check_coil_sets(ctrl::ForcingTermsControl) hint = isempty(cfg.machine) && sc.source in ("", "file") && isempty(sc.dat_file) && isempty(sc.h5_file) ? " (a named set is read from _.dat: pass machine=\"d3d\", say, or give a dat_file)" : "" - throw(ArgumentError("RMPField: $label: $(sprint(showerror, err))$hint")) + throw(ArgumentError("ForcingField: $label: $(sprint(showerror, err))$hint")) end length(sc.currents) > cs.ncoil && - throw(ArgumentError("RMPField: $label has $(length(sc.currents)) currents for its $(cs.ncoil) conductors")) + throw(ArgumentError("ForcingField: $label has $(length(sc.currents)) currents for its $(cs.ncoil) conductors")) end return nothing end -function RMPField(ctrl::ForcingTermsControl; scale::Number=1.0) - fmt = ctrl.forcing_data_format - fmt == "coil" ? _check_coil_sets(ctrl) : - fmt in ("ascii", "hdf5") ? _check_forcing_file(ctrl) : - throw(ArgumentError("RMPField: forcing_data_format \"$fmt\" is not \"ascii\", \"hdf5\" or \"coil\"")) - return RMPSource(ctrl, ComplexF64(scale)) +function ForcingField(cfg::CoilConfig; scale::Number=1.0) + _check_coil_sets(cfg) + return ForcingCoils(cfg, ComplexF64(scale)) end -RMPField(path::AbstractString; format::String=_infer_format(path), scale::Number=1.0, kwargs...) = - RMPField(ForcingTermsControl(; forcing_data_file=abspath(path), forcing_data_format=format, kwargs...); scale) +function ForcingField(path::AbstractString; format::String=_infer_format(path), scale::Number=1.0) + format in ("ascii", "hdf5") || throw(ArgumentError("ForcingField: forcing_data_format \"$format\" is not \"ascii\", \"hdf5\" or \"coil\"")) + _check_forcing_file(abspath(path), format) + return ForcingFile(abspath(path), format, ComplexF64(scale)) +end + +function ForcingField(ctrl::ForcingTermsControl; scale::Number=1.0) + ctrl.forcing_data_format == "coil" && return ForcingField(CoilConfig(ctrl); scale) + return ForcingField(ctrl.forcing_data_file; format=ctrl.forcing_data_format, scale) +end # A coil set as the [[ForcingTerms.coil_set]] table: a Dict or NamedTuple, keys as strings or symbols. _coil_dict(d::AbstractDict) = Dict{String,Any}(string(k) => v for (k, v) in d) _coil_dict(nt::NamedTuple) = Dict{String,Any}(string(k) => v for (k, v) in pairs(nt)) -_coil_dict(x) = throw(ArgumentError("RMPField: a coil set is a Dict or NamedTuple of [[ForcingTerms.coil_set]] keys, not a $(typeof(x))")) +_coil_dict(x) = throw(ArgumentError("ForcingField: a coil set is a Dict or NamedTuple of [[ForcingTerms.coil_set]] keys, not a $(typeof(x))")) -RMPField(coil_sets::AbstractVector; scale::Number=1.0, kwargs...) = - RMPField(ForcingTermsControl(; forcing_data_format="coil", coil_sets_raw=_coil_dict.(coil_sets), kwargs...); scale) -RMPField(coil_set::Union{AbstractDict,NamedTuple}; kwargs...) = RMPField([coil_set]; kwargs...) +ForcingField(coil_sets::AbstractVector; scale::Number=1.0, kwargs...) = + ForcingField(ForcingTermsControl(; forcing_data_format="coil", coil_sets_raw=_coil_dict.(coil_sets), kwargs...); scale) +ForcingField(coil_set::Union{AbstractDict,NamedTuple}; kwargs...) = ForcingField([coil_set]; kwargs...) -function RMPField(modes::Vector{ForcingMode}; scale::Number=1.0) - isempty(modes) && throw(ArgumentError("RMPField: a spectrum needs at least one forcing mode")) - all(m -> isfinite(m.amplitude), modes) || throw(ArgumentError("RMPField: a spectrum amplitude is not finite")) - return RMPSpectrum([ForcingMode(; n=mode.n, m=mode.m, amplitude=mode.amplitude) for mode in modes], ComplexF64(scale)) +function ForcingField(modes::Vector{ForcingMode}; scale::Number=1.0) + isempty(modes) && throw(ArgumentError("ForcingField: a spectrum needs at least one forcing mode")) + all(m -> isfinite(m.amplitude), modes) || throw(ArgumentError("ForcingField: a spectrum amplitude is not finite")) + return ForcingSpectrum([ForcingMode(; n=mode.n, m=mode.m, amplitude=mode.amplitude) for mode in modes], ComplexF64(scale)) end """ -Flattened weighted-leaf list of any [`RMPField`](@ref). +Flattened weighted-leaf list of any [`ForcingField`](@ref). """ -_rmp_terms(f::Union{RMPSource,RMPSpectrum}) = Union{RMPSource,RMPSpectrum}[f] -_rmp_terms(f::RMPFieldSum) = f.terms - -Base.:+(a::RMPField, b::RMPField) = RMPFieldSum(vcat(_rmp_terms(a), _rmp_terms(b))) -Base.:-(a::RMPField) = -1 * a -Base.:-(a::RMPField, b::RMPField) = a + (-1 * b) -Base.:*(c::Number, f::RMPSource) = RMPSource(f.ctrl, ComplexF64(c) * f.scale) -Base.:*(c::Number, f::RMPSpectrum) = RMPSpectrum(f.modes, ComplexF64(c) * f.scale) -Base.:*(c::Number, f::RMPFieldSum) = RMPFieldSum([c * t for t in f.terms]) -Base.:*(f::RMPField, c::Number) = c * f +_field_terms(f::ForcingLeaf) = ForcingLeaf[f] +_field_terms(f::ForcingFieldSum) = f.terms + +Base.:+(a::ForcingField, b::ForcingField) = ForcingFieldSum(vcat(_field_terms(a), _field_terms(b))) +Base.:-(a::ForcingField) = -1 * a +Base.:-(a::ForcingField, b::ForcingField) = a + (-1 * b) +Base.:*(c::Number, f::ForcingFile) = ForcingFile(f.path, f.format, ComplexF64(c) * f.scale) +Base.:*(c::Number, f::ForcingCoils) = ForcingCoils(f.config, ComplexF64(c) * f.scale) +Base.:*(c::Number, f::ForcingSpectrum) = ForcingSpectrum(f.modes, ComplexF64(c) * f.scale) +Base.:*(c::Number, f::ForcingFieldSum) = ForcingFieldSum([c * t for t in f.terms]) +Base.:*(f::ForcingField, c::Number) = c * f """ - snapshot_forcing_modes(ft_ctrl; verbose=false) -> Union{Nothing,Tuple{Vector{ForcingMode},String}} + snapshot_forcing_modes(f::ForcingFile; verbose=false) -> Tuple{Vector{ForcingMode},String} -The raw file-based forcing modes of `ft_ctrl` and their normalization tag, so they land in -`Input/RawInputs/ForcingTerms/` alongside the TOML blob; `nothing` for coil forcing, which is -recomputed from the `[[coil_set]]` TOML blob on replay. They feed only the writer: +The raw forcing modes of the file `f` and their normalization tag, so they land in +`Input/RawInputs/ForcingTerms/` alongside the TOML blob. They feed only the writer: the PE stage reads and normalizes the forcing file itself. """ -function snapshot_forcing_modes(ft_ctrl::ForcingTermsControl; verbose::Bool=false) - ft_ctrl.forcing_data_format in ("ascii", "hdf5") || return nothing +function snapshot_forcing_modes(f::ForcingFile; verbose::Bool=false) modes = ForcingMode[] - normalization = load_forcing_data!(modes, ft_ctrl.forcing_data_file, ft_ctrl.forcing_data_format, verbose) + normalization = load_forcing_data!(modes, f.path, f.format, verbose) return modes, normalization end -export ForcingTermsControl, ForcingMode, RMPField, load_forcing_data!, save_forcing_to_h5, load_forcing_from_h5_group! +export ForcingTermsControl, ForcingMode, ForcingField, load_forcing_data!, save_forcing_to_h5, load_forcing_from_h5_group! end # module ForcingTerms diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 50ecbd2550..6c907a065e 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -77,7 +77,7 @@ using .InnerLayer: GGJModel, SLAYERModel using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult using .Equilibrium: PlasmaEquilibrium, KineticProfileData, IonSpecies, attach_kinetic_profiles! -using .ForcingTerms: RMPField +using .ForcingTerms: ForcingField using .PerturbedEquilibrium: plasma_response, perturbed_equilibrium, PlasmaResponse using .KineticForces: kinetic_torque, KineticForcesControl, KineticForcesResult using .ErrorFields: error_fields, ErrorFieldsResult @@ -199,8 +199,8 @@ function main_from_inputs( if deck.perturbed_equilibrium !== nothing pe = _timed_stage(runtimes, "perturbed_equilibrium", "Perturbed Equilibrium") do # A replay from stored coil geometry skips the source checks: the geometry files they read may be gone. - rmp = preloaded_coil_sets === nothing ? RMPField(deck.forcing) : ForcingTerms.RMPSource(deck.forcing, 1) - perturbed_equilibrium(plasma_response(ffs; verbose, deck.perturbed_equilibrium.response...), rmp; + field = preloaded_coil_sets === nothing ? ForcingField(deck.forcing) : ForcingTerms.ForcingCoils(ForcingTerms.CoilConfig(deck.forcing), 1) + perturbed_equilibrium(plasma_response(ffs; verbose, deck.perturbed_equilibrium.response...), field; coil_sets=preloaded_coil_sets, verbose, deck.perturbed_equilibrium.drive...) end resp = pe.response @@ -312,7 +312,8 @@ end export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! -export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse +export PlasmaEquilibrium, + ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, ForcingField, PlasmaResponse export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 2bd8d63382..ab868a4b8f 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -458,7 +458,7 @@ kinetic profiles attached to its equilibrium. A multi-species plasma adds each s Call `save_to_hdf5(target, kf)` to save the result. ```julia -pe = perturbed_equilibrium(plasma_response(ffs), RMPField("forcing.dat")) +pe = perturbed_equilibrium(plasma_response(ffs), ForcingField("forcing.dat")) kf = kinetic_torque(pe, KineticForcesControl()) ``` """ diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 2beaa59163..eef980750c 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -48,7 +48,7 @@ control-surface response operators (inductances, permeability, reluctance; Park Plasmas 14, 052110 (2007)) and, with `compute_singular_coupling`, the singular-coupling matrices and their dominant modes (Park, Phys. Plasmas 16, 056115 (2009)), and with `compute_torque_response` the ψ-resolved torque response `T_xe` of a single-`n` forward kinetic solve. Apply it to an external field -with `perturbed_equilibrium(resp, rmp)` as many times as needed. +with `perturbed_equilibrium(resp, field)` as many times as needed. Products the producing integrator could not supply (free-boundary energies, a ξ solution) leave the corresponding operators empty: the step warns and is skipped. Galerkin and `:matched` results publish @@ -88,11 +88,11 @@ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling:: end """ - perturbed_equilibrium(resp, rmp; regularization_width=0.0, verbose=false, coil_sets=nothing) + perturbed_equilibrium(resp, field; regularization_width=0.0, verbose=false, coil_sets=nothing) -> PerturbedEquilibriumState -Drive the plasma response `resp` (from [`plasma_response`](@ref)) with the external field `rmp`, an -[`RMPField`](@ref GeneralizedPerturbedEquilibrium.ForcingTerms.RMPField), and return the perturbed-equilibrium result. One response serves any number of +Drive the plasma response `resp` (from [`plasma_response`](@ref)) with the external field `field`, an +[`ForcingField`](@ref GeneralizedPerturbedEquilibrium.ForcingTerms.ForcingField), and return the perturbed-equilibrium result. One response serves any number of fields, so a scan rebuilds no operator. `regularization_width` smooths the ideal 1/(m−nq) divergence in the profile reconstruction (Fortran @@ -101,19 +101,19 @@ Call `save_to_hdf5(target, pe)` to save the result. `coil_sets` injects already- ```julia resp = plasma_response(ffs) -pe = perturbed_equilibrium(resp, RMPField("forcing.dat")) +pe = perturbed_equilibrium(resp, ForcingField("forcing.dat")) ``` """ function perturbed_equilibrium( resp::PlasmaResponse, - rmp::ForcingTerms.RMPField; + field::ForcingTerms.ForcingField; regularization_width::Real=0.0, verbose::Bool=false, coil_sets::Union{Nothing,Vector{CoilSet}}=nothing ) # Injected coil geometry (gpec.h5 replay with `--coil-source coils`) lets the coil field # be recomputed from stored geometry without the .dat/.h5 files. - pe_state = compute_perturbed_equilibrium(resp, rmp; regularization_width, verbose, + pe_state = compute_perturbed_equilibrium(resp, field; regularization_width, verbose, preloaded_coil_sets=coil_sets === nothing ? CoilSet[] : copy(coil_sets)) return pe_state @@ -123,8 +123,8 @@ end compute_perturbed_equilibrium(resp, forcing; regularization_width=0.0, verbose=true, preloaded_coil_sets=CoilSet[])::PerturbedEquilibriumState -Drive the operators of `resp` with one external field: materialize `forcing` (a `ForcingTermsControl` -or any [`ForcingTerms.RMPField`](@ref)) on the control surface, then apply the response and the +Drive the operators of `resp` with one external field: materialize `forcing` (any +[`ForcingTerms.ForcingField`](@ref)) on the control surface, then apply the response and the singular coupling to it, reconstruct the perturbed profiles, and contract any torque response with the forcing sources. Steps whose operators `resp` does not carry are skipped, leaving only the forcing spectrum. @@ -137,15 +137,14 @@ forcing spectrum. """ function compute_perturbed_equilibrium( resp::PlasmaResponse, - forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}; + forcing::ForcingTerms.ForcingField; regularization_width::Real=0.0, verbose::Bool=true, preloaded_coil_sets::Vector{CoilSet}=CoilSet[] )::PerturbedEquilibriumState ffs = resp.ffs modes, coil_sets = materialize_forcing_modes(ffs, forcing; preloaded_coil_sets, verbose) - state = PerturbedEquilibriumState(; response=resp, forcing=forcing isa ForcingTerms.RMPField ? forcing : ForcingTerms.RMPField(forcing), - forcing_modes=modes, + state = PerturbedEquilibriumState(; response=resp, forcing, forcing_modes=modes, coil_sets=isempty(coil_sets) ? preloaded_coil_sets : coil_sets, regularization_width=Float64(regularization_width)) isempty(resp.permeability_flux) || apply_response!(state; verbose) @@ -158,92 +157,76 @@ function compute_perturbed_equilibrium( end """ - materialize_forcing_modes(ffs, forcing; preloaded_coil_sets=CoilSet[], verbose=false) + materialize_forcing_modes(ffs, field::ForcingTerms.ForcingCoils; preloaded_coil_sets=CoilSet[], verbose=false) -> (modes::Vector{ForcingMode}, coil_sets::Vector{CoilSet}) -Pure computation of the control-surface forcing spectrum: turn a forcing description into -fresh `ForcingMode`s, for every toroidal mode of the solve and in the unit-norm convention -the response step consumes. Nothing is mutated — the caller owns where the result lands. - -Coil formats build the geometry (or reuse `preloaded_coil_sets` from the gpec.h5 replay -path) and integrate the Biot-Savart field on the control surface; file formats read the -modes from `forcing_data_file` and convert their normalization onto the run's `mlow:mhigh`. Both branches evaluate on -`ffs.psilim`, the integration limit. The second return value is the coil geometry actually -used (empty for file formats), which the writer snapshots for replay. +Pure computation of the control-surface forcing spectrum of a coil leaf: build the geometry (or reuse +`preloaded_coil_sets` from the gpec.h5 replay path), integrate the Biot-Savart field on the control surface at +`ffs.psilim`, the integration limit, for every toroidal mode of the solve, and apply the leaf's weight. Modes are in +the unit-norm convention the response step consumes. The second return value is the coil geometry used, which the +writer snapshots for replay. Every `materialize_forcing_modes` method returns fresh modes and mutates nothing. """ function materialize_forcing_modes( ffs::ForceFreeStatesResult, - forcing::ForcingTerms.ForcingTermsControl; + field::ForcingTerms.ForcingCoils; preloaded_coil_sets::Vector{ForcingTerms.CoilSet}=ForcingTerms.CoilSet[], verbose::Bool=false ) - equil = ffs.equil + equil, cfg = ffs.equil, field.config + coil_sets = isempty(preloaded_coil_sets) ? ForcingTerms.load_coil_sets(cfg, ffs.nlow; equil=equil) : preloaded_coil_sets modes = ForcingMode[] - coil_sets = ForcingTerms.CoilSet[] - - if forcing.forcing_data_format == "coil" - cfg = ForcingTerms.CoilConfig(forcing) - coil_sets = isempty(preloaded_coil_sets) ? - ForcingTerms.load_coil_sets(cfg, ffs.nlow; equil=equil) : preloaded_coil_sets - for n in ffs.nlow:ffs.nhigh - modes_n = ForcingMode[] - ForcingTerms.compute_coil_forcing_modes!( - modes_n, coil_sets, equil, cfg, n, ffs.mlow, ffs.mhigh; - psi=ffs.psilim, verbose=verbose - ) - append!(modes, modes_n) - end - else - file_modes = ForcingMode[] - norm_tag = load_forcing_data!(file_modes, forcing.forcing_data_file, forcing.forcing_data_format, verbose) - for n in ffs.nlow:ffs.nhigh - modes_n = filter(m -> m.n == n, file_modes) - isempty(modes_n) && continue - # Normalize on the same control surface as the coil branch: psilim, the integration limit. - append!(modes, convert_forcing_normalization(modes_n, norm_tag, equil, n, ffs.mlow:ffs.mhigh; psi=ffs.psilim)) - end + for n in ffs.nlow:ffs.nhigh + modes_n = ForcingMode[] + ForcingTerms.compute_coil_forcing_modes!(modes_n, coil_sets, equil, cfg, n, ffs.mlow, ffs.mhigh; psi=ffs.psilim, verbose=verbose) + append!(modes, modes_n) end - + field.scale == 1 || (modes = [ForcingMode(; n=mode.n, m=mode.m, amplitude=field.scale * mode.amplitude) for mode in modes]) return modes, coil_sets end """ - materialize_forcing_modes(ffs, rmp::ForcingTerms.RMPSource; kwargs...) -> (modes, coil_sets) + materialize_forcing_modes(ffs, field::ForcingTerms.ForcingFile; verbose=false, kwargs...) -> (modes, coil_sets) -Materialize a single-source [`ForcingTerms.RMPField`](@ref) leaf and apply its weight. The -returned modes are fresh objects, so weighting never reaches back into any caller state. +Materialize a forcing-mode file: read its modes, convert their normalization onto the run's `mlow:mhigh` on the +control surface `ffs.psilim`, and apply the leaf's weight. No coil geometry is produced. """ -function materialize_forcing_modes(ffs::ForceFreeStatesResult, rmp::ForcingTerms.RMPSource; kwargs...) - modes, coil_sets = materialize_forcing_modes(ffs, rmp.ctrl; kwargs...) - if rmp.scale != 1 - modes = [ForcingMode(; n=mode.n, m=mode.m, amplitude=rmp.scale * mode.amplitude) for mode in modes] +function materialize_forcing_modes(ffs::ForceFreeStatesResult, field::ForcingTerms.ForcingFile; verbose::Bool=false, kwargs...) + file_modes = ForcingMode[] + norm_tag = load_forcing_data!(file_modes, field.path, field.format, verbose) + modes = ForcingMode[] + for n in ffs.nlow:ffs.nhigh + modes_n = filter(m -> m.n == n, file_modes) + isempty(modes_n) && continue + # Normalize on the same control surface as the coil field: psilim, the integration limit. + append!(modes, convert_forcing_normalization(modes_n, norm_tag, ffs.equil, n, ffs.mlow:ffs.mhigh; psi=ffs.psilim)) end - return modes, coil_sets + field.scale == 1 || (modes = [ForcingMode(; n=mode.n, m=mode.m, amplitude=field.scale * mode.amplitude) for mode in modes]) + return modes, ForcingTerms.CoilSet[] end """ - materialize_forcing_modes(ffs, rmp::ForcingTerms.RMPSpectrum; kwargs...) -> (modes, coil_sets) + materialize_forcing_modes(ffs, field::ForcingTerms.ForcingSpectrum; kwargs...) -> (modes, coil_sets) Materialize a spectrum-literal leaf: fresh copies of its unit-norm modes, weighted. Nothing is read or integrated, and no coil geometry is produced. """ -function materialize_forcing_modes(::ForceFreeStatesResult, rmp::ForcingTerms.RMPSpectrum; kwargs...) - return [ForcingMode(; n=mode.n, m=mode.m, amplitude=rmp.scale * mode.amplitude) for mode in rmp.modes], ForcingTerms.CoilSet[] +function materialize_forcing_modes(::ForceFreeStatesResult, field::ForcingTerms.ForcingSpectrum; kwargs...) + return [ForcingMode(; n=mode.n, m=mode.m, amplitude=field.scale * mode.amplitude) for mode in field.modes], ForcingTerms.CoilSet[] end """ - materialize_forcing_modes(ffs, rmp::ForcingTerms.RMPFieldSum; kwargs...) -> (modes, coil_sets) + materialize_forcing_modes(ffs, field::ForcingTerms.ForcingFieldSum; kwargs...) -> (modes, coil_sets) Materialize a lazy linear combination of forcing sources: each weighted leaf is evaluated independently and the mode amplitudes are summed per (n, m), sorted by (n, m). Coil geometry from every coil-format leaf is concatenated in the second return value; `preloaded_coil_sets` is ignored for sums. """ -function materialize_forcing_modes(ffs::ForceFreeStatesResult, rmp::ForcingTerms.RMPFieldSum; +function materialize_forcing_modes(ffs::ForceFreeStatesResult, field::ForcingTerms.ForcingFieldSum; preloaded_coil_sets::Vector{ForcingTerms.CoilSet}=ForcingTerms.CoilSet[], verbose::Bool=false) amplitudes = Dict{Tuple{Int,Int},ComplexF64}() coil_sets = ForcingTerms.CoilSet[] - for term in rmp.terms + for term in field.terms term_modes, term_coils = materialize_forcing_modes(ffs, term; verbose=verbose) for mode in term_modes key = (mode.n, mode.m) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index b31ed49d3a..77248c769b 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -39,7 +39,7 @@ end The source-independent half of a perturbed equilibrium: the control-surface response and singular-coupling operators of one force-free-states result, built once by [`plasma_response`](@ref) and applied to any number of external fields by -`perturbed_equilibrium(resp, rmp)` (Park et al., Phys. Plasmas 14, 052110 (2007); Park, +`perturbed_equilibrium(resp, field)` (Park et al., Phys. Plasmas 14, 052110 (2007); Park, Phys. Plasmas 16, 056115 (2009)). Operators the result could not supply (no free-boundary energies or no ξ solution) are left empty, and forcing such a response yields only the forcing spectra. @@ -155,7 +155,7 @@ The perturbed equilibrium driven by one external field: the forcing materialized ## Fields - `response::PlasmaResponse` - The operators this field was applied to. - - `forcing` - The [`ForcingTerms.RMPField`](@ref) applied. + - `forcing` - The [`ForcingTerms.ForcingField`](@ref) applied. - `forcing_modes::Vector{ForcingMode}` - The materialized control-surface forcing Φ_x. - `coil_sets::Vector{CoilSet}` - Coil geometry used (coil-format forcing); captured for the gpec.h5 rerun snapshot. - `regularization_width::Float64` - Width the profiles were reconstructed with (Fortran `reg_spot`): 0, no @@ -210,7 +210,7 @@ Torque response (kinetic forward solves with `compute_torque_response`): """ @kwdef mutable struct PerturbedEquilibriumState response::PlasmaResponse - forcing::Union{Nothing,ForcingTerms.RMPField} = nothing + forcing::Union{Nothing,ForcingTerms.ForcingField} = nothing forcing_modes::Vector{ForcingMode} = ForcingMode[] coil_sets::Vector{CoilSet} = CoilSet[] regularization_width::Float64 = 0.0 diff --git a/src/PerturbedEquilibrium/TorqueResponse.jl b/src/PerturbedEquilibrium/TorqueResponse.jl index 09cb18fc76..62935ba2ba 100644 --- a/src/PerturbedEquilibrium/TorqueResponse.jl +++ b/src/PerturbedEquilibrium/TorqueResponse.jl @@ -126,28 +126,22 @@ end Contract the torque response of `state.response` with the forcing of `state` into `state.torque_response`. Every forcing source is one column of `T_source`: each coil set -(`state.coil_sets`, on the grid `forcing` carries) when the forcing is coils, otherwise the -applied spectrum as a single source named after its forcing file. A source's Φ_x column is +(`state.coil_sets`, on the grid `forcing` carries) when the forcing is a coil leaf, otherwise the +applied spectrum as a single source, named after its forcing file when it is one. A source's Φ_x column is referred to b̃ through `R = S·A`, so `T_source = B̃†·T_xe·B̃` with `B̃ = R⁻¹·Φ_source` (the Fortran `gcoil`). """ -function apply_torque_response!( - state::PerturbedEquilibriumState, - forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}; - verbose::Bool=false -) +function apply_torque_response!(state::PerturbedEquilibriumState, forcing::ForcingTerms.ForcingField; verbose::Bool=false) resp = state.response ffs = resp.ffs - ft_ctrl = forcing isa ForcingTerms.ForcingTermsControl ? forcing : - forcing isa ForcingTerms.RMPSource ? forcing.ctrl : nothing source_b = zeros(ComplexF64, ffs.numpert_total, 0) source_names = String[] - if !isempty(state.coil_sets) && ft_ctrl !== nothing && ft_ctrl.forcing_data_format == "coil" + if forcing isa ForcingTerms.ForcingCoils && !isempty(state.coil_sets) flux_conform = resp.rootarea_to_area_weight .* resp.surface_area - source_b = flux_conform \ coil_flux_spectra(ffs, state.coil_sets, ForcingTerms.CoilConfig(ft_ctrl)) + source_b = flux_conform \ coil_flux_spectra(ffs, state.coil_sets, forcing.config) source_names = [cs.name for cs in state.coil_sets] elseif !isempty(state.forcing_b_rootarea) source_b = reshape(state.forcing_b_rootarea, :, 1) - source_names = [ft_ctrl === nothing ? "forcing" : basename(ft_ctrl.forcing_data_file)] + source_names = [forcing isa ForcingTerms.ForcingFile ? basename(forcing.path) : "forcing"] end T_xe = resp.T_xe diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 1185379dd9..572c4f01a1 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -195,9 +195,9 @@ end function _record_forcing_snapshot!(file, state::PerturbedEquilibriumState) forcing = state.forcing forcing === nothing && return nothing - snapshot = forcing isa ForcingTerms.RMPSource && forcing.scale == 1 ? ForcingTerms.snapshot_forcing_modes(forcing.ctrl) : + forcing isa ForcingTerms.ForcingCoils && forcing.scale == 1 && return nothing + snapshot = forcing isa ForcingTerms.ForcingFile && forcing.scale == 1 ? ForcingTerms.snapshot_forcing_modes(forcing) : (state.forcing_modes, "unit_flux") - snapshot === nothing && return nothing modes, normalization = snapshot ForcingTerms.save_forcing_to_h5(modes, replace_group!(file, "Input/RawInputs/ForcingTerms"); normalization) return nothing diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 6951228387..b3d799edc5 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -109,10 +109,10 @@ using TOML bare, solved = A.Equilibrium._surfaces(ret.equil), A.Equilibrium._surfaces(ret.ffs) @test all(q -> any(≈(q), bare.sing_q), solved.sing_q) && !isempty(bare.sing_q) @test plot_equilibrium_summary(ret.equil; n=2) isa Figure - # An RMPField draws its coil sources in 3D before any plasma; a mode table has no geometry to draw. - hoop = RMPField((name="hoop", source="pf_hoop", radius=12.0, height=1.2, currents=[2.0e3], tiltx=[3.0])) + # A ForcingField draws its coil sources in 3D before any plasma; a mode table has no geometry to draw. + hoop = ForcingField((name="hoop", source="pf_hoop", radius=12.0, height=1.2, currents=[2.0e3], tiltx=[3.0])) @test plot_sources(hoop) isa Figure && plot_sources(2hoop - 0.5im * hoop; color=:set) isa Figure - @test_throws A.REASONS[:coil_source] plot_sources(RMPField([GPEC.ForcingTerms.ForcingMode(; n=1, m=2, amplitude=1e-4)])) + @test_throws A.REASONS[:coil_source] plot_sources(ForcingField([GPEC.ForcingTerms.ForcingMode(; n=1, m=2, amplitude=1e-4)])) # A summary none of whose panels this run supplies is listed as unavailable, and says so when called. @test occursin("✗ plot_error_field_summary", repr(MIME"text/plain"(), available_plots(h5))) @test_throws A.REASONS[:ef_sensitivities] plot_error_field_summary(h5) diff --git a/test/runtests_deck_roundtrip.jl b/test/runtests_deck_roundtrip.jl index d0df7e00a3..92e0d1779e 100644 --- a/test/runtests_deck_roundtrip.jl +++ b/test/runtests_deck_roundtrip.jl @@ -62,7 +62,7 @@ end ffs = solve(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm) ls = solve(LocalStabilityProblem(ffs.equil; deck.local_stability.problem...), deck.local_stability.algorithm) resp = plasma_response(ffs; deck.perturbed_equilibrium.response...) - pe = perturbed_equilibrium(resp, RMPField(deck.forcing); deck.perturbed_equilibrium.drive...) + pe = perturbed_equilibrium(resp, ForcingField(deck.forcing); deck.perturbed_equilibrium.drive...) script_h5 = joinpath(dir, "script.h5") foreach(r -> save_to_hdf5(script_h5, r), (ffs, ls, pe)) @test deck_same(stages(recorded(script_h5), dir), stages(recorded(deck_h5), dir)) @@ -82,7 +82,7 @@ end @test replay.pe.forcing_b == pe.forcing_b # A weighted field has no [ForcingTerms]: the replay drives the recorded control-surface modes. - doubled = perturbed_equilibrium(resp, 2 * RMPField(deck.forcing); deck.perturbed_equilibrium.drive...) + doubled = perturbed_equilibrium(resp, 2 * ForcingField(deck.forcing); deck.perturbed_equilibrium.drive...) weighted_h5 = joinpath(dir, "weighted.h5") foreach(r -> save_to_hdf5(weighted_h5, r), (ffs, doubled)) @test recorded(weighted_h5)["ForcingTerms"]["forcing_data_file"] == "Input/RawInputs/ForcingTerms" diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 41883afea2..448c053254 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -232,7 +232,7 @@ include("display_check.jl") ffs = solve(equil, forward; nn=1, wall=wall, ffs_kwargs...) resp = plasma_response(ffs; verbose=false) - pe = perturbed_equilibrium(resp, RMPField(forcing); pe_kw...) + pe = perturbed_equilibrium(resp, ForcingField(forcing); pe_kw...) # The response matrices depend on the equilibrium and the solve, not on the drive. @test resp isa PlasmaResponse && pe.response === resp @@ -241,38 +241,38 @@ include("display_check.jl") @test size(resp.C_delta_prime) == size(reference.response.C_delta_prime) @test !isempty(pe.resonant_area_weighted_field) - # The deck reads and normalizes its forcing file exactly as a fresh RMPField does. + # The deck reads and normalizes its forcing file exactly as a fresh ForcingField does. @test pe.resonant_area_weighted_field ≈ reference.resonant_area_weighted_field rtol = 1e-10 @test pe.forcing_b ≈ reference.forcing_b rtol = 1e-10 # One response serves any number of fields: forcing it again reproduces the result. - again = perturbed_equilibrium(resp, RMPField(forcing); pe_kw...) + again = perturbed_equilibrium(resp, ForcingField(forcing); pe_kw...) @test again.forcing_b == pe.forcing_b @test again.resonant_area_weighted_field == pe.resonant_area_weighted_field # `scale` is a uniform multiplier on the materialized forcing, and the response is # linear in it. - scaled = perturbed_equilibrium(resp, RMPField(forcing; scale=2.0); pe_kw...) + scaled = perturbed_equilibrium(resp, ForcingField(forcing; scale=2.0); pe_kw...) @test scaled.forcing_b ≈ 2 .* pe.forcing_b rtol = 1e-10 @test scaled.resonant_area_weighted_field ≈ 2 .* pe.resonant_area_weighted_field rtol = 1e-10 # Lazy source algebra materializes to the combined field: 3A - A == 2A drives # the same perturbed equilibrium as scale=2 (exercises +, -, * and the merge). - combo = perturbed_equilibrium(resp, 3 * RMPField(forcing) - RMPField(forcing); pe_kw...) + combo = perturbed_equilibrium(resp, 3 * ForcingField(forcing) - ForcingField(forcing); pe_kw...) @test combo.forcing_b ≈ scaled.forcing_b rtol = 1e-10 @test combo.resonant_area_weighted_field ≈ scaled.resonant_area_weighted_field rtol = 1e-10 # A spectrum-literal leaf of the materialized modes drives the same perturbed equilibrium, # takes a weight, and sums with a file leaf. - literal = perturbed_equilibrium(resp, RMPField(pe.forcing_modes); pe_kw...) + literal = perturbed_equilibrium(resp, ForcingField(pe.forcing_modes); pe_kw...) @test literal.forcing_b == pe.forcing_b @test literal.resonant_area_weighted_field == pe.resonant_area_weighted_field - @test perturbed_equilibrium(resp, 2 * RMPField(pe.forcing_modes); pe_kw...).forcing_b ≈ scaled.forcing_b rtol = 1e-10 - @test perturbed_equilibrium(resp, RMPField(pe.forcing_modes) + RMPField(forcing); pe_kw...).forcing_b ≈ scaled.forcing_b rtol = 1e-10 + @test perturbed_equilibrium(resp, 2 * ForcingField(pe.forcing_modes); pe_kw...).forcing_b ≈ scaled.forcing_b rtol = 1e-10 + @test perturbed_equilibrium(resp, ForcingField(pe.forcing_modes) + ForcingField(forcing); pe_kw...).forcing_b ≈ scaled.forcing_b rtol = 1e-10 # No regularization by default; an explicit width changes the regularized profiles, not the control-surface quantities. @test pe.regularization_width == 0.0 - reg = perturbed_equilibrium(resp, RMPField(forcing); regularization_width=0.05, pe_kw...) + reg = perturbed_equilibrium(resp, ForcingField(forcing); regularization_width=0.05, pe_kw...) @test reg.regularization_width == 0.05 @test reg.forcing_b == pe.forcing_b @test !(reg.xi_modes.clebsch_alpha ≈ pe.xi_modes.clebsch_alpha) @@ -288,12 +288,12 @@ include("display_check.jl") f["n"], f["m"], f["amplitude_real"], f["amplitude_imag"] = [1, 1], [2, 3], [1e-4, 5e-5], [0.0, 0.0] end - @testset "RMPField algebra is lazy and flattens" begin - a = RMPField(table) - b = RMPField(table; scale=0.5) - c = RMPField(table_h5) + @testset "ForcingField algebra is lazy and flattens" begin + a = ForcingField(table) + b = ForcingField(table; scale=0.5) + c = ForcingField(table_h5) s = a + b - @test s isa GPEC.ForcingTerms.RMPFieldSum + @test s isa GPEC.ForcingTerms.ForcingFieldSum @test length(s.terms) == 2 @test length((a + b + c).terms) == 3 d = 2.0 * s @@ -303,37 +303,48 @@ include("display_check.jl") @test (a * 3).scale == 3.0 + 0.0im @test (a-b).terms[2].scale == -0.5 + 0.0im @test (-a).scale == -1.0 + 0.0im - @test repr(a) == "RMPField(forcing.dat (ascii))" && repr(3 * a) == "RMPField(3 × forcing.dat (ascii))" + @test repr(a) == "ForcingField(forcing.dat (ascii))" && repr(3 * a) == "ForcingField(3 × forcing.dat (ascii))" @test occursin("-0.5im × forcing.h5 (hdf5)", repr(MIME"text/plain"(), a - 0.5im * c)) end - @testset "RMPField infers its format and carries its scale" begin - @test RMPField(table).ctrl.forcing_data_format == "ascii" - @test RMPField(table_h5).ctrl.forcing_data_format == "hdf5" - @test isabspath(RMPField(table).ctrl.forcing_data_file) - @test RMPField(table; scale=3.0).scale == 3.0 + @testset "ForcingField infers its format and carries its scale" begin + @test ForcingField(table) isa GPEC.ForcingTerms.ForcingFile + @test ForcingField(table).format == "ascii" + @test ForcingField(table_h5).format == "hdf5" + @test isabspath(ForcingField(table).path) + @test ForcingField(table; scale=3.0).scale == 3.0 # One coil set or several, as NamedTuples or Dicts with the [[ForcingTerms.coil_set]] keys. c_coil = (name="c", currents=[1e3, 500, -500, -1e3, -500, 500]) - coil_field = RMPField(c_coil; machine="d3d") - @test coil_field.ctrl.forcing_data_format == "coil" - @test coil_field.ctrl.machine == "d3d" - @test coil_field.ctrl.coil_sets_raw == [Dict{String,Any}("name" => "c", "currents" => c_coil.currents)] - @test length(RMPField([c_coil, Dict("name" => "iu", "currents" => fill(1e3, 6))]; machine="d3d").ctrl.coil_sets_raw) == 2 - @test repr(coil_field) == "RMPField(d3d coils c)" + coil_field = ForcingField(c_coil; machine="d3d") + @test coil_field isa GPEC.ForcingTerms.ForcingCoils + @test coil_field.config.machine == "d3d" + @test [sc.name for sc in coil_field.config.coil_sets] == ["c"] + @test coil_field.config.coil_sets[1].currents == c_coil.currents + @test length(ForcingField([c_coil, Dict("name" => "iu", "currents" => fill(1e3, 6))]; machine="d3d").config.coil_sets) == 2 + @test ForcingField(coil_field.config) isa GPEC.ForcingTerms.ForcingCoils + @test repr(coil_field) == "ForcingField(d3d coils c)" + # A source's [ForcingTerms] table builds the same source back; the bundled geometry directory stays implicit. + coil_table = GPEC._forcing_table(coil_field) + @test !haskey(coil_table, "dat_dir") + rebuilt = ForcingField(GPEC.forcing_terms_control(Dict{String,Any}("ForcingTerms" => coil_table), forcing_dir)) + @test rebuilt.config.machine == "d3d" && rebuilt.config.dat_dir == coil_field.config.dat_dir + @test [(sc.name, sc.currents) for sc in rebuilt.config.coil_sets] == [(sc.name, sc.currents) for sc in coil_field.config.coil_sets] + file_table = GPEC._forcing_table(ForcingField(table_h5)) + @test ForcingField(GPEC.forcing_terms_control(Dict{String,Any}("ForcingTerms" => file_table), forcing_dir)).path == table_h5 end - @testset "RMPField checks its source when built" begin - @test_throws "no forcing file" RMPField(joinpath(forcing_dir, "missing.dat")) + @testset "ForcingField checks its source when built" begin + @test_throws "no forcing file" ForcingField(joinpath(forcing_dir, "missing.dat")) geometry = joinpath(pkgdir(GPEC), "src", "ForcingTerms", "coil_geometries", "d3d_c.dat") - @test_throws "is not a forcing-mode table" RMPField(geometry) - @test_throws "is not a forcing-mode table" RMPField(table; format="hdf5") - @test_throws "is not \"ascii\", \"hdf5\" or \"coil\"" RMPField(table; format="csv") - @test_throws "has no currents" RMPField((name="c",); machine="d3d") - @test_throws "currents for its 6 conductors" RMPField((name="c", currents=ones(7)); machine="d3d") - @test_throws "pass machine" RMPField((name="c", currents=ones(6))) - @test_throws "unknown keys" RMPField((name="c", curents=ones(6)); machine="d3d") - @test_throws "a coil set is a Dict or NamedTuple" RMPField([1, 2]) - @test_throws "at least one forcing mode" RMPField(GPEC.ForcingTerms.ForcingMode[]) + @test_throws "is not a forcing-mode table" ForcingField(geometry) + @test_throws "is not a forcing-mode table" ForcingField(table; format="hdf5") + @test_throws "is not \"ascii\", \"hdf5\" or \"coil\"" ForcingField(table; format="csv") + @test_throws "has no currents" ForcingField((name="c",); machine="d3d") + @test_throws "currents for its 6 conductors" ForcingField((name="c", currents=ones(7)); machine="d3d") + @test_throws "pass machine" ForcingField((name="c", currents=ones(6))) + @test_throws "unknown keys" ForcingField((name="c", curents=ones(6)); machine="d3d") + @test_throws "a coil set is a Dict or NamedTuple" ForcingField([1, 2]) + @test_throws "at least one forcing mode" ForcingField(GPEC.ForcingTerms.ForcingMode[]) end @testset "rejected keyword combinations" begin @@ -512,10 +523,10 @@ include("display_check.jl") "PerturbedEquilibrium" => Dict{String,Any}(String(k) => v for (k, v) in pairs(pe_kwargs)))) reference = GPEC.main([dir]) resp = plasma_response(reference.ffs; compute_singular_coupling=true, verbose=false) - pe = perturbed_equilibrium(resp, RMPField(joinpath(dir, "forcing.h5")); verbose=false) + pe = perturbed_equilibrium(resp, ForcingField(joinpath(dir, "forcing.h5")); verbose=false) @test pe.response.permeability ≈ reference.pe.response.permeability rtol = 1e-10 - # The deck reads and normalizes its forcing file exactly as a fresh RMPField does. + # The deck reads and normalizes its forcing file exactly as a fresh ForcingField does. @test pe.forcing_b ≈ reference.pe.forcing_b rtol = 1e-10 end end diff --git a/test/runtests_torque_response.jl b/test/runtests_torque_response.jl index 8ceb87036f..542fd978f6 100644 --- a/test/runtests_torque_response.jl +++ b/test/runtests_torque_response.jl @@ -146,8 +146,8 @@ compute_torque_response = true # Build the ψ-resolved torque response @testset "forcing-source matrix: a forcing file is one source" begin file_state = deepcopy(pe) - file_ctrl = _GPE_TR.ForcingTerms.ForcingTermsControl(; forcing_data_file="forcing.dat") - _PE_TR.apply_torque_response!(file_state, file_ctrl) + file_source = _GPE_TR.ForcingTerms.ForcingFile(abspath("forcing.dat"), "ascii", 1) + _PE_TR.apply_torque_response!(file_state, file_source) ts = file_state.torque_response @test ts.source_names == ["forcing.dat"] @test size(ts.T_source) == (1, 1, npsi) @@ -157,7 +157,7 @@ compute_torque_response = true # Build the ψ-resolved torque response @testset "regularization on a kinetic result" begin # No regularization by default; an explicit width still reaches the regularized profiles. @test pe.regularization_width == 0.0 - reg = perturbed_equilibrium(resp, _GPE_TR.ForcingTerms.RMPField(ft_ctrl); regularization_width=0.05, + reg = perturbed_equilibrium(resp, _GPE_TR.ForcingTerms.ForcingField(ft_ctrl); regularization_width=0.05, coil_sets=pe.coil_sets, verbose=false) @test reg.regularization_width == 0.05 @test reg.forcing_b ≈ pe.forcing_b rtol = 1e-12 @@ -169,7 +169,7 @@ compute_torque_response = true # Build the ψ-resolved torque response # A response without its operators cannot carry a torque response. @test_throws ArgumentError _PE_TR.build_torque_response!(_PE_TR.PlasmaResponse(; ffs)) # The torque response is opt-in: a response built without it contracts nothing. - plain = perturbed_equilibrium(plasma_response(ffs; compute_singular_coupling=false, verbose=false), _GPE_TR.ForcingTerms.RMPField(ft_ctrl); + plain = perturbed_equilibrium(plasma_response(ffs; compute_singular_coupling=false, verbose=false), _GPE_TR.ForcingTerms.ForcingField(ft_ctrl); coil_sets=pe.coil_sets, verbose=false) @test isempty(plain.response.T_xe) && plain.torque_response === nothing end