From e06487dfddadd1be8eafaa8fb7d488659166b223 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Wed, 12 Aug 2026 15:58:45 -0400 Subject: [PATCH 1/3] VAC - REFACTOR - Give the Vacuum module its own output struct VacuumData lived in ForceFreeStates and was passed into Vacuum as an untyped output buffer, so the Vacuum module could not name the type it wrote to. Its 14 fields mixed vacuum outputs, free-boundary energies computed in Free.jl, and sizing scratch. Split into two focused immutable structs: - Vacuum.VacuumResponse (src/Vacuum/DataTypes.jl) holds what the Vacuum module computes: wv, grri, grre, plasma_pts, wall_pts, mtheta, nzeta. compute_vacuum_response returns it instead of a positional 5-tuple that every caller partly discarded, and compute_vacuum_response! is typed rather than duck-typed on field names. - ForceFreeStates.FreeBoundaryResult holds the energy decomposition, built once at the end of free_run instead of being mutated in place across 40 lines. grri/grre are no longer retained past free_run: nothing read them there, and downstream consumers recompute their own. Drops ~2 MB of live allocation per run at mthvac=512. The dead numpoints/numpert_total/mthvac fields are gone too. Consumers narrowed to what they actually read (issue #139 principle 1): galerkin_solve takes wv, compute_perturbed_equilibrium takes wt0 and mthvac and no longer imports a ForceFreeStates struct at all, and build_flux_matrix loses a parameter it never used. No intended numerical change. Co-Authored-By: Claude Opus 5 --- benchmarks/benchmark_threads.jl | 2 +- .../pfac_study/gal_pfac_worker.jl | 2 +- .../tokamaker_beta/gal_tok_worker.jl | 2 +- benchmarks/plot_xi_eigenmode.jl | 2 +- docs/src/stability.md | 4 +- docs/src/vacuum.md | 4 +- .../run_imas_example.jl | 4 +- src/ForceFreeStates/EulerLagrange.jl | 2 +- src/ForceFreeStates/ForceFreeStatesStructs.jl | 61 ++++----- src/ForceFreeStates/Free.jl | 119 +++++++++--------- src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 14 +-- src/ForceFreeStates/Riccati.jl | 7 +- src/GeneralizedPerturbedEquilibrium.jl | 63 +++++----- .../PerturbedEquilibrium.jl | 18 +-- src/PerturbedEquilibrium/Response.jl | 17 +-- src/PerturbedEquilibrium/ResponseMatrices.jl | 17 ++- src/PerturbedEquilibrium/SingularCoupling.jl | 12 +- src/Vacuum/DataTypes.jl | 42 +++++++ src/Vacuum/Vacuum.jl | 63 +++------- test/runtests_imas.jl | 10 +- test/runtests_parallel_integration.jl | 6 +- test/runtests_riccati.jl | 6 +- test/runtests_vacuum.jl | 38 +++--- 23 files changed, 260 insertions(+), 255 deletions(-) diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 4831f6988..048c64c66 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -32,7 +32,7 @@ function run_ffs(ex; use_parallel, use_riccati=false) metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) return real(vac.et[1]), intr.numpert_total end diff --git a/benchmarks/gal_validation/pfac_study/gal_pfac_worker.jl b/benchmarks/gal_validation/pfac_study/gal_pfac_worker.jl index f0a344875..2aa197680 100644 --- a/benchmarks/gal_validation/pfac_study/gal_pfac_worker.jl +++ b/benchmarks/gal_validation/pfac_study/gal_pfac_worker.jl @@ -37,7 +37,7 @@ intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) ffit = FFS.make_matrix(equil, intr, metric) odet, fm_propagators, fm_chunks, fm_S_left = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) -vac_data = FFS.free_run!(odet, ctrl, equil, ffit, intr) +vac_data = FFS.free_run(odet, ctrl, equil, ffit, intr) if intr.msing > 0 && fm_propagators !== nothing FFS.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; wv=vac_data.wv, psio=equil.psio, debug=false, diff --git a/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl b/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl index f461601e8..2c72d9d33 100644 --- a/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl +++ b/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl @@ -73,7 +73,7 @@ function run_pipeline(psihigh) metric = FFS.make_metric(equil, intr.mpert) ffit = FFS.make_matrix(equil, intr, metric) odet, fm_propagators, fm_chunks, fm_S_left = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac_data = FFS.free_run!(odet, ctrl, equil, ffit, intr) + vac_data = FFS.free_run(odet, ctrl, equil, ffit, intr) if intr.msing > 0 && fm_propagators !== nothing FFS.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; wv=vac_data.wv, psio=equil.psio, debug=false, diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl index e7c675207..f39682caa 100644 --- a/benchmarks/plot_xi_eigenmode.jl +++ b/benchmarks/plot_xi_eigenmode.jl @@ -1,6 +1,6 @@ # Plot the ξ(ψ) profile of the eigenmode with the highest total-energy eigenvalue. # -# The total energy operator W = W_plasma + W_vacuum (free_run!, Free.jl); its eigenvectors are the +# The total energy operator W = W_plasma + W_vacuum (free_run, Free.jl); its eigenvectors are the # free-boundary edge displacement patterns and the eigenvalues are δW. This picks the eigenmode with the # largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix # (integration/xi_psi) onto that edge eigenvector: c = U_edge \ w, ξ(ψ) = U(ψ)·c. diff --git a/docs/src/stability.md b/docs/src/stability.md index 9e2d7f8d4..f5fb4b0c1 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -311,7 +311,7 @@ ffit = FFS.make_matrix(equil, intr, metric) # and always returns a 4-tuple (odet, propagators, chunks, S_at_surface_left). odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) -vac = FFS.free_run!(odet, ctrl, equil, ffit, intr) +vac = FFS.free_run(odet, ctrl, equil, ffit, intr) println("Energy eigenvalue et[1] = ", real(vac.et[1])) ``` @@ -356,5 +356,5 @@ end - `docs/src/galerkin.md` — RDCON outer-region Galerkin Δ′ solver (part of this module) - `docs/src/equilibrium.md` — build the `PlasmaEquilibrium` object required by this module -- `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run!` +- `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run` - `docs/src/perturbed_equilibrium.md` — downstream singular coupling analysis using Δ' diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 46ca5f44d..40b7c08fb 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -73,8 +73,8 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings( equal_arc_wall = true # Use equal arc length spacing ) -# Compute vacuum response matrix -wv, grri, xzpts = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) +# Compute vacuum response; returns a VacuumResponse with wv, grri, grre, plasma_pts, wall_pts +vac = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) ``` ### Vacuum Field Calculation at Observation Points diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl index df9641531..f7f6cfb9e 100644 --- a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl +++ b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl @@ -61,7 +61,7 @@ try TOML.print(io, config_imas); end result_imas = GPEC.main([tmpdir_imas]; dd=dd) - global et_imas = real(result_imas.vac_data.et[1]) + global et_imas = real(result_imas.free_energies.et[1]) global mpert_imas = result_imas.intr.mpert GPEC.write_imas(dd, result_imas) @assert dd.mhd_linear.time_slice[1].toroidal_mode[1].energy_perturbed ≈ et_imas @@ -80,7 +80,7 @@ try TOML.print(io, config_gfile); end result_gfile = GPEC.main([tmpdir_gfile]) - global et_gfile = real(result_gfile.vac_data.et[1]) + global et_gfile = real(result_gfile.free_energies.et[1]) global mpert_gfile = result_gfile.intr.mpert finally rm(tmpdir_gfile; recursive=true) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 631fea55a..1ae442098 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -252,7 +252,7 @@ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::E end odet.nzero = evaluate_stability_criterion!(odet, equil.profiles) - # Undo Gaussian reduction to get true solution vectors (for free_run! eigenvector use) + # Undo Gaussian reduction to get true solution vectors transform_u!(odet, intr) return (odet, nothing, nothing, nothing) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 5c854c5d4..046d4474c 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -248,7 +248,7 @@ gpec.toml. - `use_riccati::Bool` - Use the dual Riccati reformulation S = U₁·U₂⁻¹ instead of the standard U₁/U₂ ODE. Reduces stiffness for faster integration. See Glasser (2018) Phys. Plasmas 25, 032507. - `use_parallel::Bool` - Parallel fundamental matrix (propagator) integration using `Threads.@threads`. Each chunk is integrated independently from identity IC and assembled serially. Requires `singfac_min != 0`. Uses the same chunk bounds as the standard path but sub-divides chunks for load balancing. Crossings use the Riccati-style algorithm (no Gaussian reduction). - `parallel_threads::Int` - Cap on the number of threads the parallel BVP uses. **Default `2`** parallelises the FM chunks across two threads (the BVP has ~10 chunks; 2 threads is enough to amortize them — speedup saturates here, raising to 4 adds scheduling overhead). Set `parallel_threads = 1` to run the FM chunks SERIALLY (no `Threads.@threads`), which is bit-deterministic and immune to the thread-schedule sensitivity that can cause intermittent BVP divergence on numerically delicate equilibria. The parallel path produces bit-identical Δ′ across thread counts; `parallel_threads = 2` is about 20% faster than serial and saturates the speedup. If a parallel run diverges, drop to `parallel_threads = 1` rather than switching `use_parallel = false` — the latter is silently wrong. Capped at `Threads.nthreads()`. - - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline. This populates `u_store` / `ud_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly. Without it the parallel path stores only chunk-endpoint Riccati S matrices and zeros for `ud_store` (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable. Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`. Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run!` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases). **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage). Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case. Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N). + - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline. This populates `u_store` / `ud_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly. Without it the parallel path stores only chunk-endpoint Riccati S matrices and zeros for `ud_store` (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable. Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`. Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases). **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage). Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case. Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N). - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64. """ @kwdef struct ForceFreeStatesControl @@ -394,52 +394,39 @@ end FourFitVars(mpert::Int, numpert_total::Int) = FourFitVars(; mpert, numpert_total) """ - VacuumData + FreeBoundaryResult -A struct containing relevant data from the vacuum calculation. -Populated in `Free.jl`. +Result of the free-boundary calculation, returned by `free_run`. All matrices are in the ξ Fourier +basis and are `numpert_total × numpert_total`; the energies are generalized (W, N) pencil values, +power-normalized and invariant to the working (Jacobian) coordinate. ## Fields - - `numpoints::Int` - Total number of points in the vacuum calculation (mthvac * nzvac) - - `numpert_total::Int` - Total number of modes (mpert × npert) - - `mthvac::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput) - only needed for GPEC functionality currently - - `wt::Array{ComplexF64, 2}` - Free-boundary eigenvector matrix of the generalized eigenproblem W·v = λ·N·v (numpert_total × numpert_total). Columns are the eigenmodes sorted most-unstable first, normalized to unit power norm v†·N·v = 1. - - `wt0::Array{ComplexF64, 2}` - Free-boundary total-energy matrix W = wp + wv before diagonalisation (numpert_total × numpert_total). ξ Fourier basis. - - `wp::Array{ComplexF64, 2}` - Plasma energy matrix (numpert_total × numpert_total). ξ Fourier basis. - - `wv::Array{ComplexF64, 2}` - Vacuum energy matrix (numpert_total × numpert_total). ξ Fourier basis. + - `wt::Matrix{ComplexF64}` - Eigenvector matrix of W·v = λ·N·v. Columns are eigenmodes sorted most-unstable first, normalized to unit power norm v†·N·v = 1. + - `wt0::Matrix{ComplexF64}` - Total-energy matrix W = wp + wv before diagonalisation + - `wp::Matrix{ComplexF64}` - Plasma energy matrix + - `wv::Matrix{ComplexF64}` - Vacuum energy matrix, singfac-scaled at `qlim` - `ep::Vector{ComplexF64}` - Plasma energy per eigenmode (power quotient v†·wp·v with v†·N·v = 1) - `ev::Vector{ComplexF64}` - Vacuum energy per eigenmode (power quotient v†·wv·v with v†·N·v = 1) - - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N): power-normalized and invariant to the working (Jacobian) coordinate; et = ep + ev per mode - - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas` + - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N); et = ep + ev per mode + - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total) - - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total) - - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z) - - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z) + - `plasma_pts`, `wall_pts::Matrix{Float64}` - Cartesian (x, y, z) surface coordinates, `numpoints × 3`, retained for HDF5 output """ -@kwdef mutable struct VacuumData - numpoints::Int - numpert_total::Int - mthvac::Int # this is only needed to not break GPEC functionality currently - - wt::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total) - wt0::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total) - wp::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total) - wv::Array{ComplexF64,2} = Array{ComplexF64}(undef, numpert_total, numpert_total) - ep::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) - ev::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) - et::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) - n_tor_idx::Vector{Int} = zeros(Int, numpert_total) - vacuum_eigenvalue::Float64 = NaN - grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) - grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) - plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) - wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) +struct FreeBoundaryResult + wt::Matrix{ComplexF64} + wt0::Matrix{ComplexF64} + wp::Matrix{ComplexF64} + wv::Matrix{ComplexF64} + ep::Vector{ComplexF64} + ev::Vector{ComplexF64} + et::Vector{ComplexF64} + n_tor_idx::Vector{Int} + vacuum_eigenvalue::Float64 + plasma_pts::Matrix{Float64} + wall_pts::Matrix{Float64} end -VacuumData(numpoints::Int, numpert_total::Int, mthvac::Int) = VacuumData(; numpoints, numpert_total, mthvac) - """ EdgeScanState diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 52acb9ba2..96f3c37c0 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -9,7 +9,7 @@ Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.jmats` spline is the flux-surface average of the squared boundary displacement — the DCON power normalization as a quadratic form. N is Hermitian Toeplitz within each n-block (block-diagonal over n since the Jacobian is axisymmetric) and positive definite (J > 0). It is the metric of the generalized -eigenproblem W·v = λ·N·v solved in `free_run!` and `free_compute_total`: because W and N +eigenproblem W·v = λ·N·v solved in `free_run` and `free_compute_total`: because W and N transform by the same congruence under a change of working (Jacobian) coordinate, the pencil eigenvalues are power-normalized mode energies that are invariant to that coordinate choice. """ @@ -26,21 +26,41 @@ function power_norm_matrix!(Nmat::AbstractMatrix{ComplexF64}, jmat::AbstractVect end """ - free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> VacuumData + normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64) -> OdeState + +Rescale the stored EL solution vectors in `odet.u_store` / `odet.ud_store` so the edge +displacement matches the free-boundary eigenvectors `wt` (scaled by `2π·psio·1e-3`). +Modifies `odet` in place. Call after `free_run` when downstream code consumes the stored ξ profiles. +""" +@with_pool pool function normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64) + N = size(wt, 1) + tmp_mat = zeros!(pool, ComplexF64, N, N) + coeffs = odet.u[:, :, 1, end] \ (wt .* (2π * psio * 1e-3)) + @views for istep in 1:odet.step + mul!(tmp_mat, odet.u_store[:, :, 1, istep], coeffs) + odet.u_store[:, :, 1, istep] .= tmp_mat + mul!(tmp_mat, odet.u_store[:, :, 2, istep], coeffs) + odet.u_store[:, :, 2, istep] .= tmp_mat + mul!(tmp_mat, odet.ud_store[:, :, 1, istep], coeffs) + odet.ud_store[:, :, 1, istep] .= tmp_mat + mul!(tmp_mat, odet.ud_store[:, :, 2, istep], coeffs) + odet.ud_store[:, :, 2, istep] .= tmp_mat + end +end + +""" + free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> FreeBoundaryResult Compute the free boundary energies using the Julia port of the VACUUM code. Performs the same function as `free_run` -in the Fortran code, except now all data is passed in memory instead of via files. This -modifies `odet` in place to normalize the eigenfunctions stored in `u_store` and `ud_store`, -and returns a `VacuumData` struct containing the data needed for perturbed equilibrium calculations -and data dumping. +in the Fortran code. + +Returns a `FreeBoundaryResult` struct containing the data needed for perturbed equilibrium +calculations and data dumping. """ -@with_pool pool function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) +@with_pool pool function free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) # Initializations and allocations (; mpert, mlow, mhigh, numpert_total, psilim, qlim, npert, nlow, nhigh, wall_settings) = intr - vac_data = VacuumData(ctrl.mthvac * ctrl.nzvac, intr.numpert_total, ctrl.mthvac) - etemp = zeros!(pool, ComplexF64, numpert_total) - wp = zeros!(pool, ComplexF64, numpert_total, numpert_total) wpt = zeros!(pool, ComplexF64, numpert_total, numpert_total) wvt = zeros!(pool, ComplexF64, numpert_total, numpert_total) tmp_mat = zeros!(pool, ComplexF64, numpert_total, numpert_total) @@ -49,15 +69,17 @@ and data dumping. dV_dpsi = equil.profiles.dVdpsi_spline(psilim) # Compute plasma response matrix W = U₂ * U₁⁻¹ + wp = zeros(ComplexF64, numpert_total, numpert_total) @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 - # Compute vacuum response matrix in-place (handles 2D single-n, 2D multi-n block-diagonal, and 3D) + # Compute vacuum response (handles 2D single-n, 2D multi-n block-diagonal, and 3D) vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mlow:mhigh, nlow:nhigh) - Vacuum.compute_vacuum_response!(vac_data, vac_inputs, wall_settings) + vac = Vacuum.compute_vacuum_response(vac_inputs, wall_settings) + wv = vac.wv # Scale by (m - n*q)(m' - n'*q) [Chance Phys. Plasmas 1997 2161 eq. 126] singfac = vec((mlow:mhigh) .- qlim .* (nlow:nhigh)') - vac_data.wv .*= singfac .* singfac' + wv .*= singfac .* singfac' # Power-normalization matrix N at the plasma edge: ξ†·N·ξ = ⟨|ξ|²⟩ (see power_norm_matrix!). # The Jacobian band is evaluated at psilim (same surface as W), not at the last grid surface. @@ -68,79 +90,64 @@ and data dumping. # Least stable eigenvalue of the vacuum matrix alone, power-normalized via the pencil # (wv, N) so it shares the units of the mode energies (should be PSD; clamp noise to zero) - vac_data.vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(vac_data.wv), Hermitian(Nmat))))) - - # Preserve wp in vac_data so it can be written to HDF5 as W_plasma - vac_data.wp .= wp + vacuum_eigenvalue = max(0.0, minimum(real.(eigvals(Hermitian(wv), Hermitian(Nmat))))) # Complex energy eigenvalues and vectors of the generalized eigenproblem W·v = λ·N·v. # The eigenvalues are stationary values of the power quotient ξ†Wξ/ξ†Nξ — power-normalized # mode energies invariant to the working (Jacobian) coordinate, since W and N transform by # the same congruence under a poloidal coordinate change. - vac_data.wt .= wp .+ vac_data.wv - vac_data.wt0 .= vac_data.wt - Ev = eigen(vac_data.wt, Nmat) - vac_data.et .= Ev.values - eindex = sortperm(real.(vac_data.et); rev=true) + wt = wp .+ wv + wt0 = copy(wt) + Ev = eigen(wt, Nmat) + et = Ev.values + eindex = sortperm(real.(et); rev=true) - etemp .= vac_data.et # Rearrange wt columns for ascending real eigenvalues (most unstable first) + etemp = et + n_tor_idx = zeros(Int, numpert_total) for ipert in 1:numpert_total orig = eindex[numpert_total+1-ipert] - vac_data.wt[:, ipert] .= Ev.vectors[:, orig] - vac_data.et[ipert] = etemp[orig] + wt[:, ipert] .= Ev.vectors[:, orig] + et[ipert] = etemp[orig] # Store which n this eigenvector corresponds to (needed to write IMAS data) # This relies on the block diagonal matrix structure due to n decoupling in tokamaks imax = argmax(abs.(Ev.vectors[:, orig])) - vac_data.n_tor_idx[ipert] = (imax - 1) ÷ mpert + n_tor_idx[ipert] = (imax - 1) ÷ mpert end # Normalize eigenvectors to unit power norm v†·N·v = 1. The generalized eigenvalues are # already the power-normalized energies, so et is not rescaled here. for isol in 1:numpert_total - v = @view vac_data.wt[:, isol] + v = @view wt[:, isol] v ./= sqrt(real(dot(v, Nmat, v))) end # Normalize phase imax = 0 for isol in 1:numpert_total - imax = argmax(abs.(vac_data.wt[:, isol])) - phase = abs(vac_data.wt[imax, isol]) / vac_data.wt[imax, isol] - vac_data.wt[:, isol] .*= phase + imax = argmax(abs.(wt[:, isol])) + phase = abs(wt[imax, isol]) / wt[imax, isol] + wt[:, isol] .*= phase end # Project W_p and W_v into the eigenmode basis. Diagonal entries give # the plasma/vacuum energy split for each mode: et[i] = ep[i] + ev[i] - mul!(tmp_mat, wp, vac_data.wt) - mul!(wpt, vac_data.wt', tmp_mat) - mul!(tmp_mat, vac_data.wv, vac_data.wt) - mul!(wvt, vac_data.wt', tmp_mat) - vac_data.ep .= diag(wpt) - vac_data.ev .= diag(wvt) - - # Normalize eigenvectors based on scaled wt - coeffs = odet.u[:, :, 1, end] \ (vac_data.wt .* (2π * equil.psio * 1e-3)) - @views for istep in 1:odet.step - mul!(tmp_mat, odet.u_store[:, :, 1, istep], coeffs) - odet.u_store[:, :, 1, istep] .= tmp_mat - mul!(tmp_mat, odet.u_store[:, :, 2, istep], coeffs) - odet.u_store[:, :, 2, istep] .= tmp_mat - mul!(tmp_mat, odet.ud_store[:, :, 1, istep], coeffs) - odet.ud_store[:, :, 1, istep] .= tmp_mat - mul!(tmp_mat, odet.ud_store[:, :, 2, istep], coeffs) - odet.ud_store[:, :, 2, istep] .= tmp_mat - end + mul!(tmp_mat, wp, wt) + mul!(wpt, wt', tmp_mat) + mul!(tmp_mat, wv, wt) + mul!(wvt, wt', tmp_mat) + ep = diag(wpt) + ev = diag(wvt) # Write energies to screen if ctrl.verbose @info "Least Stable Eigenmode Energies:\n" * - " Plasma = $((@sprintf "%+.3e %+.3ei" real(vac_data.ep[1]) imag(vac_data.ep[1])))\n" * - " Vacuum = $((@sprintf "%+.3e %+.3ei" real(vac_data.ev[1]) imag(vac_data.ev[1])))\n" * - " Total = $((@sprintf "%+.3e %+.3ei" real(vac_data.et[1]) imag(vac_data.et[1])))" + " Plasma = $((@sprintf "%+.3e %+.3ei" real(ep[1]) imag(ep[1])))\n" * + " Vacuum = $((@sprintf "%+.3e %+.3ei" real(ev[1]) imag(ev[1])))\n" * + " Total = $((@sprintf "%+.3e %+.3ei" real(et[1]) imag(et[1])))" end - return vac_data + return FreeBoundaryResult(wt, wt0, wp, wv, ep, ev, et, n_tor_idx, vacuum_eigenvalue, vac.plasma_pts, vac.wall_pts) end """ @@ -175,8 +182,8 @@ q-window minimum. # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically) vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh) - wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) - @views wv_array[i, :, :] .= wv + vac = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) + @views wv_array[i, :, :] .= vac.wv end # Flatten 3D array to (npsi+1 × numpert_total^2) for series interpolant @@ -232,7 +239,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma power_norm_matrix!(Nmat, jmat_local, intr.mpert, intr.npert, dV_dpsi) # Total energy matrix and generalized eigen-decomposition of the pencil (W, N) — the - # eigenvalues are power-normalized, Jacobian-invariant mode energies (see free_run!) + # eigenvalues are power-normalized, Jacobian-invariant mode energies wt .= wp .+ wv Ev = eigen(wt, Nmat) diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index 0b7dd658e..7d534cd3f 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -48,17 +48,17 @@ end """ galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, - intr::ForceFreeStatesInternal; vac_data=nothing) -> GalerkinResult + intr::ForceFreeStatesInternal; wv=nothing) -> GalerkinResult Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of `gal_solve` (gal.f). Single toroidal mode only (`intr.npert == 1`). Returns a `GalerkinResult`; if there are no resonant surfaces in the domain it returns an empty result. -`vac_data` (a `VacuumData` from `free_run!`) supplies the free-boundary edge term -`wv_edge = vac_data.wv · psio²`; pass `nothing` for a fixed-boundary edge. +`wv` is the vacuum energy matrix from `free_run`, which supplies the free-boundary edge term +`wv_edge = wv · psio²`; pass `nothing` for a fixed-boundary edge. """ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, - intr::ForceFreeStatesInternal; vac_data=nothing) + intr::ForceFreeStatesInternal; wv=nothing) intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported") ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"") @@ -122,12 +122,12 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, ws.rhs = zeros(ComplexF64, ws.ndim, nsol) ws.sol = zeros(ComplexF64, ws.ndim, nsol) - # Free-boundary edge term wvac·psio² (vac_data.wv is already singfac-scaled at qlim in free_run!). + # Free-boundary edge term wvac·psio² (wv is already singfac-scaled at qlim in free_run). # rpec passes wv_edge=nothing; gal_set_boundary! then applies the identity edge AND injects the coil # unit sources (its three-way branch). The vacuum block is only built for the non-rpec free case. wv_edge = nothing - if ctrl.vac_flag && vac_data !== nothing && ncoil == 0 - wv_edge = Matrix{ComplexF64}(vac_data.wv .* equil.psio^2) + if ctrl.vac_flag && wv !== nothing && ncoil == 0 + wv_edge = Matrix{ComplexF64}(wv .* equil.psio^2) end gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge) diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index cb7cfeee8..d266d4dd2 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -1648,7 +1648,7 @@ function parallel_eulerlagrange_integration( chunks, propagators = _handle_edge_dW_scan!(odet, chunks, propagators, ctrl, equil, ffit, intr) - # compute_delta_prime_matrix! is called from the main pipeline (after free_run!) so + # compute_delta_prime_matrix! is called from the main pipeline (after free_run) so # that vacuum response wv is available for the edge BC. With self-consistent truncation, # the propagators/chunks returned here match intr.psilim exactly, so Δ' is well-defined # for both truncate_at_dW_peak=false (full domain) and =true (peak). @@ -1883,8 +1883,9 @@ end Replace the propagator-BVP's `odet` with a fresh serial-EL `odet` that has dense `u_store` / `ud_store` populated in axis basis (the PerturbedEquilibrium convention). The caller's `odet` is fully replaced by the fresh one because -`free_run!` downstream uses `odet.u[:,:,1,end]` to normalize `odet.u_store`, -so both must be in the same basis. The parallel BVP results that survive +`free_run` / `normalize_eigenfunctions!` downstream use `odet.u[:,:,1,end]` to +normalize `odet.u_store`, so both must be in the same basis. The parallel BVP +results that survive downstream are stored in `intr` (psilim/qlim, sing[*].delta_prime, …) and in the externally-returned `propagators` / `chunks` / `S_at_surface_left` — none of those live on `odet`, so replacing `odet` is safe. diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 8c45231dd..7c1a158d4 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -62,12 +62,12 @@ export Analysis include("Rerun.jl") # Import ForceFreeStates types and functions needed for main -using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, VacuumData, OdeState, FourFitVars +using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, FreeBoundaryResult, OdeState, FourFitVars using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix using .ForceFreeStates: find_kinetic_singular_surfaces! -using .ForceFreeStates: eulerlagrange_integration, free_run! +using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag") @@ -449,8 +449,9 @@ function main_from_inputs( wall_desc = intr.wall_settings.shape == "nowall" ? "no wall" : intr.wall_settings.shape @info "Computing free boundary energies ($wall_desc)" end - vac_data = free_run!(odet, ctrl, equil, ffit, intr) - if real(vac_data.et[1]) < 0 + free_energies = free_run(odet, ctrl, equil, ffit, intr) + normalize_eigenfunctions!(odet, free_energies.wt, equil.psio) + if real(free_energies.et[1]) < 0 if ctrl.verbose @warn "Free-boundary mode unstable for n = $nstring" end @@ -461,13 +462,13 @@ function main_from_inputs( end # Compute inter-surface Δ' matrix (STRIDE BVP) using vacuum edge BC. - # Requires propagators from parallel FM path and wv from free_run!. + # Requires propagators from parallel FM path and wv from free_run. if ctrl.kinetic_factor == 0 && intr.msing > 0 && fm_propagators !== nothing if ctrl.verbose @info "Computing Δ' matrix (STRIDE BVP with vacuum coupling)" end ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; - wv=vac_data.wv, psio=equil.psio, debug=ctrl.verbose, + wv=free_energies.wv, psio=equil.psio, debug=ctrl.verbose, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) end @@ -477,7 +478,7 @@ function main_from_inputs( gal_data = nothing if ctrl.gal_flag gal_start = time() - gal_data = galerkin_solve(ctrl, equil, ffit, intr; vac_data=ctrl.vac_flag ? vac_data : nothing) + gal_data = galerkin_solve(ctrl, equil, ffit, intr; wv=ctrl.vac_flag ? free_energies.wv : nothing) @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s" end @@ -487,7 +488,7 @@ function main_from_inputs( equil, intr, odet, - ctrl.vac_flag ? vac_data : nothing, + ctrl.vac_flag ? free_energies : nothing, ffit, git_version, inputs, @@ -540,7 +541,7 @@ function main_from_inputs( slayer_result = _run_slayer_stage(nothing) @info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER" return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, - vac_data=ctrl.vac_flag ? vac_data : nothing, + free_energies=ctrl.vac_flag ? free_energies : nothing, slayer=slayer_result) end @@ -597,9 +598,9 @@ function main_from_inputs( end # Run perturbed equilibrium calculations - # Pass vac_data and intr for response matrix calculations + # Free-boundary wt0 drives the plasma inductance; mthvac sizes the Green's-function solves pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium( - equil, pe_odet, ctrl.vac_flag ? vac_data : nothing, intr, ft_ctrl, pe_ctrl, pe_intr, + equil, pe_odet, ctrl.vac_flag ? free_energies.wt0 : nothing, ctrl.mthvac, intr, ft_ctrl, pe_ctrl, pe_intr, metric, ffit ) @@ -670,7 +671,7 @@ function main_from_inputs( # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, - vac_data=ctrl.vac_flag ? vac_data : nothing, + free_energies=ctrl.vac_flag ? free_energies : nothing, slayer=slayer_result) end @@ -693,7 +694,7 @@ function write_outputs_to_HDF5( equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, - vac_data::Union{VacuumData,Nothing}, + free_energies::Union{FreeBoundaryResult,Nothing}, ffit::Union{FourFitVars,Nothing}=nothing, git_version::String="unknown", inputs::Union{Nothing,Dict{String,Any}}=nothing, @@ -891,23 +892,23 @@ function write_outputs_to_HDF5( # working (Jacobian) coordinate. W_freeboundary_eigenmodes holds the generalized # eigenvectors, columns sorted most-unstable first, normalized to unit power norm with # the largest-magnitude entry made real-positive. - out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? vac_data.wt0 : ComplexF64[] - out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? vac_data.wp : ComplexF64[] - out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? vac_data.wv : ComplexF64[] - out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.wt : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? vac_data.et : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.ep : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.ev : ComplexF64[] - out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = ctrl.vac_flag ? vac_data.vacuum_eigenvalue : NaN + out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? free_energies.wt0 : ComplexF64[] + out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? free_energies.wp : ComplexF64[] + out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? free_energies.wv : ComplexF64[] + out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? free_energies.wt : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? free_energies.et : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? free_energies.ep : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? free_energies.ev : ComplexF64[] + out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = ctrl.vac_flag ? free_energies.vacuum_eigenvalue : NaN # Cartesian surface point clouds used downstream for visualisation and # perturbed-equilibrium plotting. - out_h5["SurfaceGeometries/Plasma/x"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 1] : Float64[] - out_h5["SurfaceGeometries/Plasma/y"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 2] : Float64[] - out_h5["SurfaceGeometries/Plasma/z"] = ctrl.vac_flag ? vac_data.plasma_pts[:, 3] : Float64[] - out_h5["SurfaceGeometries/Wall/x"] = ctrl.vac_flag ? vac_data.wall_pts[:, 1] : Float64[] - out_h5["SurfaceGeometries/Wall/y"] = ctrl.vac_flag ? vac_data.wall_pts[:, 2] : Float64[] - out_h5["SurfaceGeometries/Wall/z"] = ctrl.vac_flag ? vac_data.wall_pts[:, 3] : Float64[] + out_h5["SurfaceGeometries/Plasma/x"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 1] : Float64[] + out_h5["SurfaceGeometries/Plasma/y"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 2] : Float64[] + out_h5["SurfaceGeometries/Plasma/z"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 3] : Float64[] + out_h5["SurfaceGeometries/Wall/x"] = ctrl.vac_flag ? free_energies.wall_pts[:, 1] : Float64[] + out_h5["SurfaceGeometries/Wall/y"] = ctrl.vac_flag ? free_energies.wall_pts[:, 2] : Float64[] + out_h5["SurfaceGeometries/Wall/z"] = ctrl.vac_flag ? free_energies.wall_pts[:, 3] : Float64[] # Write kinetic parameters when kinetic mode is enabled if ctrl.kinetic_factor > 0 @@ -997,9 +998,9 @@ receives the correct least-stable δW regardless of how modes are interleaved in The `result` argument is the named tuple returned by `main`. """ function write_imas(dd, result) - result.vac_data === nothing && return + result.free_energies === nothing && return - vac_data = result.vac_data + free_energies = result.free_energies intr = result.intr # Top-level metadata @@ -1014,10 +1015,10 @@ function write_imas(dd, result) # n_tor_idx[i] (0-based) identifies which n-block eigenvalue i belongs to. resize!(ts.toroidal_mode, intr.npert) for j in 0:(intr.npert-1) - n_indices = findall(==(j), vac_data.n_tor_idx) # indices of eigenvalues in the j-th n-block + n_indices = findall(==(j), free_energies.n_tor_idx) # indices of eigenvalues in the j-th n-block mode = ts.toroidal_mode[j+1] mode.n_tor = intr.nlow + j - mode.energy_perturbed = minimum(real.(vac_data.et[n_indices])) # least-stable energy for this n-toroidal mode + mode.energy_perturbed = minimum(real.(free_energies.et[n_indices])) # least-stable energy for this n-toroidal mode end return dd diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 5f1fe83b5..3775d8d97 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -11,7 +11,7 @@ using FastInterpolations # Import parent modules import ..Equilibrium import ..ForceFreeStates -import ..ForceFreeStates: OdeState, VacuumData, ForceFreeStatesInternal, FourFitVars, MetricData +import ..ForceFreeStates: OdeState, ForceFreeStatesInternal, FourFitVars, MetricData import ..Vacuum import ..ForcingTerms import ..ForcingTerms: ForcingMode, CoilSet, load_forcing_data!, convert_forcing_normalization! @@ -39,7 +39,7 @@ export write_outputs_to_HDF5 """ compute_perturbed_equilibrium( - equil, ForceFreeStates_results, vac_data, ffs_intr, + equil, ForceFreeStates_results, wt0, mthvac, ffs_intr, ft_ctrl, ctrl, intr, metric, ffit )::PerturbedEquilibriumState @@ -52,7 +52,8 @@ coupling metrics. - `equil`: Equilibrium solution from Equilibrium module - `ForceFreeStates_results`: Stability calculation results from ForceFreeStates module - - `vac_data`: Vacuum response data from ForceFreeStates free boundary calculation + - `wt0`: Free-boundary total-energy matrix W = wp + wv from `free_run`, or `nothing` when the free-boundary calculation was not run (response and singular coupling are then skipped) + - `mthvac`: Vacuum poloidal grid resolution to reuse for the Green's-function solves - `ffs_intr`: ForceFreeStates internal state with mode information - `ft_ctrl`: Forcing terms control parameters from [ForcingTerms] section - `ctrl`: Control parameters from [PerturbedEquilibrium] section @@ -67,7 +68,8 @@ coupling metrics. function compute_perturbed_equilibrium( equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::Union{VacuumData,Nothing}, + wt0::Union{Matrix{ComplexF64},Nothing}, + mthvac::Int, ffs_intr::ForceFreeStates.ForceFreeStatesInternal, ft_ctrl::ForcingTerms.ForcingTermsControl, ctrl::PerturbedEquilibriumControl, @@ -115,19 +117,19 @@ function compute_perturbed_equilibrium( # Step 2: Compute plasma response if ctrl.compute_response - if vac_data === nothing + if wt0 === nothing @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section." else - compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit) + compute_plasma_response!(state, equil, ForceFreeStates_results, wt0, mthvac, ffs_intr, intr, ctrl, metric, ffit) end end # Step 3: Compute singular coupling metrics if ctrl.compute_singular_coupling - if vac_data === nothing + if wt0 === nothing @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section." else - compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl) + compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, mthvac, ffs_intr, intr, ctrl) end end diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index ba2827246..e501bad29 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -1,6 +1,6 @@ """ compute_plasma_response!( - state, equil, ForceFreeStates_results, vac_data, ffs_intr, + state, equil, ForceFreeStates_results, wt0, mthvac, ffs_intr, intr, ctrl, metric, ffit ) @@ -18,7 +18,8 @@ function compute_plasma_response!( state::PerturbedEquilibriumState, equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::VacuumData, + wt0::Matrix{ComplexF64}, + mthvac::Int, ffs_intr::ForceFreeStatesInternal, intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl, @@ -30,19 +31,19 @@ function compute_plasma_response!( end # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode] - flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, vac_data, ffs_intr) + flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr) # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default) - plasma_inductance = calc_plasma_inductance(vac_data, ffs_intr, equil.psio) + plasma_inductance = calc_plasma_inductance(wt0, ffs_intr, equil.psio) # Surface inductance L from Green's functions at psilim. # Requires a 2D (nzvac=1) vacuum response so rows are theta points only, nn = ffs_intr.nlow - vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn]) + vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn]) wall_nowall = Vacuum.WallShapeSettings(; shape="nowall") - _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall) - grri_2d = Matrix{ComplexF64}(grri_2d_raw) - grre_2d = Matrix{ComplexF64}(grre_2d_raw) + vac_2d = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall) + grri_2d = Matrix{ComplexF64}(vac_2d.grri) + grre_2d = Matrix{ComplexF64}(vac_2d.grre) ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac) permeability = calc_permeability(plasma_inductance, surface_inductance) diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index a4b798a6d..791ed9c57 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -161,7 +161,6 @@ end build_flux_matrix( equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::VacuumData, intr::ForceFreeStatesInternal )::Matrix{ComplexF64} @@ -181,7 +180,6 @@ The flux matrix relates eigenmode displacements to vacuum poloidal flux: - `equil`: Equilibrium solution containing flux surfaces and q-profile - `ForceFreeStates_results`: ForceFreeStates ODE integration results containing eigenmodes - - `vac_data`: Vacuum response data from free boundary calculation - `intr`: ForceFreeStates internal state with mode information ## Returns @@ -192,7 +190,6 @@ The flux matrix relates eigenmode displacements to vacuum poloidal flux: function build_flux_matrix( equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::VacuumData, intr::ForceFreeStatesInternal )::Matrix{ComplexF64} @@ -209,7 +206,7 @@ end """ calc_plasma_inductance( - vac_data::VacuumData, + wt0::Matrix{ComplexF64}, ffs_intr::ForceFreeStatesInternal, psio::Float64 )::Matrix{ComplexF64} @@ -221,13 +218,13 @@ Calculate plasma inductance matrix Λ using the wt0-based energy formula where t₁ = im/(χ₁·s_i·2π), t₂ = -im/(χ₁·s_j·2π), s_i = m_i - n·q_lim. -Note: `vac_data.wt0` already contains singfac² factors (s_i·s_j) baked into the -vacuum term via the scaling in `free_run!`. The t₁/t₂ factors divide by s_i·s_j, +Note: `wt0` already contains singfac² factors (s_i·s_j) baked into the +vacuum term via the scaling in `free_run`. The t₁/t₂ factors divide by s_i·s_j, correctly recovering the properly-normalized inductance. ## Arguments - - `vac_data`: Vacuum data containing wt0 (total energy matrix before eigenvector sorting) + - `wt0`: Total energy matrix W = wp + wv, before eigenvector sorting - `ffs_intr`: ForceFreeStates internal state with mode info (mlow, mpert, nlow, qlim) - `psio`: Total toroidal flux [Wb/rad] from equilibrium (equil.psio) @@ -236,7 +233,7 @@ correctly recovering the properly-normalized inductance. - Plasma inductance matrix Lambda [numpert_total × numpert_total] """ function calc_plasma_inductance( - vac_data::VacuumData, + wt0::Matrix{ComplexF64}, ffs_intr::ForceFreeStatesInternal, psio::Float64 )::Matrix{ComplexF64} @@ -250,9 +247,9 @@ function calc_plasma_inductance( s = [((i-1) % ffs_intr.mpert + ffs_intr.mlow) - n * qlim for i in 1:mpert] # Fortran idcon_norm: wt0 = wt0/(mu0*2)*psio^2 - # Julia's vac_data.wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2) + # Julia's wt0 is raw wp+wv; Fortran additionally scales by psio^2/(mu0*2) mu0 = 4π * 1e-7 - wt0_norm = vac_data.wt0 .* (psio^2 / (mu0 * 2)) + wt0_norm = wt0 .* (psio^2 / (mu0 * 2)) # Build temp2[i,j] = 2·t1_i·wt0[i,j]·t2_j (matches Fortran gpresp_induct) temp2 = Matrix{ComplexF64}(undef, mpert, mpert) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index d9cd1ebc9..d219e9466 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -48,7 +48,7 @@ _reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) state::PerturbedEquilibriumState, equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::VacuumData, + mthvac::Int, ffs_intr::ForceFreeStatesInternal, intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl @@ -79,7 +79,7 @@ function compute_singular_coupling_metrics!( state::PerturbedEquilibriumState, equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::VacuumData, + mthvac::Int, ffs_intr::ForceFreeStatesInternal, intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl @@ -100,7 +100,7 @@ function compute_singular_coupling_metrics!( chi1 = 2π * equil.psio twopi = 2π - mtheta = vac_data.mthvac + mtheta = mthvac wall_settings = Vacuum.WallShapeSettings(; shape="nowall") # Phase 1: Collect all resonant (surface, n) pairs in psi order @@ -175,9 +175,9 @@ function compute_singular_coupling_metrics!( # Compute Green's functions at this surface for this n (once per pair) vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings) - grri = Matrix{ComplexF64}(grri_raw) - grre = Matrix{ComplexF64}(grre_raw) + vac = Vacuum.compute_vacuum_response(vac_input, wall_settings) + grri = Matrix{ComplexF64}(vac.grri) + grre = Matrix{ComplexF64}(vac.grre) # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation) ν_vac = Vacuum.PlasmaGeometry(vac_input).ν diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 06ee4a6a0..c1f0786e6 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -143,6 +143,48 @@ function expand_field_periods(inputs::VacuumInput) ) end +""" + VacuumResponse + +Output of [`compute_vacuum_response`](@ref): the vacuum energy matrix and the surface data the +boundary-integral solve produces along the way. + +## Fields + + - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D + - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on the 3D nowall path + - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`) + - `mtheta`, `nzeta::Int`: Grid resolution the response was computed on +""" +struct VacuumResponse + wv::Matrix{ComplexF64} + grri::Matrix{ComplexF64} + grre::Matrix{ComplexF64} + plasma_pts::Matrix{Float64} + wall_pts::Matrix{Float64} + mtheta::Int + nzeta::Int +end + +""" + VacuumResponse(inputs::VacuumInput) -> VacuumResponse + +Allocate zeroed output arrays sized for `inputs` (full torus, so `nfp` field periods). +""" +function VacuumResponse(inputs::VacuumInput) + num_points = inputs.mtheta * inputs.nzeta * inputs.nfp + num_modes = length(inputs.m_modes) * length(inputs.n_modes) + return VacuumResponse( + zeros(ComplexF64, num_modes, num_modes), + zeros(ComplexF64, 2 * num_points, num_modes), + zeros(ComplexF64, 2 * num_points, num_modes), + zeros(num_points, 3), + zeros(num_points, 3), + inputs.mtheta, + inputs.nzeta + ) +end + """ WallShapeSettings diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b060f59fc..f44dff518 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -18,7 +18,7 @@ include("Kernel2D.jl") include("Kernel3D.jl") include("Field.jl") -export VacuumInput, WallShapeSettings +export VacuumInput, VacuumResponse, WallShapeSettings export compute_vacuum_response, compute_vacuum_response!, compute_vacuum_field export extract_plasma_surface_at_psi export PlasmaGeometry @@ -53,7 +53,7 @@ function _symmetrize_vacuum_energy!(wv::AbstractMatrix) end """ - _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) 2D (axisymmetric) vacuum response calculation. @@ -62,7 +62,7 @@ building the double-/single-layer operators, solving the exterior and interior s filling the corresponding diagonal block of the response matrix and the matching column block of the Green's functions. """ -@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) mpert = length(inputs.m_modes) num_points_surf = inputs.mtheta @@ -147,7 +147,7 @@ of the Green's functions. end """ - _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) 3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1` the block-circulant assembly and residue-class loop are skipped in favour of a more efficient @@ -169,7 +169,7 @@ Only `wv` is produced currently; `grri`/`grre` are returned zeroed and are not y Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) @@ -257,57 +257,26 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) -> VacuumResponse -Allocate and return the vacuum response matrix and Green's functions for the given vacuum -inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the -output arrays for the full torus and forwards to the same 2D/3D workers. For performance-critical -paths that already own preallocated storage (e.g. `ForceFreeStates.VacuumData`), prefer the -in-place method to avoid extra heap allocations. - -# Returns - - - `wv`: complex vacuum response matrix (`num_modes × num_modes`). - - `grri`, `grre`: interior/exterior Green's functions (zeroed on the 3D nowall path). - - `plasma_pts`, `wall_pts`: surface coordinate arrays. +Compute the vacuum response for the given inputs. Allocating wrapper around +[`compute_vacuum_response!`](@ref); pass a preallocated [`VacuumResponse`](@ref) to that method +instead when reusing storage across calls. """ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) - - num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D - num_modes = length(inputs.m_modes) * length(inputs.n_modes) - - vac = ( - wv=zeros(ComplexF64, num_modes, num_modes), - grri=zeros(ComplexF64, 2 * num_points, num_modes), - grre=zeros(ComplexF64, 2 * num_points, num_modes), - plasma_pts=zeros(num_points, 3), - wall_pts=zeros(num_points, 3) - ) + vac = VacuumResponse(inputs) compute_vacuum_response!(vac, inputs, wall_settings) - - return vac.wv, vac.grri, vac.grre, vac.plasma_pts, vac.wall_pts + return vac end """ - compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) - -In-place variant that computes the vacuum response and directly populates the arrays stored in -`vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to -[`_compute_vacuum_response_2d!`], 3D to [`_compute_vacuum_response_3d!`]. - -The `vac_data` argument is expected to provide the following writable fields with compatible -sizes: - - - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - - `grri::AbstractMatrix{ComplexF64}` – interior Green's functions - - `grre::AbstractMatrix{ComplexF64}` – exterior Green's functions - - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates + compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) -This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete -type (duck-typed on field names only). +In-place variant that populates the arrays of an existing [`VacuumResponse`](@ref). Dispatches on +dimensionality only: 2D (`inputs.nzeta == 1`) routes to [`_compute_vacuum_response_2d!`], 3D to +[`_compute_vacuum_response_3d!`]. """ -function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) if inputs.nzeta == 1 _compute_vacuum_response_2d!(vac_data, inputs, wall_settings) else diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl index d2afd3091..774fdbfe8 100644 --- a/test/runtests_imas.jl +++ b/test/runtests_imas.jl @@ -115,7 +115,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium mock_et = [0.4+0im, 0.7+0im, 1.1+0im] mock_n_idx = [0, 0, 0] # all belong to n=1 (j=0) mock_result = ( - vac_data = (et=mock_et, n_tor_idx=mock_n_idx), + free_energies = (et=mock_et, n_tor_idx=mock_n_idx), intr = (numpert_total=3, npert=1, nlow=1), ) @@ -145,7 +145,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium mock_et = [0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im] mock_n_idx = [0, 1, 0, 1] mock_result = ( - vac_data = (et=mock_et, n_tor_idx=mock_n_idx), + free_energies = (et=mock_et, n_tor_idx=mock_n_idx), intr = (numpert_total=4, npert=2, nlow=1), ) @@ -169,7 +169,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium # Single n=1 run dd_single = IMASdd.dd() result_single = ( - vac_data = (et=[0.3+0im, 0.6+0im], n_tor_idx=[0, 0]), + free_energies = (et=[0.3+0im, 0.6+0im], n_tor_idx=[0, 0]), intr = (numpert_total=2, npert=1, nlow=1), ) GeneralizedPerturbedEquilibrium.write_imas(dd_single, result_single) @@ -177,7 +177,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium # Single n=2 run dd_single2 = IMASdd.dd() result_single2 = ( - vac_data = (et=[0.5+0im, 0.7+0im], n_tor_idx=[0, 0]), + free_energies = (et=[0.5+0im, 0.7+0im], n_tor_idx=[0, 0]), intr = (numpert_total=2, npert=1, nlow=2), ) GeneralizedPerturbedEquilibrium.write_imas(dd_single2, result_single2) @@ -185,7 +185,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium # Combined n=1,2 run (same eigenvalues, now interleaved) dd_multi = IMASdd.dd() result_multi = ( - vac_data = (et=[0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im], n_tor_idx=[0, 1, 0, 1]), + free_energies = (et=[0.3+0im, 0.5+0im, 0.6+0im, 0.7+0im], n_tor_idx=[0, 1, 0, 1]), intr = (numpert_total=4, npert=2, nlow=1), ) GeneralizedPerturbedEquilibrium.write_imas(dd_multi, result_multi) diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index ad95f6415..1a2bab276 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -254,7 +254,7 @@ using TOML metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) return real(vac.et[1]), intr end @@ -313,7 +313,7 @@ using TOML metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) return real(vac.et[1]), intr end @@ -521,7 +521,7 @@ using TOML ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet, fm_propagators, fm_chunks, fm_S_left = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_delta_prime_matrix!( intr, fm_propagators, fm_chunks; wv=vac.wv, psio=equil.psio, diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index e156a5be0..a08867499 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -118,14 +118,14 @@ end # (which overwrites intr_ric.sing[s].delta_prime) delta_prime_inline = [copy(intr_ric.sing[s].delta_prime) for s in 1:intr_ric.msing] - vac_ric = FFS.free_run!(odet_ric, ctrl, equil, ffit, intr_ric) + vac_ric = FFS.free_run(odet_ric, ctrl, equil, ffit, intr_ric) et_ric = real(vac_ric.et[1]) # Standard integration (needed only for energy comparison). eulerlagrange_integration # returns (odet, propagators, chunks, S_at_surface_left); only odet is used here. intr_std = make_solovev_intr(inputs, ctrl, equil, ex) odet_std, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr_std) - vac_std = FFS.free_run!(odet_std, ctrl, equil, ffit, intr_std) + vac_std = FFS.free_run(odet_std, ctrl, equil, ffit, intr_std) et_std = real(vac_std.et[1]) # ───────────────────────────────────────────────────────────────────────── @@ -215,7 +215,7 @@ end # result must be bit-for-bit identical (not just approximately equal). # # Note: this call overwrites intr_ric.sing[s].delta_prime; delta_prime_inline was - # saved before free_run! above so it holds the original inline values. + # saved before free_run above so it holds the original inline values. # # See benchmarks/benchmark_delta_prime_methods.jl for the extended version. FFS.compute_delta_prime_from_ca!(odet_ric, intr_ric, equil) diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 7e8b51b46..ff07f8be5 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -368,7 +368,7 @@ @testset "nowall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -391,7 +391,7 @@ @testset "conformal wall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="conformal", a=0.5) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -408,7 +408,7 @@ @testset "edge: single poloidal mode mpert=1" begin inputs = _make_inputs(m_modes=[1], n_modes=[1]) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) @test size(grri, 2) == 1 @@ -418,7 +418,7 @@ # Keep mtheta_eq=17 so boundary has enough points for periodic spline inputs = _make_inputs(mtheta=16, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) @test size(grri) == (32, 2) # 2*mtheta, num_modes=2 @test size(plasma_pts) == (16, 3) @@ -426,22 +426,19 @@ @testset "in-place compute_vacuum_response! matches wrapper" begin # The allocating wrapper is a thin caller of the in-place routine; verify the - # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically. + # in-place entry populates caller-owned storage identically. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - wv, grri, grre, pp, wp = compute_vacuum_response(inputs, wall_settings) + ref = compute_vacuum_response(inputs, wall_settings) - numpoints = inputs.mtheta * inputs.nzeta - num_modes = length(inputs.m_modes) * length(inputs.n_modes) - vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_modes), - plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) + vac = VacuumResponse(inputs) compute_vacuum_response!(vac, inputs, wall_settings) - @test vac.wv ≈ wv - @test vac.grri ≈ grri - @test vac.grre ≈ grre - @test vac.plasma_pts ≈ pp - @test vac.wall_pts ≈ wp + @test vac.wv ≈ ref.wv + @test vac.grri ≈ ref.grri + @test vac.grre ≈ ref.grre + @test vac.plasma_pts ≈ ref.plasma_pts + @test vac.wall_pts ≈ ref.wall_pts end end end @@ -581,7 +578,7 @@ @testset "compute_vacuum_response 3D nowall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -604,7 +601,7 @@ @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -625,7 +622,7 @@ @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="conformal", a=0.3) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + (; wv, grri, grre, plasma_pts, wall_pts) = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -676,13 +673,14 @@ wall_settings = WallShapeSettings(shape="nowall") # Reduced (block-circulant) path - wv_red, _, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) + vac_red = compute_vacuum_response(inputs_red, wall_settings) + wv_red, plasma_pts_red = vac_red.wv, vac_red.plasma_pts # Full-torus reference: pre-expand so nfp=1 forces the dense path inputs_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs_red) @test inputs_full.nfp == 1 @test inputs_full.nzeta == nzeta_p * nfp - wv_full, _, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) + wv_full = compute_vacuum_response(inputs_full, wall_settings).wv num_modes = length(m_modes) * length(n_modes) @test size(wv_red) == (num_modes, num_modes) From d60fd31f86201f27dc667917244b373b3f5c51a8 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 14 Aug 2026 12:28:05 -0400 Subject: [PATCH 2/3] GPEC - BUGFIX - Address review findings on the vacuum response struct refactor --- src/ForceFreeStates/Free.jl | 3 +- src/GeneralizedPerturbedEquilibrium.jl | 41 +++++++++++++------------- src/Vacuum/DataTypes.jl | 7 +---- src/Vacuum/Vacuum.jl | 4 ++- 4 files changed, 27 insertions(+), 28 deletions(-) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 9785a12e6..51efa7610 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -48,6 +48,7 @@ Modifies `odet` in place. Call after `free_run` when downstream code consumes th mul!(tmp_mat, odet.xi_s_store[:, :, istep], coeffs) odet.xi_s_store[:, :, istep] .= tmp_mat end + return odet end """ @@ -105,7 +106,7 @@ calculations and data dumping. eindex = sortperm(real.(et); rev=true) # Rearrange wt columns for ascending real eigenvalues (most unstable first) - etemp = et + etemp = copy(et) n_tor_idx = zeros(Int, numpert_total) for ipert in 1:numpert_total orig = eindex[numpert_total+1-ipert] diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index dad33a57f..9e6cd0aa2 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -446,7 +446,8 @@ function main_from_inputs( @warn "Fixed-boundary mode unstable for n = $nstring" end - # Compute free boundary energies + # Compute free boundary energies. + free_energies = nothing if ctrl.vac_flag && !(ctrl.ksing > 0 && ctrl.ksing <= intr.msing + 1) if ctrl.verbose wall_desc = intr.wall_settings.shape == "nowall" ? "no wall" : intr.wall_settings.shape @@ -481,7 +482,7 @@ function main_from_inputs( gal_data = nothing if ctrl.gal_flag gal_start = time() - gal_data = galerkin_solve(ctrl, equil, ffit, intr; wv=ctrl.vac_flag ? free_energies.wv : nothing) + gal_data = galerkin_solve(ctrl, equil, ffit, intr; wv=free_energies !== nothing ? free_energies.wv : nothing) @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s" end @@ -491,7 +492,7 @@ function main_from_inputs( equil, intr, odet, - ctrl.vac_flag ? free_energies : nothing, + free_energies, ffit, git_version, inputs, @@ -544,7 +545,7 @@ function main_from_inputs( slayer_result = _run_slayer_stage(nothing) @info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER" return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, - free_energies=ctrl.vac_flag ? free_energies : nothing, + free_energies=free_energies, slayer=slayer_result) end @@ -607,7 +608,7 @@ function main_from_inputs( # Run perturbed equilibrium calculations # Free-boundary wt0 drives the plasma inductance; mthvac sizes the Green's-function solves pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium( - equil, pe_odet, ctrl.vac_flag ? free_energies.wt0 : nothing, ctrl.mthvac, intr, ft_ctrl, pe_ctrl, pe_intr, + equil, pe_odet, free_energies !== nothing ? free_energies.wt0 : nothing, ctrl.mthvac, intr, ft_ctrl, pe_ctrl, pe_intr, metric, ffit ) @@ -678,7 +679,7 @@ function main_from_inputs( # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, - free_energies=ctrl.vac_flag ? free_energies : nothing, + free_energies=free_energies, slayer=slayer_result) end @@ -907,23 +908,23 @@ function write_outputs_to_HDF5( # working (Jacobian) coordinate. W_freeboundary_eigenmodes holds the generalized # eigenvectors, columns sorted most-unstable first, normalized to unit power norm with # the largest-magnitude entry made real-positive. - out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? free_energies.wt0 : ComplexF64[] - out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? free_energies.wp : ComplexF64[] - out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? free_energies.wv : ComplexF64[] - out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? free_energies.wt : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? free_energies.et : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? free_energies.ep : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? free_energies.ev : ComplexF64[] - out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = ctrl.vac_flag ? free_energies.vacuum_eigenvalue : NaN + out_h5["FreeBoundaryStability/W_freeboundary"] = free_energies !== nothing ? free_energies.wt0 : ComplexF64[] + out_h5["FreeBoundaryStability/W_plasma"] = free_energies !== nothing ? free_energies.wp : ComplexF64[] + out_h5["FreeBoundaryStability/W_vacuum"] = free_energies !== nothing ? free_energies.wv : ComplexF64[] + out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = free_energies !== nothing ? free_energies.wt : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_energies"] = free_energies !== nothing ? free_energies.et : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = free_energies !== nothing ? free_energies.ep : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = free_energies !== nothing ? free_energies.ev : ComplexF64[] + out_h5["FreeBoundaryStability/vacuum_eigenvalue"] = free_energies !== nothing ? free_energies.vacuum_eigenvalue : NaN # Cartesian surface point clouds used downstream for visualisation and # perturbed-equilibrium plotting. - out_h5["SurfaceGeometries/Plasma/x"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 1] : Float64[] - out_h5["SurfaceGeometries/Plasma/y"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 2] : Float64[] - out_h5["SurfaceGeometries/Plasma/z"] = ctrl.vac_flag ? free_energies.plasma_pts[:, 3] : Float64[] - out_h5["SurfaceGeometries/Wall/x"] = ctrl.vac_flag ? free_energies.wall_pts[:, 1] : Float64[] - out_h5["SurfaceGeometries/Wall/y"] = ctrl.vac_flag ? free_energies.wall_pts[:, 2] : Float64[] - out_h5["SurfaceGeometries/Wall/z"] = ctrl.vac_flag ? free_energies.wall_pts[:, 3] : Float64[] + out_h5["SurfaceGeometries/Plasma/x"] = free_energies !== nothing ? free_energies.plasma_pts[:, 1] : Float64[] + out_h5["SurfaceGeometries/Plasma/y"] = free_energies !== nothing ? free_energies.plasma_pts[:, 2] : Float64[] + out_h5["SurfaceGeometries/Plasma/z"] = free_energies !== nothing ? free_energies.plasma_pts[:, 3] : Float64[] + out_h5["SurfaceGeometries/Wall/x"] = free_energies !== nothing ? free_energies.wall_pts[:, 1] : Float64[] + out_h5["SurfaceGeometries/Wall/y"] = free_energies !== nothing ? free_energies.wall_pts[:, 2] : Float64[] + out_h5["SurfaceGeometries/Wall/z"] = free_energies !== nothing ? free_energies.wall_pts[:, 3] : Float64[] # Write kinetic parameters when kinetic mode is enabled if ctrl.kinetic_factor > 0 diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index c1f0786e6..297273818 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -154,7 +154,6 @@ boundary-integral solve produces along the way. - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on the 3D nowall path - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`) - - `mtheta`, `nzeta::Int`: Grid resolution the response was computed on """ struct VacuumResponse wv::Matrix{ComplexF64} @@ -162,8 +161,6 @@ struct VacuumResponse grre::Matrix{ComplexF64} plasma_pts::Matrix{Float64} wall_pts::Matrix{Float64} - mtheta::Int - nzeta::Int end """ @@ -179,9 +176,7 @@ function VacuumResponse(inputs::VacuumInput) zeros(ComplexF64, 2 * num_points, num_modes), zeros(ComplexF64, 2 * num_points, num_modes), zeros(num_points, 3), - zeros(num_points, 3), - inputs.mtheta, - inputs.nzeta + zeros(num_points, 3) ) end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index f44dff518..30af03b9a 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -67,7 +67,9 @@ of the Green's functions. mpert = length(inputs.m_modes) num_points_surf = inputs.mtheta - vac_data.wv .= 0 + fill!(vac_data.wv, 0) + fill!(vac_data.grri, 0) + fill!(vac_data.grre, 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) From 162faa5d95ca99124606fc9e286c1e72088612a2 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 14 Aug 2026 13:53:49 -0400 Subject: [PATCH 3/3] VACUUM - MINOR - Update Green's functions description and add tests for vacuum response buffer clearing --- src/Vacuum/DataTypes.jl | 2 +- test/runtests_vacuum.jl | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 297273818..63dec0b9e 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -152,7 +152,7 @@ boundary-integral solve produces along the way. ## Fields - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D - - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on the 3D nowall path + - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on every 3D path - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`) """ struct VacuumResponse diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index ff07f8be5..630c60436 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -441,6 +441,26 @@ @test vac.wall_pts ≈ ref.wall_pts end end + + @testset "in-place compute_vacuum_response! clears a reused buffer" begin + # The nowall path writes only the plasma rows of grri/grre, so a buffer left + # over from a wall run must not leak its wall rows into the next result. + inputs = _make_inputs() + wall_rows = (inputs.mtheta+1):(2*inputs.mtheta) + + vac = VacuumResponse(inputs) + compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="conformal", a=0.5)) + @test any(!iszero, view(vac.grre, wall_rows, :)) + + compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall")) + fresh = compute_vacuum_response(inputs, WallShapeSettings(; shape="nowall")) + + @test all(iszero, view(vac.grri, wall_rows, :)) + @test all(iszero, view(vac.grre, wall_rows, :)) + @test vac.wv ≈ fresh.wv + @test vac.grri ≈ fresh.grri + @test vac.grre ≈ fresh.grre + end end @testset "extract_plasma_surface_at_psi" begin