ForcingTerms - FEATURE! - Refactor ForcingTerms naming and add capability to model magnetic materials at saturation magnetization in the direction of equilibrium fields

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

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
 

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
 

src/Analysis/CoilForcing.jl

@@ -232,34 +232,34 @@ 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
     return (; sets, skipped, complex_weight)
 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

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}()

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

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

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

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]