diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0147177ab..d82cfe6c1 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|local_stability_flag)\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|local_stability_flag)\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 a77d13d02..2b51564bb 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/benchmarks/benchmark_against_fortran_run.jl @@ -440,7 +440,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") @@ -492,7 +491,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/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index ade40501f..f2b4c1611 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -265,6 +265,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/development/architecture.md b/docs/development/architecture.md index 150a9d554..832dc8634 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -215,12 +215,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 ec71fc00c..160b9debf 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 take their parameters from their config object, with the @@ -42,8 +43,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 @@ -57,13 +59,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 @@ -154,6 +158,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/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/kinetic_forces.md b/docs/src/kinetic_forces.md index e1a9c06d3..ced2a590d 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -153,7 +153,7 @@ The torque is three nested adaptive Gauss–Kronrod integrals: over ψ (paneled QuadGK bounds the error of the whole integral, not of each step as the LSODE integrators of Fortran PENTRC did, so a level tighter than its margin only multiplies evaluations: the ψ quadrature cannot see λ-level noise below `rtol_psi/10`. The calculated kinetic matrix build (`kinetic_source = "calculated"`) has no ψ quadrature; it uses `rtol_pitch` when set and 1e-3 otherwise. Tightening `rtol_psi` tightens the whole ladder. -## 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 @@ -172,7 +172,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 **forces `regularization_width = 0` whenever `kinetic_factor > 0`**, logging the override. +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/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index 4bbdd6081..942b65adb 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! ``` ## Torque response matrices @@ -86,8 +89,9 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_dominant_coupling! A kinetic forward solve (`integrator = "forward"`, `kinetic_factor > 0`) carries the drift-kinetic response inside its Euler-Lagrange energy: the complex plasma energy inside a flux surface, `δW(ψ) = ξ†·U₂(ψ)·U₁(ψ)⁻¹·ξ/(2μ₀)`, has an anti-Hermitian part that is the -neoclassical toroidal viscosity torque (Logan & Park 2013, eq. 19). With -`compute_torque_response = true` the response stage refers that quadratic form back to the +neoclassical toroidal viscosity torque (Logan & Park 2013, eq. 19). On request, +`plasma_response(ffs; compute_torque_response=true)` (deck: `compute_torque_response = true` +in `[PerturbedEquilibrium]`) refers that quadratic form back to the applied root-area-weighted control-surface field b̃ — through the permeability, the boundary flux-to-displacement factor `χ₁(m − n·q_lim)·2πi` and the b̃ → Φ conform `R` — and stores the result at every stored node as the **torque response matrix** `T_xe(ψ)`, cumulative @@ -105,8 +109,9 @@ Two conventions are worth stating explicitly: its eigenvectors are the applied spectra of extremal torque and its eigenvalues their torque per unit `‖b̃‖²/2`. The anti-Hermitian part carries the energy. -The forcing-source form `T_source(ψ) = M†·T·M` is written alongside, with one column of `M` per -forcing source of the run's `[ForcingTerms]`, each as built: every coil set at its deck currents, +Each `perturbed_equilibrium(resp, rmp)` call contracts it with that field into a `TorqueResponse`: +the applied torque `T_applied` and the forcing-source form `T_source(ψ) = M†·T·M`, with one column +of `M` per forcing source, each as built: every coil set at its deck currents, or the forcing file's spectrum as a single source. Driving the sources at scale factors `s` gives `s†·T_source·s/2`. @@ -124,10 +129,11 @@ T_h = Hermitian((tr.T_xe[:, :, end] + tr.T_xe[:, :, end]') / 2) ```@docs GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.TorqueResponse -GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_response_matrices +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_response_matrix GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_profile GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coil_flux_spectra -GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_torque_response! +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.build_torque_response! +GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.apply_torque_response! ``` ## Plotting per-surface results against ψ or q 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_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index ef1800049..8282d285f 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/examples/DIIID-like_error_field_example/gpec.toml @@ -206,7 +206,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 5121cb5ef..a2c4a6f34 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 @@ -103,8 +103,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) +regularization_width = 0.0 # 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 d6bd8ef24..d706bbed7 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -75,7 +75,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 d2f04d6a3..e4cb5c343 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -77,8 +77,7 @@ 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 -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 a81491fb7..eded9dce8 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -36,11 +36,10 @@ 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 -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/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 1ee11181b..1ef0d3f5c 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -67,7 +67,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 8b74f409d..024ad168a 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 f395d34d6..60a9fe8b6 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -34,12 +34,10 @@ 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 compute_torque_response = true # Build the ψ-resolved torque response matrices T_xe and T_source 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) # Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords [LocalStability] diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index 803622ba9..679ecd4a9 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 eb66ad588..283bfdc7c 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). """ @@ -404,7 +404,7 @@ end convert_forcing_normalization(modes, from_tag, equil, n, m_out; psi, mtheta) -> 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. diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index 313fa0871..f4e00fab8 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -268,12 +268,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 @@ -287,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 @@ -319,20 +321,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 """ @@ -354,14 +373,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 9510681b1..886d551c8 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -92,6 +92,7 @@ using .Tearing.Runner: TearingProblem using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult 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", @@ -99,7 +100,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") const _DEPRECATED_KEYS = ("Equilibrium" => _DEPRECATED_EQUIL_KEYS, "ForceFreeStates" => _DEPRECATED_FFS_KEYS, "KineticForces" => _DEPRECATED_KF_KEYS, "PerturbedEquilibrium" => _DEPRECATED_PE_KEYS) @@ -994,70 +995,66 @@ 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. + 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, - (Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])...) + (Symbol(k) => v for (k, v) in section if !(k in response_keys))...) end return pe_state end """ - perturbed_equilibrium(ffs, rmp; coil_sets=nothing, kwargs...) -> PerturbedEquilibriumState + perturbed_equilibrium(resp, rmp; regularization_width=resp.regularization_width, verbose=true, + write_outputs_to_HDF5=true, output_filename="", coil_sets=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. +`regularization_width` smooths the ideal 1/(m−nq) divergence in the profile reconstruction (Fortran +`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. ```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; - coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing, - kwargs... + 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 ) + 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) + pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(resp, rmp; regularization_width, verbose, + preloaded_coil_sets=coil_sets === nothing ? ForcingTerms.CoilSet[] : copy(coil_sets)) # 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 @@ -1197,7 +1194,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) @@ -1252,7 +1249,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" @@ -1275,22 +1273,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, @@ -1305,15 +1304,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. @@ -1335,10 +1335,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 @@ -1878,8 +1875,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 export LocalStabilityProblem, Ballooning, LocalStabilityResult diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index c4f073bc1..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): @@ -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 aa58cafa8..ebb845221 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -29,99 +29,92 @@ 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 -export TorqueResponse, torque_response_matrices, torque_profile +export TorqueResponse, torque_response_matrix, torque_profile """ - compute_perturbed_equilibrium(ffs, forcing, ctrl, intr)::PerturbedEquilibriumState - -Main entry point for perturbed equilibrium calculations. - -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 - -## Returns - - - `PerturbedEquilibriumState`: Calculation results + plasma_response(ffs; compute_singular_coupling=true, compute_torque_response=false, verbose=true) -> PlasmaResponse + +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)), and with `compute_torque_response` the +ψ-resolved torque response `T_xe` of a single-`n` forward kinetic solve. Apply it to an external field +with `perturbed_equilibrium(resp, rmp)` as many times as needed. + +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 -)::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, compute_torque_response::Bool=false, + 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) - modes, coil_sets = materialize_forcing_modes(ffs, forcing; - dir_path=intr.dir_path, preloaded_coil_sets=intr.coil_sets, verbose=ctrl.verbose) - intr.forcing_modes = modes - isempty(coil_sets) || (intr.coil_sets = coil_sets) - end + resp = PlasmaResponse(; ffs, m_modes, n_modes, + 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 + compute_torque_response && build_torque_response!(resp; verbose) + return resp +end - # Torque response matrices of a kinetic forward solve, on the plasma response just computed - if ctrl.compute_torque_response - ctrl.compute_response || throw(ArgumentError("compute_torque_response needs compute_response = true")) - ForceFreeStates.require_solution(ffs, "torque response matrices") || - throw(ArgumentError("compute_torque_response needs the dense forward-integrated ξ solution")) - compute_torque_response!(state, ffs, forcing, intr, ctrl) - end +""" + compute_perturbed_equilibrium(resp, forcing; regularization_width=resp.regularization_width, verbose=true, + preloaded_coil_sets=CoilSet[])::PerturbedEquilibriumState - # 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 +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, reconstruct the perturbed profiles, and contract any torque response with +the forcing sources. Steps whose operators `resp` does not carry are skipped, leaving only the +forcing spectrum. +## 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) +""" +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[] +)::PerturbedEquilibriumState + ffs = resp.ffs + modes, coil_sets = materialize_forcing_modes(ffs, forcing; dir_path=ffs.dir_path, preloaded_coil_sets, verbose) + 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) + isempty(resp.T_xe) || apply_torque_response!(state, forcing; 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 @@ -193,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 7ddbb64d1..29efcea89 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -34,124 +34,152 @@ struct TorqueResponse end """ - 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, 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]. + - `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) - - `compute_torque_response::Bool` - Build the ψ-resolved torque response matrices `T_xe` and `T_source` (default: false). Needs a forward-integrated kinetic solve (`integrator = "forward"`, `kinetic_factor > 0`) and `compute_response`; errors otherwise - - `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 - compute_torque_response::Bool = false - 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] + - `T_xe` - ψ-resolved torque response on b̃ of a single-`n` forward kinetic solve, cumulative in ψ + [numpert_total × numpert_total × npsi]; built only when `plasma_response` is asked for it, else empty. + Each `perturbed_equilibrium` call contracts it with its forcing into a [`TorqueResponse`](@ref). """ -@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 = 0.0 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[] + + # Torque response of a kinetic forward solve, opt-in [numpert_total × numpert_total × npsi] + T_xe::Array{ComplexF64,3} = zeros(ComplexF64, 0, 0, 0) end """ 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): @@ -160,19 +188,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, Λ: @@ -192,6 +207,11 @@ Torque response (kinetic forward solves with `compute_torque_response`): applied b̃ and on the coil sets, with the run's own forcing contracted; see [`TorqueResponse`](@ref) """ @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[] @@ -201,13 +221,6 @@ Torque response (kinetic forward solves with `compute_torque_response`): 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[] @@ -215,25 +228,11 @@ Torque response (kinetic forward solves with `compute_torque_response`): 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) @@ -243,14 +242,6 @@ Torque response (kinetic forward solves with `compute_torque_response`): 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 @@ -260,25 +251,3 @@ Torque response (kinetic forward solves with `compute_torque_response`): # ψ-resolved torque response matrices of a kinetic forward solve torque_response::Union{Nothing,TorqueResponse} = nothing 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/TorqueResponse.jl b/src/PerturbedEquilibrium/TorqueResponse.jl index fd6cc7d70..a7adbd2be 100644 --- a/src/PerturbedEquilibrium/TorqueResponse.jl +++ b/src/PerturbedEquilibrium/TorqueResponse.jl @@ -10,13 +10,11 @@ eigenvectors are the applied spectra of extremal torque [Logan 2015, ch. 7]. """ """ - torque_response_matrices(ffs, permeability, flux_conform; - forcing_b_rootarea=ComplexF64[], source_flux=zeros(ComplexF64, N, 0), source_names=String[]) -> TorqueResponse + torque_response_matrix(ffs, permeability, flux_conform) -> Array{ComplexF64,3} -Build the ψ-resolved torque response of a forward kinetic solve from its stored solution. +ψ-resolved torque response `T_xe` of a forward kinetic solve from its stored solution. `permeability` is the flux-space plasma response operator `P` (Φ_tot = P·Φ_x) and `flux_conform` -the b̃ → Φ operator `R = S·A`, both at the control surface; `source_flux` holds one unit-norm -Φ_x column per forcing source and `forcing_b_rootarea` the run's applied b̃ for `T_applied`. +the b̃ → Φ operator `R = S·A`, both at the control surface. Per stored node `k`, with `U₁`, `U₂` the displacement and momentum blocks of `u_store` and `U₁ₗ` the displacement block at the control surface: @@ -26,20 +24,12 @@ Per stored node `k`, with `U₁`, `U₂` the displacement and momentum blocks of `K` maps an applied b̃ to the coefficients of the solution basis (b̃ → Φ_x → Φ_tot → boundary ξ → coefficients), so `T_xe` is the Fortran `gresp` chain (`dwks`, `bsurfmat`, `gind`, `gres`, `ptof`) with the single edge inverse the stored, globally consistent basis allows. -`T_source = Ks†·[U₁†U₂]·Ks · 2n·i/(2μ₀)` with `Ks = U₁ₗ⁻¹·D·P·source_flux` is the Fortran `gcoil`. Errors unless the solve is a single-`n` forward integration with kinetic matrices; the ideal Euler-Lagrange energy is Hermitian and carries no torque. Warns when the kinetic matrices are the synthetic `"fixed"` test matrices. """ -function torque_response_matrices( - ffs::ForceFreeStatesResult, - permeability::AbstractMatrix{ComplexF64}, - flux_conform::AbstractMatrix{ComplexF64}; - forcing_b_rootarea::AbstractVector{ComplexF64}=ComplexF64[], - source_flux::AbstractMatrix{ComplexF64}=zeros(ComplexF64, ffs.numpert_total, 0), - source_names::Vector{String}=String[] -)::TorqueResponse +function torque_response_matrix(ffs::ForceFreeStatesResult, permeability::AbstractMatrix{ComplexF64}, flux_conform::AbstractMatrix{ComplexF64}) ffs.integrator === :forward || throw(ArgumentError("torque response matrices need the dense forward-integrated solution; this result came from the $(ffs.integrator) integrator")) ffs.solution === nothing && throw(ArgumentError("torque response matrices need a ξ solution; the forward integrator produced none")) @@ -48,8 +38,6 @@ function torque_response_matrices( ffs.npert == 1 || throw(ArgumentError("torque response matrices are defined per toroidal harmonic; this solve stacks $(ffs.npert) (nn_low != nn_high)")) ffs.control.kinetic_source == "fixed" && @warn "kinetic_source = \"fixed\": the torque response is built from synthetic test matrices and has no physical meaning" - size(source_flux, 2) == length(source_names) || - throw(DimensionMismatch("source_flux has $(size(source_flux, 2)) columns but $(length(source_names)) source names were given")) solution = ffs.solution N = ffs.numpert_total @@ -58,29 +46,21 @@ function torque_response_matrices( μ₀ = 4π * 1e-7 chi1 = 2π * ffs.equil.psio m_modes = [(i - 1) % ffs.mpert + ffs.mlow for i in 1:N] - n_modes = fill(nn, N) # Boundary flux → boundary displacement, Φ_tot = χ₁·(m − n·q_lim)·2π·i·ξ (the same factor Λ carries). t = [-im / (chi1 * (m_modes[i] - nn * ffs.qlim) * 2π) for i in 1:N] U1_lim = solution.u_store[:, :, 1, npsi] - boundary_xi = Diagonal(t) * permeability - K = U1_lim \ (boundary_xi * flux_conform) - Ks = U1_lim \ (boundary_xi * source_flux) + K = U1_lim \ (Diagonal(t) * permeability * flux_conform) - nsource = size(source_flux, 2) T_xe = Array{ComplexF64}(undef, N, N, npsi) - T_source = Array{ComplexF64}(undef, nsource, nsource, npsi) prefactor = 2 * nn * im / (2 * μ₀) W = Matrix{ComplexF64}(undef, N, N) for k in 1:npsi mul!(W, adjoint(@view(solution.u_store[:, :, 1, k])), @view(solution.u_store[:, :, 2, k])) W .*= prefactor T_xe[:, :, k] = K' * W * K - nsource > 0 && (T_source[:, :, k] = Ks' * W * Ks) end - - T_applied = isempty(forcing_b_rootarea) ? ComplexF64[] : torque_profile(T_xe, forcing_b_rootarea) - return TorqueResponse(solution.psi_store[1:npsi], m_modes, n_modes, T_xe, source_names, T_source, T_applied) + return T_xe end """ @@ -128,43 +108,59 @@ function coil_flux_spectra(ffs::ForceFreeStatesResult, coil_sets::Vector{CoilSet end """ - compute_torque_response!(state, ffs, forcing, intr, ctrl) + build_torque_response!(resp; verbose=false) -Build the torque response of the solve from the plasma response already on `state` and -`intr` and store it in `state.torque_response`. Every forcing source is one column of -`T_source`: each coil set (`intr.coil_sets`, on the grid `forcing` carries) when the forcing is -coils, otherwise the run's applied spectrum as a single source named after its forcing file. -Errors when the plasma response was not computed. +Store the torque response `T_xe` of the solve on `resp`, from its flux-space permeability. +Errors when `resp` carries no plasma response operators. """ -function compute_torque_response!( +function build_torque_response!(resp::PlasmaResponse; verbose::Bool=false) + isempty(resp.permeability_flux) && + throw(ArgumentError("compute_torque_response needs the plasma response operators, which this result could not supply")) + verbose && @info "Computing torque response matrices" + resp.T_xe = torque_response_matrix(resp.ffs, resp.permeability_flux, resp.rootarea_to_area_weight .* resp.surface_area) + return resp +end + +""" + apply_torque_response!(state, forcing; verbose=false) + +Contract the torque response of `state.response` with the forcing of `state` into +`state.torque_response`. Every forcing source is one column of `T_source`: each coil set +(`state.coil_sets`, on the grid `forcing` carries) when the forcing is coils, otherwise the +applied spectrum as a single source named after its forcing file. A source's Φ_x column is +referred to b̃ through `R = S·A`, so `T_source = B̃†·T_xe·B̃` with `B̃ = R⁻¹·Φ_source` (the Fortran `gcoil`). +""" +function apply_torque_response!( state::PerturbedEquilibriumState, - ffs::ForceFreeStatesResult, - forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}, - intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl + forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}; + verbose::Bool=false ) - isempty(intr.plasma_response) && - throw(ArgumentError("compute_torque_response needs the plasma response (compute_response = true)")) - ctrl.verbose && @info "Computing torque response matrices" - flux_conform = state.rootarea_to_area_weight .* state.surface_area - + resp = state.response + ffs = resp.ffs ft_ctrl = forcing isa ForcingTerms.ForcingTermsControl ? forcing : forcing isa ForcingTerms.RMPSource ? forcing.ctrl : nothing - source_flux = zeros(ComplexF64, ffs.numpert_total, 0) + source_b = zeros(ComplexF64, ffs.numpert_total, 0) source_names = String[] - if !isempty(intr.coil_sets) && ft_ctrl !== nothing && ft_ctrl.forcing_data_format == "coil" - source_flux = coil_flux_spectra(ffs, intr.coil_sets, ForcingTerms.CoilConfig(ft_ctrl)) - source_names = [cs.name for cs in intr.coil_sets] + if !isempty(state.coil_sets) && ft_ctrl !== nothing && ft_ctrl.forcing_data_format == "coil" + flux_conform = resp.rootarea_to_area_weight .* resp.surface_area + source_b = flux_conform \ coil_flux_spectra(ffs, state.coil_sets, ForcingTerms.CoilConfig(ft_ctrl)) + source_names = [cs.name for cs in state.coil_sets] elseif !isempty(state.forcing_b_rootarea) - source_flux = reshape(flux_conform * state.forcing_b_rootarea, :, 1) + source_b = reshape(state.forcing_b_rootarea, :, 1) source_names = [ft_ctrl === nothing ? "forcing" : basename(ft_ctrl.forcing_data_file)] end - state.torque_response = torque_response_matrices(ffs, intr.plasma_response, flux_conform; - forcing_b_rootarea=state.forcing_b_rootarea, source_flux=source_flux, source_names=source_names) + T_xe = resp.T_xe + npsi = size(T_xe, 3) + T_source = Array{ComplexF64}(undef, length(source_names), length(source_names), npsi) + for k in 1:npsi + T_source[:, :, k] = source_b' * @view(T_xe[:, :, k]) * source_b + end + T_applied = isempty(state.forcing_b_rootarea) ? ComplexF64[] : torque_profile(T_xe, state.forcing_b_rootarea) + state.torque_response = TorqueResponse(ffs.solution.psi_store[1:npsi], copy(resp.m_modes), copy(resp.n_modes), T_xe, source_names, T_source, T_applied) - if ctrl.verbose - total = state.torque_response.T_applied[end] + if verbose && !isempty(T_applied) + total = T_applied[end] @info "Torque response complete: applied-forcing torque $(@sprintf("%.4e", real(total))) N·m, 2n·δW $(@sprintf("%.4e", imag(total))) J" end return state diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index e82a62d39..e9d6f6d73 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). @@ -116,23 +82,20 @@ PerturbedEquilibrium/ └── total_torque # Re T_applied at the control surface [N·m] ``` """ -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). @@ -148,12 +111,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") @@ -215,11 +178,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) @@ -228,23 +191,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_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, diff --git a/test/runtests_dominant_coupling.jl b/test/runtests_dominant_coupling.jl index fe7cab712..edeb23f8f 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 + # 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 + @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 216ca9647..da15fe93f 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -173,7 +173,7 @@ using TOML @test_throws r"`moment`, `save_records`, which are no longer keys" 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`, which is no longer a key" 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 # The stop runs through main, before any equilibrium work. mktempdir() do deck_dir @@ -187,7 +187,8 @@ using TOML # 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 57201abfb..860c9ff9f 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. @@ -131,7 +131,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 @@ -533,7 +533,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, @@ -634,7 +634,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..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,24 +89,4 @@ @test psi_c_spiked ≈ psi_c @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 - end end diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 076b62730..a0b51b271 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 56e9f34fe..5d33ded18 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -121,9 +121,8 @@ using HDF5 @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) reference = GPEC.main([_stage_deck(dir; ffs_overrides=Dict{String,Any}("integrator" => "galerkin", "gal_nx" => 32))]).ffs @test reference.galerkin.msing == msing @@ -135,10 +134,10 @@ using HDF5 @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; @@ -151,32 +150,51 @@ using HDF5 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 + + # 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 @@ -312,11 +330,8 @@ using HDF5 @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 @@ -353,7 +368,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) @@ -363,9 +378,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 diff --git a/test/runtests_torque_response.jl b/test/runtests_torque_response.jl index 2e827efbf..54e5edb15 100644 --- a/test/runtests_torque_response.jl +++ b/test/runtests_torque_response.jl @@ -35,7 +35,6 @@ gap_fraction = 0.15 # Toroidal gap between neighbouring fram currents = [1000.0, 0.0, -1000.0, -1000.0, 0.0, 1000.0] # n=1 phasing shifted by 60 degrees, 1 kA peak [A] [PerturbedEquilibrium] -compute_response = true # Compute plasma response to forcing compute_singular_coupling = false # Compute singular layer coupling metrics compute_torque_response = true # Build the ψ-resolved torque response matrices T_xe and T_source verbose = false # Enable verbose logging @@ -90,8 +89,9 @@ write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to # The applied forcing, referred back to the solution basis without the function under test: # Φ_x = R·b̃, boundary ξ = D·P·Φ_x with the flux-space permeability P recovered from the stored b̃ form. - R = pe.rootarea_to_area_weight .* pe.surface_area - P = R * pe.permeability / R + resp = pe.response + R = resp.rootarea_to_area_weight .* resp.surface_area + P = R * resp.permeability / R chi1 = 2π * ffs.equil.psio D = Diagonal([-im / (chi1 * (tr.m_modes[i] - nn * ffs.qlim) * 2π) for i in 1:N]) b̃ = pe.forcing_b_rootarea @@ -145,22 +145,32 @@ write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to @testset "forcing-source matrix: a forcing file is one source" begin file_state = deepcopy(pe) - file_intr = _PE_TR.PerturbedEquilibriumInternal(; plasma_response=P) file_ctrl = _GPE_TR.ForcingTerms.ForcingTermsControl(; forcing_data_file="forcing.dat") - _PE_TR.compute_torque_response!(file_state, ffs, file_ctrl, file_intr, _PE_TR.PerturbedEquilibriumControl()) + _PE_TR.apply_torque_response!(file_state, file_ctrl) ts = file_state.torque_response @test ts.source_names == ["forcing.dat"] @test size(ts.T_source) == (1, 1, npsi) @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 - flux_conform = pe.rootarea_to_area_weight .* pe.surface_area - @test_throws DimensionMismatch _PE_TR.torque_response_matrices(ffs, P, flux_conform; source_flux=zeros(ComplexF64, N, 1)) @test_throws DimensionMismatch _PE_TR.torque_profile(tr, ComplexF64[1.0]) - empty_state = _PE_TR.PerturbedEquilibriumState() - empty_intr = _PE_TR.PerturbedEquilibriumInternal() - @test_throws ArgumentError _PE_TR.compute_torque_response!(empty_state, ffs, ft_ctrl, empty_intr, _PE_TR.PerturbedEquilibriumControl()) + # A response without its operators cannot carry a torque response. + @test_throws ArgumentError _PE_TR.build_torque_response!(_PE_TR.PlasmaResponse(; ffs)) + # The torque response is opt-in: a response built without it contracts nothing. + plain = perturbed_equilibrium(plasma_response(ffs; compute_singular_coupling=false, verbose=false), _GPE_TR.ForcingTerms.RMPField(ft_ctrl); + coil_sets=pe.coil_sets, verbose=false, write_outputs_to_HDF5=false) + @test isempty(plain.response.T_xe) && plain.torque_response === nothing end end end