From 8663910e7580a3f4219cbe429db2bb9796df2e36 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 20:31:05 +0200 Subject: [PATCH 1/6] PerturbedEquilibrium - API! - Build the plasma response once and force it per source --- .pre-commit-config.yaml | 2 +- benchmarks/benchmark_against_fortran_run.jl | 2 - docs/development/architecture.md | 10 +- docs/development/toml-conventions.md | 2 +- docs/src/api.md | 16 +- docs/src/error_fields.md | 2 +- docs/src/perturbed_equilibrium.md | 13 +- .../DIIID-like_error_field_example/gpec.toml | 1 - .../gpec.toml | 13 +- examples/DIIID-like_ideal_example/gpec.toml | 1 - .../gpec.toml | 1 - examples/Solovev_ideal_example/gpec.toml | 1 - examples/Solovev_ideal_example_3D/gpec.toml | 1 - .../Solovev_kinetic_NTV_example/gpec.toml | 1 - .../gpec.toml | 2 - .../gpec.toml | 1 - .../cases/gal_resistive_pe.toml | 8 +- src/ForcingTerms/CoilFourier.jl | 2 +- src/ForcingTerms/ForcingTerms.jl | 35 ++- src/GeneralizedPerturbedEquilibrium.jl | 106 ++++---- .../FieldReconstruction.jl | 21 +- .../PerturbedEquilibrium.jl | 140 +++++----- .../PerturbedEquilibriumStructs.jl | 256 ++++++++---------- src/PerturbedEquilibrium/ResonantCoupling.jl | 54 ++-- src/PerturbedEquilibrium/Response.jl | 90 +++--- src/PerturbedEquilibrium/SingularCoupling.jl | 230 ++++++++-------- src/PerturbedEquilibrium/Utils.jl | 107 +++----- test/runtests_dominant_coupling.jl | 35 +-- test/runtests_equil.jl | 5 +- test/runtests_error_fields.jl | 8 +- test/runtests_kinetic_validity.jl | 22 +- test/runtests_rerun_from_h5.jl | 2 +- test/runtests_result_struct.jl | 5 +- test/runtests_solve_api.jl | 51 ++-- 34 files changed, 589 insertions(+), 657 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 697ae70a2..c41090f0b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -77,5 +77,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index a4de7bdcd..e55f8f3b3 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/benchmarks/benchmark_against_fortran_run.jl @@ -439,7 +439,6 @@ function write_gpec_toml_coil( end println(io) println(io, "[PerturbedEquilibrium]") - println(io, "compute_response = true") println(io, "compute_singular_coupling = true") println(io, "verbose = true") println(io, "write_outputs_to_HDF5 = true") @@ -490,7 +489,6 @@ function write_gpec_toml_file( println(io, "forcing_data_format = \"ascii\"") println(io) println(io, "[PerturbedEquilibrium]") - println(io, "compute_response = true") println(io, "compute_singular_coupling = true") println(io, "verbose = true") println(io, "write_outputs_to_HDF5 = true") diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 059539a8b..d38a044f7 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -214,12 +214,12 @@ The complete GPEC analysis pipeline: is built ### Perturbed Equilibrium -- `PerturbedEquilibriumControl` - User-facing TOML configuration parameters -- `PerturbedEquilibriumInternal` - Internal state with mode arrays -- `PerturbedEquilibriumState` - Results including: +- `PlasmaResponse` - Source-independent operators of one FFS result (`plasma_response`): inductances, + permeability and reluctance, the singular coupling matrices [n_rational × numpert_total] and their + dominant modes; applied to any number of external fields +- `PerturbedEquilibriumState` - The perturbed equilibrium driven by one field (`perturbed_equilibrium`), including: - Response fields (ξ_modes, b_modes) in mode space - - Singular coupling matrices [msing × numpert_total] - - Island diagnostics (half-widths, Chirikov parameters) + - Applied resonant vectors and island diagnostics (half-widths, Chirikov parameters) - `ForcingMode` - External forcing specification (m, n, amplitude, phase) ## Module Dependencies diff --git a/docs/development/toml-conventions.md b/docs/development/toml-conventions.md index 97139921a..43bd03a08 100644 --- a/docs/development/toml-conventions.md +++ b/docs/development/toml-conventions.md @@ -11,7 +11,7 @@ modify them. **Do not invent a new convention** — match what the existing exam 2. **Inline annotation on every variable line.** Every `key = value` gets a trailing `# description`. Source the wording from the matching config struct's docstring (`EquilibriumConfig`, `WallShapeSettings`, `ForceFreeStatesControl`, `ForcingTermsControl`, - `PerturbedEquilibriumControl`, `SolovevConfig`, `TJAnalyticConfig`, `KineticForcesControl`) + the `plasma_response`/`perturbed_equilibrium` docstrings, `SolovevConfig`, `TJAnalyticConfig`, `KineticForcesControl`) and keep it to one terse line. **Use the same description for the same variable across all files** — `examples/DIIID-like_ideal_example/gpec.toml` is the canonical reference for the common sections. Exception: in `regression-harness/cases/*.toml` the repeated diff --git a/docs/src/api.md b/docs/src/api.md index bf734f762..33e0eba1b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -11,7 +11,8 @@ using GeneralizedPerturbedEquilibrium eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") ffs = solve(eq, Riccati(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) rmp = RMPField("coils.dat") -pe = perturbed_equilibrium(ffs, rmp) +resp = plasma_response(ffs) +pe = perturbed_equilibrium(resp, rmp) ``` `solve(eq, alg; kwargs...)` is sugar for the canonical problem form: a @@ -25,7 +26,7 @@ prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true ffs = solve(prob, Riccati()) ``` -The four objects map one-to-one onto the pipeline stages: +The objects map one-to-one onto the pipeline stages: - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]` section. Analytic equilibria (`sol`, `lar`, `tj_analytic`) take their parameters from a @@ -35,8 +36,9 @@ The four objects map one-to-one onto the pipeline stages: a `ForceFreeStatesResult`, exactly the object the TOML driver publishes. - [`RMPField`](@ref) describes the external field — the `[ForcingTerms]` section — without reading anything from disk until it is applied. - - `perturbed_equilibrium` runs the plasma-response stage — the `[PerturbedEquilibrium]` - section. + - [`plasma_response`](@ref) builds the plasma-response and singular-coupling operators of a + `ForceFreeStatesResult` once, and `perturbed_equilibrium(resp, rmp)` drives them with one field — + together the `[PerturbedEquilibrium]` section. ## Combining forcing sources @@ -50,13 +52,15 @@ phase-rotates a source. ```julia nominal = RMPField("nominal_efc.dat") weld_field = RMPField("weld_fields.h5") -pe = perturbed_equilibrium(ffs, 2.0 * nominal + 0.5 * weld_field) +pe = perturbed_equilibrium(resp, 2.0 * nominal + 0.5 * weld_field) ``` This is the substrate for error-field workflows: build per-unit sources independently (a coil set per ampere, a displacement field per millimeter), then assemble physical cases by weighting and summing — without recomputing anything per combination until the final -`perturbed_equilibrium` call. +`perturbed_equilibrium` call, and with one `plasma_response` serving every combination. +`RMPField(modes)` wraps a control-surface spectrum computed elsewhere, already in the unit-norm +`Φ_x` convention, as a leaf that combines like any other. The weights are linear-combination coefficients, not physical amplitudes: sources that are not scalar multiples of each other get their own description and the algebra. A coil set diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index 1aa1e337f..66862d8c0 100644 --- a/docs/src/error_fields.md +++ b/docs/src/error_fields.md @@ -439,7 +439,7 @@ table.abs_delta_rim_per_mm # the same for tilt, as rim displacement at t table.tilt_sensitivity_per_deg[1, :] # ∂δ/∂θx per degree, per coil set # In memory, from the run's returned state (no I/O) -rc = PerturbedEquilibrium.ResonantCoupling(run.pe, run.ffs) +rc = PerturbedEquilibrium.ResonantCoupling(run.pe.response) dom = PerturbedEquilibrium.dominant_coupling(rc; psi_low=0.5) table = EF.sensitivity_table(run.coil_sensitivities, dom; mode=1) ``` diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index 5e079cfee..cd3bea745 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/docs/src/perturbed_equilibrium.md @@ -1,18 +1,21 @@ # Perturbed Equilibrium -The `PerturbedEquilibrium` module computes the plasma response to external magnetic perturbations. +The `PerturbedEquilibrium` module computes the plasma response to external magnetic perturbations, +in two stages: `plasma_response(ffs)` builds the source-independent response and singular-coupling +operators of a force-free-states result once, and `perturbed_equilibrium(resp, rmp)` drives them with +one external field. A scan over fields reuses one response. ## Types ```@docs -GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumControl -GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumInternal +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PlasmaResponse GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumState ``` ## Functions ```@docs +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_perturbed_equilibrium GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.write_outputs_to_HDF5 ``` @@ -53,7 +56,7 @@ the same way from a finished run in memory or from its `gpec.h5`: ```julia using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium -rc = ResonantCoupling("gpec.h5") # post hoc; or ResonantCoupling(pe_state, ffs) in memory +rc = ResonantCoupling("gpec.h5") # post hoc; or ResonantCoupling(resp) in memory dom = dominant_coupling(rc; psi_low=0.0, psi_high=0.9) # SVD over the surfaces in the window b̃ = rootarea_field(rc, coil_modes) # unit-norm forcing modes → root-area-weighted field @@ -78,7 +81,7 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.CORE_PSI_HIGH GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.rootarea_field GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coupling_overlap GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.check_mode_basis -GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_dominant_coupling! +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.build_dominant_coupling! ``` ## Plotting per-surface results against ψ or q diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index eee1c3e35..7a66f832d 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/examples/DIIID-like_error_field_example/gpec.toml @@ -196,7 +196,6 @@ height = -1.5780 # Hoop height: centroid Z of the pack cr currents = [55000.0] # 1 kA per turn × 55 turns; sensitivities scale linearly with current [A] [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = false # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index f291a598c..e27abe72d 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -1,9 +1,9 @@ # DIII-D-like H-mode DRIVEN resistive example: gal Δ′ + inner-layer matching, with DIII-D -# C-coil forcing (n=1 cosine phasing) staged for PerturbedEquilibrium. The response and the -# singular coupling are switched off: they require the free-boundary δW, which the Galerkin -# formalism does not yet produce, and PE refuses them on a Galerkin or matched result. The -# gal-matched resistive ξ and penetrated field are carried on the result for when gal-side -# δW lands. Same equilibrium and Galerkin settings as the sibling +# C-coil forcing (n=1 cosine phasing) staged for PerturbedEquilibrium. The run stops at the +# PerturbedEquilibrium stage for now: the plasma response needs the free-boundary δW, which the +# Galerkin formalism does not yet produce, and a follow-up makes it run. The gal-matched +# resistive ξ and penetrated field are carried on the result for when gal-side δW lands. +# Same equilibrium and Galerkin settings as the sibling # DIIID-like_gal_resistive_example; differences start at [ForcingTerms]. [Equilibrium] eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file @@ -93,8 +93,7 @@ name = "c" # DIII-D C-coil (6 coils); resolves to c currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing, 1 kA peak [A] [PerturbedEquilibrium] -compute_response = false # Compute plasma response to forcing -compute_singular_coupling = false # Compute singular layer coupling metrics +compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index cb9378287..b40d8a190 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -65,7 +65,6 @@ name = "c" # Coil set; resolves to coil_geometries/ currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing, 1 kA peak [A] [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index 898db615e..e3ab0ea88 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -67,7 +67,6 @@ name = "c" # Coil set; resolves to coil_geometries/ currents = [1800.0, -50.0, -450.0, -800.0, -450.0, -50.0] # Sum of n=1, 2, 3 cosine phasings with 1000, 500, 300 A peaks; n=3 lies outside nn_high [A] [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = false # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 87fb2ad36..5e299b068 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -36,7 +36,6 @@ forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_rea forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 411b33d98..c45ba99bd 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -57,7 +57,6 @@ equal_arc_wall = false # Equal arc length distribution of nodes # [PerturbedEquilibrium] # forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) # forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) -# compute_response = true # Compute plasma response to forcing # compute_singular_coupling = true # Compute singular layer coupling metrics # verbose = true # Enable verbose logging # write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index afd921754..ed53207a4 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -38,7 +38,6 @@ forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_rea forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 6e5bef877..21f81ef9a 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -33,11 +33,9 @@ forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_rea forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = false # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) [ForceFreeStates] local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index f66eeac15..00fa39ecd 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -38,7 +38,6 @@ forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_rea forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 diff --git a/regression-harness/cases/gal_resistive_pe.toml b/regression-harness/cases/gal_resistive_pe.toml index 3bc8b1ed1..04a8b6dd7 100644 --- a/regression-harness/cases/gal_resistive_pe.toml +++ b/regression-harness/cases/gal_resistive_pe.toml @@ -2,8 +2,8 @@ # fed into PerturbedEquilibrium with DIII-D C-coil forcing. Tracks the finite # per-rational-surface singular-coupling scalars (island half-widths, Chirikov # overlap, resonant flux, penetrated field, Δ′). -# Expected N/A: the example switches the PE response and coupling off until Galerkin and -# matched results publish the free-boundary δW, so none of these datasets are written yet. +# Expected N/A: the example's PerturbedEquilibrium stage stops until Galerkin and matched results +# publish the free-boundary δW, so the run ends after the matched solve (see [overrides]). # Each [quantities.*] block names an HDF5 path in the run output, how to extract it, # and the noise floor below which a difference is treated as zero. [case] @@ -11,6 +11,10 @@ name = "gal_resistive_pe" description = "DIII-D-like, n=1, DRIVEN/RPEC: gal-matched resistive ξ → PerturbedEquilibrium (coil-driven singular coupling)" example_dir = "examples/DIIID-like_gal_resistive_pe_example" +# Stop after force-free states so the run completes; the follow-up Galerkin-PE PR removes this. +[overrides] +"ForceFreeStates.force_termination" = true + # Full-pipeline coverage for gal_match_flag → gal_matched_odestate → PerturbedEquilibrium. Tracks the # FINITE, per-rational-surface singular-coupling scalars (the physics-meaningful driven quantities) — not # any whole-profile ξ/ξ′ norm (those diverge at the rationals and are ill-posed as benchmarks). diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 75b967640..6dc5eb462 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -374,7 +374,7 @@ coil sets on the plasma boundary — a [`CoilForcingGrid`](@ref) at `psi` follow [`coil_forcing_modes`](@ref) over the whole assembly in one Biot-Savart pass. Output amplitudes are in unit-norm convention (= Fortran `Phi_x`). -No normalization conversion is needed when using these modes with `compute_plasma_response!`. +No normalization conversion is needed when using these modes to drive a plasma response. Result is appended to `forcing_modes` (existing content is cleared first). """ diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index d5294c6ce..50d0eb3c4 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -267,12 +267,14 @@ several solves. RMPField(path; format=..., scale=1.0, kwargs...) RMPField(coil_sets::Vector{Dict{String,Any}}; scale=1.0, kwargs...) RMPField(ctrl::ForcingTermsControl; scale=1.0) + RMPField(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 TOML-shaped coil blocks (the `[[ForcingTerms.coil_set]]` layout) and selects the coil format. Remaining `kwargs` are -[`ForcingTermsControl`](@ref) fields. All three return a single-source leaf -([`RMPSource`](@ref)). +[`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)). ## Algebra @@ -318,20 +320,37 @@ struct RMPSource <: RMPField scale::ComplexF64 end +""" + RMPSpectrum + +A spectrum-literal [`RMPField`](@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. + +## 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). +""" +struct RMPSpectrum <: RMPField + modes::Vector{ForcingMode} + scale::ComplexF64 +end + """ RMPFieldSum -A lazy linear combination of [`RMPSource`](@ref) leaves, built by `+`/`-` on +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 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{RMPSource}` - The flattened weighted sources. + - `terms::Vector{Union{RMPSource,RMPSpectrum}}` - The flattened weighted leaves. """ struct RMPFieldSum <: RMPField - terms::Vector{RMPSource} + terms::Vector{Union{RMPSource,RMPSpectrum}} end """ @@ -353,14 +372,18 @@ function RMPField(coil_sets::Vector{Dict{String,Any}}; scale::Number=1.0, kwargs return RMPSource(ctrl, ComplexF64(scale)) end +RMPField(modes::Vector{ForcingMode}; scale::Number=1.0) = + RMPSpectrum([ForcingMode(; n=mode.n, m=mode.m, amplitude=mode.amplitude) for mode in modes], ComplexF64(scale)) + "Flattened weighted-leaf list of any [`RMPField`](@ref)." -_rmp_terms(f::RMPSource) = [f] +_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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 236741200..af4277b16 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -91,6 +91,7 @@ using .InnerLayer: GGJModel, SLAYERModel using .Tearing.Runner: TearingProblem using .Equilibrium: PlasmaEquilibrium using .ForcingTerms: RMPField +using .PerturbedEquilibrium: plasma_response, PlasmaResponse const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim", "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi", @@ -98,7 +99,7 @@ const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_ const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc", "use_galgrid") const _DEPRECATED_KF_KEYS = ("moment", "eq_out", "theta_out", "xlmda_out", "eqpsi_out", "save_records", "gpec_file", "fnml_flag", "ellip_flag", "diag_flag", "force_xialpha") -const _DEPRECATED_PE_KEYS = ("fixed_boundary", "output_eigenmodes", "filter_modes", "singular_point_method") +const _DEPRECATED_PE_KEYS = ("fixed_boundary", "output_eigenmodes", "filter_modes", "singular_point_method", "compute_response") # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated. function _stop_on_deprecated_keys(table, deprecated_keys, section::String) @@ -987,12 +988,14 @@ function run_perturbed_equilibrium( pe_state = nothing if "PerturbedEquilibrium" in keys(inputs) ft_ctrl = forcing_terms_control(inputs) + section = inputs["PerturbedEquilibrium"] - # 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); + # The deck's section feeds the two scripting calls, and its forcing block is an unscaled + # RMPField, so the TOML path and the scripting API share one stage. + resp = plasma_response(result; (Symbol(k) => v for (k, v) in section if k in ("compute_singular_coupling", "verbose"))...) + pe_state = perturbed_equilibrium(resp, ForcingTerms.RMPField(ft_ctrl); coil_sets=preloaded_coil_sets, runtimes=runtimes, - (Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])...) + (Symbol(k) => v for (k, v) in section if k != "compute_singular_coupling")...) end @info "Perturbed Equilibrium completed in $(@sprintf("%.3f", time() - pe_start)) s" @@ -1001,59 +1004,54 @@ function run_perturbed_equilibrium( end """ - perturbed_equilibrium(ffs, rmp; coil_sets=nothing, runtimes=nothing, kwargs...) -> PerturbedEquilibriumState + perturbed_equilibrium(resp, rmp; regularization_width=resp.regularization_width, verbose=true, + write_outputs_to_HDF5=true, output_filename="", coil_sets=nothing, runtimes=nothing) + -> 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 -arguments are `PerturbedEquilibrium.PerturbedEquilibriumControl` fields. +Drive the plasma response `resp` (from [`plasma_response`](@ref)) with the external field `rmp`, an +[`RMPField`](@ref), and write the perturbed-equilibrium outputs. One response serves any number of +fields, so a scan rebuilds no operator. -Products the producing integrator could not supply gate the corresponding calculation: the -step warns and is skipped rather than erroring, so a Riccati-fed result still flows through. -Galerkin and `:matched` results publish no free-boundary δW yet and are refused when the -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`. +`regularization_width` smooths the ideal 1/(m−nq) divergence in the profile reconstruction (Fortran +`reg_spot`); its default is 0 for a kinetic result and 0.05 otherwise. `write_outputs_to_HDF5` appends +the outputs to `output_filename`, or to the force-free-states file when that is empty. `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`. ```julia -pe = perturbed_equilibrium(ffs, RMPField("forcing.dat")) +resp = plasma_response(ffs) +pe = perturbed_equilibrium(resp, RMPField("forcing.dat")) ``` """ function perturbed_equilibrium( - ffs::ForceFreeStatesResult, + resp::PerturbedEquilibrium.PlasmaResponse, rmp::ForcingTerms.RMPField; + regularization_width::Real=resp.regularization_width, + verbose::Bool=true, + write_outputs_to_HDF5::Bool=true, + output_filename::String="", coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing, - runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing, - kwargs... + runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing ) + ffs = resp.ffs ctrl = ffs.control - pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; PerturbedEquilibrium.kinetic_regularization_kwargs(ffs, kwargs)...) - pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=ffs.dir_path) - - # Inner-layer penetrated resonant field; zeros under ideal closure. - pe_intr.inner_bpen = ffs.bpen - # Injected coil geometry (gpec.h5 replay with `--coil-source coils`) lets the coil field # be recomputed from stored geometry without the .dat/.h5 files. - if coil_sets !== nothing - 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(resp, rmp; regularization_width, verbose, + preloaded_coil_sets=coil_sets === nothing ? ForcingTerms.CoilSet[] : copy(coil_sets), runtimes) # Write perturbed equilibrium outputs to same HDF5 file - if pe_ctrl.write_outputs_to_HDF5 - output_file = isempty(pe_ctrl.output_filename) ? ctrl.HDF5_filename : pe_ctrl.output_filename - PerturbedEquilibrium.write_outputs_to_HDF5( - pe_state, pe_intr, joinpath(ffs.dir_path, output_file) - ) + if write_outputs_to_HDF5 + output_file = isempty(output_filename) ? ctrl.HDF5_filename : output_filename + PerturbedEquilibrium.write_outputs_to_HDF5(pe_state, joinpath(ffs.dir_path, output_file)) @info "Results written to $output_file" end # Snapshot the coil geometry actually used into the gpec.h5 output so the run # is replayable from the output file alone (see `main_from_h5 --coil-source coils`). - if ctrl.write_outputs_to_HDF5 && !isempty(pe_intr.coil_sets) - _write_coil_snapshot!(joinpath(ffs.dir_path, ctrl.HDF5_filename), pe_intr.coil_sets) + if ctrl.write_outputs_to_HDF5 && !isempty(pe_state.coil_sets) + _write_coil_snapshot!(joinpath(ffs.dir_path, ctrl.HDF5_filename), pe_state.coil_sets) end return pe_state @@ -1201,7 +1199,7 @@ function run_error_fields( cfg = ForcingTerms.CoilConfig(ft_ctrl) coil_sets = preloaded_coil_sets === nothing ? ForcingTerms.load_coil_sets(cfg, result.nlow; equil=result.equil) : preloaded_coil_sets - rc = PerturbedEquilibrium.ResonantCoupling(pe_state, result) + rc = PerturbedEquilibrium.ResonantCoupling(pe_state.response) sens = ErrorFields.compute_coil_sensitivities(coil_sets, rc, result.equil, cfg, ef_ctrl; psi=result.psilim, b_t0=result.equil.params.bt0) @@ -1256,7 +1254,8 @@ function run_error_fields( missing = setdiff(ntv_ctrl.efc_coils, [cs.name for cs in efc_sets]) isempty(missing) || error("[ErrorFields.NTV] efc_coils not among the run's coil sets: $(join(missing, ", "))") ntv_start = time() - couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, + couplings = efc_couplings(pe_state.response, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, + regularization_width=pe_state.regularization_width, rotation_scan=ntv_ctrl.rotation_scan, scan_points=ntv_ctrl.rotation_scan_points, scan_max_points=ntv_ctrl.rotation_scan_max_points, scan_tolerance=ntv_ctrl.rotation_scan_tolerance, span_factor=ntv_ctrl.rotation_span_factor, offset_factor=ntv_ctrl.rotation_offset_factor) @info "NTV couplings of $(length(couplings)) correction arrays in $(@sprintf("%.1f", time() - ntv_start)) s" @@ -1281,22 +1280,23 @@ function run_error_fields( end """ - efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} + efc_couplings(resp, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; mode=1, method="fgar", + regularization_width=resp.regularization_width, verbose=false) -> Vector{ErrorFields.EFCCoupling} The couplings of each correction coil array in `coil_sets` per kilo-ampere-turn: its dominant-mode overlap and resonant fraction from its spectrum on the control surface, and its NTV torque for the whole field and for the field with mode `mode` of `dom` projected out. The -torque is a quadratic form of the applied spectrum, so each is one plasma-response evaluation -of the unit-current spectrum (injected as forcing modes, nothing written) followed by the -kinetic torque of `method`; both scale exactly with the square of the current. `cfg` is the -run's `[ForcingTerms]` coil configuration, so the boundary grid matches the one the run's -sensitivities were evaluated on; `rc` must be the coupling of `ffs`'s own solve and +torque is a quadratic form of the applied spectrum, so each is one evaluation of the run's plasma +response `resp` on the unit-current spectrum (a spectrum-literal `RMPField`, nothing written, with the +run's `regularization_width`) followed by the kinetic torque of `method`; both scale exactly with the +square of the current. `cfg` is the run's `[ForcingTerms]` coil configuration, so the boundary grid +matches the one the run's sensitivities were evaluated on; `rc` must be the coupling of `resp` and `kinetic_profiles` its kinetic context. The spectrum is normalized by the magnitude of the array's ampere-turns, `|nw| × max|I|`, keeping its winding sense and current pattern as the phase reference; the torque is stored with its sign. """ function efc_couplings( - ffs::ForceFreeStatesResult, + resp::PerturbedEquilibrium.PlasmaResponse, coil_sets::Vector{ForcingTerms.CoilSet}, rc::PerturbedEquilibrium.ResonantCoupling, dom::PerturbedEquilibrium.DominantCoupling, @@ -1311,15 +1311,16 @@ function efc_couplings( scan_tolerance::Real=0.05, span_factor::Real=1.0, offset_factor::Real=2.0, + regularization_width::Real=resp.regularization_width, verbose::Bool=false ) + ffs = resp.ffs kinetic_profiles === nothing && error("efc_couplings needs kinetic profiles: add a [KineticForces] section with a kinetic_file") getfield(kf_ctrl, Symbol(method * "_flag")) || error("efc_couplings: KineticForces method \"$method\" is not enabled in the control") equil = ffs.equil grids = [(n, ForcingTerms.CoilForcingGrid(equil, cfg, n; psi=ffs.psilim)) for n in sort(unique(rc.n_modes))] m_low, m_high = extrema(rc.m_modes) v = dom.right_singular_vectors[:, mode] - quiet = PerturbedEquilibrium.PerturbedEquilibriumControl(; compute_response=true, compute_singular_coupling=false, verbose=verbose, write_outputs_to_HDF5=false) # Reference rotation and scan span from the kinetic profiles, with the torque kernel's own # diamagnetic frequencies (Logan & Park 2013 Eq. 7): ω_φ = ω_E + ω_*n + ω_*T. @@ -1341,10 +1342,7 @@ function efc_couplings( "scan ±$(@sprintf("%.3e", span)) rad/s on $(length(grid0)) initial points (max $(scan_max_points))" function response_of(modes::Vector{ForcingTerms.ForcingMode}) - pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=ffs.dir_path) - pe_intr.inner_bpen = ffs.bpen - pe_intr.forcing_modes = modes - pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(ffs, ForcingTerms.RMPField(ForcingTerms.ForcingTermsControl()), quiet, pe_intr) + pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(resp, ForcingTerms.RMPField(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 @@ -1903,8 +1901,8 @@ function write_imas(dd, result) end export main, write_imas -export solve, perturbed_equilibrium -export PlasmaEquilibrium, EulerLagrangeProblem, Forward, Riccati, Galerkin, ForceFreeStatesResult, RMPField +export solve, plasma_response, perturbed_equilibrium +export PlasmaEquilibrium, EulerLagrangeProblem, Forward, Riccati, Galerkin, ForceFreeStatesResult, RMPField, PlasmaResponse export MatchProblem, TearingProblem, GGJModel, SLAYERModel include("Precompile.jl") diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index c4f073bc1..18e634ed6 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -45,7 +45,7 @@ _mode_mn(ffs::ForceFreeStatesResult, j::Int) = (ffs.mlow + (j - 1) % ffs.mpert, """ reconstruct_physical_fields( response_vector, flux_matrix, solution, - equil, ffs, intr, metric, mats, ctrl + equil, ffs, metric, mats, regularization_width ) -> (xi_modes, b_modes) Reconstruct displacement and perturbed magnetic field from eigenmode response. @@ -80,10 +80,9 @@ function reconstruct_physical_fields( solution::SolutionProfiles, equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, - intr::PerturbedEquilibriumInternal, metric::MetricData, mats::MatrixSplines, - ctrl::PerturbedEquilibriumControl + regularization_width::Float64 ) npsi = size(solution.u_store, 4) psi_grid = solution.psi_store[1:npsi] @@ -114,7 +113,7 @@ function reconstruct_physical_fields( # Compute Clebsch displacements with regularization (matches Fortran gpeq_sol + gpout_xclebsch) clebsch_psi, clebsch_psi1, clebsch_alpha = compute_clebsch_displacements( xi_psi_modes, xi_psi1_modes, xi_s_modes, - psi_grid, equil, ffs, mats, ctrl + psi_grid, equil, ffs, mats, regularization_width ) # Compute regularized (modified) b-field components (matches Fortran gpeq_sol bmt/bmz) @@ -126,7 +125,7 @@ function reconstruct_physical_fields( # Compute contravariant displacement via Jacobian convolution (matches Fortran gpeq_contra) xwp_modes, xwt_modes, xwz_modes, xmt_modes, xmz_modes = compute_contra_displacements( xi_psi_modes, clebsch_psi1, clebsch_alpha, - psi_grid, equil, ffs, metric, ctrl + psi_grid, equil, ffs, metric, regularization_width ) # Compute covariant components via metric tensor contraction (matches Fortran gpeq_cova) @@ -327,7 +326,7 @@ end """ compute_clebsch_displacements( xi_psi_modes, xi_psi1_modes, xi_s_modes, - psi_grid, equil, ffs, mats, ctrl + psi_grid, equil, ffs, mats, regularization_width ) -> (clebsch_psi, clebsch_psi1, clebsch_alpha) Compute Clebsch displacement components for PENTRC output. @@ -351,7 +350,7 @@ function compute_clebsch_displacements( equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, mats::MatrixSplines, - ctrl::PerturbedEquilibriumControl + regularization_width::Float64 ) npsi, mpert = size(xi_psi_modes) chi1 = 2π * equil.psio @@ -361,7 +360,7 @@ function compute_clebsch_displacements( clebsch_psi1 = copy(xi_psi1_modes) # will be regularized below clebsch_alpha = xi_s_modes ./ chi1 # ξ^α/χ₁ (will be regularized below) - reg_spot = ctrl.regularization_width # Fortran gpec reg_spot + reg_spot = regularization_width # Fortran gpec reg_spot @assert reg_spot >= 0 "regularization_width must be non-negative (got $reg_spot)" if reg_spot == 0 @@ -484,7 +483,7 @@ end """ compute_contra_displacements( xi_psi_modes, clebsch_psi1, clebsch_alpha, - psi_grid, equil, ffs, metric, ctrl + psi_grid, equil, ffs, metric, regularization_width ) -> (xwp_modes, xwt_modes, xwz_modes, xmt_modes, xmz_modes) Compute contravariant displacement via Jacobian mode coupling convolution. @@ -507,12 +506,12 @@ function compute_contra_displacements( equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, metric::MetricData, - ctrl::PerturbedEquilibriumControl + regularization_width::Float64 ) npsi, numpert_total = size(xi_psi_modes) (; mpert, mlow, nlow, npert) = ffs chi1 = 2π * equil.psio - reg_spot = ctrl.regularization_width # Fortran gpec reg_spot + reg_spot = regularization_width # Fortran gpec reg_spot fc = metric.fourier_coeffs xwp_modes = zeros(ComplexF64, npsi, numpert_total) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 9df339eac..c757f9cec 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -28,95 +28,91 @@ include("ResonantCoupling.jl") include("Utils.jl") # Export main types -export PerturbedEquilibriumControl -export PerturbedEquilibriumInternal +export PlasmaResponse export PerturbedEquilibriumState export ForcingMode # Export main functions +export plasma_response export compute_perturbed_equilibrium export write_outputs_to_HDF5 export ResonantCoupling, DominantCoupling, dominant_coupling, rootarea_field, coupling_overlap, check_mode_basis, CORE_PSI_LOW, CORE_PSI_HIGH """ - compute_perturbed_equilibrium(ffs, forcing, ctrl, intr; runtimes=nothing)::PerturbedEquilibriumState + plasma_response(ffs; compute_singular_coupling=true, verbose=true) -> PlasmaResponse -Main entry point for perturbed equilibrium calculations. +Build the source-independent half of a perturbed equilibrium from a force-free-states result: the +control-surface response operators (inductances, permeability, reluctance; Park et al., Phys. +Plasmas 14, 052110 (2007)) and, with `compute_singular_coupling`, the singular-coupling matrices and +their dominant modes (Park, Phys. Plasmas 16, 056115 (2009)). Apply it to an external field with +`perturbed_equilibrium(resp, rmp)` as many times as needed. -Computes plasma response to external forcing and calculates singular layer -coupling metrics. Every ForceFreeStates input — the equilibrium, the mode space, the -metric and matrix fits, the free-boundary energies and the ξ solution — is read off `ffs`. -Products the producing integrator could not supply gate the corresponding calculation: -the step warns and is skipped instead of erroring. Galerkin and `:matched` results publish -no free-boundary δW yet, so asking for the response or the coupling on one is an error. - -## Arguments - - - `ffs`: `ForceFreeStates.ForceFreeStatesResult` from the stability solve - - `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 - - - `PerturbedEquilibriumState`: Calculation results +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 +no free-boundary δW yet, so building a response from one is an error. """ -function compute_perturbed_equilibrium( - ffs::ForceFreeStatesResult, - forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}, - ctrl::PerturbedEquilibriumControl, - intr::PerturbedEquilibriumInternal; - runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing -)::PerturbedEquilibriumState - - state = PerturbedEquilibriumState() - equil = ffs.equil - mats = ffs.mats - mtheta_vacuum = ffs.control.mtheta_vacuum - - # Galerkin and matched results publish no plasma edge response wp, so neither step could run. - if (ctrl.compute_response || ctrl.compute_singular_coupling) && ffs.free_boundary === nothing && - (ffs.integrator === :galerkin || ffs.closure === :matched) - error("PerturbedEquilibrium needs the free-boundary δW, which a $(ffs.integrator) result with $(ffs.closure) closure " * - "does not publish yet; set compute_response = false and compute_singular_coupling = false, or use a forward solve") +function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling::Bool=true, verbose::Bool=true)::PlasmaResponse + # Galerkin and matched results publish no plasma edge response wp, so no operator could be built. + if ffs.free_boundary === nothing && (ffs.integrator === :galerkin || ffs.closure === :matched) + error("plasma_response needs the free-boundary δW, which a $(ffs.integrator) result with $(ffs.closure) closure " * + "does not publish yet; use a forward solve") end - # Mode arrays for convenient indexing - initialize_mode_arrays!(intr, ffs) - + m_modes, n_modes = mode_arrays(ffs) # A run has at most one ξ solution, already closed at the rationals and carrying populated # Ξ′ / Ξ_s stores. The gal-native basis takes the analytic Galerkin Ξ′ downstream. solution = ffs.solution - intr.odet_from_gal = solution !== nothing && solution.basis === :gal_native - - # 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 + resp = PlasmaResponse(; ffs, m_modes, n_modes, regularization_width=default_regularization_width(ffs), + odet_from_gal=solution !== nothing && solution.basis === :gal_native) - # Plasma response - if ctrl.compute_response && - ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") && + if ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") && ForceFreeStates.require_solution(ffs, "plasma response calculation") - compute_plasma_response!(state, equil, solution, ffs.free_boundary.wt0, mtheta_vacuum, ffs, intr, ctrl, ffs.metric, mats) + build_response_operators!(resp; verbose) + if compute_singular_coupling + build_coupling_operators!(resp; verbose) + build_dominant_coupling!(resp; verbose) + end end + return resp +end - # Singular coupling metrics - if ctrl.compute_singular_coupling && - ForceFreeStates.require(ffs, :free_boundary, "singular coupling calculation") && - ForceFreeStates.require_solution(ffs, "singular coupling calculation") - compute_singular_coupling_metrics!(state, equil, solution, mtheta_vacuum, ffs, intr, ctrl, mats) - compute_dominant_coupling!(state, ctrl) - end +""" + compute_perturbed_equilibrium(resp, forcing; regularization_width=resp.regularization_width, verbose=true, + preloaded_coil_sets=CoilSet[], runtimes=nothing)::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 +singular coupling to it and reconstruct the perturbed profiles. Steps whose operators `resp` does not +carry are skipped, leaving only the forcing spectrum. + +## Arguments + - `regularization_width`: width of the ideal 1/(m−nq) smoothing in the profile reconstruction + (Fortran `reg_spot`); must be ≥ 0 + - `preloaded_coil_sets`: coil geometry to integrate instead of reading it (the gpec.h5 replay path) + - `runtimes`: optional collector; receives `"forcing_terms" => dt` for the forcing-mode + materialization (coil Biot-Savart onto the plasma surface) +""" +function compute_perturbed_equilibrium( + resp::PlasmaResponse, + forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}; + regularization_width::Real=resp.regularization_width, + verbose::Bool=true, + preloaded_coil_sets::Vector{CoilSet}=CoilSet[], + runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing +)::PerturbedEquilibriumState + ffs = resp.ffs + ft_start = time() + modes, coil_sets = materialize_forcing_modes(ffs, forcing; dir_path=ffs.dir_path, preloaded_coil_sets, verbose) + runtimes === nothing || push!(runtimes, "forcing_terms" => time() - ft_start) + state = PerturbedEquilibriumState(; response=resp, 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) + isempty(resp.coupling_flux) || apply_coupling!(state; verbose) + if !isempty(resp.dominant_right_singular_vectors) && !isempty(state.forcing_b_rootarea) + state.dominant_forcing_overlap = resp.dominant_right_singular_vectors' * state.forcing_b_rootarea + end return state end @@ -190,6 +186,16 @@ function materialize_forcing_modes(ffs::ForceFreeStatesResult, rmp::ForcingTerms return modes, coil_sets end +""" + materialize_forcing_modes(ffs, rmp::ForcingTerms.RMPSpectrum; 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[] +end + """ materialize_forcing_modes(ffs, rmp::ForcingTerms.RMPFieldSum; kwargs...) -> (modes, coil_sets) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 3b0815439..a48078223 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -1,120 +1,151 @@ """ - PerturbedEquilibriumControl + PlasmaResponse -User-facing control parameters from TOML [PerturbedEquilibrium] section. +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, +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. ## Fields -Note: Forcing data file settings are now in [ForcingTerms] section. + - `ffs::ForceFreeStatesResult` - The result the operators were built from. + - `m_modes`, `n_modes` - Poloidal and toroidal mode number of each entry of a spectrum [numpert_total]. + - `regularization_width::Float64` - Default profile-reconstruction width (Fortran `reg_spot`): 0 for a + kinetic result, whose Euler-Lagrange operator has no ideal 1/(m−nq) divergence to smooth (Park & Logan, + Phys. Plasmas 24, 032505 (2017) §III D), 0.05 otherwise. + - `odet_from_gal::Bool` - The solution is the gal-native matched basis, whose `du_store` carries the + analytic Galerkin Ξ′; selects that branch of the singular-coupling Ξ′ evaluation. + - `flux_matrix::Matrix{ComplexF64}` - Control-surface flux of the free-boundary eigenmodes [mode × eigenmode]. + - `plasma_inductance_flux`, `surface_inductance_flux`, `permeability_flux` - Λ, L and P = Λ·L⁻¹ acting + on the physical control-surface flux Φ_x; the applied quantities and energies are evaluated with these. -Calculation: - - - `compute_response::Bool` - Compute plasma response (default: true) - - `compute_singular_coupling::Bool` - Compute singular coupling metrics (default: true) - - `verbose::Bool` - Enable verbose logging (default: true) - -Output Settings: +Control surface matrices [numpert_total × numpert_total], stored in the coordinate-invariant +root-area-weighted field (b̃) space (Pharr 2026). Writing `S ≡ rootarea_to_area_weight` +(b̃→b̄) and `A ≡ surface_area`, the brief internal flux-conform operator is `R = S·A` (Φ = R·b̃). +Recover the area-weighted (b̄) forms by conforming with `S` (e.g. `L_b̄ = S·L̃·S†`); recover flux +with `Φ = A·b̄`. - - `output_filename::String` - Combined output file with ForceFreeStates results (default: uses ForceFreeStates HDF5_filename) - - `write_outputs_to_HDF5::Bool` - Write perturbed equilibrium outputs to HDF5 (default: true) + - `plasma_inductance` - Λ̃ = R⁻¹·Λ·R⁻† (wt0-based plasma inductance, congruence) + - `surface_inductance` - L̃ = R⁻¹·L·R⁻† (vacuum surface inductance, congruence) + - `permeability` - P̃ = R⁻¹·P·R (plasma response operator P=Λ·L⁻¹, similarity) + - `reluctance` - ϱ̃ = R†·ϱ·R (ϱ = L⁻¹·(Λ−L)·L⁻¹, congruence) + - `rootarea_to_area_weight` - S = Σ/√A at psilim (b̃→b̄ recovery operator) + - `surface_area` - scalar A = ∫J|∇ψ|dθ at psilim (flux recovery Φ = A·b̄) -Regularization: +Coupling matrices [n_rational × numpert_total] — one row per resonant (surface, n) pair, mapping the +applied root-area-weighted field b̃_x to the resonant response at that surface. Matches Fortran +`C_f_x_out`, `C_i_x_out`, etc. (shape [mode_C, m_out]). - - `regularization_width::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 - to disable. Must be ≥ 0. Forced to 0 in kinetic runs, whose Euler-Lagrange operator has no - resonant singularity to smooth. -""" -@kwdef struct PerturbedEquilibriumControl - # Calculation - compute_response::Bool = true - compute_singular_coupling::Bool = true - verbose::Bool = true - - # Output settings - output_filename::String = "" # Empty means use ForceFreeStates filename - write_outputs_to_HDF5::Bool = true - - # Regularization width for singular surface smoothing (matches Fortran gpec.f reg_spot). - # Set to 0 to disable regularization. Must be non-negative. - regularization_width::Float64 = 5e-2 -end + - `C_resonant_area_weighted_field` - Φ_r/A^r coupling (resonant area-weighted field b^r in tesla; singcoup row 1) [Pharr 2026] + - `C_resonant_current` - Resonant current coupling (singcoup row 2) + - `C_island_width_sq` - (w/2)² coupling (singcoup row 3) + - `C_penetrated_area_weighted_field` - Penetrated area-weighted field coupling (singcoup row 4) + - `C_delta_prime` - Δ' coupling (singcoup row 5) + - `coupling_flux::NamedTuple` - The same five matrices acting on the physical flux Φ_x (keys without the + `C_` prefix), used for the applied resonant vectors; empty when the coupling was not built. + - `C_coeffs::Matrix{ComplexF64}` - Edge-to-solution coefficients `u_bnd⁻¹·edge_mn` of each unit forcing mode. -""" - PerturbedEquilibriumInternal +Metadata [n_rational] — identifies each (surface, n) row: -Internal state variables for perturbed equilibrium calculations. + - `rational_psi`, `rational_q`, `rational_m_res`, `rational_n`, `rational_surface_idx`, `rational_area` -## Fields +Dominant resonant-coupling modes — the run summary SVD `U·diag(σ)·Vᴴ` of +`C_resonant_area_weighted_field` over the core window (window post hoc with `dominant_coupling` +on a `ResonantCoupling`): - - `dir_path::String` - Working directory path - - `forcing_modes::Vector{ForcingMode}` - Loaded forcing mode data - - `coil_sets::Vector{CoilSet}` - Coil geometry used (when `forcing_data_format = "coil"`); - captured for the gpec.h5 rerun snapshot - - `plasma_response::Matrix{ComplexF64}` - Plasma response matrix - - `m_modes::Vector{Int}` - Poloidal mode numbers for each index i in 1:numpert_total - - `n_modes::Vector{Int}` - Toroidal mode numbers for each index i in 1:numpert_total + - `dominant_singular_values` - σ, descending [rank] + - `dominant_right_singular_vectors` - V, applied-b̃ spectra ranked by resonant drive [numpert_total × rank] + - `dominant_left_singular_vectors` - U, resonant-field patterns over the retained surfaces [n_retained × rank] + - `dominant_rational_index` - rows of the `rational_*` arrays retained in the window [n_retained] """ -@kwdef mutable struct PerturbedEquilibriumInternal - dir_path::String = "./" - forcing_modes::Vector{ForcingMode} = ForcingMode[] - coil_sets::Vector{CoilSet} = CoilSet[] - plasma_response::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) +@kwdef mutable struct PlasmaResponse + ffs::ForceFreeStatesResult m_modes::Vector{Int} = Int[] n_modes::Vector{Int} = Int[] - # ForceFreeStates-provided B_pen per (match surface × coil-drive column) from inner layer. - inner_bpen::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - # True when the consumed solution came from gal matching, whose du_store carries the - # analytic galerkin Ξ′; selects the gal branch of the singular-coupling Ξ′ evaluation. + regularization_width::Float64 = 5e-2 odet_from_gal::Bool = false + + # Flux-space operators [numpert_total × numpert_total] + flux_matrix::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + plasma_inductance_flux::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + surface_inductance_flux::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + permeability_flux::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + + # Control surface matrices, root-area-weighted field (b̃) space + plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + surface_area::Float64 = 0.0 + + # Coupling matrices [n_rational × numpert_total], b̃ input basis and flux input basis + C_resonant_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_resonant_current::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_island_width_sq::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_penetrated_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_delta_prime::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + coupling_flux::NamedTuple = (;) + C_coeffs::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + + # Metadata [n_rational] + rational_psi::Vector{Float64} = Float64[] + rational_q::Vector{Float64} = Float64[] + rational_m_res::Vector{Int} = Int[] + rational_n::Vector{Int} = Int[] + rational_surface_idx::Vector{Int} = Int[] + rational_area::Vector{Float64} = Float64[] + + # Dominant resonant-coupling modes + dominant_singular_values::Vector{Float64} = Float64[] + dominant_right_singular_vectors::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + dominant_left_singular_vectors::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + dominant_rational_index::Vector{Int} = Int[] end +""" + default_regularization_width(ffs) -> Float64 + +Profile-reconstruction width for a result: 0 when it carries a kinetic fit, which removes the ideal +resonant singularity, 0.05 otherwise. `ffs` is duck-typed: only `ffs.mats.kinetic` is read. +""" +default_regularization_width(ffs) = ffs.mats.kinetic === nothing ? 5e-2 : 0.0 + """ PerturbedEquilibriumState -Results from perturbed equilibrium calculations. +The perturbed equilibrium driven by one external field: the forcing materialized against a +[`PlasmaResponse`](@ref) and the response, profiles and resonant quantities it produces. ## Fields + - `response::PlasmaResponse` - The operators this field was applied to. + - `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. + Response fields (mode space): + - `psi_grid` - FFS ODE integration ψ_n values [npsi] - `xi_modes::Union{Nothing, NamedTuple}` - Displacement (psi, theta, zeta) [npsi, numpert_total] - `b_modes::Union{Nothing, NamedTuple}` - Magnetic field; psi=b^ψ, b_psi_area_weighted=b^ψ/⟨J·|∇ψ|⟩_θ, theta/zeta=unregularized, theta_reg/zeta_reg=regularized [npsi, numpert_total] - `b_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal field b_n [npsi, numpert_total] - `xi_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal displacement xi_n [npsi, numpert_total] -Coupling matrices [n_rational × numpert_total] — one row per resonant (surface, n) pair. -Each row maps the full applied field to the resonant response at that surface. -Matches Fortran `C_f_x_out`, `C_i_x_out`, etc. (shape [mode_C, m_out]). - - - `C_resonant_area_weighted_field` - Φ_r/A^r coupling (resonant area-weighted field b^r in tesla; singcoup row 1) [Pharr 2026] - - `C_resonant_current` - Resonant current coupling (singcoup row 2) - - `C_island_width_sq` - (w/2)² coupling (singcoup row 3) - - `C_penetrated_area_weighted_field` - Penetrated area-weighted field coupling (singcoup row 4) - - `C_delta_prime` - Δ' coupling (singcoup row 5) - -Applied resonant vectors [n_rational] = C · forcing_amplitudes. +Applied resonant vectors [n_rational] = C · forcing, on the rows of `response`. Matches Fortran `Phi_res`, `w_isl`, `K_isl`, `Delta`. - `resonant_area_weighted_field`, `resonant_current`, `island_width_sq`, `penetrated_area_weighted_field`, `delta_prime` + - `forcing_solution_weights` - `C_coeffs`·Φ_x, the forcing's weights on the edge-driven solutions Diagnostics [n_rational]: - `island_half_width::Vector{Float64}` - w/2 = sqrt(|island_width_sq|) from applied forcing - `chirikov_parameter::Vector{Float64}` - Island overlap metric - -Metadata [n_rational] — identifies each (surface, n) row: - - - `rational_psi`, `rational_q`, `rational_m_res`, `rational_n`, `rational_surface_idx` - -Dominant resonant-coupling modes — the run summary SVD `U·diag(σ)·Vᴴ` of -`C_resonant_area_weighted_field` over every resonant row (window post hoc with `dominant_coupling` -on a `ResonantCoupling`): - - - `dominant_singular_values` - σ, descending [rank] - - `dominant_right_singular_vectors` - V, applied-b̃ spectra ranked by resonant drive [numpert_total × rank] - - `dominant_left_singular_vectors` - U, resonant-field patterns over the retained surfaces [n_retained × rank] - - `dominant_rational_index` - rows of the `rational_*` arrays retained in the window [n_retained] - - `dominant_forcing_overlap` - Vᴴ·b̃_x, the applied forcing's coefficient on each mode [rank] + - `dominant_forcing_overlap` - Vᴴ·b̃_x, the applied forcing's coefficient on each dominant mode [rank] Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field representations (all tesla; no flux/weber is stored): @@ -123,19 +154,6 @@ representations (all tesla; no flux/weber is stored): - `forcing_b_rootarea`/`response_b_rootarea` - root-area-weighted field b̃ (coordinate-invariant) - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄) -Control surface matrices [numpert_total × numpert_total], stored in the coordinate-invariant -root-area-weighted field (b̃) space (Pharr 2026). Writing `S ≡ rootarea_to_area_weight` -(b̃→b̄) and `A ≡ surface_area`, the brief internal flux-conform operator is `R = S·A` (Φ = R·b̃). -Recover the area-weighted (b̄) forms by conforming with `S` (e.g. `L_b̄ = S·L̃·S†`); recover flux -with `Φ = A·b̄`. - - - `plasma_inductance` - Λ̃ = R⁻¹·Λ·R⁻† (wt0-based plasma inductance, congruence) - - `surface_inductance` - L̃ = R⁻¹·L·R⁻† (vacuum surface inductance, congruence) - - `permeability` - P̃ = R⁻¹·P·R (plasma response operator P=Λ·L⁻¹, similarity) - - `reluctance` - ϱ̃ = R†·ϱ·R (ϱ = L⁻¹·(Λ−L)·L⁻¹, congruence) - - `rootarea_to_area_weight` - S = Σ/√A at psilim (b̃→b̄ recovery operator) - - `surface_area` - scalar A = ∫J|∇ψ|dθ at psilim (flux recovery Φ = A·b̄) - Energies (Joules; Fortran gpout convention). Congruence-invariant scalars (energy = Φ†·G⁻¹·Φ = b̃†·G̃⁻¹·b̃), evaluated from the brief internal flux vectors Φ_x, Φ_tot with the well-conditioned flux-space inductances L, Λ: @@ -150,6 +168,11 @@ well-conditioned flux-space inductances L, Λ: self-consistent solutions, and is a distinct construction from the KineticForces NTV torque. """ @kwdef mutable struct PerturbedEquilibriumState + response::PlasmaResponse + forcing_modes::Vector{ForcingMode} = ForcingMode[] + coil_sets::Vector{CoilSet} = CoilSet[] + regularization_width::Float64 = response.regularization_width + # Radial grid (FFS ODE integration ψ_n values) [npsi] psi_grid::Vector{Float64} = Float64[] @@ -159,13 +182,6 @@ well-conditioned flux-space inductances L, Λ: b_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing # physical normal field b_n [npsi, numpert_total] xi_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing # physical normal displacement xi_n [npsi, numpert_total] - # Coupling matrices [n_rational × numpert_total] - C_resonant_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - C_resonant_current::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - C_island_width_sq::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - C_penetrated_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - C_delta_prime::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - # Applied resonant vectors [n_rational] = C · amp_vec resonant_area_weighted_field::Vector{ComplexF64} = ComplexF64[] resonant_current::Vector{ComplexF64} = ComplexF64[] @@ -173,25 +189,11 @@ well-conditioned flux-space inductances L, Λ: penetrated_area_weighted_field::Vector{ComplexF64} = ComplexF64[] delta_prime::Vector{ComplexF64} = ComplexF64[] forcing_solution_weights::Vector{ComplexF64} = ComplexF64[] - rational_area::Vector{Float64} = Float64[] # Diagnostics [n_rational] island_half_width::Vector{Float64} = Float64[] chirikov_parameter::Vector{Float64} = Float64[] - - # Metadata [n_rational] - rational_psi::Vector{Float64} = Float64[] - rational_q::Vector{Float64} = Float64[] - rational_m_res::Vector{Int} = Int[] - rational_n::Vector{Int} = Int[] - rational_surface_idx::Vector{Int} = Int[] - - # Dominant resonant-coupling modes: SVD of C_resonant_area_weighted_field over the retained rows - dominant_singular_values::Vector{Float64} = Float64[] - dominant_right_singular_vectors::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # V [numpert_total × rank] - dominant_left_singular_vectors::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # U [n_retained × rank] - dominant_rational_index::Vector{Int} = Int[] # retained rows of rational_* - dominant_forcing_overlap::Vector{ComplexF64} = ComplexF64[] # Vᴴ·b̃_x [rank] + dominant_forcing_overlap::Vector{ComplexF64} = ComplexF64[] # Vᴴ·b̃_x [rank] # Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x) @@ -201,39 +203,9 @@ well-conditioned flux-space inductances L, Λ: response_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ - # Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space - plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space) - surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # L̃ (field space) - permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R - reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R - rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # S = Σ/√A at psilim: b̃→b̄ recovery operator - surface_area::Float64 = 0.0 # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A) - # Energies — see the struct docstring for formulas vacuum_energy::Float64 = 0.0 surface_energy::Float64 = 0.0 plasma_energy::Float64 = 0.0 toroidal_torque::Float64 = 0.0 end - -""" - kinetic_regularization_kwargs(ffs, kwargs) -> NamedTuple - -`regularization_width` smooths the **ideal** 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they drive the NTV -integrand. A kinetic solve has no such divergence — det(F̄) is complex and nonzero at the rationals -(Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D) — so regularizing there suppresses a finite -physical response, and does so inconsistently, since ξ^ψ is never regularized. Force -`regularization_width = 0` whenever a kinetic fit is present (either `kinetic_source`, both of which remove the -resonant singularity) and log the override; ideal solves keep their setting. `ffs` is duck-typed: -only `ffs.mats.kinetic` is read, so the override is unit-testable without a full run. -""" -function kinetic_regularization_kwargs(ffs, kwargs) - base = values(kwargs) - default_reg = PerturbedEquilibrium.PerturbedEquilibriumControl().regularization_width - prev = get(base, :regularization_width, default_reg) - kinetic = ffs.mats.kinetic !== nothing - kinetic && prev != 0 && - @info "Kinetic run: overriding regularization_width=$prev with 0 " * - "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" - return merge(base, (; regularization_width=kinetic ? 0.0 : prev)) -end diff --git a/src/PerturbedEquilibrium/ResonantCoupling.jl b/src/PerturbedEquilibrium/ResonantCoupling.jl index 2e129ffc1..9a786a7e0 100644 --- a/src/PerturbedEquilibrium/ResonantCoupling.jl +++ b/src/PerturbedEquilibrium/ResonantCoupling.jl @@ -3,7 +3,7 @@ Everything needed to project an applied control-surface field onto the resonant responses of a perturbed-equilibrium solve, detached from the run that produced it. Build it in memory from a -`PerturbedEquilibriumState` and its `ForceFreeStatesResult`, or post hoc from a `gpec.h5`; then +[`PlasmaResponse`](@ref), or post hoc from a `gpec.h5`; then window and decompose it with [`dominant_coupling`](@ref) and evaluate coil spectra with [`rootarea_field`](@ref) and [`coupling_overlap`](@ref) as often as needed without re-running. @@ -29,26 +29,16 @@ struct ResonantCoupling end """ - ResonantCoupling(state::PerturbedEquilibriumState, ffs::ForceFreeStatesResult) + ResonantCoupling(resp::PlasmaResponse) -In-memory construction from a solve: the coupling matrix and row labels come from `state`, the -column ordering from `ffs`, and the conform operator from `state` when the response step ran or -from the equilibrium otherwise. Errors when `state` carries no singular-coupling matrix. +In-memory construction from a [`PlasmaResponse`](@ref): the coupling matrix, row and column labels +and conform operator all come from `resp`. Errors when `resp` carries no singular-coupling matrix. """ -function ResonantCoupling(state::PerturbedEquilibriumState, ffs::ForceFreeStatesResult) - isempty(state.C_resonant_area_weighted_field) && - throw(ArgumentError("the perturbed-equilibrium state carries no singular-coupling matrix (compute_singular_coupling off, or no resonant surfaces)")) - flux_conform = if isempty(state.rootarea_to_area_weight) - S, A = build_control_surface_rootarea_to_area_weight(ffs.equil, ffs) - S .* A - else - state.rootarea_to_area_weight .* state.surface_area - end - N = ffs.numpert_total - m_modes = [(i - 1) % ffs.mpert + ffs.mlow for i in 1:N] - n_modes = [(i - 1) ÷ ffs.mpert + ffs.nlow for i in 1:N] - return ResonantCoupling(state.C_resonant_area_weighted_field, state.rational_psi, state.rational_q, - state.rational_m_res, state.rational_n, m_modes, n_modes, flux_conform) +function ResonantCoupling(resp::PlasmaResponse) + isempty(resp.C_resonant_area_weighted_field) && + throw(ArgumentError("the plasma response carries no singular-coupling matrix (compute_singular_coupling off, or no resonant surfaces)")) + return ResonantCoupling(resp.C_resonant_area_weighted_field, resp.rational_psi, resp.rational_q, + resp.rational_m_res, resp.rational_n, resp.m_modes, resp.n_modes, resp.rootarea_to_area_weight .* resp.surface_area) end """ @@ -63,7 +53,7 @@ function ResonantCoupling(h5path::AbstractString) haskey(f, "PerturbedEquilibrium/SingularCoupling/C_resonant_area_weighted_field") || throw(ArgumentError("$h5path has no singular-coupling matrix")) haskey(f, "PerturbedEquilibrium/ResponseMatrices/rootarea_to_area_weight_operator") || - throw(ArgumentError("$h5path has no control-surface conform operator (run with compute_response = true)")) + throw(ArgumentError("$h5path has no control-surface conform operator")) haskey(f, "Info/mn_index") || throw(ArgumentError("$h5path has no Info/mn_index mode labels")) sc = f["PerturbedEquilibrium/SingularCoupling"] mn = read(f["Info/mn_index"]) @@ -197,22 +187,20 @@ coupling_overlap(dom::DominantCoupling, b̃::AbstractVector{<:Number}) = dom.rig coupling_overlap(dom::DominantCoupling, rc::ResonantCoupling, modes) = coupling_overlap(dom, rootarea_field(rc, modes)) """ - compute_dominant_coupling!(state, ctrl) + build_dominant_coupling!(resp; verbose=false) -Store the dominant resonant-coupling decomposition of `state.C_resonant_area_weighted_field` -over the core window (`ψ_N ≤ CORE_PSI_HIGH`) on `state` as a run summary, with the applied -forcing's coefficients on each mode when the response spectrum is available. Other windows are +Store the dominant resonant-coupling decomposition of `resp.C_resonant_area_weighted_field` +over the core window (`ψ_N ≤ CORE_PSI_HIGH`) on `resp` as a run summary. Other windows are evaluated post hoc on a [`ResonantCoupling`](@ref). Returns without change when no coupling matrix exists. """ -function compute_dominant_coupling!(state::PerturbedEquilibriumState, ctrl::PerturbedEquilibriumControl) - isempty(state.C_resonant_area_weighted_field) && return nothing - dom = dominant_coupling(state.C_resonant_area_weighted_field, state.rational_psi) - state.dominant_singular_values = dom.singular_values - state.dominant_right_singular_vectors = dom.right_singular_vectors - state.dominant_left_singular_vectors = dom.left_singular_vectors - state.dominant_rational_index = dom.rational_index - isempty(state.forcing_b_rootarea) || (state.dominant_forcing_overlap = coupling_overlap(dom, state.forcing_b_rootarea)) - ctrl.verbose && +function build_dominant_coupling!(resp::PlasmaResponse; verbose::Bool=false) + isempty(resp.C_resonant_area_weighted_field) && return nothing + dom = dominant_coupling(resp.C_resonant_area_weighted_field, resp.rational_psi) + resp.dominant_singular_values = dom.singular_values + resp.dominant_right_singular_vectors = dom.right_singular_vectors + resp.dominant_left_singular_vectors = dom.left_singular_vectors + resp.dominant_rational_index = dom.rational_index + verbose && @info "Dominant resonant coupling: $(length(dom.singular_values)) modes over " * "$(length(dom.rational_index)) rational surfaces, σ₁ = $(@sprintf("%.3e", dom.singular_values[1]))" return nothing diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index 120847e73..1d47daec4 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -1,43 +1,28 @@ """ - compute_plasma_response!( - state, equil, solution, wt0, mtheta_vacuum, ffs, - intr, ctrl, metric, mats - ) - -Compute plasma response to external forcing using ForceFreeStates eigenmode solutions. + build_response_operators!(resp; verbose=false) -Implements resp_index=0 calculation from Fortran gpresp: +Build the control-surface response operators of `resp.ffs` onto `resp`, following Fortran gpresp's +resp_index=0 energy method: - Build flux matrix from eigenmodes - Calculate plasma inductance Lambda (wt0-based, resp_induct_flag=TRUE) - Calculate surface inductance L from Green's function - - Compute permeability P = Lambda * L^{-1} - - Apply forcing to get response: Phi_tot = P * Phi_x + - Compute permeability P = Lambda * L^{-1} and reluctance + - Conform them to the root-area-weighted field space for output """ -function compute_plasma_response!( - state::PerturbedEquilibriumState, - equil::Equilibrium.PlasmaEquilibrium, - solution::SolutionProfiles, - wt0::Matrix{ComplexF64}, - mtheta_vacuum::Int, - ffs::ForceFreeStatesResult, - intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl, - metric::MetricData, - mats::MatrixSplines -) - if ctrl.verbose - @info "Computing plasma response (wt0-based inductance)" - end +function build_response_operators!(resp::PlasmaResponse; verbose::Bool=false) + ffs = resp.ffs + equil = ffs.equil + verbose && @info "Computing plasma response (wt0-based inductance)" # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode] - flux_matrix = build_flux_matrix(equil, solution, ffs) + resp.flux_matrix = build_flux_matrix(equil, ffs.solution, ffs) # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default) - plasma_inductance = calc_plasma_inductance(wt0, ffs, equil.psio) + plasma_inductance = calc_plasma_inductance(ffs.free_boundary.wt0, ffs, equil.psio) # Surface inductance L from vacuum surface-current matrix at psilim, block-diagonal in n. - surface_inductance = calc_surface_inductance(equil, ffs.psilim, mtheta_vacuum, ffs.mlow:ffs.mhigh, ffs.nlow:ffs.nhigh) + surface_inductance = calc_surface_inductance(equil, ffs.psilim, ffs.control.mtheta_vacuum, ffs.mlow:ffs.mhigh, ffs.nlow:ffs.nhigh) permeability = calc_permeability(plasma_inductance, surface_inductance) # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats). @@ -45,10 +30,11 @@ function compute_plasma_response!( L_inv = inv(surface_inductance) reluctance = L_inv * (plasma_inductance' - surface_inductance) * L_inv - # Store permeability in internal state for singular coupling / field reconstruction. - # These consumers operate on the physical control-surface flux Φ_x, so the internal - # copy stays in flux space; only the stored/output quantities are conformed to fields below. - intr.plasma_response = permeability + # Singular coupling, field reconstruction and the energies operate on the physical control-surface + # flux Φ_x, so the flux-space operators are kept; only the stored/output ones are conformed to fields. + resp.plasma_inductance_flux = plasma_inductance + resp.surface_inductance_flux = surface_inductance + resp.permeability_flux = permeability # Conform the control-surface matrices to the coordinate-invariant root-area-weighted # field (b̃) space for output (Pharr 2026). Store the b̃→b̄ operator S = Σ/√A @@ -56,17 +42,37 @@ function compute_plasma_response!( # flux (Φ = A·b̄) — see Utils.jl output docs. rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs) field_mats = field_space_response_matrices(plasma_inductance, surface_inductance, permeability, reluctance, rootarea_to_area_weight, surface_area) - state.plasma_inductance = field_mats.plasma_inductance - state.surface_inductance = field_mats.surface_inductance - state.permeability = field_mats.permeability - state.reluctance = field_mats.reluctance - state.rootarea_to_area_weight = rootarea_to_area_weight - state.surface_area = surface_area + resp.plasma_inductance = field_mats.plasma_inductance + resp.surface_inductance = field_mats.surface_inductance + resp.permeability = field_mats.permeability + resp.reluctance = field_mats.reluctance + resp.rootarea_to_area_weight = rootarea_to_area_weight + resp.surface_area = surface_area + return resp +end + +""" + apply_response!(state; verbose=false) + +Apply the response operators of `state.response` to the forcing on `state`: Phi_tot = P * Phi_x, +the forcing and response spectra in the three field representations, the control-surface energies +and torque, and the reconstructed ξ and b profiles. +""" +function apply_response!(state::PerturbedEquilibriumState; verbose::Bool=false) + resp = state.response + ffs = resp.ffs + equil = ffs.equil + solution = ffs.solution + permeability = resp.permeability_flux + surface_inductance = resp.surface_inductance_flux + plasma_inductance = resp.plasma_inductance_flux + rootarea_to_area_weight = resp.rootarea_to_area_weight + surface_area = resp.surface_area # Forcing and response on the control surface. Flux Φ appears only as a brief internal bridge: # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are # formed via the conform operator R = S·A (Φ = R·b̃). - forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs) + forcing_flux = map_forcing_to_eigenmodes(state.forcing_modes, ffs) response_flux = compute_plasma_response_vector(permeability, forcing_flux) # Output forcing/response in the three Pharr field representations (all tesla): @@ -105,8 +111,8 @@ function compute_plasma_response!( state.plasma_energy = real(py) # Fortran's "total energy" is this pengy xi_modes, b_modes = reconstruct_physical_fields( - response_flux, flux_matrix, solution, equil, ffs, intr, - metric, mats, ctrl + response_flux, resp.flux_matrix, solution, equil, ffs, + ffs.metric, ffs.mats, state.regularization_width ) npsi = size(solution.u_store, 4) @@ -120,7 +126,7 @@ function compute_plasma_response!( state.b_n_modes = b_n_modes state.xi_n_modes = xi_n_modes - if ctrl.verbose - @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" + if verbose + @info "Response complete: $(length(state.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" end end diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 05a97049d..ad4132dcc 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -218,59 +218,31 @@ function _el_solution_at( end """ - compute_singular_coupling_metrics!( - state::PerturbedEquilibriumState, - equil::Equilibrium.PlasmaEquilibrium, - solution::SolutionProfiles, - mtheta_vacuum::Int, - ffs::ForceFreeStatesResult, - intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl, - mats::MatrixSplines - ) - -Compute singular layer coupling matrices and applied resonant vectors. - -Matching Fortran GPEC output: `C_f_x_out` (coupling matrix) and `Phi_res`, `w_isl`, etc. (applied vectors). -[Park Phys. Plasmas 2009 056115] - -## Modifies - -Populates `state` with: - -Coupling matrices `[n_rational × numpert_total]` — C[row, j] = coupling when forcing mode j has unit amplitude: - - - `C_resonant_area_weighted_field`, `C_resonant_current`, `C_island_width_sq`, `C_penetrated_area_weighted_field`, `C_delta_prime` + build_coupling_operators!(resp; verbose=false) -Applied resonant vectors `[n_rational]` = C · amp_vec: - - - `resonant_area_weighted_field`, `resonant_current`, `island_width_sq`, `penetrated_area_weighted_field`, `delta_prime` - -Diagnostics `[n_rational]`: `island_half_width`, `chirikov_parameter` - -Metadata `[n_rational]`: `rational_psi`, `rational_q`, `rational_m_res`, `rational_n`, `rational_surface_idx` +Build the singular-layer coupling matrices of `resp.ffs` onto `resp`, matching Fortran GPEC's +`C_f_x_out` [Park Phys. Plasmas 2009 056115]: one row per resonant (surface, n) pair, giving the +resonant area-weighted field, current, (w/2)², penetrated field and Δ' per unit forcing mode, plus +the rational-surface metadata and the edge-to-solution coefficients `C_coeffs`. The matrices are +stored acting on b̃_x, with flux-space copies in `resp.coupling_flux` for applying a forcing. """ -function compute_singular_coupling_metrics!( - state::PerturbedEquilibriumState, - equil::Equilibrium.PlasmaEquilibrium, - solution::SolutionProfiles, - mtheta_vacuum::Int, - ffs::ForceFreeStatesResult, - intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl, - mats::MatrixSplines -) - ctrl.verbose && @info "Computing singular coupling metrics (GPEC method)" +function build_coupling_operators!(resp::PlasmaResponse; verbose::Bool=false) + ffs = resp.ffs + equil = ffs.equil + solution = ffs.solution + mats = ffs.mats + mtheta_vacuum = ffs.control.mtheta_vacuum + verbose && @info "Computing singular coupling metrics (GPEC method)" (; numpert_total, mlow, mhigh, nlow, nhigh) = ffs msing = length(ffs.surfaces) if msing == 0 - ctrl.verbose && @info "No singular surfaces found. Skipping singular coupling calculation." + verbose && @info "No singular surfaces found. Skipping singular coupling calculation." return end - if size(intr.plasma_response, 1) == 0 + if isempty(resp.permeability_flux) @warn "Permeability matrix not computed. Skipping singular coupling." return end @@ -295,30 +267,30 @@ function compute_singular_coupling_metrics!( n_rational = length(resonant_pairs) if n_rational == 0 - ctrl.verbose && @info "No resonant surfaces found. Skipping singular coupling." + verbose && @info "No resonant surfaces found. Skipping singular coupling." return end - ctrl.verbose && @info "Found $n_rational resonant (surface, n) pairs" + verbose && @info "Found $n_rational resonant (surface, n) pairs" # Phase 2: Allocate output arrays - state.C_resonant_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) - state.C_resonant_current = zeros(ComplexF64, n_rational, numpert_total) - state.C_island_width_sq = zeros(ComplexF64, n_rational, numpert_total) - have_inner_bpen = !isempty(intr.inner_bpen) + resp.C_resonant_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) + resp.C_resonant_current = zeros(ComplexF64, n_rational, numpert_total) + resp.C_island_width_sq = zeros(ComplexF64, n_rational, numpert_total) + have_inner_bpen = !isempty(ffs.bpen) if have_inner_bpen - state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) + resp.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) else - state.C_penetrated_area_weighted_field = zeros(ComplexF64, 0, 0) + resp.C_penetrated_area_weighted_field = zeros(ComplexF64, 0, 0) @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog=1 end - state.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total) - state.rational_psi = zeros(Float64, n_rational) - state.rational_q = zeros(Float64, n_rational) - state.rational_m_res = zeros(Int, n_rational) - state.rational_n = zeros(Int, n_rational) - state.rational_surface_idx = zeros(Int, n_rational) - state.rational_area = zeros(Float64, n_rational) + resp.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total) + resp.rational_psi = zeros(Float64, n_rational) + resp.rational_q = zeros(Float64, n_rational) + resp.rational_m_res = zeros(Int, n_rational) + resp.rational_n = zeros(Int, n_rational) + resp.rational_surface_idx = zeros(Int, n_rational) + resp.rational_area = zeros(Float64, n_rational) # Precompute ODE coefficient matrix C_coeffs for all PE forcing modes. # For each forcing mode k: c_k = u_bnd⁻¹ × edge_mn_k @@ -326,15 +298,16 @@ function compute_singular_coupling_metrics!( # Matches Fortran gpout_resp: edge_mn = foutmn/(chi1·singfac·twopi·ifac) psi_lim = solution.psi_store[solution.step] q_lim = equil.profiles.q_spline(psi_lim) - singfac_lim = [intr.m_modes[j] - intr.n_modes[j] * q_lim for j in 1:numpert_total] + singfac_lim = [resp.m_modes[j] - resp.n_modes[j] * q_lim for j in 1:numpert_total] u_bnd = solution.u_store[:, :, 1, solution.step] # Divide each row j by singfac_lim[j] — reshape to column vector so Julia broadcasts row-wise, not column-wise. - edge_mn = intr.plasma_response ./ (chi1 * 2π * im .* reshape(singfac_lim, :, 1)) + edge_mn = resp.permeability_flux ./ (chi1 * 2π * im .* reshape(singfac_lim, :, 1)) C_coeffs = u_bnd \ edge_mn # mpert × numpert_total + resp.C_coeffs = C_coeffs # Phase 3: Compute full coupling matrix rows. Each rational surface is independent -- it # builds its own Green's functions (a fresh, allocation-local vacuum solve) and writes its - # own disjoint coupling-matrix row (state.C_*[row,:] / state.rational_*[row]) -- so thread + # own disjoint coupling-matrix row (resp.C_*[row,:] / resp.rational_*[row]) -- so thread # the loop over surfaces. Note a single surface `s` can appear on multiple rows in a # multi-n run (integer q resonates at several n), so per-row work must not write shared # sing[s] state. BLAS is pinned to one thread across the loop (each surface's solve is @@ -353,7 +326,7 @@ function compute_singular_coupling_metrics!( sing_surf = ffs.surfaces[s] m_res = round(Int, sing_surf.q * nn) - resnum = findfirst(j -> intr.m_modes[j] == m_res && intr.n_modes[j] == nn, 1:numpert_total) + resnum = findfirst(j -> resp.m_modes[j] == m_res && resp.n_modes[j] == nn, 1:numpert_total) if resnum === nothing @warn "Could not find index for resonant mode (m=$m_res, n=$nn)" maxlog = 1 continue @@ -383,7 +356,7 @@ function compute_singular_coupling_metrics!( # Evaluate u and dξ/dψ at lpsi and rpsi from the stored ODE solution. # Branch order mirrors Fortran gpeq_sol: gal-matched solutions take the analytic # galerkin Ξ′, ideal runs the EL relation, kinetic runs the stored Ξ′. - if intr.odet_from_gal + if resp.odet_from_gal # interpolate u and the analytic galerkin dξ/dψ carried in du_store u_l, ud_l = _gal_solution_at(lpsi, resnum, solution, nstep) u_r, ud_r = _gal_solution_at(rpsi, resnum, solution, nstep) @@ -420,9 +393,9 @@ function compute_singular_coupling_metrics!( # coil-drive columns as the outer solution, so it contracts with C_coeffs exactly like # the outer solution values above (xsp = transpose(u) * ck); /area matches the area-weighted # convention of the pointwise row. - if have_inner_bpen && s <= size(intr.inner_bpen, 1) - pen_row = (transpose(C_coeffs) * @view(intr.inner_bpen[s, :])) ./ area - state.C_penetrated_area_weighted_field[row, :] = pen_row + if have_inner_bpen && s <= size(ffs.bpen, 1) + pen_row = (transpose(C_coeffs) * @view(ffs.bpen[s, :])) ./ area + resp.C_penetrated_area_weighted_field[row, :] = pen_row end # LHS normalization audit — output scalar coordinate-invariance per row: @@ -431,25 +404,25 @@ function compute_singular_coupling_metrics!( # - resonant (shielding) current: j_c already integrates jac·|∇ψ| over the surface, so the # Jacobian weighting is carried inside j_c — no separate area factor needed. # - resonant flux → field: Φ^r/A^r [T], invariant [Park 2008; Pharr 2026]. - state.C_delta_prime[row, :] = jump_vec ./ (twopi * chi1) - state.C_resonant_current[row, :] = jump_vec .* (-j_c / (twopi * m_res)) + resp.C_delta_prime[row, :] = jump_vec ./ (twopi * chi1) + resp.C_resonant_current[row, :] = jump_vec .* (-j_c / (twopi * m_res)) # Matches Fortran gpout_resp: singflx = L·fkaxmn, resonant area-weighted field = singflx/area, # islandhwids = 4·singflx/(2π·shear·q·chi1) - singflx_pre = (L_mm / (twopi * nn)) .* state.C_resonant_current[row, :] - state.C_resonant_area_weighted_field[row, :] = singflx_pre ./ area + singflx_pre = (L_mm / (twopi * nn)) .* resp.C_resonant_current[row, :] + resp.C_resonant_area_weighted_field[row, :] = singflx_pre ./ area if abs(shear) > 1e-10 - state.C_island_width_sq[row, :] = (4.0 / (twopi * shear * sing_surf.q * chi1)) .* singflx_pre + resp.C_island_width_sq[row, :] = (4.0 / (twopi * shear * sing_surf.q * chi1)) .* singflx_pre end - state.rational_psi[row] = sing_surf.psifac - state.rational_area[row] = area - state.rational_q[row] = sing_surf.q - state.rational_m_res[row] = m_res - state.rational_n[row] = nn - state.rational_surface_idx[row] = s + resp.rational_psi[row] = sing_surf.psifac + resp.rational_area[row] = area + resp.rational_q[row] = sing_surf.q + resp.rational_m_res[row] = m_res + resp.rational_n[row] = nn + resp.rational_surface_idx[row] = s - if ctrl.verbose - dp_diag = real(state.C_delta_prime[row, resnum]) + if verbose + dp_diag = real(resp.C_delta_prime[row, resnum]) @info "Row $row: q=$(@sprintf("%.3f", sing_surf.q)), ψ=$(@sprintf("%.3f", sing_surf.psifac)), m=$m_res, n=$nn, Δ'(diag)=$(@sprintf("%.3e", dp_diag))" end end @@ -457,40 +430,54 @@ function compute_singular_coupling_metrics!( BLAS.set_num_threads(_blas_nthreads) end - # Phase 4: Apply forcing amplitudes → R = C · Φ_x. The applied resonant scalars are - # physical, coordinate-invariant quantities, so evaluate them from the flux-space C and - # the physical forcing Φ_x (bit-identical to pre-conform values). + # Conform the stored coupling-matrix input basis to the coordinate-invariant root-area-weighted + # field (b̃) space (Pharr 2026): C̃ = C·R, so each stored row acts on the applied field + # b̃_x (Φ_x = R·b̃_x) and its singular values become coordinate-invariant. The conform operator + # R = S·A (Σ·√A) is the only place flux briefly appears; it is built from the b̃→b̄ operator S and + # the scalar surface area A. The flux-space copies kept in `coupling_flux` apply a forcing without round-trip noise. + resp.coupling_flux = (; resonant_area_weighted_field=resp.C_resonant_area_weighted_field, resonant_current=resp.C_resonant_current, + island_width_sq=resp.C_island_width_sq, penetrated_area_weighted_field=resp.C_penetrated_area_weighted_field, delta_prime=resp.C_delta_prime) + flux_conform = resp.rootarea_to_area_weight .* resp.surface_area + resp.C_resonant_area_weighted_field = resp.C_resonant_area_weighted_field * flux_conform + resp.C_resonant_current = resp.C_resonant_current * flux_conform + resp.C_island_width_sq = resp.C_island_width_sq * flux_conform + resp.C_delta_prime = resp.C_delta_prime * flux_conform + have_inner_bpen && (resp.C_penetrated_area_weighted_field = resp.C_penetrated_area_weighted_field * flux_conform) + + return resp +end + +""" + apply_coupling!(state; verbose=false) + +Apply the coupling matrices of `state.response` to the forcing on `state`: the applied resonant vectors +(Fortran `Phi_res`, `w_isl`, `K_isl`, `Delta`), the forcing's solution weights, and the island +half-widths and Chirikov overlap. +""" +function apply_coupling!(state::PerturbedEquilibriumState; verbose::Bool=false) + resp = state.response + C = resp.coupling_flux + numpert_total = resp.ffs.numpert_total + + # Apply forcing amplitudes → R = C · Φ_x. The applied resonant scalars are physical, + # coordinate-invariant quantities, so evaluate them from the flux-space C and the physical forcing Φ_x. forcing_flux = zeros(ComplexF64, numpert_total) - for mode in intr.forcing_modes - j = findfirst(k -> intr.m_modes[k] == mode.m && intr.n_modes[k] == mode.n, 1:numpert_total) + for mode in state.forcing_modes + j = findfirst(k -> resp.m_modes[k] == mode.m && resp.n_modes[k] == mode.n, 1:numpert_total) isnothing(j) || (forcing_flux[j] = mode.amplitude) end - state.resonant_area_weighted_field = state.C_resonant_area_weighted_field * forcing_flux - state.resonant_current = state.C_resonant_current * forcing_flux - state.island_width_sq = state.C_island_width_sq * forcing_flux - state.delta_prime = state.C_delta_prime * forcing_flux - have_inner_bpen && (state.penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * forcing_flux) - state.forcing_solution_weights = C_coeffs * forcing_flux + state.resonant_area_weighted_field = C.resonant_area_weighted_field * forcing_flux + state.resonant_current = C.resonant_current * forcing_flux + state.island_width_sq = C.island_width_sq * forcing_flux + state.delta_prime = C.delta_prime * forcing_flux + isempty(C.penetrated_area_weighted_field) || (state.penetrated_area_weighted_field = C.penetrated_area_weighted_field * forcing_flux) + state.forcing_solution_weights = resp.C_coeffs * forcing_flux - # Conform the stored coupling-matrix input basis to the coordinate-invariant root-area-weighted - # field (b̃) space (Pharr 2026): C̃ = C·R, so each stored row acts on the applied field - # b̃_x (Φ_x = R·b̃_x) and its singular values become coordinate-invariant. The conform operator - # R = S·A (Σ·√A) is the only place flux briefly appears; it is built from the b̃→b̄ operator S and - # the scalar surface area A. Done after the applied-vector evaluation above so those physical - # scalars carry no round-trip noise. - rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs) - flux_conform = rootarea_to_area_weight .* surface_area - state.C_resonant_area_weighted_field = state.C_resonant_area_weighted_field * flux_conform - state.C_resonant_current = state.C_resonant_current * flux_conform - state.C_island_width_sq = state.C_island_width_sq * flux_conform - state.C_delta_prime = state.C_delta_prime * flux_conform - have_inner_bpen && (state.C_penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * flux_conform) - - # Phase 5: Island diagnostics from applied resonant vectors - compute_island_diagnostics!(state, n_rational) - - if ctrl.verbose + # Island diagnostics from the applied resonant vectors + state.island_half_width, state.chirikov_parameter = island_diagnostics(state.island_width_sq, resp.rational_psi, resp.rational_surface_idx) + + if verbose max_island = isempty(state.island_half_width) ? 0.0 : maximum(state.island_half_width) max_dp = isempty(state.delta_prime) ? 0.0 : maximum(abs.(real.(state.delta_prime))) @info "Singular coupling complete. Max island half-width = $(@sprintf("%.3e", max_island)), max |Δ'| = $(@sprintf("%.3e", max_dp))" @@ -645,32 +632,29 @@ function compute_surface_area( end """ - compute_island_diagnostics!(state::PerturbedEquilibriumState, n_rational::Int) + island_diagnostics(island_width_sq, rational_psi, rational_surface_idx) -> (island_half_width, chirikov_parameter) -Compute island half-width and Chirikov parameter from applied resonant vectors. +Island half-width and Chirikov parameter from the applied (w/2)² on each rational row. - `island_half_width[row]` = √|island_width_sq[row]| - `chirikov_parameter[row]` = half-width / (half-distance to the nearest other rational surface); rows on the same surface at another n (integer q) are not neighbors """ -function compute_island_diagnostics!(state::PerturbedEquilibriumState, n_rational::Int) - state.island_half_width = sqrt.(abs.(state.island_width_sq)) - state.chirikov_parameter = zeros(Float64, n_rational) - - n_rational <= 1 && return +function island_diagnostics(island_width_sq::AbstractVector, rational_psi::AbstractVector{Float64}, rational_surface_idx::AbstractVector{Int}) + island_half_width = sqrt.(abs.(island_width_sq)) + n_rational = length(rational_psi) + chirikov_parameter = zeros(Float64, n_rational) + n_rational <= 1 && return island_half_width, chirikov_parameter for row in 1:n_rational - psi_row = state.rational_psi[row] + psi_row = rational_psi[row] min_dist = Inf for row2 in 1:n_rational - state.rational_surface_idx[row2] == state.rational_surface_idx[row] && continue - dist = abs(state.rational_psi[row2] - psi_row) + rational_surface_idx[row2] == rational_surface_idx[row] && continue + dist = abs(rational_psi[row2] - psi_row) min_dist = min(min_dist, dist) end - if min_dist > 1e-10 - state.chirikov_parameter[row] = state.island_half_width[row] / (min_dist / 2.0) - else - state.chirikov_parameter[row] = Inf - end + chirikov_parameter[row] = min_dist > 1e-10 ? island_half_width[row] / (min_dist / 2.0) : Inf end + return island_half_width, chirikov_parameter end diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index d2064aee4..3f821e950 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -1,53 +1,19 @@ """ - initialize_mode_arrays!( - intr::PerturbedEquilibriumInternal, - ffs::ForceFreeStatesResult - ) - -Initialize mode number arrays for convenient indexing. - -Pre-computes m_modes[i] and n_modes[i] for each linear index i in 1:numpert_total, -avoiding repeated index arithmetic throughout the code. + mode_arrays(ffs::ForceFreeStatesResult) -> (m_modes, n_modes) -## Mode Indexing Convention - -For linear index i ∈ [1, numpert_total]: +Mode numbers of each linear index i ∈ [1, numpert_total], in the ForceFreeStates ordering +(m1,n1), (m2,n1), ..., (mpert,n1), (m1,n2), ...: - m_modes[i] = (i-1) % mpert + mlow - n_modes[i] = (i-1) ÷ mpert + nlow - -This matches the convention used in ForceFreeStates where modes are ordered as: -(m1,n1), (m2,n1), ..., (mpert,n1), (m1,n2), (m2,n2), ..., (mpert,npert) """ -function initialize_mode_arrays!( - intr::PerturbedEquilibriumInternal, - ffs::ForceFreeStatesResult -) - numpert_total = ffs.numpert_total - mpert = ffs.mpert - mlow = ffs.mlow - nlow = ffs.nlow - - # Allocate arrays - intr.m_modes = zeros(Int, numpert_total) - intr.n_modes = zeros(Int, numpert_total) - - # Fill arrays using indexing convention - for i in 1:numpert_total - m_idx = (i - 1) % mpert + 1 # 1-based index into mpert range - n_idx = (i - 1) ÷ mpert + 1 # 1-based index into npert range - - intr.m_modes[i] = m_idx - 1 + mlow - intr.n_modes[i] = n_idx - 1 + nlow - end +function mode_arrays(ffs::ForceFreeStatesResult) + idx = 1:ffs.numpert_total + return [(i - 1) % ffs.mpert + ffs.mlow for i in idx], [(i - 1) ÷ ffs.mpert + ffs.nlow for i in idx] end """ - write_outputs_to_HDF5( - state::PerturbedEquilibriumState, - intr::PerturbedEquilibriumInternal, - filename::String - ) + write_outputs_to_HDF5(state::PerturbedEquilibriumState, filename::String) Write perturbed equilibrium results to HDF5 file (appends to existing ForceFreeStates output). @@ -108,23 +74,20 @@ PerturbedEquilibrium/ └── toroidal_torque ``` """ -function write_outputs_to_HDF5( - state::PerturbedEquilibriumState, - intr::PerturbedEquilibriumInternal, - filename::String -) +function write_outputs_to_HDF5(state::PerturbedEquilibriumState, filename::String) + resp = state.response h5open(filename, "cw") do file pe_group = haskey(file, "PerturbedEquilibrium") ? file["PerturbedEquilibrium"] : create_group(file, "PerturbedEquilibrium") # Forcing modes forcing_group = haskey(pe_group, "ForcingModes") ? pe_group["ForcingModes"] : create_group(pe_group, "ForcingModes") - forcing_group["n"] = [mode.n for mode in intr.forcing_modes] - forcing_group["m"] = [mode.m for mode in intr.forcing_modes] - forcing_group["amplitude"] = [mode.amplitude for mode in intr.forcing_modes] + forcing_group["n"] = [mode.n for mode in state.forcing_modes] + forcing_group["m"] = [mode.m for mode in state.forcing_modes] + forcing_group["amplitude"] = [mode.amplitude for mode in state.forcing_modes] # Labels of the mode axis shared by the spectra, matrices and profiles below - !isempty(intr.m_modes) && (pe_group["mode_m"] = intr.m_modes) - !isempty(intr.n_modes) && (pe_group["mode_n"] = intr.n_modes) + !isempty(resp.m_modes) && (pe_group["mode_m"] = resp.m_modes) + !isempty(resp.n_modes) && (pe_group["mode_n"] = resp.n_modes) # Control-surface forcing/response spectra in the three Pharr field representations # (all tesla; flux/weber is never stored). b̃ = root-area-weighted (coordinate-invariant). @@ -140,12 +103,12 @@ function write_outputs_to_HDF5( # operator S ≡ rootarea_to_area_weight (b̄ = S·b̃): e.g. L_b̄ = S·L̃·S†; recover flux with the # scalar surface_area A: Φ = A·b̄ (internally R = S·A, Φ = R·b̃). [Pharr 2026] mat_group = haskey(pe_group, "ResponseMatrices") ? pe_group["ResponseMatrices"] : create_group(pe_group, "ResponseMatrices") - !isempty(state.plasma_inductance) && (mat_group["plasma_inductance"] = state.plasma_inductance) - !isempty(state.surface_inductance) && (mat_group["surface_inductance"] = state.surface_inductance) - !isempty(state.permeability) && (mat_group["permeability"] = state.permeability) - !isempty(state.reluctance) && (mat_group["reluctance"] = state.reluctance) - !isempty(state.rootarea_to_area_weight) && (mat_group["rootarea_to_area_weight_operator"] = state.rootarea_to_area_weight) - (state.surface_area != 0.0) && (mat_group["surface_area"] = state.surface_area) + !isempty(resp.plasma_inductance) && (mat_group["plasma_inductance"] = resp.plasma_inductance) + !isempty(resp.surface_inductance) && (mat_group["surface_inductance"] = resp.surface_inductance) + !isempty(resp.permeability) && (mat_group["permeability"] = resp.permeability) + !isempty(resp.reluctance) && (mat_group["reluctance"] = resp.reluctance) + !isempty(resp.rootarea_to_area_weight) && (mat_group["rootarea_to_area_weight_operator"] = resp.rootarea_to_area_weight) + (resp.surface_area != 0.0) && (mat_group["surface_area"] = resp.surface_area) # Response fields (ComplexF64 directly) response_group = haskey(pe_group, "Response") ? pe_group["Response"] : create_group(pe_group, "Response") @@ -207,11 +170,11 @@ function write_outputs_to_HDF5( coupling_group = haskey(pe_group, "SingularCoupling") ? pe_group["SingularCoupling"] : create_group(pe_group, "SingularCoupling") # Coupling matrices [n_rational × numpert_total] - !isempty(state.C_resonant_area_weighted_field) && (coupling_group["C_resonant_area_weighted_field"] = state.C_resonant_area_weighted_field) - !isempty(state.C_resonant_current) && (coupling_group["C_resonant_current"] = state.C_resonant_current) - !isempty(state.C_island_width_sq) && (coupling_group["C_island_width_sq"] = state.C_island_width_sq) - !isempty(state.C_penetrated_area_weighted_field) && (coupling_group["C_penetrated_area_weighted_field"] = state.C_penetrated_area_weighted_field) - !isempty(state.C_delta_prime) && (coupling_group["C_Delta_prime"] = state.C_delta_prime) + !isempty(resp.C_resonant_area_weighted_field) && (coupling_group["C_resonant_area_weighted_field"] = resp.C_resonant_area_weighted_field) + !isempty(resp.C_resonant_current) && (coupling_group["C_resonant_current"] = resp.C_resonant_current) + !isempty(resp.C_island_width_sq) && (coupling_group["C_island_width_sq"] = resp.C_island_width_sq) + !isempty(resp.C_penetrated_area_weighted_field) && (coupling_group["C_penetrated_area_weighted_field"] = resp.C_penetrated_area_weighted_field) + !isempty(resp.C_delta_prime) && (coupling_group["C_Delta_prime"] = resp.C_delta_prime) # Applied resonant vectors [n_rational] !isempty(state.resonant_area_weighted_field) && (coupling_group["resonant_area_weighted_field"] = state.resonant_area_weighted_field) @@ -220,23 +183,23 @@ function write_outputs_to_HDF5( !isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field) !isempty(state.delta_prime) && (coupling_group["Delta_prime"] = state.delta_prime) !isempty(state.forcing_solution_weights) && (coupling_group["forcing_solution_weights"] = state.forcing_solution_weights) - !isempty(state.rational_area) && (coupling_group["rational_area"] = state.rational_area) + !isempty(resp.rational_area) && (coupling_group["rational_area"] = resp.rational_area) !isempty(state.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width) !isempty(state.chirikov_parameter) && (coupling_group["chirikov_parameter"] = state.chirikov_parameter) # Metadata [n_rational] - !isempty(state.rational_psi) && (coupling_group["rational_psi"] = state.rational_psi) - !isempty(state.rational_q) && (coupling_group["rational_q"] = state.rational_q) - !isempty(state.rational_m_res) && (coupling_group["rational_m"] = state.rational_m_res) - !isempty(state.rational_n) && (coupling_group["rational_n"] = state.rational_n) + !isempty(resp.rational_psi) && (coupling_group["rational_psi"] = resp.rational_psi) + !isempty(resp.rational_q) && (coupling_group["rational_q"] = resp.rational_q) + !isempty(resp.rational_m_res) && (coupling_group["rational_m"] = resp.rational_m_res) + !isempty(resp.rational_n) && (coupling_group["rational_n"] = resp.rational_n) # Dominant resonant-coupling modes (SVD of the b̃-space coupling matrix over the retained rows) - if !isempty(state.dominant_singular_values) + if !isempty(resp.dominant_singular_values) dominant_group = haskey(coupling_group, "DominantMode") ? coupling_group["DominantMode"] : create_group(coupling_group, "DominantMode") - dominant_group["singular_values"] = state.dominant_singular_values - dominant_group["right_singular_vectors"] = state.dominant_right_singular_vectors - dominant_group["left_singular_vectors"] = state.dominant_left_singular_vectors - dominant_group["rational_index"] = state.dominant_rational_index + dominant_group["singular_values"] = resp.dominant_singular_values + dominant_group["right_singular_vectors"] = resp.dominant_right_singular_vectors + dominant_group["left_singular_vectors"] = resp.dominant_left_singular_vectors + dominant_group["rational_index"] = resp.dominant_rational_index !isempty(state.dominant_forcing_overlap) && (dominant_group["forcing_overlap"] = state.dominant_forcing_overlap) end diff --git a/test/runtests_dominant_coupling.jl b/test/runtests_dominant_coupling.jl index fe7cab712..21ed66108 100644 --- a/test/runtests_dominant_coupling.jl +++ b/test/runtests_dominant_coupling.jl @@ -74,31 +74,32 @@ include("h5_metadata_check.jl") inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = true inputs["ForcingTerms"] = Dict{String,Any}("forcing_data_file" => "forcing.dat", "forcing_data_format" => "ascii") inputs["PerturbedEquilibrium"] = Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, + "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => true) open(io -> TOML.print(io, inputs), toml_path, "w") run = GPEC.main([dir]) pe, ffs = run.pe, run.ffs - # The in-memory entry point needs a coupling matrix on the state. - @test_throws ArgumentError PE.ResonantCoupling(PE.PerturbedEquilibriumState(), ffs) + # The in-memory entry point needs a coupling matrix on the response. + @test_throws ArgumentError PE.ResonantCoupling(PE.PlasmaResponse(; ffs)) + resp = pe.response h5path = joinpath(dir, "gpec.h5") # Run summary: the core-window decomposition, rebuilt into the stored applied resonant # field (which the writer evaluates before conforming the matrix to b̃ space). - @test !isempty(pe.dominant_singular_values) - rank = length(pe.dominant_singular_values) - rows = pe.dominant_rational_index - @test rows == 1:length(pe.rational_psi) - @test issorted(pe.dominant_singular_values; rev=true) - V, σ, U = pe.dominant_right_singular_vectors, pe.dominant_singular_values, pe.dominant_left_singular_vectors + @test !isempty(resp.dominant_singular_values) + rank = length(resp.dominant_singular_values) + rows = resp.dominant_rational_index + @test rows == 1:length(resp.rational_psi) + @test issorted(resp.dominant_singular_values; rev=true) + V, σ, U = resp.dominant_right_singular_vectors, resp.dominant_singular_values, resp.dominant_left_singular_vectors @test size(V) == (length(pe.forcing_b_rootarea), rank) @test size(U) == (length(rows), rank) @test pe.dominant_forcing_overlap[1] ≈ dot(V[:, 1], pe.forcing_b_rootarea) @test U * (σ .* pe.dominant_forcing_overlap) ≈ pe.resonant_area_weighted_field[rows] rtol = 1e-8 # The same ResonantCoupling from memory and from the file. - rc = PE.ResonantCoupling(pe, ffs) + rc = PE.ResonantCoupling(resp) rc_h5 = PE.ResonantCoupling(h5path) @test rc_h5.C ≈ rc.C @test rc_h5.rational_psi == rc.rational_psi @@ -138,10 +139,6 @@ include("h5_metadata_check.jl") @test isempty(bare.m_modes) @test PE.check_mode_basis(bare, rc.m_modes .+ 1, rc.n_modes, "unlabelled") === nothing - # Without the response step the conform operator is rebuilt from the equilibrium. - pe.rootarea_to_area_weight = zeros(ComplexF64, 0, 0) - @test PE.ResonantCoupling(pe, ffs).flux_conform ≈ rc.flux_conform - h5open(h5path, "r") do h5 g = h5["PerturbedEquilibrium/SingularCoupling/DominantMode"] @test read(g["singular_values"]) ≈ σ @@ -157,12 +154,8 @@ end @testset "island diagnostics across toroidal modes" begin PE = GeneralizedPerturbedEquilibrium.PerturbedEquilibrium # q = 2 surface (index 1) resonates at n = 1 and n = 2; q = 2.5 (index 2) at n = 2 only. - state = PE.PerturbedEquilibriumState() - state.rational_psi = [0.4, 0.4, 0.6] - state.rational_surface_idx = [1, 1, 2] - state.island_width_sq = ComplexF64[0.01, 0.0004, 0.0025] - PE.compute_island_diagnostics!(state, 3) - @test state.island_half_width ≈ [0.1, 0.02, 0.05] + half_width, chirikov = PE.island_diagnostics(ComplexF64[0.01, 0.0004, 0.0025], [0.4, 0.4, 0.6], [1, 1, 2]) + @test half_width ≈ [0.1, 0.02, 0.05] # The repeated surface is not its own neighbour: every row sees the 0.2 gap to the other surface. - @test state.chirikov_parameter ≈ [0.1, 0.02, 0.05] ./ 0.1 + @test chirikov ≈ [0.1, 0.02, 0.05] ./ 0.1 end diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 7f776e1f1..a8a52c78e 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -172,12 +172,13 @@ using Statistics @test_throws r"`moment`, `save_records`" GPEC._stop_on_deprecated_keys( Dict{String,Any}("moment" => "pressure", "save_records" => false, "fgar_flag" => true), GPEC._DEPRECATED_KF_KEYS, "KineticForces") @test_throws r"\[PerturbedEquilibrium\] sets `fixed_boundary`" GPEC._stop_on_deprecated_keys( - Dict{String,Any}("fixed_boundary" => false, "compute_response" => false), GPEC._DEPRECATED_PE_KEYS, "PerturbedEquilibrium") + Dict{String,Any}("fixed_boundary" => false, "compute_singular_coupling" => false), GPEC._DEPRECATED_PE_KEYS, "PerturbedEquilibrium") @test GPEC._stop_on_deprecated_keys(Dict{String,Any}("nn_low" => 1), GPEC._DEPRECATED_FFS_KEYS, "ForceFreeStates") === nothing # A listed deck key must not still be a control field, or listing it would stop runs that set a live key. @test !any(k -> hasfield(GPEC.ForceFreeStates.ForceFreeStatesControl, Symbol(k)), GPEC._DEPRECATED_FFS_KEYS) @test !any(k -> hasfield(GPEC.KineticForces.KineticForcesControl, Symbol(k)), GPEC._DEPRECATED_KF_KEYS) - @test !any(k -> hasfield(GPEC.PerturbedEquilibrium.PerturbedEquilibriumControl, Symbol(k)), GPEC._DEPRECATED_PE_KEYS) + pe_keywords = union((Base.kwarg_decl(m) for f in (GPEC.plasma_response, GPEC.perturbed_equilibrium) for m in methods(f))...) + @test !any(k -> Symbol(k) in pe_keywords, GPEC._DEPRECATED_PE_KEYS) end @testset "Inverse grid types (original, rho)" begin diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index ccd9e098e..b2a985bed 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -75,7 +75,7 @@ include("h5_metadata_check.jl") Dict{String,Any}("name" => "hoop_efc", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.0, "currents" => [2.0e3], "tilty" => [2.0]) ]) inputs["PerturbedEquilibrium"] = Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, + "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => true) cp(joinpath(@__DIR__, "test_data", "ErrorFields", "tolerances_two_hoops.toml"), joinpath(dir, "tolerances.toml")) # Kinetic profiles consistent with this Solovev pressure, for the NTV torque couplings. @@ -129,7 +129,7 @@ include("h5_metadata_check.jl") @test Ty ≈ im .* Tx rtol = 1e-8 # Projection onto the full-window dominant mode: the sets' overlaps sum to the run's. - rc = PE.ResonantCoupling(pe, ffs) + rc = PE.ResonantCoupling(pe.response) dom = PE.dominant_coupling(rc) table = EF.sensitivity_table(sens, dom) @test table.mode == 1 @@ -531,7 +531,7 @@ include("h5_metadata_check.jl") @test length(with_profiles.subplots) == 3 # and the deposition profile its own panel @test with_profiles[3].series_list[1][:y] == c.torque_full_profile[:, 2] @test with_profiles[3].series_list[2][:y] == c.torque_residual_profile[:, 2] - @test_throws ErrorException GPEC.efc_couplings(ffs, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl(), nothing) + @test_throws ErrorException GPEC.efc_couplings(pe.response, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl(), nothing) # Central differences: doubling the step moves the derivatives at O(h²). coarse = EF.compute_coil_sensitivities(sets, rc, ffs.equil, cfg, @@ -632,7 +632,7 @@ include("h5_metadata_check.jl") "height" => 0.4, "currents" => [2.0e3], "tiltx" => [3.0]) ]) inputs["PerturbedEquilibrium"] = Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, + "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => false) inputs["ErrorFields"] = Dict{String,Any}( "verbose" => false, "output_filename" => "ef_only.h5") diff --git a/test/runtests_kinetic_validity.jl b/test/runtests_kinetic_validity.jl index cca53fa45..86044955e 100644 --- a/test/runtests_kinetic_validity.jl +++ b/test/runtests_kinetic_validity.jl @@ -91,23 +91,9 @@ @test psi_c_spiked < 0.5 end - @testset "kinetic_regularization_kwargs" begin - default_reg = PE.PerturbedEquilibriumControl().regularization_width - ideal = (; mats=(; kinetic=nothing)) - kinetic = (; mats=(; kinetic=(;))) - - # Ideal runs keep whatever the deck asked for, explicit or defaulted. - @test PE.kinetic_regularization_kwargs(ideal, pairs((; regularization_width=0.05))).regularization_width == 0.05 - @test PE.kinetic_regularization_kwargs(ideal, pairs((;))).regularization_width == default_reg - - # Kinetic runs force 0 — including when the deck omits the key entirely, which must pick up - # the struct default rather than leaving regularization_width unset. - @test PE.kinetic_regularization_kwargs(kinetic, pairs((; regularization_width=0.05))).regularization_width == 0.0 - @test PE.kinetic_regularization_kwargs(kinetic, pairs((;))).regularization_width == 0.0 - @test default_reg != 0.0 # otherwise the override above proves nothing - - # Unrelated kwargs pass through untouched. - out = PE.kinetic_regularization_kwargs(kinetic, pairs((; regularization_width=0.05, verbose=false))) - @test out.verbose == false + @testset "default_regularization_width" begin + # A kinetic fit removes the ideal resonant singularity, so its default is no regularization. + @test PE.default_regularization_width((; mats=(; kinetic=nothing))) == 0.05 + @test PE.default_regularization_width((; mats=(; kinetic=(;)))) == 0.0 end end diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 6a67801ca..e3f620868 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -68,7 +68,7 @@ end inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = true inputs["ForcingTerms"] = Dict{String,Any}("forcing_data_file" => "forcing.dat", "forcing_data_format" => "ascii") inputs["PerturbedEquilibrium"] = Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, + "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => true) open(io -> TOML.print(io, inputs), toml_path, "w") diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index 7c03a0964..8c0c3f5fc 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -118,7 +118,6 @@ include("h5_metadata_check.jl") "forcing_data_file" => "forcing.dat", "forcing_data_format" => "ascii"), "PerturbedEquilibrium" => Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => false))) @@ -129,8 +128,8 @@ include("h5_metadata_check.jl") @test_logs (:warn,) match_mode = :any (ret = GeneralizedPerturbedEquilibrium.main([dir])) # The stage runs to completion; both sub-calculations are simply not populated. @test ret.pe !== nothing - @test isempty(ret.pe.permeability) - @test isempty(ret.pe.C_delta_prime) + @test isempty(ret.pe.response.permeability) + @test isempty(ret.pe.response.C_delta_prime) end end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index bc77acf04..88be30031 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -91,19 +91,18 @@ using TOML @test dp.matrix ≈ FFS.pest3_decompose(dp.raw).Δ @test dp.A !== nothing # the parity blocks the galerkin solve persists - # With no free-boundary δW, PE refuses the response rather than returning an empty state. - forcing = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat") - @test_throws ErrorException perturbed_equilibrium(api, RMPField(forcing); compute_response=true, write_outputs_to_HDF5=false) + # With no free-boundary δW, no plasma response can be built from it. + @test_throws ErrorException plasma_response(api) end end @testset "perturbed_equilibrium round-trips a forward result" begin forcing = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat") pe_section = Dict{String,Any}( - "compute_response" => true, "compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => false) + pe_kw = (; verbose=false, write_outputs_to_HDF5=false) mktempdir() do dir _stage_deck(dir; @@ -116,32 +115,51 @@ using TOML reference = GPEC.main([dir]).pe ffs = solve(equil, Forward(); nn=1, wall=wall, dir_path=dir, ffs_kwargs...) - pe = perturbed_equilibrium(ffs, RMPField(forcing); - (Symbol(k) => v for (k, v) in pe_section)...) + resp = plasma_response(ffs; verbose=false) + pe = perturbed_equilibrium(resp, RMPField(forcing); pe_kw...) # The response matrices depend on the equilibrium and the solve, not on the drive. - @test !isempty(pe.permeability) - @test pe.permeability ≈ reference.permeability rtol = 1e-10 - @test size(pe.C_delta_prime) == size(reference.C_delta_prime) + @test resp isa PlasmaResponse && pe.response === resp + @test !isempty(resp.permeability) + @test resp.permeability ≈ reference.response.permeability rtol = 1e-10 + @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. @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...) + @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(ffs, RMPField(forcing; scale=2.0); - (Symbol(k) => v for (k, v) in pe_section)...) + scaled = perturbed_equilibrium(resp, RMPField(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(ffs, 3 * RMPField(forcing) - RMPField(forcing); - (Symbol(k) => v for (k, v) in pe_section)...) + combo = perturbed_equilibrium(resp, 3 * RMPField(forcing) - RMPField(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...) + @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 + + # An explicit width is respected: it changes the regularized profiles, not the control-surface quantities. + @test resp.regularization_width == pe.regularization_width == 0.05 + raw = perturbed_equilibrium(resp, RMPField(forcing); regularization_width=0.0, pe_kw...) + @test raw.regularization_width == 0.0 + @test raw.forcing_b == pe.forcing_b + @test !(raw.xi_modes.clebsch_alpha ≈ pe.xi_modes.clebsch_alpha) end end @@ -274,11 +292,8 @@ using TOML @test matched.match !== nothing && matched.bpen == matched.match.bpen @test matched.wp === nothing && matched.free_boundary === nothing @test matched.solution === nothing # no Galerkin basis to build a matched ξ from - # PE refuses the response and coupling on it, but still materializes the forcing alone. - forcing = joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example", "forcing.dat") - @test_throws ErrorException perturbed_equilibrium(matched, RMPField(forcing); compute_singular_coupling=true, write_outputs_to_HDF5=false) - idle = perturbed_equilibrium(matched, RMPField(forcing); compute_response=false, compute_singular_coupling=false, write_outputs_to_HDF5=false) - @test idle isa GPEC.PerturbedEquilibrium.PerturbedEquilibriumState + # It carries no ideal δW, so no plasma response can be built from it. + @test_throws ErrorException plasma_response(matched) # A matched result is a new solution, not an outer solve to match again. @test_throws ErrorException MatchProblem(matched; ideal=true) end From d09bf3c220a9df93df1e99532627eab868ae004c Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Sat, 24 Oct 2026 16:11:38 +0200 Subject: [PATCH 2/6] ForcingTerms - DOCS - Name the forced perturbed-equilibrium stage in the forcing normalization docstring --- src/ForcingTerms/CoilFourier.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 7a450912b..b90c88218 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -404,7 +404,7 @@ end convert_forcing_normalization(modes, from_tag, equil, n, m_out; psi, mtheta, nzeta) -> Vector{ForcingMode} Convert `ForcingMode` amplitudes from `from_tag` normalization to the internal -unit-norm convention (= Fortran `Phi_x`), expected by `compute_plasma_response!`. +unit-norm convention (= Fortran `Phi_x`) that `perturbed_equilibrium` applies the plasma response to. If `from_tag == "sfl_flux_Wb"`, the user has provided amplitudes in the 2π-angle SFL-flux convention. These are scaled by (2π)² to reach unit-norm, keeping their harmonics. From 85cf56210ce1ffb4d96b5db53cd8f32fd9170b0c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sun, 25 Oct 2026 20:38:32 -0400 Subject: [PATCH 3/6] Docs - DOCS - Describe the kinetic regularization width as a default Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- docs/src/kinetic_forces.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 3b2ecf710..8e0d954ab 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -160,7 +160,7 @@ for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmet Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own notation. -GPEC therefore **forces `regularization_width = 0` whenever `kinetic_factor > 0`**, logging the override. +GPEC therefore **defaults `regularization_width` to 0 whenever `kinetic_factor > 0`**; an explicit deck value is respected. Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` is never regularized, so damping the other two breaks their near-resonance cancellation in ``\delta B/B`` and leaves a spurious residue driving the NTV integrand. From 60ed66bc4d6ddec2f141dac6675257da88c44158 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Mon, 26 Oct 2026 14:52:18 +0100 Subject: [PATCH 4/6] PerturbedEquilibrium - REFACTOR! - Default the regularization width to 0 for every result --- .../benchmark_diiid_ideal_ntv_torque.jl | 1 + docs/src/api.md | 4 ++++ docs/src/kinetic_forces.md | 6 ++++-- docs/src/workflow.md | 10 ++++++---- .../gpec.toml | 2 +- .../gpec.toml | 2 +- examples/Solovev_ideal_example/gpec.toml | 2 +- src/ForcingTerms/ForcingTerms.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 15 ++++++--------- .../FieldReconstruction.jl | 2 +- .../PerturbedEquilibrium.jl | 11 +++-------- .../PerturbedEquilibriumStructs.jl | 17 +++++------------ test/runtests_dominant_coupling.jl | 4 ++-- test/runtests_kinetic_validity.jl | 7 ------- test/runtests_solve_api.jl | 19 ++++++++++--------- 15 files changed, 46 insertions(+), 58 deletions(-) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 430c353cc..2af0d2419 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -264,6 +264,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir()) # Construct PerturbedEquilibriumState with Fortran Clebsch data. # Fortran gpec_xclebsch already writes ξ^α = xms/χ₁, exactly what the KF operator consumes — feed directly. pe_state = PE.PerturbedEquilibriumState(; + response=PE.PlasmaResponse(; ffs), psi_grid=psi_grid_f, xi_modes=( clebsch_psi=clebsch_psi, diff --git a/docs/src/api.md b/docs/src/api.md index af898b589..5c99bae0d 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -137,6 +137,10 @@ GeneralizedPerturbedEquilibrium.solve GeneralizedPerturbedEquilibrium.perturbed_equilibrium ``` +[`plasma_response`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response), which builds the +[`PlasmaResponse`](@ref GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PlasmaResponse) that `perturbed_equilibrium` +applies, is documented with the [perturbed equilibrium](perturbed_equilibrium.md) module. + ## Integrator selectors ```@autodocs diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 8e0d954ab..3f08e32e6 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -141,7 +141,7 @@ terms respectively, and do not modify the stored kinetic profile splines. when computing the `toroidal_rotation_factor` back-solve. Julia uses a clean reimplementation with consistent pre-scaling derivatives throughout. -## Regularization: why kinetic runs set `regularization_width = 0` +## Regularization: off by default, and never needed by kinetic runs `[PerturbedEquilibrium] regularization_width` smooths the displacements before they drive the NTV integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by @@ -160,7 +160,9 @@ for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmet Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own notation. -GPEC therefore **defaults `regularization_width` to 0 whenever `kinetic_factor > 0`**; an explicit deck value is respected. +GPEC therefore **defaults `regularization_width` to 0** for every result, as the Fortran `reg_flag = f` +template does. Set it, typically to 0.05, only to compute the NTV torque of an ideal solve; a kinetic run +should leave it at 0. Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` is never regularized, so damping the other two breaks their near-resonance cancellation in ``\delta B/B`` and leaves a spurious residue driving the NTV integrand. diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 85d8f6692..1d729c962 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -113,14 +113,16 @@ The single `gpec.toml` file supplies user-selected options to every module. The **Module**: `PerturbedEquilibrium` -**Purpose**: Compute the self-consistent plasma response to the external perturbation specified in Stage 4. Uses the force-free eigenfunctions and vacuum matrices to project the response into mode space, then evaluates singular coupling diagnostics at each rational surface. +**Purpose**: Compute the self-consistent plasma response to the external perturbation specified in Stage 4, in two steps: +- `plasma_response(ffs)` builds the source-independent operators once from the force-free eigenfunctions and vacuum matrices: the permeability, inductance and reluctance matrices on the control surface, and the singular coupling matrices at each rational surface (`PlasmaResponse`). +- `perturbed_equilibrium(resp, rmp)` applies them to one forcing field: it projects the response into mode space and evaluates the singular coupling diagnostics for that field. One response serves any number of fields. **Configuration**: `[PerturbedEquilibrium]` section of `gpec.toml` **Inputs**: -- ForceFreeStates results from Stage 3 (eigenfunctions, singular surface data) -- Vacuum response matrices from Stage 2 -- `ForcingMode` array from Stage 4 +- ForceFreeStates results from Stage 3 (eigenfunctions, singular surface data), for `plasma_response` +- Vacuum response matrices from Stage 2, for `plasma_response` +- The forcing field from Stage 4, for `perturbed_equilibrium` **Outputs** (`PerturbedEquilibriumState`): - Mode-space normal displacement ξ(m, ψ) across the plasma diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index e27abe72d..e86e2e6af 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -96,4 +96,4 @@ currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.0 # Regularization width for singular surfaces (0 = disabled) diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index e3ab0ea88..14b1350bb 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -70,4 +70,4 @@ currents = [1800.0, -50.0, -450.0, -800.0, -450.0, -50.0] # Sum of n=1, 2, 3 co compute_singular_coupling = true # Compute singular layer coupling metrics verbose = false # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.0 # Regularization width for singular surfaces (0 = disabled) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 5e299b068..8881e9ff2 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -39,7 +39,7 @@ forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Bi compute_singular_coupling = true # Compute singular layer coupling metrics verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 -regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +regularization_width = 0.0 # Regularization width for singular surfaces (0 = disabled) # Note: SLAYER tearing analysis is not run on the Solovev analytic equilibrium # — its Δ' / inner-layer dispersion does not yield meaningful tearing roots. diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index b6f7375be..f4e00fab8 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -289,7 +289,7 @@ evaluated on the control surface and the mode amplitudes are summed. nominal = RMPField("nominal_efc.dat") weld_field = RMPField("weld_fields.h5"; scale=0.5) total = 2.0 * nominal + weld_field # lazy: records terms and weights, computes nothing -pe = perturbed_equilibrium(ffs, total) # materializes each term, sums, drives PE +pe = perturbed_equilibrium(plasma_response(ffs), total) # materializes each term, sums, drives PE ``` Amplitude lives with the forcing description itself: per-conductor currents for the coil diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 7eb8b0a64..9bf740153 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -1005,7 +1005,7 @@ function run_perturbed_equilibrium( response_keys = ("compute_singular_coupling", "compute_torque_response") resp = plasma_response(result; (Symbol(k) => v for (k, v) in section if k in response_keys || k == "verbose")...) pe_state = perturbed_equilibrium(resp, ForcingTerms.RMPField(ft_ctrl); - coil_sets=preloaded_coil_sets, runtimes=runtimes, + coil_sets=preloaded_coil_sets, (Symbol(k) => v for (k, v) in section if !(k in response_keys))...) end @@ -1014,7 +1014,7 @@ end """ perturbed_equilibrium(resp, rmp; regularization_width=resp.regularization_width, verbose=true, - write_outputs_to_HDF5=true, output_filename="", coil_sets=nothing, runtimes=nothing) + write_outputs_to_HDF5=true, output_filename="", coil_sets=nothing) -> PerturbedEquilibriumState Drive the plasma response `resp` (from [`plasma_response`](@ref)) with the external field `rmp`, an @@ -1022,11 +1022,9 @@ Drive the plasma response `resp` (from [`plasma_response`](@ref)) with the exter fields, so a scan rebuilds no operator. `regularization_width` smooths the ideal 1/(m−nq) divergence in the profile reconstruction (Fortran -`reg_spot`); its default is 0 for a kinetic result and 0.05 otherwise. `write_outputs_to_HDF5` appends +`reg_spot`); its default is 0, and 0.05 is typical for driving KineticForces with an ideal solve. `write_outputs_to_HDF5` appends the outputs to `output_filename`, or to the force-free-states file when that is empty. `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`. +injects already-built coil geometry (the gpec.h5 replay path), bypassing the geometry read. ```julia resp = plasma_response(ffs) @@ -1040,15 +1038,14 @@ function perturbed_equilibrium( verbose::Bool=true, write_outputs_to_HDF5::Bool=true, output_filename::String="", - coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing, - runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing + coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing ) ffs = resp.ffs ctrl = ffs.control # 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 = PerturbedEquilibrium.compute_perturbed_equilibrium(resp, rmp; regularization_width, verbose, - preloaded_coil_sets=coil_sets === nothing ? ForcingTerms.CoilSet[] : copy(coil_sets), runtimes) + preloaded_coil_sets=coil_sets === nothing ? ForcingTerms.CoilSet[] : copy(coil_sets)) # Write perturbed equilibrium outputs to same HDF5 file if write_outputs_to_HDF5 diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 18e634ed6..14cf6d7cd 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -19,7 +19,7 @@ where χ₁ = 2π·Ψ₀ [Park Phys. Plasmas 14, 052110 (2007) eq. 8-10]. Clebsch displacement components for PENTRC (matches Fortran gpout_xclebsch): ξ^ψ = xsp_mn (unregularized) ∂ξ^ψ/∂ψ = xmp1_mn (regularized: xsp1 * singfac²/(singfac² + reg_spot²)) - regularization_width = 0 in kinetic runs — no ideal singularity to smooth + unregularized at the default regularization_width = 0 ξ^α = xms_mn (regularized: -A⁻¹(B·xmp1 + C·xsp) ideal, ξ_s·singfac factor kinetic; divided by χ₁ in output) Contravariant displacement from Jacobian convolution (matches Fortran gpeq_contra): diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 4109eb207..ebb845221 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -66,7 +66,7 @@ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling:: # A run has at most one ξ solution, already closed at the rationals and carrying populated # Ξ′ / Ξ_s stores. The gal-native basis takes the analytic Galerkin Ξ′ downstream. solution = ffs.solution - resp = PlasmaResponse(; ffs, m_modes, n_modes, regularization_width=default_regularization_width(ffs), + resp = PlasmaResponse(; ffs, m_modes, n_modes, odet_from_gal=solution !== nothing && solution.basis === :gal_native) if ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") && @@ -83,7 +83,7 @@ end """ compute_perturbed_equilibrium(resp, forcing; regularization_width=resp.regularization_width, verbose=true, - preloaded_coil_sets=CoilSet[], runtimes=nothing)::PerturbedEquilibriumState + 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 @@ -96,21 +96,16 @@ forcing spectrum. - `regularization_width`: width of the ideal 1/(m−nq) smoothing in the profile reconstruction (Fortran `reg_spot`); must be ≥ 0 - `preloaded_coil_sets`: coil geometry to integrate instead of reading it (the gpec.h5 replay path) - - `runtimes`: optional collector; receives `"forcing_terms" => dt` for the forcing-mode - materialization (coil Biot-Savart onto the plasma surface) """ function compute_perturbed_equilibrium( resp::PlasmaResponse, forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}; regularization_width::Real=resp.regularization_width, verbose::Bool=true, - preloaded_coil_sets::Vector{CoilSet}=CoilSet[], - runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing + preloaded_coil_sets::Vector{CoilSet}=CoilSet[] )::PerturbedEquilibriumState ffs = resp.ffs - ft_start = time() modes, coil_sets = materialize_forcing_modes(ffs, forcing; dir_path=ffs.dir_path, preloaded_coil_sets, verbose) - runtimes === nothing || push!(runtimes, "forcing_terms" => time() - ft_start) state = PerturbedEquilibriumState(; response=resp, forcing_modes=modes, coil_sets=isempty(coil_sets) ? preloaded_coil_sets : coil_sets, regularization_width=Float64(regularization_width)) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 5681131a8..29efcea89 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -48,9 +48,10 @@ spectra. - `ffs::ForceFreeStatesResult` - The result the operators were built from. - `m_modes`, `n_modes` - Poloidal and toroidal mode number of each entry of a spectrum [numpert_total]. - - `regularization_width::Float64` - Default profile-reconstruction width (Fortran `reg_spot`): 0 for a - kinetic result, whose Euler-Lagrange operator has no ideal 1/(m−nq) divergence to smooth (Park & Logan, - Phys. Plasmas 24, 032505 (2017) §III D), 0.05 otherwise. + - `regularization_width::Float64` - Default profile-reconstruction width (Fortran `reg_spot`): 0, no + regularization (Fortran `reg_flag = f`). Set it, typically to 0.05, only to drive KineticForces with an + ideal solve, whose 1/(m−nq) divergence leaves the NTV torque integral unconverged; a kinetic solve has + none to smooth (Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D). - `odet_from_gal::Bool` - The solution is the gal-native matched basis, whose `du_store` carries the analytic Galerkin Ξ′; selects that branch of the singular-coupling Ξ′ evaluation. - `flux_matrix::Matrix{ComplexF64}` - Control-surface flux of the free-boundary eigenmodes [mode × eigenmode]. @@ -103,7 +104,7 @@ on a `ResonantCoupling`): ffs::ForceFreeStatesResult m_modes::Vector{Int} = Int[] n_modes::Vector{Int} = Int[] - regularization_width::Float64 = 5e-2 + regularization_width::Float64 = 0.0 odet_from_gal::Bool = false # Flux-space operators [numpert_total × numpert_total] @@ -147,14 +148,6 @@ on a `ResonantCoupling`): T_xe::Array{ComplexF64,3} = zeros(ComplexF64, 0, 0, 0) end -""" - default_regularization_width(ffs) -> Float64 - -Profile-reconstruction width for a result: 0 when it carries a kinetic fit, which removes the ideal -resonant singularity, 0.05 otherwise. `ffs` is duck-typed: only `ffs.mats.kinetic` is read. -""" -default_regularization_width(ffs) = ffs.mats.kinetic === nothing ? 5e-2 : 0.0 - """ PerturbedEquilibriumState diff --git a/test/runtests_dominant_coupling.jl b/test/runtests_dominant_coupling.jl index 21ed66108..edeb23f8f 100644 --- a/test/runtests_dominant_coupling.jl +++ b/test/runtests_dominant_coupling.jl @@ -85,8 +85,8 @@ include("h5_metadata_check.jl") resp = pe.response h5path = joinpath(dir, "gpec.h5") - # Run summary: the core-window decomposition, rebuilt into the stored applied resonant - # field (which the writer evaluates before conforming the matrix to b̃ space). + # Run summary: the core-window decomposition, rebuilt into the applied resonant field + # that perturbed_equilibrium evaluates from the flux-space coupling. @test !isempty(resp.dominant_singular_values) rank = length(resp.dominant_singular_values) rows = resp.dominant_rational_index diff --git a/test/runtests_kinetic_validity.jl b/test/runtests_kinetic_validity.jl index 86044955e..d2e9aca16 100644 --- a/test/runtests_kinetic_validity.jl +++ b/test/runtests_kinetic_validity.jl @@ -3,7 +3,6 @@ GPE = GeneralizedPerturbedEquilibrium KF = GPE.KineticForces EQ = GPE.Equilibrium - PE = GPE.PerturbedEquilibrium # Duck-typed stand-in for PlasmaEquilibrium: the boundary scan reads only these fields. mock_equil(q, r, R; ro=1.7, b0=2.0) = @@ -90,10 +89,4 @@ @test psi_c_spiked ≈ psi_c @test psi_c_spiked < 0.5 end - - @testset "default_regularization_width" begin - # A kinetic fit removes the ideal resonant singularity, so its default is no regularization. - @test PE.default_regularization_width((; mats=(; kinetic=nothing))) == 0.05 - @test PE.default_regularization_width((; mats=(; kinetic=(;)))) == 0.0 - end end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index fbbc08909..734936879 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -188,12 +188,12 @@ using HDF5 @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 - # An explicit width is respected: it changes the regularized profiles, not the control-surface quantities. - @test resp.regularization_width == pe.regularization_width == 0.05 - raw = perturbed_equilibrium(resp, RMPField(forcing); regularization_width=0.0, pe_kw...) - @test raw.regularization_width == 0.0 - @test raw.forcing_b == pe.forcing_b - @test !(raw.xi_modes.clebsch_alpha ≈ pe.xi_modes.clebsch_alpha) + # No regularization by default; an explicit width changes the regularized profiles, not the control-surface quantities. + @test resp.regularization_width == pe.regularization_width == 0.0 + reg = perturbed_equilibrium(resp, RMPField(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) end end @@ -367,7 +367,7 @@ using HDF5 # A known drift is `@test_broken`, so the change that removes it surfaces as an unexpected pass. @testset "deck and API parity" begin @testset "hdf5 file forcing" begin - pe_kwargs = (compute_response=true, compute_singular_coupling=true, verbose=false, write_outputs_to_HDF5=false) + pe_kwargs = (compute_singular_coupling=true, verbose=false, write_outputs_to_HDF5=false) mktempdir() do dir modes = GPEC.ForcingTerms.ForcingMode[] GPEC.ForcingTerms.load_forcing_data!(modes, joinpath(@__DIR__, "..", "examples", "Solovev_ideal_example"), "forcing.dat", "ascii", false) @@ -377,9 +377,10 @@ using HDF5 "ForcingTerms" => Dict{String,Any}("forcing_data_file" => "forcing.h5", "forcing_data_format" => "hdf5"), "PerturbedEquilibrium" => Dict{String,Any}(String(k) => v for (k, v) in pairs(pe_kwargs)))) reference = GPEC.main([dir]) - pe = perturbed_equilibrium(reference.ffs, RMPField(joinpath(dir, "forcing.h5")); pe_kwargs...) + resp = plasma_response(reference.ffs; compute_singular_coupling=true, verbose=false) + pe = perturbed_equilibrium(resp, RMPField(joinpath(dir, "forcing.h5")); verbose=false, write_outputs_to_HDF5=false) - @test pe.permeability ≈ reference.pe.permeability rtol = 1e-10 + @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. @test pe.forcing_b ≈ reference.pe.forcing_b rtol = 1e-10 end From 5e4981075f6c1f5a58227c356a44ede4315ec4ba Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Mon, 26 Oct 2026 15:18:35 +0100 Subject: [PATCH 5/6] PerturbedEquilibrium - TEST - Check the regularization default and an explicit width on a kinetic result --- test/runtests_torque_response.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/runtests_torque_response.jl b/test/runtests_torque_response.jl index a315e28b5..54e5edb15 100644 --- a/test/runtests_torque_response.jl +++ b/test/runtests_torque_response.jl @@ -153,6 +153,16 @@ write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to @test ts.T_source[1, 1, :] ./ 2 ≈ tr.T_applied rtol = 1e-8 end + @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, + coil_sets=pe.coil_sets, verbose=false, write_outputs_to_HDF5=false) + @test reg.regularization_width == 0.05 + @test reg.forcing_b ≈ pe.forcing_b rtol = 1e-12 + @test !(reg.xi_modes.clebsch_alpha ≈ pe.xi_modes.clebsch_alpha) + end + @testset "argument checks" begin @test_throws DimensionMismatch _PE_TR.torque_profile(tr, ComplexF64[1.0]) # A response without its operators cannot carry a torque response. From 10c7350db4e21eb59d3c0d4a1b7d1f159810aef3 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Tue, 27 Oct 2026 09:40:36 +0100 Subject: [PATCH 6/6] Test - TEST - Stage the analysis-plots deck without the retired compute_response key --- test/runtests_analysis_plots.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 63bd4b598..e9c0a50fc 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -32,7 +32,7 @@ using TOML "ForcingTerms" => Dict{String,Any}("forcing_data_format" => "coil", "mtheta_coil" => 240, "nzeta_coil" => 32, "coil_set" => [Dict{String,Any}("name" => "hoop", "source" => "pf_hoop", "radius" => 1.5, "height" => 0.4, "currents" => [2.0e3], "tiltx" => [3.0])]), - "PerturbedEquilibrium" => Dict{String,Any}("compute_response" => true, "compute_singular_coupling" => true, + "PerturbedEquilibrium" => Dict{String,Any}("compute_singular_coupling" => true, "verbose" => false, "write_outputs_to_HDF5" => true))) for f in (A.Equilibrium.plot_qprofile, A.Equilibrium.plot_pressure_profile, A.Equilibrium.plot_f_profile,