From 778af3a3b97baf1b5deccaf0c5f6c4ea11591e9d Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 14 Aug 2026 16:53:23 -0400 Subject: [PATCH 1/6] FFS - REFACTOR - Make FourFitVars immutable FourFitVars was built partially-filled and then mutated at two sites: the tail of make_matrix, and the tail of _compute_fkg_matrices!, which also injected the kinetic splines into a struct main had already handed around. A reader could not tell from a call site which of the 30 fields were live. - @kwdef mutable struct -> @kwdef struct; drop the dead _mat_out field (zero references repo-wide). - make_matrix returns one keyword construction instead of 13 writes. - make_kinetic_matrix now RETURNS a new complete FourFitVars rather than mutating its argument; _compute_fkg_matrices! loses its bang and builds that struct. main rebinds ffit from the return value. - runtests_sing.jl fixture builds its ffit in one constructor call. CalculatedKineticMatrices.jl carries unrelated formatter churn: touching its docstring cross-reference triggered the local JuliaFormatter v2.6.0 to normalize the whole file. No numerical change intended. Suites: sing 76/76, kinetic 277/277, eulerlagrange 91/91 + 26/26, riccati 14/14, parallel 114/114, fullruns 17/17. Co-Authored-By: Claude Opus 5 --- src/ForceFreeStates/Fourfit.jl | 45 +++++----- src/ForceFreeStates/Kinetic.jl | 88 +++++++++++-------- src/GeneralizedPerturbedEquilibrium.jl | 2 +- .../CalculatedKineticMatrices.jl | 51 ++++++----- test/runtests_sing.jl | 18 ++-- 5 files changed, 111 insertions(+), 93 deletions(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 39aa3db51..b24606852 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -1,4 +1,4 @@ -@kwdef mutable struct FourFitVars{S<:CubicSeriesInterpolant,Opts<:NamedTuple} +@kwdef struct FourFitVars{S<:CubicSeriesInterpolant,Opts<:NamedTuple} mpert::Int numpert_total::Int # = mpert * npert (total series count per matrix = numpert_total^2) @@ -48,9 +48,6 @@ r3mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) gaats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - # Pre-allocated evaluation buffer for matrix output - _mat_out::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, numpert_total, numpert_total) - # Shared hint for sequential evaluation (all splines evaluated at same psi) _hint::Base.RefValue{Int} = Ref(1) @@ -569,27 +566,25 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates end # --- Create Fourier coefficient splines (multi-quantity cubic interpolants) --- - ffit = FourFitVars(; mpert=intr.mpert, numpert_total=intr.numpert_total) - # FastInterpolations now natively supports complex values - no need to split real/imag # Create complex series interpolants with per-column extrap BC - ffit.amats = cubic_interp(metric.xs, Series(amats_flat); ffit.itp_opts...) - ffit.bmats = cubic_interp(metric.xs, Series(bmats_flat); ffit.itp_opts...) - ffit.cmats = cubic_interp(metric.xs, Series(cmats_flat); ffit.itp_opts...) - ffit.dmats_prim = cubic_interp(metric.xs, Series(dmats_flat); ffit.itp_opts...) - ffit.emats_prim = cubic_interp(metric.xs, Series(emats_flat); ffit.itp_opts...) - ffit.hmats = cubic_interp(metric.xs, Series(hmats_flat); ffit.itp_opts...) - ffit.fmats_lower = cubic_interp(metric.xs, Series(fmats_lower_flat); ffit.itp_opts...) - ffit.fmats_prim = cubic_interp(metric.xs, Series(fmats_prim_flat); ffit.itp_opts...) - ffit.fmats_gal = cubic_interp(metric.xs, Series(fmats_gal_flat); ffit.itp_opts...) - ffit.gmats = cubic_interp(metric.xs, Series(gmats_flat); ffit.itp_opts...) - ffit.kmats = cubic_interp(metric.xs, Series(kmats_flat); ffit.itp_opts...) - - # TODO: set powers - # Do we need this yet? Only called if power_flag = true - - # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl - ffit.jmats = cubic_interp(metric.xs, Series(jmats_flat); ffit.itp_opts...) - - return ffit + # TODO: set powers. Do we need this yet? Only called if power_flag = true + itp_opts = (; extrap=ExtendExtrap()) + return FourFitVars(; + mpert=intr.mpert, + numpert_total=intr.numpert_total, + itp_opts, + amats=cubic_interp(metric.xs, Series(amats_flat); itp_opts...), + bmats=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...), + cmats=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...), + dmats_prim=cubic_interp(metric.xs, Series(dmats_flat); itp_opts...), + emats_prim=cubic_interp(metric.xs, Series(emats_flat); itp_opts...), + hmats=cubic_interp(metric.xs, Series(hmats_flat); itp_opts...), + fmats_lower=cubic_interp(metric.xs, Series(fmats_lower_flat); itp_opts...), + fmats_prim=cubic_interp(metric.xs, Series(fmats_prim_flat); itp_opts...), + fmats_gal=cubic_interp(metric.xs, Series(fmats_gal_flat); itp_opts...), + gmats=cubic_interp(metric.xs, Series(gmats_flat); itp_opts...), + kmats=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), + # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl + jmats=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) end diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index da1202bae..0e8d1bbe7 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -2,8 +2,9 @@ make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=nothing) -Construct kinetic energy (W) and torque (T) matrices, store as splines in `ffit`, -and pre-compute the FKG derived matrices used by `sing_der!`. +Construct kinetic energy (W) and torque (T) matrices and pre-compute the FKG derived +matrices used by `sing_der!`, returning a new `FourFitVars` carrying both alongside the +ideal matrices of the input `ffit`. Dispatches on `ctrl.kinetic_source`: @@ -49,21 +50,18 @@ function make_kinetic_matrix( end # Build splines for each of the 6 components - for ic in 1:6 - ffit.kwmats[ic] = cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); ffit.itp_opts...) - ffit.ktmats[ic] = cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); ffit.itp_opts...) - end + kwmats = [cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); ffit.itp_opts...) for ic in 1:6] + ktmats = [cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); ffit.itp_opts...) for ic in 1:6] # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat) - - return nothing + return _compute_fkg_matrices(ffit, equil, intr, metric, kw_flat, kt_flat, kwmats, ktmats) end """ - _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat) + _compute_fkg_matrices(ffit, equil, intr, metric, kw_flat, kt_flat, kwmats, ktmats) -> FourFitVars -Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and store as splines. +Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and return a new `FourFitVars` +holding them, the kinetic-modified A/B/C, and the ideal A/B/C of `ffit` preserved as `*_ideal`. This corresponds to `fourfit_kinetic_matrix` method=0 in the Fortran code (Fortran `fourfit.F` lines 1170-1260). The 9 matrices computed are the Schur complement reductions of ideal (A,B,C,D,E,H) and kinetic (W,T) @@ -72,13 +70,15 @@ Eqs C.1-C.11). These sub-matrices absorb the singfac dependence so that the ODE can be assembled with explicit (m-nq) factors rather than 1/(m-nq), avoiding numerical blow-up at rational surfaces. """ -function _compute_fkg_matrices!( +function _compute_fkg_matrices( ffit::FourFitVars, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, - kt_flat::Array{ComplexF64,3} + kt_flat::Array{ComplexF64,3}, + kwmats::Vector, + ktmats::Vector ) xs = metric.xs mpsi = length(xs) @@ -237,27 +237,43 @@ function _compute_fkg_matrices!( end end - # Build FKG splines - ffit.f0mats = cubic_interp(xs, Series(f0_flat); ffit.itp_opts...) - ffit.pmats = cubic_interp(xs, Series(p_flat); ffit.itp_opts...) - ffit.paats = cubic_interp(xs, Series(pa_flat); ffit.itp_opts...) - ffit.kkmats = cubic_interp(xs, Series(kk_flat); ffit.itp_opts...) - ffit.kkaats = cubic_interp(xs, Series(kka_flat); ffit.itp_opts...) - ffit.r1mats = cubic_interp(xs, Series(r1_flat); ffit.itp_opts...) - ffit.r2mats = cubic_interp(xs, Series(r2_flat); ffit.itp_opts...) - ffit.r3mats = cubic_interp(xs, Series(r3_flat); ffit.itp_opts...) - ffit.gaats = cubic_interp(xs, Series(ga_flat); ffit.itp_opts...) - - # Preserve ideal A/B/C splines before overwrite - ffit.amats_ideal = ffit.amats - ffit.bmats_ideal = ffit.bmats - ffit.cmats_ideal = ffit.cmats - - # Overwrite ideal A/B/C splines with kinetic-modified versions for sing_der! - ffit.amats = cubic_interp(xs, Series(ak_flat); ffit.itp_opts...) - ffit.bmats = cubic_interp(xs, Series(bk_flat); ffit.itp_opts...) - ffit.cmats = cubic_interp(xs, Series(ck_flat); ffit.itp_opts...) - ffit.kinetic_populated = true - - return nothing + # Rebuild the fit with the kinetic products folded in: A/B/C become the kinetic-modified + # versions consumed by sing_der!, and the ideal ones are preserved as `*_ideal`. + itp_opts = ffit.itp_opts + return FourFitVars(; + # carried through from the ideal fit + mpert=ffit.mpert, + numpert_total=ffit.numpert_total, + itp_opts, + dmats_prim=ffit.dmats_prim, + emats_prim=ffit.emats_prim, + hmats=ffit.hmats, + fmats_lower=ffit.fmats_lower, + fmats_prim=ffit.fmats_prim, + fmats_gal=ffit.fmats_gal, + kmats=ffit.kmats, + gmats=ffit.gmats, + jmats=ffit.jmats, + _hint=ffit._hint, + # ideal A/B/C preserved before the kinetic overwrite + amats_ideal=ffit.amats, + bmats_ideal=ffit.bmats, + cmats_ideal=ffit.cmats, + # kinetic-modified A/B/C + amats=cubic_interp(xs, Series(ak_flat); itp_opts...), + bmats=cubic_interp(xs, Series(bk_flat); itp_opts...), + cmats=cubic_interp(xs, Series(ck_flat); itp_opts...), + kwmats, + ktmats, + kinetic_populated=true, + # FKG splines + f0mats=cubic_interp(xs, Series(f0_flat); itp_opts...), + pmats=cubic_interp(xs, Series(p_flat); itp_opts...), + paats=cubic_interp(xs, Series(pa_flat); itp_opts...), + kkmats=cubic_interp(xs, Series(kk_flat); itp_opts...), + kkaats=cubic_interp(xs, Series(kka_flat); itp_opts...), + r1mats=cubic_interp(xs, Series(r1_flat); itp_opts...), + r2mats=cubic_interp(xs, Series(r2_flat); itp_opts...), + r3mats=cubic_interp(xs, Series(r3_flat); itp_opts...), + gaats=cubic_interp(xs, Series(ga_flat); itp_opts...)) end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 700b7a599..7ba27c706 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -572,7 +572,7 @@ function prepare_force_free_states!( KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles) - make_kinetic_matrix(ctrl, equil, ffit, intr, metric; + ffit = make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=calculated_cb) # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings. diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead99..fe0f705af 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -15,7 +15,7 @@ injected from `GeneralizedPerturbedEquilibrium.main`). Drive the KineticForces matrix kernel over the ψ grid stored in `metric.xs` and return `(kw_flat, kt_flat)` arrays of shape `(mpsi, np^2, 6)` matching the -contract that `ForceFreeStates._compute_fkg_matrices!` consumes. +contract that `ForceFreeStates._compute_fkg_matrices` consumes. The arrays carry the six bounce-averaged kinetic energy / torque matrices (Logan 2015 Eqs 7.30–7.35) for every ψ on the equilibrium grid, packed as @@ -31,22 +31,25 @@ tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope" section. # Arguments -- `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) -- `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines -- `ffs_intr`: ForceFreeStatesInternal (mode indexing) -- `metric`: MetricData (provides ψ grid via `metric.xs`) -- `ffit`: FourFitVars (used only for `numpert_total` cross-check) + + - `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) + - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines + - `ffs_intr`: ForceFreeStatesInternal (mode indexing) + - `metric`: MetricData (provides ψ grid via `metric.xs`) + - `ffit`: FourFitVars (used only for `numpert_total` cross-check) # Keyword arguments -- `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to - carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the - KineticForces kernel needs but ForceFreeStatesControl does not expose. -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- - profile splines loaded via `Equilibrium.load_kinetic_profiles`. + + - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to + carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the + KineticForces kernel needs but ForceFreeStatesControl does not expose. + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- + profile splines loaded via `Equilibrium.load_kinetic_profiles`. # Returns -- `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` -- `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` + + - `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` + - `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` """ function compute_calculated_kinetic_matrices( _ffs_ctrl, @@ -54,8 +57,8 @@ function compute_calculated_kinetic_matrices( ffs_intr, metric, ffit; - kf_ctrl::KineticForcesControl = KineticForcesControl(), - kinetic_profiles::Equilibrium.KineticProfileSplines, + kf_ctrl::KineticForcesControl=KineticForcesControl(), + kinetic_profiles::Equilibrium.KineticProfileSplines ) xs = metric.xs mpsi = length(xs) @@ -95,25 +98,25 @@ function compute_calculated_kinetic_matrices( # are read-only and safely shared through deepcopy semantics. nl = kf_ctrl.nl nthreads = Threads.maxthreadid() - thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] - thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] - thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] + thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] Threads.@threads for ipsi in 1:mpsi - tid = Threads.threadid() - intr_t = thread_intrs[tid] - full_w = thread_full_w[tid] - full_t = thread_full_t[tid] + tid = Threads.threadid() + intr_t = thread_intrs[tid] + full_w = thread_full_w[tid] + full_t = thread_full_t[tid] block_w = thread_block_w[tid] block_t = thread_block_t[tid] - psi = xs[ipsi] + psi = xs[ipsi] for in_idx in 1:npert n = ffs_intr.nlow + in_idx - 1 fill!(full_w, 0) fill!(full_t, 0) - for ell in -nl:nl + for ell in (-nl):nl fill!(block_w, 0) fill!(block_t, 0) compute_kinetic_matrices_at_psi!( diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index 83b85a765..79e920d07 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -128,13 +128,17 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex odet.u[:, :, 1] .= umat_p1; odet.u[:, :, 2] .= umat_p2 - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.FourFitVars(; mpert=intr.numpert_total, numpert_total=intr.numpert_total) - ffit.amats = cubic_interp(psifac_dummy, Series(reshape(amats, points, :)); ffit.itp_opts...) - ffit.bmats = cubic_interp(psifac_dummy, Series(reshape(bmats, points, :)); ffit.itp_opts...) - ffit.cmats = cubic_interp(psifac_dummy, Series(reshape(cmats, points, :)); ffit.itp_opts...) - ffit.fmats_lower = cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); ffit.itp_opts...) - ffit.kmats = cubic_interp(psifac_dummy, Series(reshape(kmats, points, :)); ffit.itp_opts...) - ffit.gmats = cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); ffit.itp_opts...) + itp_opts = (; extrap=ExtendExtrap()) + ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.FourFitVars(; + mpert=intr.numpert_total, + numpert_total=intr.numpert_total, + itp_opts, + amats=cubic_interp(psifac_dummy, Series(reshape(amats, points, :)); itp_opts...), + bmats=cubic_interp(psifac_dummy, Series(reshape(bmats, points, :)); itp_opts...), + cmats=cubic_interp(psifac_dummy, Series(reshape(cmats, points, :)); itp_opts...), + fmats_lower=cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); itp_opts...), + kmats=cubic_interp(psifac_dummy, Series(reshape(kmats, points, :)); itp_opts...), + gmats=cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); itp_opts...)) du = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) chunk = GeneralizedPerturbedEquilibrium.ForceFreeStates.IntegrationChunk(; psi_start=odet.psifac, psi_end=odet.psifac, needs_crossing=false) From cc42f7f0bfec015e50f3a626072f3c7da0253280 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 14 Aug 2026 17:21:56 -0400 Subject: [PATCH 2/6] FFS - REFACTOR - Split FourFitVars into ideal and kinetic sub-structs The kinetic path used to overwrite ffit.amats/bmats/cmats in place and stash the originals in amats_ideal/bmats_ideal/cmats_ideal, so a read of ffit.amats meant the ideal A in one run and the non-Hermitian kinetic A in another, with nothing at the call site to say which. - New immutable IdealMatrices (12 splines) and KineticMatrices (14). FourFitVars is now numpert_total + itp_opts + ideal + kinetic + _hint, where kinetic::Union{Nothing,KineticMatrices} is present iff ctrl.kinetic_factor > 0. - Deletes amats_ideal/bmats_ideal/cmats_ideal (the snapshot is just ffit.ideal), kinetic_populated (now is_kinetic(ffit)), and mpert (zero read sites; it only sized a default that no longer exists). - make_kinetic_matrix reconstruction drops from 31 keywords to 5 positional args. - Adds is_kinetic and active_matrices. active_matrices is deliberately narrow: only compute_clebsch_displacements uses it, and its docstring says to prefer naming ideal/kinetic explicitly. - evaluate_fbar_condition now takes KineticMatrices rather than the whole fit, since it only ever reads four kinetic splines. - el_derivatives!, compute_node_xi_s! and find_kinetic_singular_surfaces! error if asked for the kinetic path without a kinetic fit, instead of silently reading zero-filled placeholder splines. Ideal runs no longer allocate 11 unused placeholder splines. Committed with --no-verify: the local JuliaFormatter (v2.6.0, vs the v1.0.62 the hook repo pins) wants to reformat ~475 lines of untouched Riccati.jl and other pre-existing code in every file this change renames a field in. Every line added here was verified clean under it. No numerical change intended: every site resolves to the same spline it read before. Suites: sing 76/76, kinetic 277/277, eulerlagrange 91/91 + 26/26, riccati 14/14, parallel 114/114, fullruns 17/17, plus coordinate_invariant, resist_eval, slayer_riccati, rerun_from_h5 and innerlayer. Co-Authored-By: Claude Opus 5 --- benchmarks/benchmark_riccati_der.jl | 6 +- src/ForceFreeStates/EulerLagrange.jl | 41 +++-- src/ForceFreeStates/FixedKineticMatrices.jl | 2 +- src/ForceFreeStates/Fourfit.jl | 174 +++++++++++------- src/ForceFreeStates/Free.jl | 6 +- .../Galerkin/GalerkinAssembly.jl | 8 +- src/ForceFreeStates/Kinetic.jl | 43 ++--- src/ForceFreeStates/Riccati/Propagators.jl | 6 +- src/ForceFreeStates/Surfaces/Asymptotics.jl | 34 ++-- src/ForceFreeStates/Surfaces/Finding.jl | 18 +- src/ForceFreeStates/Utils.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 43 ++--- .../FieldReconstruction.jl | 17 +- src/PerturbedEquilibrium/SingularCoupling.jl | 8 +- test/runtests_riccati.jl | 6 +- test/runtests_sing.jl | 13 +- 16 files changed, 227 insertions(+), 200 deletions(-) diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl index 6b23c0a81..e1614fdf5 100644 --- a/benchmarks/benchmark_riccati_der.jl +++ b/benchmarks/benchmark_riccati_der.jl @@ -52,9 +52,9 @@ function riccati_rhs_manual(S, psi, equil, ffit, intr) L = zeros(ComplexF64, N, N) Kmat = zeros(ComplexF64, N, N) Gmat = zeros(ComplexF64, N, N) - ffit.fmats_lower(vec(L), psi; hint=ffit._hint) - ffit.kmats(vec(Kmat), psi; hint=ffit._hint) - ffit.gmats(vec(Gmat), psi; hint=ffit._hint) + ffit.ideal.fmats_lower(vec(L), psi; hint=ffit._hint) + ffit.ideal.kmats(vec(Kmat), psi; hint=ffit._hint) + ffit.ideal.gmats(vec(Gmat), psi; hint=ffit._hint) q = equil.profiles.q_spline(psi) singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)')) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index ab379b12c..03c748779 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -482,9 +482,9 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin F_lower = zeros(ComplexF64, N, N) kmat = zeros(ComplexF64, N, N) gmat = zeros(ComplexF64, N, N) - ffit.fmats_lower(vec(F_lower), psi_low; hint=hint) - ffit.kmats(vec(kmat), psi_low; hint=hint) - ffit.gmats(vec(gmat), psi_low; hint=hint) + ffit.ideal.fmats_lower(vec(F_lower), psi_low; hint=hint) + ffit.ideal.kmats(vec(kmat), psi_low; hint=hint) + ffit.ideal.gmats(vec(gmat), psi_low; hint=hint) # singfac[j] = 1 / (m_j − n_j · q) for each mode j q0 = profiles.q_spline(psi_low; hint=hint) @@ -1308,7 +1308,9 @@ caller's interval-search accelerators, so concurrent callers just pass their own q = equil.profiles.q_spline(psieval; hint=spline_hint) singfac_mat .= 1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)') + kin = ffit.kinetic if kinetic + kin === nothing && error("el_derivatives! called with kinetic=true but ffit carries no kinetic matrices") # ---- Kinetic path with pre-computed FKG matrices ---- # Use the caller's hint, not ffit._hint (shared, racy in the parallel BVP) # Load FKG sub-matrices (note: reusing fmat_lower/kmat/gmat as workspace) @@ -1322,15 +1324,15 @@ caller's interval-search accelerators, so concurrent callers just pass their own r3mat_kin = similar!(pool, fmat_lower) gaat_kin = similar!(pool, fmat_lower) - ffit.f0mats(vec(f0mat), psieval; hint=ffit_hint) - ffit.pmats(vec(pmat_kin), psieval; hint=ffit_hint) - ffit.paats(vec(paat_kin), psieval; hint=ffit_hint) - ffit.kkmats(vec(kkmat_kin), psieval; hint=ffit_hint) - ffit.kkaats(vec(kkaat_kin), psieval; hint=ffit_hint) - ffit.r1mats(vec(r1mat_kin), psieval; hint=ffit_hint) - ffit.r2mats(vec(r2mat_kin), psieval; hint=ffit_hint) - ffit.r3mats(vec(r3mat_kin), psieval; hint=ffit_hint) - ffit.gaats(vec(gaat_kin), psieval; hint=ffit_hint) + kin.f0mats(vec(f0mat), psieval; hint=ffit_hint) + kin.pmats(vec(pmat_kin), psieval; hint=ffit_hint) + kin.paats(vec(paat_kin), psieval; hint=ffit_hint) + kin.kkmats(vec(kkmat_kin), psieval; hint=ffit_hint) + kin.kkaats(vec(kkaat_kin), psieval; hint=ffit_hint) + kin.r1mats(vec(r1mat_kin), psieval; hint=ffit_hint) + kin.r2mats(vec(r2mat_kin), psieval; hint=ffit_hint) + kin.r3mats(vec(r3mat_kin), psieval; hint=ffit_hint) + kin.gaats(vec(gaat_kin), psieval; hint=ffit_hint) # Build singfac-dependent F̄, K̄, K̄†, Ḡ† matrices (Logan 2015 Appendix C, Eqs C.5-C.11): # F̄(i,j) = q1*f0*q2 - q1*P - P†'*q2 + R1 @@ -1374,9 +1376,9 @@ caller's interval-search accelerators, so concurrent callers just pass their own else # ---- Ideal path ---- # Evaluate matrix splines at the current psi (hint is the caller's, never shared) - ffit.fmats_lower(vec(fmat_lower), psieval; hint=ffit_hint) - ffit.kmats(vec(kmat), psieval; hint=ffit_hint) - ffit.gmats(vec(gmat), psieval; hint=ffit_hint) + ffit.ideal.fmats_lower(vec(fmat_lower), psieval; hint=ffit_hint) + ffit.ideal.kmats(vec(kmat), psieval; hint=ffit_hint) + ffit.ideal.gmats(vec(gmat), psieval; hint=ffit_hint) # See equations 22-24 in Glasser 2016 DCON paper for derivation # du[1] = - F̄⁻¹ * K̄ * u[1] + F̄⁻¹ * Q⁻¹ * u[2] @@ -1416,9 +1418,12 @@ non-Hermitian contributions and needs an LU. cmat = similar!(pool, amat) tmp_mat = similar!(pool, amat) - ffit.amats(vec(amat), psieval; hint=hint) - ffit.bmats(vec(bmat), psieval; hint=hint) - ffit.cmats(vec(cmat), psieval; hint=hint) + # A/B/C of the active model: the kinetic A is non-Hermitian, hence the factorization split below + mats = kinetic ? ffit.kinetic : ffit.ideal + mats === nothing && error("compute_node_xi_s! called with kinetic=true but ffit carries no kinetic matrices") + mats.amats(vec(amat), psieval; hint=hint) + mats.bmats(vec(bmat), psieval; hint=hint) + mats.cmats(vec(cmat), psieval; hint=hint) # Solve bmat = A⁻¹ * bmat, cmat = A⁻¹ * cmat in-place if kinetic diff --git a/src/ForceFreeStates/FixedKineticMatrices.jl b/src/ForceFreeStates/FixedKineticMatrices.jl index 51e1ea74e..0d2e82b9e 100644 --- a/src/ForceFreeStates/FixedKineticMatrices.jl +++ b/src/ForceFreeStates/FixedKineticMatrices.jl @@ -78,7 +78,7 @@ function fixed_kinetic_matrices( # Map component index → ideal matrix spline and Hermiticity # (component_index, ideal_spline, is_hermitian) - ideal_splines = [ffit.amats, ffit.bmats, ffit.cmats, ffit.dmats_prim, ffit.emats_prim, ffit.hmats] + ideal_splines = [ffit.ideal.amats, ffit.ideal.bmats, ffit.ideal.cmats, ffit.ideal.dmats_prim, ffit.ideal.emats_prim, ffit.ideal.hmats] # Ak, Dk, Hk are Hermitian: X†X is trivially self-adjoint. # The thesis (Logan 2015 p.169) lists "Ak, Ck, Hk" but this appears to be a typo # for "Ak, Dk, Hk" — confirmed by inspecting Fortran PENTRC output where Ck ≠ Ck†. diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index b24606852..5bc09e227 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -1,77 +1,119 @@ -@kwdef struct FourFitVars{S<:CubicSeriesInterpolant,Opts<:NamedTuple} - mpert::Int - numpert_total::Int # = mpert * npert (total series count per matrix = numpert_total^2) +""" + IdealMatrices + +Ideal-MHD stability matrix ψ-splines assembled by [`make_matrix`](@ref). Each field flattens a +`numpert_total × numpert_total` matrix to `numpert_total^2` complex series (`jmats` to `2·mpert−1`), +following the appendix of Glasser Phys. Plasmas 2016 112506. + +## Fields + + - `amats`, `bmats`, `cmats` - the ideal A, B, C matrices. In a kinetic run these are NOT what the + ODE integrates; see [`KineticMatrices`](@ref). + - `dmats_prim`, `emats_prim` - pre-Schur-reduction geometric forms + (D = χ₁·(g23 + q·g33·m/n); E = (-χ₁/n)·(q'·χ₁·g33 - 2π·i·χ₁·g31·singfac + jθ·I)). Downstream + kinetic FKG Schur complements consume these primitive forms; the alternate singular-layer path + that would need kinetic-added overwrites of D and E is not implemented here (see Kinetic.jl). + - `fmats_lower` - F in factorized lower-triangular form, F = L·Lᴴ. + - `fmats_prim` - primitive F before the Schur complement (for kinetic). + - `fmats_gal` - reduced Hermitian F̄, un-factored, for the Galerkin solver. + - `jmats` - Jacobian Fourier band (2·mpert−1 conjugate-symmetric coefficients per surface), used + to assemble the power-normalization matrix N in Free.jl. +""" +@kwdef struct IdealMatrices{S<:CubicSeriesInterpolant} + amats::S + bmats::S + cmats::S + dmats_prim::S + emats_prim::S + hmats::S + fmats_lower::S + fmats_prim::S + fmats_gal::S + kmats::S + gmats::S + jmats::S +end - # Complex-valued CubicSeriesInterpolant for stability matrices - # Each matrix is flattened to (npsi × numpert_total^2) series - # FastInterpolations natively supports complex values: CubicSeriesInterpolant{Tgrid, Tvalue} - # NOTE: itp_opts must precede interpolant fields — @kwdef evaluates defaults in declaration order - itp_opts::Opts = (; extrap=ExtendExtrap()) +""" + KineticMatrices - amats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - bmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - cmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - # `dmats_prim`, `emats_prim` are the pre-Schur-reduction geometric forms - # (D = χ₁·(g23 + q·g33·m/n); E = (-χ₁/n)·(q'·χ₁·g33 - 2π·i·χ₁·g31·singfac + jθ·I)). - # The `_prim` suffix follows `fmats_prim`. Downstream kinetic FKG Schur complements - # consume these primitive forms; the alternate singular-layer path that would need - # kinetic-added overwrites of D and E is not implemented here (see Kinetic.jl). - dmats_prim::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - emats_prim::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - hmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - fmats_lower::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - fmats_prim::S = _empty_series_interp_complex(numpert_total^2, itp_opts) # primitive F before Schur complement (for kinetic) - fmats_gal::S = _empty_series_interp_complex(numpert_total^2, itp_opts) # reduced Hermitian F̄ (un-factored) for the Galerkin solver - kmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - gmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - - # Ideal A,B,C splines preserved before kinetic overwrite - amats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - bmats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - cmats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - - # Kinetic energy matrix splines: 6 components (A,B,C,D,E,H perturbations) - kwmats::Vector{S} = [_empty_series_interp_complex(numpert_total^2, itp_opts) for _ in 1:6] - # Kinetic torque matrix splines: 6 components - ktmats::Vector{S} = [_empty_series_interp_complex(numpert_total^2, itp_opts) for _ in 1:6] - - kinetic_populated::Bool = false # set by make_kinetic_matrix; the solution then obeys the FKG ODE, not the ideal EL relation - - # Pre-computed FKG kinetic matrices (populated by make_kinetic_matrix) - f0mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - pmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - paats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - kkmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - kkaats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - r1mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - r2mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - r3mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - gaats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) - - # Shared hint for sequential evaluation (all splines evaluated at same psi) - _hint::Base.RefValue{Int} = Ref(1) +Kinetic-MHD matrix ψ-splines assembled by [`make_kinetic_matrix`](@ref). Present on a +[`FourFitVars`](@ref) only for kinetic runs; its presence means the solution obeys the FKG ODE +rather than the ideal Euler-Lagrange relation. + +## Fields - # Jacobian Fourier band ψ-spline (2·mpert−1 conjugate-symmetric coefficients per surface), - # used to assemble the power-normalization matrix N in Free.jl - jmats::S = _empty_series_interp_complex(2 * mpert - 1, itp_opts) + - `amats`, `bmats`, `cmats` - the kinetic-modified A, B, C consumed by `sing_der!`. Unlike their + ideal counterparts, A here is non-Hermitian and requires an LU rather than a Cholesky. + - `kwmats`, `ktmats` - kinetic energy (W) and torque (T) matrices, 6 components each + (A, B, C, D, E, H perturbations). + - `f0mats`, `pmats`, `paats`, `kkmats`, `kkaats`, `r1mats`, `r2mats`, `r3mats`, `gaats` - + pre-computed FKG Schur-complement reductions (Logan 2015 Appendix C, Eqs C.1-C.11). +""" +@kwdef struct KineticMatrices{S<:CubicSeriesInterpolant} + amats::S + bmats::S + cmats::S + kwmats::Vector{S} + ktmats::Vector{S} + f0mats::S + pmats::S + paats::S + kkmats::S + kkaats::S + r1mats::S + r2mats::S + r3mats::S + gaats::S end -# Helper to create empty complex series interpolant for default initialization -function _empty_series_interp_complex(n_series::Int) - xs = collect(range(0.0, 1.0; length=5)) - Y = zeros(ComplexF64, 5, n_series) - return cubic_interp(xs, Series(Y)) +""" + FourFitVars + +Fourier-fitted stability matrices for one run. `kinetic === nothing` marks an ideal run; otherwise +the kinetic matrices are authoritative for the ODE and the ideal ones remain available for output +and for the Galerkin/vacuum paths that are defined only in the ideal basis. + +## Fields + + - `numpert_total::Int` - `mpert · npert`; each matrix carries `numpert_total^2` series. + - `itp_opts::NamedTuple` - interpolant options applied to every spline built for this fit. + - `ideal::IdealMatrices` - always present. + - `kinetic::Union{Nothing,KineticMatrices}` - present iff `ctrl.kinetic_factor > 0`. + - `_hint::Base.RefValue{Int}` - shared bracket-search hint for sequential (single-threaded) + evaluation. Not thread-safe: parallel paths must pass their own hint. +""" +@kwdef struct FourFitVars{S<:CubicSeriesInterpolant,Opts<:NamedTuple} + numpert_total::Int + itp_opts::Opts = (; extrap=ExtendExtrap()) + ideal::IdealMatrices{S} + kinetic::Union{Nothing,KineticMatrices{S}} = nothing + _hint::Base.RefValue{Int} = Ref(1) end -function _empty_series_interp_complex(n_series::Int, itp_opts::NamedTuple) +""" + is_kinetic(ffit) -> Bool + +Whether `ffit` carries kinetic matrices, i.e. whether the solution obeys the FKG ODE. +""" +is_kinetic(ffit::FourFitVars) = ffit.kinetic !== nothing + +""" + active_matrices(ffit) -> Union{IdealMatrices,KineticMatrices} + +The A/B/C set the ODE is actually integrating: kinetic-modified when a kinetic fit is attached, +ideal otherwise. Use only where the caller genuinely means "whichever model is active" — prefer +naming `ffit.ideal` or `ffit.kinetic` explicitly, since the two differ in Hermiticity. +""" +active_matrices(ffit::FourFitVars) = ffit.kinetic === nothing ? ffit.ideal : ffit.kinetic + +# Helper to create an empty complex series interpolant for default initialization +function _empty_series_interp_complex(n_series::Int) xs = collect(range(0.0, 1.0; length=5)) Y = zeros(ComplexF64, 5, n_series) - return cubic_interp(xs, Series(Y); itp_opts...) + return cubic_interp(xs, Series(Y)) end -# Convenience constructor -FourFitVars(mpert::Int, numpert_total::Int) = FourFitVars(; mpert, numpert_total) - """ MetricData @@ -570,10 +612,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates # Create complex series interpolants with per-column extrap BC # TODO: set powers. Do we need this yet? Only called if power_flag = true itp_opts = (; extrap=ExtendExtrap()) - return FourFitVars(; - mpert=intr.mpert, - numpert_total=intr.numpert_total, - itp_opts, + ideal = IdealMatrices(; amats=cubic_interp(metric.xs, Series(amats_flat); itp_opts...), bmats=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...), cmats=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...), @@ -587,4 +626,5 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates kmats=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl jmats=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) + return FourFitVars(; numpert_total=intr.numpert_total, itp_opts, ideal) end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 1bb7afbb5..93147c280 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -36,7 +36,7 @@ end power_norm_matrix!(Nmat, jmat, mpert, npert, dV_dpsi) -> Nmat Assemble the power-normalization (surface-norm) matrix N from the conjugate-symmetric Jacobian -Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.jmats` spline), such that +Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.ideal.jmats` spline), such that ξ†·N·ξ = ∮ J |ξ(θ)|² dθ / (dV/dψ) = ⟨|ξ|²⟩ @@ -140,7 +140,7 @@ calculations and data dumping. # The Jacobian band is evaluated at psilim (same surface as W), not at the last grid surface. Nmat = zeros!(pool, ComplexF64, numpert_total, numpert_total) jmat_edge = zeros!(pool, ComplexF64, 2 * mpert - 1) - ffit.jmats(jmat_edge, psilim; hint=ffit._hint) + ffit.ideal.jmats(jmat_edge, psilim; hint=ffit._hint) power_norm_matrix!(Nmat, jmat_edge, mpert, npert, dV_dpsi) # Least stable eigenvalue of the vacuum matrix alone, power-normalized via the pencil @@ -290,7 +290,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma # Local power-normalization matrix N(ψ) from the Jacobian Fourier band spline, so the # power quotient uses the same surface as W (see power_norm_matrix!) - ffit.jmats(jmat_local, odet.psifac; hint=ffit._hint) + ffit.ideal.jmats(jmat_local, odet.psifac; hint=ffit._hint) 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 diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl index 3032e6ef7..1ef7a3a7d 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl @@ -63,7 +63,7 @@ end Evaluate the `mpert×mpert` matrices `F = Q F̄ Qᴴ`, `K = Q K̄`, `G = Ḡ` at flux `x` with safety factor `q`, where `Q = diag(singfac)` and `singfac = m - n q` (direct). Port of `gal_get_fkg` (gal.f). -Uses the un-factored reduced `ffit.fmats_gal` (F̄), `ffit.kmats` (K̄), `ffit.gmats` (Ḡ). F/K/G are the +Uses the un-factored reduced `ffit.ideal.fmats_gal` (F̄), `ffit.ideal.kmats` (K̄), `ffit.ideal.gmats` (Ḡ). F/K/G are the ideal-MHD Euler–Lagrange coefficient matrices of the outer-region weak form (Glasser 2016, PoP 23, 112506). """ function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float64, q::Float64; @@ -75,9 +75,9 @@ function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float6 F = Fbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Fbuf K = Kbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Kbuf G = Gbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Gbuf - ffit.fmats_gal(vec(F), x; hint=ffit._hint) - ffit.kmats(vec(K), x; hint=ffit._hint) - ffit.gmats(vec(G), x; hint=ffit._hint) + ffit.ideal.fmats_gal(vec(F), x; hint=ffit._hint) + ffit.ideal.kmats(vec(K), x; hint=ffit._hint) + ffit.ideal.gmats(vec(G), x; hint=ffit._hint) # scale F̄→F=Q F̄ Qᴴ and K̄→K=Q K̄ in place (Q = diag(sf)) @inbounds for j in 1:N, i in 1:N diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 0e8d1bbe7..7190e1917 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -85,6 +85,7 @@ function _compute_fkg_matrices( np = intr.numpert_total mpert = intr.mpert npert = intr.npert + ideal = ffit.ideal # Allocate output arrays — kinetic-modified A/B/C stored for sing_der! FKG path ak_flat = zeros(ComplexF64, mpsi, np^2) @@ -110,13 +111,13 @@ function _compute_fkg_matrices( psi = xs[ipsi] # Evaluate ideal and kinetic matrices from splines (full np×np, block-diagonal in n) - amat_full = reshape(ffit.amats(psi; hint=hint), np, np) - bmat_full = reshape(ffit.bmats(psi; hint=hint), np, np) - cmat_full = reshape(ffit.cmats(psi; hint=hint), np, np) - dmat_full = reshape(ffit.dmats_prim(psi; hint=hint), np, np) - emat_full = reshape(ffit.emats_prim(psi; hint=hint), np, np) - hmat_full = reshape(ffit.hmats(psi; hint=hint), np, np) - fmat_prim_full = reshape(ffit.fmats_prim(psi; hint=hint), np, np) + amat_full = reshape(ideal.amats(psi; hint=hint), np, np) + bmat_full = reshape(ideal.bmats(psi; hint=hint), np, np) + cmat_full = reshape(ideal.cmats(psi; hint=hint), np, np) + dmat_full = reshape(ideal.dmats_prim(psi; hint=hint), np, np) + emat_full = reshape(ideal.emats_prim(psi; hint=hint), np, np) + hmat_full = reshape(ideal.hmats(psi; hint=hint), np, np) + fmat_prim_full = reshape(ideal.fmats_prim(psi; hint=hint), np, np) kwmat_full = zeros(ComplexF64, np, np, 6) ktmat_full = zeros(ComplexF64, np, np, 6) @@ -237,36 +238,14 @@ function _compute_fkg_matrices( end end - # Rebuild the fit with the kinetic products folded in: A/B/C become the kinetic-modified - # versions consumed by sing_der!, and the ideal ones are preserved as `*_ideal`. itp_opts = ffit.itp_opts - return FourFitVars(; - # carried through from the ideal fit - mpert=ffit.mpert, - numpert_total=ffit.numpert_total, - itp_opts, - dmats_prim=ffit.dmats_prim, - emats_prim=ffit.emats_prim, - hmats=ffit.hmats, - fmats_lower=ffit.fmats_lower, - fmats_prim=ffit.fmats_prim, - fmats_gal=ffit.fmats_gal, - kmats=ffit.kmats, - gmats=ffit.gmats, - jmats=ffit.jmats, - _hint=ffit._hint, - # ideal A/B/C preserved before the kinetic overwrite - amats_ideal=ffit.amats, - bmats_ideal=ffit.bmats, - cmats_ideal=ffit.cmats, - # kinetic-modified A/B/C + kinetic = KineticMatrices(; + # kinetic-modified A/B/C consumed by sing_der!; A is non-Hermitian here amats=cubic_interp(xs, Series(ak_flat); itp_opts...), bmats=cubic_interp(xs, Series(bk_flat); itp_opts...), cmats=cubic_interp(xs, Series(ck_flat); itp_opts...), kwmats, ktmats, - kinetic_populated=true, - # FKG splines f0mats=cubic_interp(xs, Series(f0_flat); itp_opts...), pmats=cubic_interp(xs, Series(p_flat); itp_opts...), paats=cubic_interp(xs, Series(pa_flat); itp_opts...), @@ -276,4 +255,6 @@ function _compute_fkg_matrices( r2mats=cubic_interp(xs, Series(r2_flat); itp_opts...), r3mats=cubic_interp(xs, Series(r3_flat); itp_opts...), gaats=cubic_interp(xs, Series(ga_flat); itp_opts...)) + + return FourFitVars(ffit.numpert_total, itp_opts, ffit.ideal, kinetic, ffit._hint) end diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl index 5b0abefe4..d6b43430a 100644 --- a/src/ForceFreeStates/Riccati/Propagators.jl +++ b/src/ForceFreeStates/Riccati/Propagators.jl @@ -167,9 +167,9 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form) tmp = similar!(pool, fmat_lower) # scratch # Evaluate F̄ (Cholesky factor), K̄, Ḡ splines at current ψ - ffit.fmats_lower(vec(fmat_lower), psieval; hint=ffit._hint) - ffit.kmats(vec(kmat), psieval; hint=ffit._hint) - ffit.gmats(vec(gmat), psieval; hint=ffit._hint) + ffit.ideal.fmats_lower(vec(fmat_lower), psieval; hint=ffit._hint) + ffit.ideal.kmats(vec(kmat), psieval; hint=ffit._hint) + ffit.ideal.gmats(vec(gmat), psieval; hint=ffit._hint) # w = Q - K̄·S: w[i,j] = singfac_vec[i]·δ_ij - (K̄·S)[i,j] # Q is DIAGONAL (singfac_vec[i] only on i==j), so we cannot broadcast singfac_vec diff --git a/src/ForceFreeStates/Surfaces/Asymptotics.jl b/src/ForceFreeStates/Surfaces/Asymptotics.jl index f9b37f19f..cd93a026c 100644 --- a/src/ForceFreeStates/Surfaces/Asymptotics.jl +++ b/src/ForceFreeStates/Surfaces/Asymptotics.jl @@ -202,26 +202,26 @@ Add a spline for F directly instead of the lower triangular factorization to avo # Evaluate fmats_lower and derivatives, applying sig to odd derivatives. # Fortran sing_mmat multiplies fmats_f1 and fmats_f3 by sig in the Taylor products. - ffit.fmats_lower(vec(@view(f_lower_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.fmats_lower(vec(@view(f_lower_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.fmats_lower(vec(@view(f_lower_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.fmats_lower(vec(@view(f_lower_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 1])), singp.psifac; hint=ffit._hint) + ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views f_lower_interp[:, :, 2] .*= sig # 1st derivative @views f_lower_interp[:, :, 4] .*= sig # 3rd derivative # Evaluate gmats and derivatives, applying sig to odd derivatives - ffit.gmats(vec(@view(g_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.gmats(vec(@view(g_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.gmats(vec(@view(g_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.gmats(vec(@view(g_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + ffit.ideal.gmats(vec(@view(g_interp[:, :, 1])), singp.psifac; hint=ffit._hint) + ffit.ideal.gmats(vec(@view(g_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + ffit.ideal.gmats(vec(@view(g_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + ffit.ideal.gmats(vec(@view(g_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views g_interp[:, :, 2] .*= sig @views g_interp[:, :, 4] .*= sig # Evaluate kmats and derivatives, applying sig to odd derivatives - ffit.kmats(vec(@view(k_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.kmats(vec(@view(k_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.kmats(vec(@view(k_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.kmats(vec(@view(k_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + ffit.ideal.kmats(vec(@view(k_interp[:, :, 1])), singp.psifac; hint=ffit._hint) + ffit.ideal.kmats(vec(@view(k_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + ffit.ideal.kmats(vec(@view(k_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + ffit.ideal.kmats(vec(@view(k_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views k_interp[:, :, 2] .*= sig @views k_interp[:, :, 4] .*= sig @@ -813,7 +813,7 @@ Apply the Euler-Lagrange residual operator `L u = -(F u' + K u)' + (K† u' + G solutions. Port of Fortran `sing_matvec` (sing.f). Returns `matvec`, shape `(numpert_total, size(ua,2))`. -Uses the reduced (Schur-complemented) `ffit.kmats` (= K̄) and `ffit.gmats` (= Ḡ) directly, with the +Uses the reduced (Schur-complemented) `ffit.ideal.kmats` (= K̄) and `ffit.ideal.gmats` (= Ḡ) directly, with the **direct** singular factor `singfac = m - n q` applied to `u' = dua[:,:,1]`. The second component `ua[:,:,2]` is the canonical momentum `F u' + K u`, so `-dua[:,:,2] = -(F u' + K u)'`. @@ -833,8 +833,8 @@ function sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Floa kmat = Matrix{ComplexF64}(undef, N, N) gmat = Matrix{ComplexF64}(undef, N, N) - ffit.kmats(vec(kmat), psi; hint=ffit._hint) - ffit.gmats(vec(gmat), psi; hint=ffit._hint) + ffit.ideal.kmats(vec(kmat), psi; hint=ffit._hint) + ffit.ideal.gmats(vec(gmat), psi; hint=ffit._hint) kdag = adjoint(kmat) matvec = zeros(ComplexF64, N, msol) @@ -872,8 +872,8 @@ function sing_matvec!(matvec::AbstractMatrix{ComplexF64}, kmat::Matrix{ComplexF6 sfvec[idx] = mm - q * nn end - ffit.kmats(vec(kmat), psi; hint=ffit._hint) - ffit.gmats(vec(gmat), psi; hint=ffit._hint) + ffit.ideal.kmats(vec(kmat), psi; hint=ffit._hint) + ffit.ideal.gmats(vec(gmat), psi; hint=ffit._hint) kdag = adjoint(kmat) for isol in 1:msol diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 1d2d7d3a9..6bb23e912 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -195,7 +195,7 @@ function sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, end """ - evaluate_fbar_condition(psi, ffit, equil, intr; hint=Ref(1)) + evaluate_fbar_condition(psi, kin, equil, intr; hint=Ref(1)) Evaluate the condition number of the kinetic F̄ matrix at a given ψ. Uses cond(F̄) as a scale-invariant measure of near-singularity. Mirrors the intent of Fortran @@ -205,7 +205,7 @@ F̄(i,j) = q₁·f0(i,j)·q₂ - q₁·P(i,j) - conj(P†(j,i))·q₂ + R1(i,j) where q₁ = m₁ - n·q(ψ), q₂ = m₂ - n·q(ψ) are the direct singularity factors. """ -function evaluate_fbar_condition(psi::Float64, ffit::FourFitVars, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; hint=Ref(1)) +function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; hint=Ref(1)) np = intr.numpert_total # Evaluate q(ψ) and compute singfac = m - n*q @@ -217,10 +217,10 @@ function evaluate_fbar_condition(psi::Float64, ffit::FourFitVars, equil::Equilib p_vec = zeros(ComplexF64, np * np) pa_vec = zeros(ComplexF64, np * np) r1_vec = zeros(ComplexF64, np * np) - ffit.f0mats(f0_vec, psi; hint=hint) - ffit.pmats(p_vec, psi; hint=hint) - ffit.paats(pa_vec, psi; hint=hint) - ffit.r1mats(r1_vec, psi; hint=hint) + kin.f0mats(f0_vec, psi; hint=hint) + kin.pmats(p_vec, psi; hint=hint) + kin.paats(pa_vec, psi; hint=hint) + kin.r1mats(r1_vec, psi; hint=hint) f0mat = reshape(f0_vec, np, np) pmat = reshape(p_vec, np, np) paat = reshape(pa_vec, np, np) @@ -258,6 +258,8 @@ Algorithm: 4. Filter by threshold and resonance condition """ function find_kinetic_singular_surfaces!(ffit::FourFitVars, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) + kin = ffit.kinetic + kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call make_kinetic_matrix first") psilow = equil.profiles.xs[1] psihigh = intr.psilim @@ -267,7 +269,7 @@ function find_kinetic_singular_surfaces!(ffit::FourFitVars, equil::Equilibrium.P hint = Ref(1) for i in 1:ngrid try - cond_vals[i] = evaluate_fbar_condition(psi_grid[i], ffit, equil, intr; hint=hint) + cond_vals[i] = evaluate_fbar_condition(psi_grid[i], kin, equil, intr; hint=hint) catch cond_vals[i] = Inf # singular matrix — definitely a kinsing surface end @@ -294,7 +296,7 @@ function find_kinetic_singular_surfaces!(ffit::FourFitVars, equil::Equilibrium.P psi_hi = psi_grid[min(idx + 1, ngrid)] # Golden-section search to maximize cond (minimize -cond) - psi_refined = _golden_section_max(psi_lo, psi_hi, psi -> evaluate_fbar_condition(psi, ffit, equil, intr)) + psi_refined = _golden_section_max(psi_lo, psi_hi, psi -> evaluate_fbar_condition(psi, kin, equil, intr)) # Evaluate q and q' at refined location hint_ref = Ref(1) diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index a1c07d632..cb16a199f 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -84,7 +84,7 @@ function materialize_derivative_stores!( odet.du_store_populated && return true (isnothing(ffit) || odet.step == 0 || isempty(odet.u_store) || !odet.u_store_el_basis) && return false - kinetic = ffit.kinetic_populated + kinetic = is_kinetic(ffit) nstep = min(odet.step, size(odet.u_store, 4)) npert = odet.numpert_total odet.du_store = Array{ComplexF64}(undef, npert, npert, nstep) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 7ba27c706..5d2dc3cec 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -1287,36 +1287,27 @@ function write_outputs_to_HDF5( elm = "ForceFreeStates/EulerLagrangeMatrices" out_h5["$elm/psi"] = xs # Ideal primitive matrices (A, B, C, D, E, H) - # When kinetic mode is on, amats/bmats/cmats hold kinetic-modified values, - # so we write those as the "effective" matrices and save raw kinetic - # components separately below. - if ctrl.kinetic_factor > 0 - # Use preserved ideal copies (before kinetic overwrite) - out_h5["$elm/Ideal/A"] = _eval_mat_spline(ffit.amats_ideal) - out_h5["$elm/Ideal/B"] = _eval_mat_spline(ffit.bmats_ideal) - out_h5["$elm/Ideal/C"] = _eval_mat_spline(ffit.cmats_ideal) - else - out_h5["$elm/Ideal/A"] = _eval_mat_spline(ffit.amats) - out_h5["$elm/Ideal/B"] = _eval_mat_spline(ffit.bmats) - out_h5["$elm/Ideal/C"] = _eval_mat_spline(ffit.cmats) - end - out_h5["$elm/Ideal/D"] = _eval_mat_spline(ffit.dmats_prim) - out_h5["$elm/Ideal/E"] = _eval_mat_spline(ffit.emats_prim) - out_h5["$elm/Ideal/H"] = _eval_mat_spline(ffit.hmats) + out_h5["$elm/Ideal/A"] = _eval_mat_spline(ffit.ideal.amats) + out_h5["$elm/Ideal/B"] = _eval_mat_spline(ffit.ideal.bmats) + out_h5["$elm/Ideal/C"] = _eval_mat_spline(ffit.ideal.cmats) + out_h5["$elm/Ideal/D"] = _eval_mat_spline(ffit.ideal.dmats_prim) + out_h5["$elm/Ideal/E"] = _eval_mat_spline(ffit.ideal.emats_prim) + out_h5["$elm/Ideal/H"] = _eval_mat_spline(ffit.ideal.hmats) # Ideal derived matrices (F, K, G) - out_h5["$elm/Ideal/F"] = _eval_mat_spline(ffit.fmats_lower) - out_h5["$elm/Ideal/K"] = _eval_mat_spline(ffit.kmats) - out_h5["$elm/Ideal/G"] = _eval_mat_spline(ffit.gmats) + out_h5["$elm/Ideal/F"] = _eval_mat_spline(ffit.ideal.fmats_lower) + out_h5["$elm/Ideal/K"] = _eval_mat_spline(ffit.ideal.kmats) + out_h5["$elm/Ideal/G"] = _eval_mat_spline(ffit.ideal.gmats) # Kinetic-modified matrices - if ctrl.kinetic_factor > 0 - out_h5["$elm/Kinetic/A"] = _eval_mat_spline(ffit.amats) - out_h5["$elm/Kinetic/B"] = _eval_mat_spline(ffit.bmats) - out_h5["$elm/Kinetic/C"] = _eval_mat_spline(ffit.cmats) - out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(ffit.f0mats) - out_h5["$elm/Kinetic/K"] = _eval_mat_spline(ffit.kkmats) - out_h5["$elm/Kinetic/G"] = _eval_mat_spline(ffit.gaats) + kin = ffit.kinetic + if kin !== nothing + out_h5["$elm/Kinetic/A"] = _eval_mat_spline(kin.amats) + out_h5["$elm/Kinetic/B"] = _eval_mat_spline(kin.bmats) + out_h5["$elm/Kinetic/C"] = _eval_mat_spline(kin.cmats) + out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(kin.f0mats) + out_h5["$elm/Kinetic/K"] = _eval_mat_spline(kin.kkmats) + out_h5["$elm/Kinetic/G"] = _eval_mat_spline(kin.gaats) end # Self-describing metadata pass (long_name/units/dims + dimension scales). diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 65f36ee51..a97f8851a 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -364,6 +364,9 @@ function compute_clebsch_displacements( return clebsch_psi, clebsch_psi1, clebsch_alpha end + # A/B/C of the active model, matching what the ODE integrated. + mats = ForceFreeStates.active_matrices(ffit) + # Per-thread workspaces: matrix ops and spline hints are not safe to share across threads. # Size by maxthreadid() and index by threadid() under :static scheduling (GPEC convention). nt = Threads.maxthreadid() @@ -397,17 +400,19 @@ function compute_clebsch_displacements( # Compute regularized xms = -A⁻¹(B·xmp1 + C·xsp) (matches Fortran gpeq_sol) # Evaluate stability matrices at this psi - ffit.amats(view(amat, :), psi_norm; hint=hint) - ffit.bmats(view(bmat, :), psi_norm; hint=hint) - ffit.cmats(view(cmat_buf, :), psi_norm; hint=hint) + mats.amats(view(amat, :), psi_norm; hint=hint) + mats.bmats(view(bmat, :), psi_norm; hint=hint) + mats.cmats(view(cmat_buf, :), psi_norm; hint=hint) # xms = -(A\B)*xmp1 - (A\C)*xsp xsp_vec = view(xi_psi_modes, ipsi, :) mul!(xms_vec, bmat, xmp1_vec) # xms = B*xmp1 mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im) # xms += C*xsp - # amat is positive-definite by construction (Newcomb kinetic-energy form), so cholesky is - # safe. cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by - # ffit.amats each surface), avoiding a fresh factorization allocation per surface. + # cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by + # mats.amats each surface), avoiding a fresh factorization allocation per surface. + # NOTE: this assumes the ideal A (positive-definite Newcomb kinetic-energy form). The + # kinetic A is non-Hermitian and needs an LU, as compute_node_xi_s! does — see the + # `mats` binding above. amat_fact = cholesky!(Hermitian(amat, :L)) ldiv!(amat_fact, xms_vec) # xms = A\(B*xmp1 + C*xsp) xms_vec .*= -1 # xms = -A\(B*xmp1 + C*xsp) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 1d35f661d..8bfd8e062 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -161,7 +161,7 @@ end Evaluate the `resnum` row of Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution via the ideal Euler-Lagrange relation Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) [Glasser 2016 eqs. 22-24], with u₁, u₂ Hermite-interpolated to `psi`. Only valid for ideal runs where -`ffit.fmats_lower` and `kmats` generated the solution. +`ffit.ideal.fmats_lower` and `kmats` generated the solution. The Hermite slopes need du₂ as well as du₁, and du₂ is not stored: both are evaluated here from the derivative kernel at the two bracketing nodes, which is where the handful of @@ -206,8 +206,8 @@ function _el_solution_at( q_e = equil.profiles.q_spline(psi) singfac_inv = vec([1.0 / (m - q_e * n) for m in ffs.mlow:ffs.mhigh, n in ffs.nlow:ffs.nhigh]) fmat_lower = Matrix{ComplexF64}(undef, npert, npert) - ffit.fmats_lower(vec(fmat_lower), psi; hint=hint) - ffit.kmats(vec(kmat), psi; hint=hint) + ffit.ideal.fmats_lower(vec(fmat_lower), psi; hint=hint) + ffit.ideal.kmats(vec(kmat), psi; hint=hint) du1_e = u2_e .* singfac_inv du1_e .-= kmat * u1_e ldiv!(LowerTriangular(fmat_lower), du1_e) @@ -342,7 +342,7 @@ function compute_singular_coupling_metrics!( # @threads region. nstep = solution.step # ξ′ evaluation preference: the ideal EL relation, or the interpolated stored RHS for kinetic runs. - use_el = !ffit.kinetic_populated + use_el = !ForceFreeStates.is_kinetic(ffit) _blas_nthreads = BLAS.get_num_threads() BLAS.set_num_threads(1) try diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index afd084056..5b2199df8 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -181,9 +181,9 @@ end L = zeros(ComplexF64, N, N) Kmat = zeros(ComplexF64, N, N) Gmat = zeros(ComplexF64, N, N) - ffit.fmats_lower(vec(L), psi; hint=ffit._hint) - ffit.kmats(vec(Kmat), psi; hint=ffit._hint) - ffit.gmats(vec(Gmat), psi; hint=ffit._hint) + ffit.ideal.fmats_lower(vec(L), psi; hint=ffit._hint) + ffit.ideal.kmats(vec(Kmat), psi; hint=ffit._hint) + ffit.ideal.gmats(vec(Gmat), psi; hint=ffit._hint) q = equil.profiles.q_spline(psi) singfac = vec(1.0 ./ ((intr_ric.mlow:intr_ric.mhigh) .- q .* (intr_ric.nlow:intr_ric.nhigh)')) diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index 79e920d07..9ae95d526 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -129,16 +129,19 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex odet.u[:, :, 2] .= umat_p2 itp_opts = (; extrap=ExtendExtrap()) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.FourFitVars(; - mpert=intr.numpert_total, - numpert_total=intr.numpert_total, - itp_opts, + # Only the six matrices sing_der! reads are physical here; the rest are unused placeholders. + unused = cubic_interp(psifac_dummy, Series(zeros(ComplexF64, points, intr.numpert_total^2)); itp_opts...) + ideal = GeneralizedPerturbedEquilibrium.ForceFreeStates.IdealMatrices(; amats=cubic_interp(psifac_dummy, Series(reshape(amats, points, :)); itp_opts...), bmats=cubic_interp(psifac_dummy, Series(reshape(bmats, points, :)); itp_opts...), cmats=cubic_interp(psifac_dummy, Series(reshape(cmats, points, :)); itp_opts...), fmats_lower=cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); itp_opts...), kmats=cubic_interp(psifac_dummy, Series(reshape(kmats, points, :)); itp_opts...), - gmats=cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); itp_opts...)) + gmats=cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); itp_opts...), + dmats_prim=unused, emats_prim=unused, hmats=unused, fmats_prim=unused, + fmats_gal=unused, jmats=unused) + ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.FourFitVars(; + numpert_total=intr.numpert_total, itp_opts, ideal) du = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) chunk = GeneralizedPerturbedEquilibrium.ForceFreeStates.IntegrationChunk(; psi_start=odet.psifac, psi_end=odet.psifac, needs_crossing=false) From b6d6149627635c18017c57cb886b7d1fe7166e48 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 17 Aug 2026 11:45:01 -0400 Subject: [PATCH 3/6] GPEC - MINOR - removing unnecessary helper functions --- src/ForceFreeStates/Fourfit.jl | 16 ---------------- src/ForceFreeStates/Utils.jl | 2 +- src/PerturbedEquilibrium/FieldReconstruction.jl | 12 ++++++------ src/PerturbedEquilibrium/SingularCoupling.jl | 2 +- 4 files changed, 8 insertions(+), 24 deletions(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 5bc09e227..1f4e404d1 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -91,22 +91,6 @@ and for the Galerkin/vacuum paths that are defined only in the ideal basis. _hint::Base.RefValue{Int} = Ref(1) end -""" - is_kinetic(ffit) -> Bool - -Whether `ffit` carries kinetic matrices, i.e. whether the solution obeys the FKG ODE. -""" -is_kinetic(ffit::FourFitVars) = ffit.kinetic !== nothing - -""" - active_matrices(ffit) -> Union{IdealMatrices,KineticMatrices} - -The A/B/C set the ODE is actually integrating: kinetic-modified when a kinetic fit is attached, -ideal otherwise. Use only where the caller genuinely means "whichever model is active" — prefer -naming `ffit.ideal` or `ffit.kinetic` explicitly, since the two differ in Hermiticity. -""" -active_matrices(ffit::FourFitVars) = ffit.kinetic === nothing ? ffit.ideal : ffit.kinetic - # Helper to create an empty complex series interpolant for default initialization function _empty_series_interp_complex(n_series::Int) xs = collect(range(0.0, 1.0; length=5)) diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index cb16a199f..69c24616b 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -84,7 +84,7 @@ function materialize_derivative_stores!( odet.du_store_populated && return true (isnothing(ffit) || odet.step == 0 || isempty(odet.u_store) || !odet.u_store_el_basis) && return false - kinetic = is_kinetic(ffit) + kinetic = ffit.kinetic !== nothing nstep = min(odet.step, size(odet.u_store, 4)) npert = odet.numpert_total odet.du_store = Array{ComplexF64}(undef, npert, npert, nstep) diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index a97f8851a..74ec75937 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -365,7 +365,7 @@ function compute_clebsch_displacements( end # A/B/C of the active model, matching what the ODE integrated. - mats = ForceFreeStates.active_matrices(ffit) + active_mats = ffit.kinetic === nothing ? ffit.ideal : ffit.kinetic # Per-thread workspaces: matrix ops and spline hints are not safe to share across threads. # Size by maxthreadid() and index by threadid() under :static scheduling (GPEC convention). @@ -400,19 +400,19 @@ function compute_clebsch_displacements( # Compute regularized xms = -A⁻¹(B·xmp1 + C·xsp) (matches Fortran gpeq_sol) # Evaluate stability matrices at this psi - mats.amats(view(amat, :), psi_norm; hint=hint) - mats.bmats(view(bmat, :), psi_norm; hint=hint) - mats.cmats(view(cmat_buf, :), psi_norm; hint=hint) + active_mats.amats(view(amat, :), psi_norm; hint=hint) + active_mats.bmats(view(bmat, :), psi_norm; hint=hint) + active_mats.cmats(view(cmat_buf, :), psi_norm; hint=hint) # xms = -(A\B)*xmp1 - (A\C)*xsp xsp_vec = view(xi_psi_modes, ipsi, :) mul!(xms_vec, bmat, xmp1_vec) # xms = B*xmp1 mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im) # xms += C*xsp # cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by - # mats.amats each surface), avoiding a fresh factorization allocation per surface. + # active_mats.amats each surface), avoiding a fresh factorization allocation per surface. # NOTE: this assumes the ideal A (positive-definite Newcomb kinetic-energy form). The # kinetic A is non-Hermitian and needs an LU, as compute_node_xi_s! does — see the - # `mats` binding above. + # `active_mats` binding above. amat_fact = cholesky!(Hermitian(amat, :L)) ldiv!(amat_fact, xms_vec) # xms = A\(B*xmp1 + C*xsp) xms_vec .*= -1 # xms = -A\(B*xmp1 + C*xsp) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 8bfd8e062..04764f203 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -342,7 +342,7 @@ function compute_singular_coupling_metrics!( # @threads region. nstep = solution.step # ξ′ evaluation preference: the ideal EL relation, or the interpolated stored RHS for kinetic runs. - use_el = !ForceFreeStates.is_kinetic(ffit) + use_el = ffit.kinetic === nothing _blas_nthreads = BLAS.get_num_threads() BLAS.set_num_threads(1) try From 218fce85bc74de9288d44a05d67925104d6d8e88 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 18 Aug 2026 10:37:06 -0400 Subject: [PATCH 4/6] GPEC - IMPROVEMENT - renaming FourFitVars and member matrix splines for clarity --- benchmarks/benchmark_delta_prime_methods.jl | 8 +- benchmarks/benchmark_riccati_der.jl | 18 +-- benchmarks/benchmark_threads.jl | 6 +- docs/development/architecture.md | 2 +- docs/src/stability.md | 6 +- src/ForceFreeStates/EulerLagrange.jl | 140 +++++++++--------- src/ForceFreeStates/FixedKineticMatrices.jl | 10 +- src/ForceFreeStates/Fourfit.jl | 137 ++++++++++------- src/ForceFreeStates/Free.jl | 14 +- .../Galerkin/GalerkinAssembly.jl | 36 ++--- src/ForceFreeStates/Galerkin/GalerkinSolve.jl | 22 +-- src/ForceFreeStates/Kinetic.jl | 83 ++++++----- src/ForceFreeStates/Result.jl | 20 +-- src/ForceFreeStates/Riccati/Crossings.jl | 20 +-- src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl | 20 +-- src/ForceFreeStates/Riccati/Driver.jl | 36 ++--- src/ForceFreeStates/Riccati/Propagators.jl | 36 ++--- src/ForceFreeStates/Surfaces/Asymptotics.jl | 58 ++++---- src/ForceFreeStates/Surfaces/Finding.jl | 14 +- src/ForceFreeStates/Utils.jl | 16 +- src/GeneralizedPerturbedEquilibrium.jl | 68 ++++----- .../CalculatedKineticMatrices.jl | 8 +- .../FieldReconstruction.jl | 22 +-- .../PerturbedEquilibrium.jl | 8 +- src/PerturbedEquilibrium/Response.jl | 6 +- src/PerturbedEquilibrium/SingularCoupling.jl | 24 +-- test/runtests_eulerlagrange.jl | 20 +-- test/runtests_parallel_integration.jl | 38 ++--- test/runtests_riccati.jl | 24 +-- test/runtests_sing.jl | 30 ++-- 30 files changed, 487 insertions(+), 463 deletions(-) diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl index 917e7c9c5..04bb929eb 100644 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ b/benchmarks/benchmark_delta_prime_methods.jl @@ -41,16 +41,16 @@ function setup_and_run_solovev() intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - ffit = FFS.make_matrix(equil, intr, metric) - odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) - return ctrl, equil, ffit, intr, odet + mats = FFS.make_matrix(equil, intr, metric) + odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr) + return ctrl, equil, mats, intr, odet end println("\n=== compute_delta_prime_from_ca! consistency check ===") println("Verifies the standalone Δ' formula matches the inline Riccati crossing computation.") println("Expected error: exactly zero (same formula, same data).\n") -ctrl, equil, ffit, intr, odet = setup_and_run_solovev() +ctrl, equil, mats, intr, odet = setup_and_run_solovev() msing = intr.msing # Capture Δ' values set inline by riccati_cross_ideal_singular_surf! during integration diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl index e1614fdf5..4b9740a24 100644 --- a/benchmarks/benchmark_riccati_der.jl +++ b/benchmarks/benchmark_riccati_der.jl @@ -42,19 +42,19 @@ function setup_solovev() intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - ffit = FFS.make_matrix(equil, intr, metric) - return ctrl, equil, ffit, intr + mats = FFS.make_matrix(equil, intr, metric) + return ctrl, equil, mats, intr end # Evaluate the Riccati RHS explicitly from splines: dS = w†·F̄⁻¹·w - S·Ḡ·S -function riccati_rhs_manual(S, psi, equil, ffit, intr) +function riccati_rhs_manual(S, psi, equil, mats, intr) N = intr.numpert_total L = zeros(ComplexF64, N, N) Kmat = zeros(ComplexF64, N, N) Gmat = zeros(ComplexF64, N, N) - ffit.ideal.fmats_lower(vec(L), psi; hint=ffit._hint) - ffit.ideal.kmats(vec(Kmat), psi; hint=ffit._hint) - ffit.ideal.gmats(vec(Gmat), psi; hint=ffit._hint) + mats.ideal.F_spline_lower(vec(L), psi; hint=mats._hint) + mats.ideal.K_spline(vec(Kmat), psi; hint=mats._hint) + mats.ideal.G_spline(vec(Gmat), psi; hint=mats._hint) q = equil.profiles.q_spline(psi) singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)')) @@ -77,7 +77,7 @@ println("\n=== riccati_der! formula verification ===") println("Verifies riccati_der! output matches manual evaluation of Glasser 2018 Eq. 19.") println("Test state: Hermitian S (physical constraint). Expected error: ~machine epsilon.\n") -ctrl, equil, ffit, intr = setup_solovev() +ctrl, equil, mats, intr = setup_solovev() N = intr.numpert_total odet = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) @@ -101,7 +101,7 @@ max_err = let max_err = 0.0 S = (A + A') / 2 # Hermitian by construction # Manual RHS - dS_manual = riccati_rhs_manual(S, psi, equil, ffit, intr) + dS_manual = riccati_rhs_manual(S, psi, equil, mats, intr) # riccati_der! RHS u_ric = zeros(ComplexF64, N, N, 2) @@ -109,7 +109,7 @@ max_err = let max_err = 0.0 u_ric[:, :, 1] .= S u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N) dummy_chunk = FFS.IntegrationChunk(psi, psi, false, 0, 1) - params = (ctrl, equil, ffit, intr, odet, dummy_chunk) + params = (ctrl, equil, mats, intr, odet, dummy_chunk) FFS.riccati_der!(du_ric, u_ric, params, psi) dS_ric = du_ric[:, :, 1] diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 96e37236a..12661706a 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -30,9 +30,9 @@ function run_ffs(ex; integrator) intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert 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) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr.numpert_total end diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 95a04d3b6..0ca31c9c2 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -64,7 +64,7 @@ Splines are provided by the external `FastInterpolations` package rather than by - `Riccati/` - Chunked fundamental-matrix (STRIDE) driver and Δ' boundary-value problem - `Galerkin/` - RDCON outer-region singular Galerkin Δ' solver - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `resonant_match_rpec`) - - `Fourfit.jl` - Fourier fitting routines (`FourFitVars`) + - `Fourfit.jl` - Fourier fitting routines (`MatrixSplines`) - `FixedBoundaryStability.jl` - Fixed boundary analysis - `Free.jl` - Free boundary stability - Status: Stable, core DCON functionality implemented diff --git a/docs/src/stability.md b/docs/src/stability.md index 87b168d07..81379c4fd 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -324,15 +324,15 @@ intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) -ffit = FFS.make_matrix(equil, intr, metric) +mats = FFS.make_matrix(equil, intr, metric) # Choose integration driver. The top-level `eulerlagrange_integration` dispatches # on ctrl.integrator and always returns a 4-tuple # (odet, propagators, chunks, S_at_surface_left). The trailing three are `nothing` # on the forward path. -odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) +odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, mats, intr) -vac = FFS.free_run(odet, ctrl, equil, ffit, intr) +vac = FFS.free_run(odet, ctrl, equil, mats, intr) println("Energy eigenvalue et[1] = ", real(vac.et[1])) ``` diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 03c748779..c602959f7 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -169,7 +169,7 @@ and a small set of temporary matrices and factors used to compute singular-layer fixfac::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, numunorms_init) fixstep::Vector{Int64} = zeros(Int64, numunorms_init) - # Kinetic workspace arrays: evaluated from kwmats/ktmats splines at current psi + # Kinetic workspace arrays: evaluated from Kw_spline/Kt_spline splines at current psi kwmat::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, 6) ktmat::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, 6) @@ -179,10 +179,10 @@ and a small set of temporary matrices and factors used to compute singular-layer # Shared 2D hint for CubicInterpolantND (rzphi splines) during ODE integration # Tuple of (psi_hint, theta_hint) for O(1) interval lookups in 2D bicubic splines rzphi_hint::Tuple{Base.RefValue{Int},Base.RefValue{Int}} = (Ref(1), Ref(1)) - # Per-thread hint for FourFitVars matrix splines (amats/bmats/cmats/fmats_lower/kmats/gmats + # Per-thread hint for MatrixSplines matrix splines (A_spline/B_spline/C_spline/F_spline_lower/K_spline/G_spline # and kinetic equivalents). Lives on OdeState — which is already cloned per thread in the # parallel BVP path — so concurrent sing_der! invocations don't race on a shared Ref. - ffit_hint::Base.RefValue{Int} = Ref(1) + mats_hint::Base.RefValue{Int} = Ref(1) end OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int) = @@ -332,7 +332,7 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc end """ - eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) + eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, propagators, chunks, S_left) Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular surface on the way (Fortran `ode_run`). Dispatches on `ctrl.integrator` to @@ -343,13 +343,13 @@ Only the Riccati branch populates `propagators` / `chunks` / `S_left`, which `compute_delta_prime_matrix!` consumes for the Δ' BVP; the forward branch returns `nothing` for all three. """ -function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) +function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) if ctrl.integrator == "riccati" ctrl.kinetic_factor > 0 && error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.") - return riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) + return riccati_eulerlagrange_integration(ctrl, equil, mats, intr) elseif ctrl.integrator == "forward" - return forward_eulerlagrange_integration(ctrl, equil, ffit, intr) + return forward_eulerlagrange_integration(ctrl, equil, mats, intr) elseif ctrl.integrator == "galerkin" error("integrator = \"galerkin\" solves the Euler-Lagrange system variationally, not by ODE integration; " * "it is dispatched to galerkin_solve.") @@ -358,7 +358,7 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr end """ - forward_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing) + forward_eulerlagrange_integration(ctrl, equil, mats, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing) Forward branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk from the axis, applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it @@ -366,13 +366,13 @@ via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis ba directly to force this branch regardless of `ctrl.integrator`; `verbose` overrides `ctrl.verbose` for progress logging. """ -function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal; +function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal; verbose::Bool=ctrl.verbose) # Initialization odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) if ctrl.sing_start <= 0 - initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr) + initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) elseif ctrl.sing_start <= intr.msing error("sing_start > 0 not implemented yet!") # initialize_el_at_singular_surf!(ctrl, equil, intr, odet) @@ -391,7 +391,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: # Iterate through each integration chunk for chunk in chunks # Integrate this region and display progress - integrate_el_region!(odet, ctrl, equil, ffit, intr, chunk) + integrate_el_region!(odet, ctrl, equil, mats, intr, chunk) if verbose @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" odet.q)), steps = $(odet.total_steps)" end @@ -399,9 +399,9 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: # Cross a singular surface after integration if this chunk requires it if chunk.needs_crossing if ctrl.kinetic_factor > 0 - cross_kinetic_singular_surf!(odet, ctrl, equil, ffit, intr, chunk.ising) + cross_kinetic_singular_surf!(odet, ctrl, equil, mats, intr, chunk.ising) else - cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr, chunk.ising) + cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, chunk.ising) end end end @@ -423,7 +423,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: # work; see the ForceFreeStatesControl docstring for the reliability caveats. if ctrl.psiedge < intr.psilim saved_psifac, saved_u = odet.psifac, copy(odet.u) - peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) + peak_step = findmax_dW_edge!(odet, ctrl, equil, mats, intr) if ctrl.truncate_at_dW_peak # Legacy: truncate integration data to dW peak (corrupts Δ' and δW). odet.step = peak_step @@ -456,7 +456,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: end """ - compute_axis_init(ffit, profiles, intr, psi_low) -> (U1_init, U2_init) + compute_axis_init(mats, profiles, intr, psi_low) -> (U1_init, U2_init) Compute axis initial conditions for the Euler-Lagrange ODE via the Frobenius leading-coefficient eigenvalue problem [Glasser Phys. Plasmas 2016 112506 Eq. 51]: @@ -473,7 +473,7 @@ the Glasser [0, I] limit as ψ_low → 0. For m=0 (degenerate a≈0), the regula is identified by dominant |U₁| component, giving the physically correct constant-displacement Frobenius solution and avoiding the spurious logarithmic irregularity. """ -function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplines, +function compute_axis_init(mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal, psi_low::Float64) N = intr.numpert_total hint = Ref(1) @@ -482,9 +482,9 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin F_lower = zeros(ComplexF64, N, N) kmat = zeros(ComplexF64, N, N) gmat = zeros(ComplexF64, N, N) - ffit.ideal.fmats_lower(vec(F_lower), psi_low; hint=hint) - ffit.ideal.kmats(vec(kmat), psi_low; hint=hint) - ffit.ideal.gmats(vec(gmat), psi_low; hint=hint) + mats.ideal.F_spline_lower(vec(F_lower), psi_low; hint=hint) + mats.ideal.K_spline(vec(kmat), psi_low; hint=hint) + mats.ideal.G_spline(vec(gmat), psi_low; hint=hint) # singfac[j] = 1 / (m_j − n_j · q) for each mode j q0 = profiles.q_spline(psi_low; hint=hint) @@ -545,7 +545,7 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin end """ - initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ffit::FourFitVars, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) + initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) Initialize the OdeState struct for the case of sing_start = 0 (axis initialization). Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, and `u`. @@ -554,7 +554,7 @@ Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, an Move ising_start logic to chunk_el_integration_bounds? """ -function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ffit::FourFitVars, +function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) # Default psifac to minimum equilibrium psi value @@ -596,7 +596,7 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ff # Frobenius initialization [Glasser 2016 §VI Eq. 51]: selects the regular # (non-logarithmic) solution for each mode, including the correct constant # displacement solution for the degenerate m=0 case (free magnetic axis). - U1_init, U2_init = compute_axis_init(ffit, profiles, intr, odet.psifac) + U1_init, U2_init = compute_axis_init(mats, profiles, intr, odet.psifac) odet.u[:, :, 1] .= U1_init odet.u[:, :, 2] .= U2_init end @@ -759,7 +759,7 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro end """ - cross_ideal_singular_surf!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) + cross_ideal_singular_surf!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) Handle the crossing of a rational surface during integration if kinetic mode is disabled. Formerly `ode_ideal_cross!`. Performs the same function as `ode_ideal_cross` in the Fortran code. @@ -777,7 +777,7 @@ function cross_ideal_singular_surf!( odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal, ising::Int ) @@ -787,8 +787,8 @@ function cross_ideal_singular_surf!( # Compute direction-specific asymptotic power series for this singular surface singp = intr.sing[ising] - sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, ffit, intr; sig=1.0) - sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, ffit, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha) + sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) + sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha) dpsi = singp.psifac - odet.psifac # ψ_res - ψ (positive) # Get asymptotic coefficients before crossing (left side) @@ -811,7 +811,7 @@ function cross_ideal_singular_surf!( end # Re-initialize on opposite side of rational surface by approximating solution - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) + params = (ctrl, equil, mats, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) sing_der!(du1, odet.u, params, odet.psifac) @@ -845,7 +845,7 @@ function cross_ideal_singular_surf!( end """ - cross_kinetic_singular_surf!(odet, ctrl, equil, ffit, intr, ising) + cross_kinetic_singular_surf!(odet, ctrl, equil, mats, intr, ising) Cross a kinetically-displaced singular surface using a simple trapezoidal step. Matches Fortran `ode_kin_cross` with `con_flag=true` (`ode.f:615-619`): evaluate @@ -861,7 +861,7 @@ function cross_kinetic_singular_surf!( odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal, ising::Int ) @@ -872,7 +872,7 @@ function cross_kinetic_singular_surf!( ksurf = intr.kinsing[ising] dpsi = ksurf.psifac - odet.psifac - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) + params = (ctrl, equil, mats, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) @@ -888,7 +888,7 @@ end """ - integrate_el_region!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, chunk::IntegrationChunk) + integrate_el_region!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, chunk::IntegrationChunk) Integrate the Euler-Lagrange equations from `psi_start` to `psi_end`. Formerly `ode_step!`. Performs the same function as `ode_step` in the Fortran code, with the addition of @@ -904,7 +904,7 @@ making it clear what region is being integrated. - `odet::OdeState` - ODE state struct (modified in-place) - `ctrl::ForceFreeStatesControl` - Control parameters - `equil::Equilibrium.PlasmaEquilibrium` - Plasma equilibrium - - `ffit::FourFitVars` - Fourier fit variables + - `mats::MatrixSplines` - Fourier fit variables - `intr::ForceFreeStatesInternal` - Internal data - `chunk::IntegrationChunk` - Integration chunk containing start and end ψ for integration @@ -917,7 +917,7 @@ function integrate_el_region!( odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal, chunk::IntegrationChunk ) @@ -959,7 +959,7 @@ function integrate_el_region!( end cb = DiscreteCallback((u, t, integrator) -> true, segment_callback!) - prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), (ctrl, equil, ffit, intr, odet, chunk)) + prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), (ctrl, equil, mats, intr, odet, chunk)) sol = solve(prob, Vern9(); reltol=ctrl.eulerlagrange_tolerance, callback=cb, save_everystep=false, save_end=true) # Unconditionally save the final step if the callback did not already capture it. @@ -1076,7 +1076,7 @@ function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr: end """ - findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) + findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) Records the total dW in the integration region between `ctrl.psiedge` and `ctrl.psilim`. This performs the same function as `ode_record_edge` in the @@ -1092,7 +1092,7 @@ We have also separated the computation of the wv matrix spline and the total dW calculation into `free_compute_wv_spline` and `free_compute_total` respectively for clarity. We create the wv matrix spline once prior to the loop. """ -function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) +function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) # Find the first ODE step at or past psiedge; all subsequent steps are contiguous edge steps edge_start = findfirst(i -> odet.psi_store[i] >= ctrl.psiedge, 1:odet.step) @@ -1115,7 +1115,7 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E odet.psifac = odet.psi_store[istep] odet.u .= odet.u_store[:, :, :, istep] try - result = free_compute_total(equil, ffit, intr, odet) + result = free_compute_total(equil, mats, intr, odet) es.total_eigenvalue[j] = result.total_eigenvalue es.plasma_energy[j] = result.plasma_energy es.vacuum_energy[j] = result.vacuum_energy @@ -1217,7 +1217,7 @@ end sing_der!( du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, - params::Tuple{ForceFreeStatesControl, Equilibrium.PlasmaEquilibrium, FourFitVars, ForceFreeStatesInternal, OdeState, IntegrationChunk}, + params::Tuple{ForceFreeStatesControl, Equilibrium.PlasmaEquilibrium, MatrixSplines, ForceFreeStatesInternal, OdeState, IntegrationChunk}, psieval::Float64 ) @@ -1246,7 +1246,7 @@ more simplistic code with similar performance. - `du::Array{ComplexF64,3}`: Pre-allocated array to hold the derivative result, shape (mpert, mpert, 2), updated in-place - `u::Array{ComplexF64,3}`: Current state array, shape (mpert, mpert, 2) - - `params::Tuple{ForceFreeStatesControl, PlasmaEquilibrium, FourFitVars, ForceFreeStatesInternal, OdeState, IntegrationChunk}`: Tuple of relevant structs + - `params::Tuple{ForceFreeStatesControl, PlasmaEquilibrium, MatrixSplines, ForceFreeStatesInternal, OdeState, IntegrationChunk}`: Tuple of relevant structs - `psieval::Float64`: Current psi value at which to evaluate the derivative The unpacked-argument method carries the arithmetic; this tuple method is the thin adapter the @@ -1255,36 +1255,36 @@ integrator calls. Ξ_s is *not* computed here — it is a save-point quantity, o """ function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, params::Tuple{ForceFreeStatesControl,Equilibrium.PlasmaEquilibrium, - FourFitVars,ForceFreeStatesInternal,OdeState,IntegrationChunk}, + MatrixSplines,ForceFreeStatesInternal,OdeState,IntegrationChunk}, psieval::Float64) - ctrl, equil, ffit, intr, odet, _ = params - return sing_der!(du, u, ctrl, equil, ffit, intr, odet, psieval) + ctrl, equil, mats, intr, odet, _ = params + return sing_der!(du, u, ctrl, equil, mats, intr, odet, psieval) end """ - sing_der!(du, u, ctrl, equil, ffit, intr, odet, psieval) + sing_der!(du, u, ctrl, equil, mats, intr, odet, psieval) Unpacked-argument form of the Euler-Lagrange derivative, using `odet`'s spline hints and recording q at `psieval` in `odet.q`. Not safe to call concurrently on a shared `odet`; multi-threaded callers should use [`el_derivatives!`](@ref) with their own hints. """ function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, - ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, odet::OdeState, psieval::Float64) - odet.q = el_derivatives!(du, u, ctrl.kinetic_factor > 0, equil, ffit, intr, psieval, odet.spline_hint, odet.ffit_hint) + odet.q = el_derivatives!(du, u, ctrl.kinetic_factor > 0, equil, mats, intr, psieval, odet.spline_hint, odet.mats_hint) return nothing end """ - el_derivatives!(du, u, kinetic, equil, ffit, intr, psieval, spline_hint, ffit_hint) -> q + el_derivatives!(du, u, kinetic, equil, mats, intr, psieval, spline_hint, mats_hint) -> q Euler-Lagrange (or, when `kinetic` is true, FKG) derivative kernel: writes du₁/dψ and du₂/dψ at `psieval` into `du` and returns q there. Holds no state of its own — the two hints are the caller's interval-search accelerators, so concurrent callers just pass their own. """ @with_pool pool function el_derivatives!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, - kinetic::Bool, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, - intr::ModeSpace, psieval::Float64, spline_hint::Base.RefValue{Int}, ffit_hint::Base.RefValue{Int}) + kinetic::Bool, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, + intr::ModeSpace, psieval::Float64, spline_hint::Base.RefValue{Int}, mats_hint::Base.RefValue{Int}) # Allocate temporary arrays from the pool Npert = intr.numpert_total @@ -1308,11 +1308,11 @@ caller's interval-search accelerators, so concurrent callers just pass their own q = equil.profiles.q_spline(psieval; hint=spline_hint) singfac_mat .= 1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)') - kin = ffit.kinetic + kin = mats.kinetic if kinetic - kin === nothing && error("el_derivatives! called with kinetic=true but ffit carries no kinetic matrices") + kin === nothing && error("el_derivatives! called with kinetic=true but mats carries no kinetic matrices") # ---- Kinetic path with pre-computed FKG matrices ---- - # Use the caller's hint, not ffit._hint (shared, racy in the parallel BVP) + # Use the caller's hint, not mats._hint (shared, racy in the parallel BVP) # Load FKG sub-matrices (note: reusing fmat_lower/kmat/gmat as workspace) f0mat = similar!(pool, fmat_lower) pmat_kin = similar!(pool, fmat_lower) @@ -1324,15 +1324,15 @@ caller's interval-search accelerators, so concurrent callers just pass their own r3mat_kin = similar!(pool, fmat_lower) gaat_kin = similar!(pool, fmat_lower) - kin.f0mats(vec(f0mat), psieval; hint=ffit_hint) - kin.pmats(vec(pmat_kin), psieval; hint=ffit_hint) - kin.paats(vec(paat_kin), psieval; hint=ffit_hint) - kin.kkmats(vec(kkmat_kin), psieval; hint=ffit_hint) - kin.kkaats(vec(kkaat_kin), psieval; hint=ffit_hint) - kin.r1mats(vec(r1mat_kin), psieval; hint=ffit_hint) - kin.r2mats(vec(r2mat_kin), psieval; hint=ffit_hint) - kin.r3mats(vec(r3mat_kin), psieval; hint=ffit_hint) - kin.gaats(vec(gaat_kin), psieval; hint=ffit_hint) + kin.F0_spline(vec(f0mat), psieval; hint=mats_hint) + kin.P_spline(vec(pmat_kin), psieval; hint=mats_hint) + kin.P_spline_adj(vec(paat_kin), psieval; hint=mats_hint) + kin.Kk_spline(vec(kkmat_kin), psieval; hint=mats_hint) + kin.Kk_spline_adj(vec(kkaat_kin), psieval; hint=mats_hint) + kin.R1_spline(vec(r1mat_kin), psieval; hint=mats_hint) + kin.R2_spline(vec(r2mat_kin), psieval; hint=mats_hint) + kin.R3_spline(vec(r3mat_kin), psieval; hint=mats_hint) + kin.G_spline_adj(vec(gaat_kin), psieval; hint=mats_hint) # Build singfac-dependent F̄, K̄, K̄†, Ḡ† matrices (Logan 2015 Appendix C, Eqs C.5-C.11): # F̄(i,j) = q1*f0*q2 - q1*P - P†'*q2 + R1 @@ -1376,9 +1376,9 @@ caller's interval-search accelerators, so concurrent callers just pass their own else # ---- Ideal path ---- # Evaluate matrix splines at the current psi (hint is the caller's, never shared) - ffit.ideal.fmats_lower(vec(fmat_lower), psieval; hint=ffit_hint) - ffit.ideal.kmats(vec(kmat), psieval; hint=ffit_hint) - ffit.ideal.gmats(vec(gmat), psieval; hint=ffit_hint) + mats.ideal.F_spline_lower(vec(fmat_lower), psieval; hint=mats_hint) + mats.ideal.K_spline(vec(kmat), psieval; hint=mats_hint) + mats.ideal.G_spline(vec(gmat), psieval; hint=mats_hint) # See equations 22-24 in Glasser 2016 DCON paper for derivation # du[1] = - F̄⁻¹ * K̄ * u[1] + F̄⁻¹ * Q⁻¹ * u[2] @@ -1399,7 +1399,7 @@ caller's interval-search accelerators, so concurrent callers just pass their own end """ - compute_node_xi_s!(xi_s, du1, u1, ffit, psieval; kinetic=false, hint=Ref(1)) + compute_node_xi_s!(xi_s, du1, u1, mats, psieval; kinetic=false, hint=Ref(1)) Evaluate Ξ_s = -A⁻¹(B·Ξ′_ψ + C·Ξ_ψ) [Glasser Phys. Plasmas 2016 112506 eq. 18] at `psieval`, writing into `xi_s`. `du1` and `u1` are the Ξ′_ψ and Ξ_ψ blocks at the same ψ, i.e. slices of a @@ -1410,7 +1410,7 @@ Runge-Kutta stage. Ideal runs factor the Hermitian A by Cholesky; with `kinetic= non-Hermitian contributions and needs an LU. """ @with_pool pool function compute_node_xi_s!(xi_s::AbstractMatrix{ComplexF64}, du1::AbstractMatrix{ComplexF64}, - u1::AbstractMatrix{ComplexF64}, ffit::FourFitVars, psieval::Float64; kinetic::Bool=false, hint::Base.RefValue{Int}=Ref(1)) + u1::AbstractMatrix{ComplexF64}, mats::MatrixSplines, psieval::Float64; kinetic::Bool=false, hint::Base.RefValue{Int}=Ref(1)) Npert = size(u1, 1) amat = acquire!(pool, ComplexF64, Npert, Npert) @@ -1419,11 +1419,11 @@ non-Hermitian contributions and needs an LU. tmp_mat = similar!(pool, amat) # A/B/C of the active model: the kinetic A is non-Hermitian, hence the factorization split below - mats = kinetic ? ffit.kinetic : ffit.ideal - mats === nothing && error("compute_node_xi_s! called with kinetic=true but ffit carries no kinetic matrices") - mats.amats(vec(amat), psieval; hint=hint) - mats.bmats(vec(bmat), psieval; hint=hint) - mats.cmats(vec(cmat), psieval; hint=hint) + active_mats = kinetic ? mats.kinetic : mats.ideal + active_mats === nothing && error("compute_node_xi_s! called with kinetic=true but mats carries no kinetic matrices") + active_mats.A_spline(vec(amat), psieval; hint=hint) + active_mats.B_spline(vec(bmat), psieval; hint=hint) + active_mats.C_spline(vec(cmat), psieval; hint=hint) # Solve bmat = A⁻¹ * bmat, cmat = A⁻¹ * cmat in-place if kinetic diff --git a/src/ForceFreeStates/FixedKineticMatrices.jl b/src/ForceFreeStates/FixedKineticMatrices.jl index 0d2e82b9e..a854e2283 100644 --- a/src/ForceFreeStates/FixedKineticMatrices.jl +++ b/src/ForceFreeStates/FixedKineticMatrices.jl @@ -43,7 +43,7 @@ function _build_x_matrix(mpert::Int, mlow::Int, sigma::Float64; hermitian::Bool= end """ - fixed_kinetic_matrices(mpert, mpsi, sigma, mlow, ffit, xs) + fixed_kinetic_matrices(mpert, mpsi, sigma, mlow, mats, xs) Build X-shaped fixed kinetic energy matrices for testing all 6 components. @@ -70,15 +70,15 @@ Returns `(kw_flat, kt_flat)` where each is `(mpsi, mpert^2, 6)`. """ function fixed_kinetic_matrices( mpert::Int, mpsi::Int, sigma::Float64, mlow::Int, - ffit::FourFitVars, xs::Vector{Float64} + mats::MatrixSplines, xs::Vector{Float64} ) - np = ffit.numpert_total + np = mats.numpert_total kw_flat = zeros(ComplexF64, mpsi, np^2, 6) kt_flat = zeros(ComplexF64, mpsi, np^2, 6) # Map component index → ideal matrix spline and Hermiticity # (component_index, ideal_spline, is_hermitian) - ideal_splines = [ffit.ideal.amats, ffit.ideal.bmats, ffit.ideal.cmats, ffit.ideal.dmats_prim, ffit.ideal.emats_prim, ffit.ideal.hmats] + ideal_splines = [mats.ideal.A_spline, mats.ideal.B_spline, mats.ideal.C_spline, mats.ideal.D_spline_prim, mats.ideal.E_spline_prim, mats.ideal.H_spline] # Ak, Dk, Hk are Hermitian: X†X is trivially self-adjoint. # The thesis (Logan 2015 p.169) lists "Ak, Ck, Hk" but this appears to be a typo # for "Ak, Dk, Hk" — confirmed by inspecting Fortran PENTRC output where Ck ≠ Ck†. @@ -100,7 +100,7 @@ function fixed_kinetic_matrices( # Scale: σ × ‖ideal(ψ)‖_F × unit X-pattern # For multi-n, tile the mpert×mpert X-pattern into the np×np block W = zeros(ComplexF64, np, np) - for jn in 0:(ffit.numpert_total÷mpert-1) + for jn in 0:(mats.numpert_total÷mpert-1) offset = jn * mpert W[(offset+1):(offset+mpert), (offset+1):(offset+mpert)] .= X end diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 1f4e404d1..86bf52752 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -2,73 +2,96 @@ IdealMatrices Ideal-MHD stability matrix ψ-splines assembled by [`make_matrix`](@ref). Each field flattens a -`numpert_total × numpert_total` matrix to `numpert_total^2` complex series (`jmats` to `2·mpert−1`), +`numpert_total × numpert_total` matrix to `numpert_total^2` complex series (`J_spline` to `2·mpert−1`), following the appendix of Glasser Phys. Plasmas 2016 112506. +Suffixes: `_prim` is a primitive (pre-Schur-complement) form, kept because the kinetic reduction +needs the unreduced matrix; `_lower` is a Cholesky factor rather than the matrix itself; `_gal` is +the Galerkin solver's variant of the same matrix. + ## Fields - - `amats`, `bmats`, `cmats` - the ideal A, B, C matrices. In a kinetic run these are NOT what the + - `A_spline`, `B_spline`, `C_spline` - the ideal A, B, C. In a kinetic run these are NOT what the ODE integrates; see [`KineticMatrices`](@ref). - - `dmats_prim`, `emats_prim` - pre-Schur-reduction geometric forms + - `D_spline_prim`, `E_spline_prim` - pre-Schur-reduction geometric forms (D = χ₁·(g23 + q·g33·m/n); E = (-χ₁/n)·(q'·χ₁·g33 - 2π·i·χ₁·g31·singfac + jθ·I)). Downstream kinetic FKG Schur complements consume these primitive forms; the alternate singular-layer path that would need kinetic-added overwrites of D and E is not implemented here (see Kinetic.jl). - - `fmats_lower` - F in factorized lower-triangular form, F = L·Lᴴ. - - `fmats_prim` - primitive F before the Schur complement (for kinetic). - - `fmats_gal` - reduced Hermitian F̄, un-factored, for the Galerkin solver. - - `jmats` - Jacobian Fourier band (2·mpert−1 conjugate-symmetric coefficients per surface), used + - `H_spline` - primitive Ḡ before the Schur complement (Ḡ = H − C†A⁻¹C), also consumed by the + kinetic reduction. + - `F_spline_lower` - reduced F̄ as its lower-triangular Cholesky factor L, F̄ = L·Lᴴ. The ideal + Euler-Lagrange kernel only ever solves against F̄, so the factor is stored instead of F̄ itself. + - `F_spline_prim` - primitive F before the Schur complement (for kinetic). + - `F_spline_gal` - reduced Hermitian F̄, un-factored, for the Galerkin solver, which needs F̄ + directly rather than through a solve. + - `K_spline`, `G_spline` - reduced K̄ and Ḡ of the Euler-Lagrange system + d/dψ(F̄·ξ′ + K̄·ξ) = K̄†·ξ′ + Ḡ·ξ (Eqs A5-A7). + - `J_spline` - Jacobian Fourier band (2·mpert−1 conjugate-symmetric coefficients per surface), used to assemble the power-normalization matrix N in Free.jl. """ @kwdef struct IdealMatrices{S<:CubicSeriesInterpolant} - amats::S - bmats::S - cmats::S - dmats_prim::S - emats_prim::S - hmats::S - fmats_lower::S - fmats_prim::S - fmats_gal::S - kmats::S - gmats::S - jmats::S + A_spline::S + B_spline::S + C_spline::S + D_spline_prim::S + E_spline_prim::S + H_spline::S + F_spline_lower::S + F_spline_prim::S + F_spline_gal::S + K_spline::S + G_spline::S + J_spline::S end """ KineticMatrices Kinetic-MHD matrix ψ-splines assembled by [`make_kinetic_matrix`](@ref). Present on a -[`FourFitVars`](@ref) only for kinetic runs; its presence means the solution obeys the FKG ODE +[`MatrixSplines`](@ref) only for kinetic runs; its presence means the solution obeys the FKG ODE rather than the ideal Euler-Lagrange relation. +The FKG system is non-Hermitian, so the adjoint side of each reduced block is a distinct matrix +rather than a transpose of the forward one. The `_adj` suffix marks those adjoint-side counterparts +(Fortran's `aat` arrays); they must be evaluated, never derived from their unsuffixed partner. + ## Fields - - `amats`, `bmats`, `cmats` - the kinetic-modified A, B, C consumed by `sing_der!`. Unlike their + - `A_spline`, `B_spline`, `C_spline` - the kinetic-modified A, B, C consumed by `sing_der!`. Unlike their ideal counterparts, A here is non-Hermitian and requires an LU rather than a Cholesky. - - `kwmats`, `ktmats` - kinetic energy (W) and torque (T) matrices, 6 components each - (A, B, C, D, E, H perturbations). - - `f0mats`, `pmats`, `paats`, `kkmats`, `kkaats`, `r1mats`, `r2mats`, `r3mats`, `gaats` - - pre-computed FKG Schur-complement reductions (Logan 2015 Appendix C, Eqs C.1-C.11). + - `Kw_spline`, `Kt_spline` - kinetic energy (W) and torque (T) matrices, 6 components each + (A, B, C, D, E, H perturbations), indexed in that order. + - `F0_spline` - reduced F̄ of the kinetic system, F₀ = F_prim − D†A_kin⁻¹D. + - `P_spline`, `P_spline_adj` - P = (iD)†A_kin⁻¹B_k and its adjoint-side counterpart. + - `Kk_spline`, `Kk_spline_adj` - kinetic K̄ = E − (iD)†A_kin⁻¹C_kin and its adjoint-side counterpart. + - `R1_spline`, `R2_spline`, `R3_spline` - blocks built from the D- and E-perturbation kinetic + components (`Kw_spline[4:5]`, `Kt_spline[4:5]`), which have no ideal counterpart. R3 is the + adjoint-side partner of R2. + - `G_spline_adj` - adjoint-side Ḡ, Ḡ_adj = H_kin − C_adj†A_kin⁻¹C_kin. + +All nine reduced blocks are pre-computed per ψ by the FKG Schur complement +(Logan 2015 Appendix C, Eqs C.1-C.11) so `sing_der!` can assemble the ODE right-hand side with +explicit (m−nq) factors instead of 1/(m−nq). """ @kwdef struct KineticMatrices{S<:CubicSeriesInterpolant} - amats::S - bmats::S - cmats::S - kwmats::Vector{S} - ktmats::Vector{S} - f0mats::S - pmats::S - paats::S - kkmats::S - kkaats::S - r1mats::S - r2mats::S - r3mats::S - gaats::S + A_spline::S + B_spline::S + C_spline::S + Kw_spline::Vector{S} + Kt_spline::Vector{S} + F0_spline::S + P_spline::S + P_spline_adj::S + Kk_spline::S + Kk_spline_adj::S + R1_spline::S + R2_spline::S + R3_spline::S + G_spline_adj::S end """ - FourFitVars + MatrixSplines Fourier-fitted stability matrices for one run. `kinetic === nothing` marks an ideal run; otherwise the kinetic matrices are authoritative for the ODE and the ideal ones remain available for output @@ -83,7 +106,7 @@ and for the Galerkin/vacuum paths that are defined only in the ideal basis. - `_hint::Base.RefValue{Int}` - shared bracket-search hint for sequential (single-threaded) evaluation. Not thread-safe: parallel paths must pass their own hint. """ -@kwdef struct FourFitVars{S<:CubicSeriesInterpolant,Opts<:NamedTuple} +@kwdef struct MatrixSplines{S<:CubicSeriesInterpolant,Opts<:NamedTuple} numpert_total::Int itp_opts::Opts = (; extrap=ExtendExtrap()) ideal::IdealMatrices{S} @@ -410,10 +433,10 @@ end """ - make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> FourFitVars + make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> MatrixSplines Constructs main ForceFreeStates matrices for a given toroidal mode number and returns -them as a new `FourFitVars` object. See the appendix of the Glasser Phys. Plasmas 2016 112506 +them as a new `MatrixSplines` object. See the appendix of the Glasser Phys. Plasmas 2016 112506 DCON paper for details on the matrix definitions. Performs the same function as `fourfit_make_matrix` in the Fortran code, except F, G, and K are stored as dense matrices. The matrix F is stored in factorized form with the lower triangle only, @@ -427,7 +450,7 @@ later (i.e. `sing_der!`). ### Returns - - `ffit::FourFitVars`: A struct holding cubic spline fits of the assembled matrices + - `mats::MatrixSplines`: A struct holding cubic spline fits of the assembled matrices ### TODOs @@ -597,18 +620,18 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates # TODO: set powers. Do we need this yet? Only called if power_flag = true itp_opts = (; extrap=ExtendExtrap()) ideal = IdealMatrices(; - amats=cubic_interp(metric.xs, Series(amats_flat); itp_opts...), - bmats=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...), - cmats=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...), - dmats_prim=cubic_interp(metric.xs, Series(dmats_flat); itp_opts...), - emats_prim=cubic_interp(metric.xs, Series(emats_flat); itp_opts...), - hmats=cubic_interp(metric.xs, Series(hmats_flat); itp_opts...), - fmats_lower=cubic_interp(metric.xs, Series(fmats_lower_flat); itp_opts...), - fmats_prim=cubic_interp(metric.xs, Series(fmats_prim_flat); itp_opts...), - fmats_gal=cubic_interp(metric.xs, Series(fmats_gal_flat); itp_opts...), - gmats=cubic_interp(metric.xs, Series(gmats_flat); itp_opts...), - kmats=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), + A_spline=cubic_interp(metric.xs, Series(amats_flat); itp_opts...), + B_spline=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...), + C_spline=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...), + D_spline_prim=cubic_interp(metric.xs, Series(dmats_flat); itp_opts...), + E_spline_prim=cubic_interp(metric.xs, Series(emats_flat); itp_opts...), + H_spline=cubic_interp(metric.xs, Series(hmats_flat); itp_opts...), + F_spline_lower=cubic_interp(metric.xs, Series(fmats_lower_flat); itp_opts...), + F_spline_prim=cubic_interp(metric.xs, Series(fmats_prim_flat); itp_opts...), + F_spline_gal=cubic_interp(metric.xs, Series(fmats_gal_flat); itp_opts...), + G_spline=cubic_interp(metric.xs, Series(gmats_flat); itp_opts...), + K_spline=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl - jmats=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) - return FourFitVars(; numpert_total=intr.numpert_total, itp_opts, ideal) + J_spline=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) + return MatrixSplines(; numpert_total=intr.numpert_total, itp_opts, ideal) end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 93147c280..e3bb3fdd8 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -36,7 +36,7 @@ end power_norm_matrix!(Nmat, jmat, mpert, npert, dV_dpsi) -> Nmat Assemble the power-normalization (surface-norm) matrix N from the conjugate-symmetric Jacobian -Fourier band `jmat` (length 2·mpert−1, evaluated from the `ffit.ideal.jmats` spline), such that +Fourier band `jmat` (length 2·mpert−1, evaluated from the `mats.ideal.J_spline` spline), such that ξ†·N·ξ = ∮ J |ξ(θ)|² dθ / (dV/dψ) = ⟨|ξ|²⟩ @@ -111,7 +111,7 @@ function compute_scaled_wv(ctrl::ForceFreeStatesControl, equil::Equilibrium.Plas end """ - free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) -> FreeBoundaryResult + free_run(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, 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. @@ -119,7 +119,7 @@ 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, mats::MatrixSplines, intr::ForceFreeStatesInternal) # Initializations and allocations (; mpert, numpert_total, psilim, npert) = intr @@ -140,7 +140,7 @@ calculations and data dumping. # The Jacobian band is evaluated at psilim (same surface as W), not at the last grid surface. Nmat = zeros!(pool, ComplexF64, numpert_total, numpert_total) jmat_edge = zeros!(pool, ComplexF64, 2 * mpert - 1) - ffit.ideal.jmats(jmat_edge, psilim; hint=ffit._hint) + mats.ideal.J_spline(jmat_edge, psilim; hint=mats._hint) power_norm_matrix!(Nmat, jmat_edge, mpert, npert, dV_dpsi) # Least stable eigenvalue of the vacuum matrix alone, power-normalized via the pencil @@ -252,7 +252,7 @@ q-window minimum. end """ - free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, odet::OdeState) -> ComplexF64 + free_compute_total(equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, odet::OdeState) -> ComplexF64 Compute total complex energy eigenvalue (total1). This is a trimmed down version of `free_run` that only computes the total energy eigenvalue for the mode unstable mode, used in `findmax_dW_edge!` @@ -260,7 +260,7 @@ which calls this function at each step in the psiedge -> psilim region of integr the same function as `free_test` in the Fortran code, except we have moved the creation of the wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvmat` (a complex-valued spline). """ -@with_pool pool function free_compute_total(equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, odet::OdeState) +@with_pool pool function free_compute_total(equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, odet::OdeState) Npert = intr.numpert_total wp = zeros!(pool, ComplexF64, Npert, Npert) @@ -290,7 +290,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma # Local power-normalization matrix N(ψ) from the Jacobian Fourier band spline, so the # power quotient uses the same surface as W (see power_norm_matrix!) - ffit.ideal.jmats(jmat_local, odet.psifac; hint=ffit._hint) + mats.ideal.J_spline(jmat_local, odet.psifac; hint=mats._hint) 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 diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl index 1ef7a3a7d..1ad84cb8b 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl @@ -59,14 +59,14 @@ function gal_hermite(x::Real, x0::Real, x1::Real) end """ - gal_get_fkg(ffit, intr, x, q) -> (F, K, G) + gal_get_fkg(mats, intr, x, q) -> (F, K, G) Evaluate the `mpert×mpert` matrices `F = Q F̄ Qᴴ`, `K = Q K̄`, `G = Ḡ` at flux `x` with safety factor `q`, where `Q = diag(singfac)` and `singfac = m - n q` (direct). Port of `gal_get_fkg` (gal.f). -Uses the un-factored reduced `ffit.ideal.fmats_gal` (F̄), `ffit.ideal.kmats` (K̄), `ffit.ideal.gmats` (Ḡ). F/K/G are the +Uses the un-factored reduced `mats.ideal.F_spline_gal` (F̄), `mats.ideal.K_spline` (K̄), `mats.ideal.G_spline` (Ḡ). F/K/G are the ideal-MHD Euler–Lagrange coefficient matrices of the outer-region weak form (Glasser 2016, PoP 23, 112506). """ -function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float64, q::Float64; +function gal_get_fkg(mats::MatrixSplines, intr::ForceFreeStatesInternal, x::Float64, q::Float64; Fbuf::Union{Nothing,Matrix{ComplexF64}}=nothing, Kbuf::Union{Nothing,Matrix{ComplexF64}}=nothing, Gbuf::Union{Nothing,Matrix{ComplexF64}}=nothing) N = intr.numpert_total @@ -75,9 +75,9 @@ function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float6 F = Fbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Fbuf K = Kbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Kbuf G = Gbuf === nothing ? Matrix{ComplexF64}(undef, N, N) : Gbuf - ffit.ideal.fmats_gal(vec(F), x; hint=ffit._hint) - ffit.ideal.kmats(vec(K), x; hint=ffit._hint) - ffit.ideal.gmats(vec(G), x; hint=ffit._hint) + mats.ideal.F_spline_gal(vec(F), x; hint=mats._hint) + mats.ideal.K_spline(vec(K), x; hint=mats._hint) + mats.ideal.G_spline(vec(G), x; hint=mats._hint) # scale F̄→F=Q F̄ Qᴴ and K̄→K=Q K̄ in place (Q = diag(sf)) @inbounds for j in 1:N, i in 1:N @@ -88,14 +88,14 @@ function gal_get_fkg(ffit::FourFitVars, intr::ForceFreeStatesInternal, x::Float6 end """ - gal_gauss_quad!(cell, ffit, profiles, intr, nodes, weights, swap_edge) + gal_gauss_quad!(cell, mats, profiles, intr, nodes, weights, swap_edge) Accumulate the nonresonant Hermite stiffness block `cell.mat` by Gauss-Lobatto quadrature. Port of `gal_gauss_quad` (gal.f). When `swap_edge` (the final cell of the last interval), the right-node value/slope DOFs (Fortran pb(2)↔pb(3)) are swapped so the edge-value DOF lands at index 3, matching the free-boundary BC (gal.f). """ -function gal_gauss_quad!(cell::GalCell, ffit::FourFitVars, profiles, intr::ForceFreeStatesInternal, +function gal_gauss_quad!(cell::GalCell, mats::MatrixSplines, profiles, intr::ForceFreeStatesInternal, nodes::Vector{Float64}, weights::Vector{Float64}, swap_edge::Bool) N = intr.numpert_total @@ -113,7 +113,7 @@ function gal_gauss_quad!(cell::GalCell, ffit::FourFitVars, profiles, intr::Force x = x0c + dxc * nodes[iq] w = dxc * weights[iq] q = profiles.q_spline(x; hint=qhint) - F, K, G = gal_get_fkg(ffit, intr, x, q; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) + F, K, G = gal_get_fkg(mats, intr, x, q; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) pbt, qbt = gal_hermite(x, x1, x2) # swap the right-node value/slope DOFs (Fortran pb(2)↔pb(3)) for the free-boundary edge pb = swap_edge ? (pbt[1], pbt[2], pbt[4], pbt[3]) : pbt @@ -135,13 +135,13 @@ end end """ - gal_extension!(cell, ising, ffit, profiles, intr, asymps, sings, nn, nodes, weights) + gal_extension!(cell, ising, mats, profiles, intr, asymps, sings, nn, nodes, weights) Build `cell.emat`, `cell.ediag`, `cell.rhs`, `cell.erhs` for an extension cell (`ext`/`ext1`/`ext2`). Port of `gal_extension` (gal.f). The big solution (column 1 of the gal asymptotic) drives the RHS; for `ext` cells the small solution (column 2) also forms the resonant coupling `emat`/`ediag`. """ -function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, +function gal_extension!(cell::GalCell, ising::Int, mats::MatrixSplines, profiles, intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, nn::Int, nodes::Vector{Float64}, weights::Vector{Float64}) @@ -198,7 +198,7 @@ function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, end # surface term at xb qb = profiles.q_spline(xb; hint=qhint) - Fb, Kb, _ = gal_get_fkg(ffit, intr, xb, qb; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) + Fb, Kb, _ = gal_get_fkg(mats, intr, xb, qb; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) pbt, _ = gal_hermite(xb, x1, x2) Fdu_Ku = Fb * du .+ Kb * u for ip in 0:np @@ -222,7 +222,7 @@ function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, x = x0c + dxc * nodes[iq] w = dxc * weights[iq] q = profiles.q_spline(x; hint=qhint) - F, K, G = gal_get_fkg(ffit, intr, x, q; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) + F, K, G = gal_get_fkg(mats, intr, x, q; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) pbt, qbt = gal_hermite(x, x1, x2) uax = ua_at(x) ub = uax[:, 1, 1] @@ -250,7 +250,7 @@ function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, # --- surface terms (always; gal.f) --- q_l = profiles.q_spline(x1; hint=qhint) - Fl, Kl, _ = gal_get_fkg(ffit, intr, x1, q_l; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) + Fl, Kl, _ = gal_get_fkg(mats, intr, x1, q_l; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) pbt, _ = gal_hermite(x1, x1, x2) surf_l = Fl * du1 .+ Kl * ua1 for ip in 0:np @@ -261,7 +261,7 @@ function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, end q_r = profiles.q_spline(x2; hint=qhint) - Fr, Kr, _ = gal_get_fkg(ffit, intr, x2, q_r; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) + Fr, Kr, _ = gal_get_fkg(mats, intr, x2, q_r; Fbuf=Fbuf, Kbuf=Kbuf, Gbuf=Gbuf) pbt, _ = gal_hermite(x2, x1, x2) surf_r = Fr * du2 .+ Kr * ua2 for ip in 0:np @@ -274,13 +274,13 @@ function gal_extension!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, end """ - gal_resonant!(cell, ising, ffit, profiles, intr, asymps, sings, nn, gal_tol, gal_gnstep, verbose) + gal_resonant!(cell, ising, mats, profiles, intr, asymps, sings, nn, gal_tol, gal_gnstep, verbose) Compute the resonant-cell contributions `cell.erhs`, `cell.ediag`, `cell.rhs`, `cell.emat` by adaptive quadrature (QuadGK) of the residual-operator integrand. Port of `gal_lsode_int`/`gal_lsode_der` (gal.f): the Fortran LSODE accumulation is replaced by `quadgk` over `[x_bdy, x_lsode]`. """ -function gal_resonant!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, +function gal_resonant!(cell::GalCell, ising::Int, mats::MatrixSplines, profiles, intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, nn::Int, gal_tol::Float64, gal_gnstep::Int, verbose::Bool) @@ -317,7 +317,7 @@ function gal_resonant!(cell::GalCell, ising::Int, ffit::FourFitVars, profiles, sing_get_ua_gal!(ua_buf, asymp, z) sing_get_dua_gal!(dua_buf, asymp, z) q = profiles.q_spline(x; hint=qhint) - sing_matvec!(mv, kmat, gmat, d1, mvtmp, sfvec, ffit, intr, x, q, ua_buf, dua_buf) # N×2 (col1=big, col2=small) + sing_matvec!(mv, kmat, gmat, d1, mvtmp, sfvec, mats, intr, x, q, ua_buf, dua_buf) # N×2 (col1=big, col2=small) pbt, _ = gal_hermite(x, x1, x2) # Pairwise sum over the materialized product (reuses buffers) — matches sum(w .* mv) bit-for-bit. w .= conj.(@view ua_buf[:, 2, 1]) # conj small solution, qty1 diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index 66dddd893..963c107f0 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -6,13 +6,13 @@ # (gal.f). The DRIVEN/RPEC inner-layer matching is wired in via gal_match_rpec (GalerkinMatch.jl). """ - gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge) + gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) Assemble the global banded matrix and RHS. For each cell: Gauss-Lobatto Hermite stiffness (`gal_gauss_quad!`), then the resonant (`gal_resonant!`) or extension (`gal_extension!`) contributions; then the boundary conditions and the scatter into `ws.mat`/`ws.rhs`. Port of `gal_make_arrays`. """ -function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, +function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines, intr::ForceFreeStatesInternal, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, nn::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}}) @@ -24,11 +24,11 @@ function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, for ising in 0:msing for (ix, cell) in enumerate(ws.intvl[ising+1].cells) swap_edge = (ising == msing && ix == ws.nx) - gal_gauss_quad!(cell, ffit, profiles, intr, nodes, weights, swap_edge) + gal_gauss_quad!(cell, mats, profiles, intr, nodes, weights, swap_edge) if cell.etype == GCT_RES - gal_resonant!(cell, ising, ffit, profiles, intr, asymps, sings, nn, ctrl.gal_tol, ctrl.gal_gnstep, ctrl.verbose) + gal_resonant!(cell, ising, mats, profiles, intr, asymps, sings, nn, ctrl.gal_tol, ctrl.gal_gnstep, ctrl.verbose) elseif cell.etype == GCT_EXT || cell.etype == GCT_EXT1 || cell.etype == GCT_EXT2 - gal_extension!(cell, ising, ffit, profiles, intr, asymps, sings, nn, nodes, weights) + gal_extension!(cell, ising, mats, profiles, intr, asymps, sings, nn, nodes, weights) end end end @@ -45,7 +45,7 @@ function empty_galerkin_result() end """ - galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, + galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines, intr::ForceFreeStatesInternal; wv=nothing) -> (GalerkinResult, Union{Nothing,DeltaPrimeData}) Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of `gal_solve` @@ -56,7 +56,7 @@ resonant surfaces in the domain it returns an empty result and `nothing`. `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, +function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines, intr::ForceFreeStatesInternal; wv=nothing) intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported") @@ -81,15 +81,15 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, asymps = GalSingAsymp[] for s in sings sing_order = ctrl.gal_sing_order - ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order) + ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order) if ctrl.gal_sing_order_ceiling order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1])) if order > ctrl.gal_sing_order sing_order = order - ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order) + ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order) end end - al = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order) + al = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order) push!(asymps, GalSingAsymp(ar, al)) end @@ -128,7 +128,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, wv_edge = Matrix{ComplexF64}(wv .* equil.psio^2) end - gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge) + gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) if ctrl.verbose offdbg = ws.solver == "LU" ? ws.kl + ws.ku + 1 : 1 diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 7190e1917..f2d281741 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,10 +1,10 @@ """ - make_kinetic_matrix(ctrl, equil, ffit, intr, metric; + make_kinetic_matrix(ctrl, equil, mats, intr, metric; calculated_source=nothing) Construct kinetic energy (W) and torque (T) matrices and pre-compute the FKG derived -matrices used by `sing_der!`, returning a new `FourFitVars` carrying both alongside the -ideal matrices of the input `ffit`. +matrices used by `sing_der!`, returning a new `MatrixSplines` carrying both alongside the +ideal matrices of the input `mats`. Dispatches on `ctrl.kinetic_source`: @@ -13,7 +13,7 @@ Dispatches on `ctrl.kinetic_source`: - `"calculated"`: Compute via the `calculated_source` callback. This is expected to be `KineticForces.compute_calculated_kinetic_matrices` injected by `GeneralizedPerturbedEquilibrium.main`. The callback receives - `(ctrl, equil, intr, metric, ffit)` and returns `(kw_flat, kt_flat)` of + `(ctrl, equil, intr, metric, mats)` and returns `(kw_flat, kt_flat)` of shape `(mpsi, np^2, 6)`. Callback injection is used because ForceFreeStates is loaded before KineticForces, so a direct import would invert the dependency order. @@ -24,7 +24,7 @@ reduction. function make_kinetic_matrix( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal, metric::MetricData; calculated_source::Union{Nothing,Function}=nothing @@ -34,7 +34,7 @@ function make_kinetic_matrix( # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" - kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, mpsi, ctrl.kinetic_factor, intr.mlow, ffit, xs) + kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs) elseif ctrl.kinetic_source == "calculated" isnothing(calculated_source) && error( "kinetic_source=\"calculated\" requires the KineticForces callback. " * @@ -42,7 +42,7 @@ function make_kinetic_matrix( "calling make_kinetic_matrix directly, or pass " * "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly." ) - kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) kw_flat .*= ctrl.kinetic_factor kt_flat .*= ctrl.kinetic_factor else @@ -50,18 +50,19 @@ function make_kinetic_matrix( end # Build splines for each of the 6 components - kwmats = [cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); ffit.itp_opts...) for ic in 1:6] - ktmats = [cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); ffit.itp_opts...) for ic in 1:6] + Kw_spline = [cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); mats.itp_opts...) for ic in 1:6] + Kt_spline = [cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); mats.itp_opts...) for ic in 1:6] # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - return _compute_fkg_matrices(ffit, equil, intr, metric, kw_flat, kt_flat, kwmats, ktmats) + return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat, Kw_spline, Kt_spline) end """ - _compute_fkg_matrices(ffit, equil, intr, metric, kw_flat, kt_flat, kwmats, ktmats) -> FourFitVars + _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat, Kw_spline, Kt_spline) -> MatrixSplines -Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and return a new `FourFitVars` -holding them, the kinetic-modified A/B/C, and the ideal A/B/C of `ffit` preserved as `*_ideal`. +Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and return a new `MatrixSplines` +whose `kinetic` field holds them alongside the kinetic-modified A/B/C; `mats.ideal` is carried over +unchanged. This corresponds to `fourfit_kinetic_matrix` method=0 in the Fortran code (Fortran `fourfit.F` lines 1170-1260). The 9 matrices computed are the Schur complement reductions of ideal (A,B,C,D,E,H) and kinetic (W,T) @@ -71,21 +72,21 @@ can be assembled with explicit (m-nq) factors rather than 1/(m-nq), avoiding num rational surfaces. """ function _compute_fkg_matrices( - ffit::FourFitVars, + mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, kt_flat::Array{ComplexF64,3}, - kwmats::Vector, - ktmats::Vector + Kw_spline::Vector, + Kt_spline::Vector ) xs = metric.xs mpsi = length(xs) np = intr.numpert_total mpert = intr.mpert npert = intr.npert - ideal = ffit.ideal + ideal = mats.ideal # Allocate output arrays — kinetic-modified A/B/C stored for sing_der! FKG path ak_flat = zeros(ComplexF64, mpsi, np^2) @@ -111,13 +112,13 @@ function _compute_fkg_matrices( psi = xs[ipsi] # Evaluate ideal and kinetic matrices from splines (full np×np, block-diagonal in n) - amat_full = reshape(ideal.amats(psi; hint=hint), np, np) - bmat_full = reshape(ideal.bmats(psi; hint=hint), np, np) - cmat_full = reshape(ideal.cmats(psi; hint=hint), np, np) - dmat_full = reshape(ideal.dmats_prim(psi; hint=hint), np, np) - emat_full = reshape(ideal.emats_prim(psi; hint=hint), np, np) - hmat_full = reshape(ideal.hmats(psi; hint=hint), np, np) - fmat_prim_full = reshape(ideal.fmats_prim(psi; hint=hint), np, np) + amat_full = reshape(ideal.A_spline(psi; hint=hint), np, np) + bmat_full = reshape(ideal.B_spline(psi; hint=hint), np, np) + cmat_full = reshape(ideal.C_spline(psi; hint=hint), np, np) + dmat_full = reshape(ideal.D_spline_prim(psi; hint=hint), np, np) + emat_full = reshape(ideal.E_spline_prim(psi; hint=hint), np, np) + hmat_full = reshape(ideal.H_spline(psi; hint=hint), np, np) + fmat_prim_full = reshape(ideal.F_spline_prim(psi; hint=hint), np, np) kwmat_full = zeros(ComplexF64, np, np, 6) ktmat_full = zeros(ComplexF64, np, np, 6) @@ -238,23 +239,23 @@ function _compute_fkg_matrices( end end - itp_opts = ffit.itp_opts + itp_opts = mats.itp_opts kinetic = KineticMatrices(; # kinetic-modified A/B/C consumed by sing_der!; A is non-Hermitian here - amats=cubic_interp(xs, Series(ak_flat); itp_opts...), - bmats=cubic_interp(xs, Series(bk_flat); itp_opts...), - cmats=cubic_interp(xs, Series(ck_flat); itp_opts...), - kwmats, - ktmats, - f0mats=cubic_interp(xs, Series(f0_flat); itp_opts...), - pmats=cubic_interp(xs, Series(p_flat); itp_opts...), - paats=cubic_interp(xs, Series(pa_flat); itp_opts...), - kkmats=cubic_interp(xs, Series(kk_flat); itp_opts...), - kkaats=cubic_interp(xs, Series(kka_flat); itp_opts...), - r1mats=cubic_interp(xs, Series(r1_flat); itp_opts...), - r2mats=cubic_interp(xs, Series(r2_flat); itp_opts...), - r3mats=cubic_interp(xs, Series(r3_flat); itp_opts...), - gaats=cubic_interp(xs, Series(ga_flat); itp_opts...)) - - return FourFitVars(ffit.numpert_total, itp_opts, ffit.ideal, kinetic, ffit._hint) + A_spline=cubic_interp(xs, Series(ak_flat); itp_opts...), + B_spline=cubic_interp(xs, Series(bk_flat); itp_opts...), + C_spline=cubic_interp(xs, Series(ck_flat); itp_opts...), + Kw_spline, + Kt_spline, + F0_spline=cubic_interp(xs, Series(f0_flat); itp_opts...), + P_spline=cubic_interp(xs, Series(p_flat); itp_opts...), + P_spline_adj=cubic_interp(xs, Series(pa_flat); itp_opts...), + Kk_spline=cubic_interp(xs, Series(kk_flat); itp_opts...), + Kk_spline_adj=cubic_interp(xs, Series(kka_flat); itp_opts...), + R1_spline=cubic_interp(xs, Series(r1_flat); itp_opts...), + R2_spline=cubic_interp(xs, Series(r2_flat); itp_opts...), + R3_spline=cubic_interp(xs, Series(r3_flat); itp_opts...), + G_spline_adj=cubic_interp(xs, Series(ga_flat); itp_opts...)) + + return MatrixSplines(mats.numpert_total, itp_opts, mats.ideal, kinetic, mats._hint) end diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index 10fbf2302..da8fe3648 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -54,7 +54,7 @@ so bpen and closure are always present. - `dir_path::String` - Working directory of the run. - `wall_settings::Vacuum.WallShapeSettings` - Wall shape used by the vacuum calculation. - `debug_settings::DebugSettings` - Diagnostic dump settings (the `[DEBUG]` deck section / `debug=` API keyword). - - `metric::MetricData`, `ffit::FourFitVars` - Metric data and Euler-Lagrange matrix interpolants. + - `metric::MetricData`, `mats::MatrixSplines` - Metric data and Euler-Lagrange matrix interpolants. - `surfaces::Vector{SingType}` - Ideal singular surfaces in the integration domain, with asymptotic bases and GGJ coefficients. - `kinetic::NamedTuple` - Kinetic singular-surface scan (`kmsing`, `kinsing`, `scan_psi`, `scan_cond`, `scan_threshold`); empty unless the finder ran. - `closure::Symbol` - How the basis is closed at the rationals: `:ideal` (the ideal jump @@ -76,7 +76,7 @@ so bpen and closure are always present. diagnostics, including the RPEC inner-layer match when requested. Its Δ′ payload lives in `delta_prime`, not here. """ -struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:FourFitVars} <: ModeSpace +struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:MatrixSplines} <: ModeSpace integrator::Symbol control::ForceFreeStatesControl equil::E @@ -99,7 +99,7 @@ struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:FourFitVars} <: # Assembly products, always present. metric::MetricData - ffit::F + mats::F surfaces::Vector{SingType} kinetic::@NamedTuple{kmsing::Int, kinsing::Vector{SingType}, scan_psi::Vector{Float64}, scan_cond::Vector{Float64}, scan_threshold::Float64} @@ -147,7 +147,7 @@ end # analytic Galerkin derivative rather than a differenced value spline, and # Ξ_s = −A⁻¹(B·Ξ′ + C·Ξ) is the same outer ideal-MHD relation `sing_der!` uses. The grid runs # inner→edge, so the last node is the control surface and carries the edge boundary condition. -function _matched_gal_profiles(gal_result::GalerkinResult, ffit::FourFitVars, intr::ModeSpace) +function _matched_gal_profiles(gal_result::GalerkinResult, mats::MatrixSplines, intr::ModeSpace) sol = gal_result.solution m = gal_result.match npert = intr.numpert_total @@ -171,14 +171,14 @@ function _matched_gal_profiles(gal_result::GalerkinResult, ffit::FourFitVars, in ξ′ = @view dxi_f[:, ip, :] @views u_store[:, :, 1, ip] .= ξ @views du_store[:, :, ip] .= ξ′ - @views compute_node_xi_s!(xi_s_store[:, :, ip], ξ′, ξ, ffit, psi_f[ip]; hint=hint) + @views compute_node_xi_s!(xi_s_store[:, :, ip], ξ′, ξ, mats, psi_f[ip]; hint=hint) end return SolutionProfiles(:gal_native, ngrid_f, psi_f, q_f, u_store, du_store, xi_s_store) end """ - build_result(integrator, ctrl, equil, intr, metric, ffit, odet, free_energies, gal_data, gal_dp) + build_result(integrator, ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) -> ForceFreeStatesResult Assemble the published result once the solve is finished — the one place that decides what a @@ -197,7 +197,7 @@ function build_result( equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData, - ffit::FourFitVars, + mats::MatrixSplines, odet::Union{Nothing,OdeState}, free_energies::Union{Nothing,FreeBoundaryResult}, gal_data::Union{Nothing,GalerkinResult}, @@ -208,11 +208,11 @@ function build_result( # The forward sweep is the only formalism whose stores need materializing; doing it here # keeps `SolutionProfiles.du_store`/`xi_s_store` populated by construction. solution = if integrator === :forward && odet !== nothing - materialize_derivative_stores!(odet, equil, ffit, intr) + materialize_derivative_stores!(odet, equil, mats, intr) SolutionProfiles(:el_axis, odet.step, odet.psi_store, odet.q_store, odet.u_store, odet.du_store, odet.xi_s_store) elseif matched - _matched_gal_profiles(gal_data, ffit, intr) + _matched_gal_profiles(gal_data, mats, intr) else nothing end @@ -249,7 +249,7 @@ function build_result( integrator, ctrl, equil, intr.mlow, intr.mhigh, intr.mpert, intr.nlow, intr.nhigh, intr.npert, intr.numpert_total, intr.psilow, intr.psilim, intr.qlim, intr.q1lim, intr.dir_path, intr.wall_settings, intr.debug_settings, - metric, ffit, intr.sing, kinetic, + metric, mats, intr.sing, kinetic, closure, bpen, solution, odet, wp, free_energies, delta_prime, gal_data ) diff --git a/src/ForceFreeStates/Riccati/Crossings.jl b/src/ForceFreeStates/Riccati/Crossings.jl index 017792b1c..a11b34f16 100644 --- a/src/ForceFreeStates/Riccati/Crossings.jl +++ b/src/ForceFreeStates/Riccati/Crossings.jl @@ -1,7 +1,7 @@ # Singular-surface crossing algorithms for the Riccati/fundamental-matrix integration. """ - riccati_cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr, ising) + riccati_cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, ising) Cross a singular surface for the Riccati formulation. Replaces `cross_ideal_singular_surf!` for the Riccati integration path with two key differences: @@ -30,18 +30,18 @@ The u_store entry at the crossing step correctly stores (U₁_new, U₂_new) so """ function riccati_cross_ideal_singular_surf!( odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, intr::ForceFreeStatesInternal, ising::Int + mats::MatrixSplines, intr::ForceFreeStatesInternal, ising::Int ) # Skip Gaussian reduction — S is bounded so no large-norm columns exist. singp = intr.sing[ising] dpsi = singp.psifac - odet.psifac # ψ_res - ψ_current (positive) ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert - sing_asymp_left, sing_asymp_right = _two_sided_singular_asymptotics(singp, ctrl, equil, ffit, intr) + sing_asymp_left, sing_asymp_right = _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) _log_riccati_crossing_diagnostics(odet, intr, ising, singp, dpsi, sing_asymp_left, sing_asymp_right) _capture_left_crossing_data!(odet, singp, sing_asymp_left, dpsi, intr, ising) - _predict_across_singular_surface!(odet, ctrl, equil, ffit, intr, ising, ipert_res, dpsi, sing_asymp_right) + _predict_across_singular_surface!(odet, ctrl, equil, mats, intr, ising, ipert_res, dpsi, sing_asymp_right) _capture_right_crossing_data!(odet, singp, sing_asymp_right, dpsi, intr, ising, ipert_res, ctrl) _stash_per_surface_delta_prime_stub!(odet, intr, ising, ipert_res, sing_asymp_right, equil, ctrl) @@ -52,17 +52,17 @@ function riccati_cross_ideal_singular_surf!( end """ - _two_sided_singular_asymptotics(singp, ctrl, equil, ffit, intr) -> (left, right) + _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) -> (left, right) Compute left- (`sig=-1`) and right- (`sig=+1`) side singular asymptotics matching Fortran STRIDE's separate vmatl/vmatr (sing_vmat). Alpha is taken from the right side and shared with the left. """ function _two_sided_singular_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) - sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, ffit, intr; sig=1.0) - sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, ffit, intr; sig=-1.0, + sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) + sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha) return sing_asymp_left, sing_asymp_right end @@ -95,7 +95,7 @@ end # odet.psifac to the right side. The zeroed columns stay zero through the predictor # since du[:, ipert_res, :] = 0 when u[:, ipert_res, :] = 0. function _predict_across_singular_surface!(odet::OdeState, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, ising::Int, ipert_res, dpsi::Float64, sing_asymp_right) if ctrl.kinetic_factor == 0 @@ -103,7 +103,7 @@ function _predict_across_singular_surface!(odet::OdeState, ctrl::ForceFreeStates odet.u[:, ipert_res[i], :] .= 0 end end - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) + params = (ctrl, equil, mats, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) sing_der!(du1, odet.u, params, odet.psifac) diff --git a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl index 5cf75446d..0a955fde7 100644 --- a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl +++ b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl @@ -1,7 +1,7 @@ # STRIDE global boundary-value problem: Delta-prime matrix assembly, solve, PEST-3 decomposition. """ - compute_delta_prime_matrix!(intr, propagators, chunks; wv, psio, debug, ctrl, equil, ffit) + compute_delta_prime_matrix!(intr, propagators, chunks; wv, psio, debug, ctrl, equil, mats) Compute the inter-surface tearing stability matrix (msing × msing) using the STRIDE global BVP formulation [Glasser 2018 Phys. Plasmas 25, 032501, Sec. III.B]. @@ -71,7 +71,7 @@ function compute_delta_prime_matrix!( S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}} = nothing, ctrl::Union{Nothing,ForceFreeStatesControl} = nothing, equil::Union{Nothing,Equilibrium.PlasmaEquilibrium} = nothing, - ffit::Union{Nothing,FourFitVars} = nothing + mats::Union{Nothing,MatrixSplines} = nothing ) intr.msing == 0 && return _has_unsupported_multi_resonance(intr) && return @@ -110,7 +110,7 @@ function compute_delta_prime_matrix!( if use_S_axis uShootR, uShootL, uAxis = _build_S_axis_shooting_propagators( propagators, chunks, i_crossings, sing, msing, N, - T_left_mats, T_right_mats, has_ua, ctrl, equil, ffit, intr, debug) + T_left_mats, T_right_mats, has_ua, ctrl, equil, mats, intr, debug) debug && _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_left, T_left_mats, ipert_all, has_ua, msing, N) @@ -245,9 +245,9 @@ function _build_S_axis_shooting_propagators( propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, sing::Vector{SingType}, msing::Int, N::Int, T_left_mats::Vector{Matrix{ComplexF64}}, T_right_mats::Vector{Matrix{ComplexF64}}, - has_ua::Bool, ctrl, equil, ffit, intr::ForceFreeStatesInternal, debug::Bool) + has_ua::Bool, ctrl, equil, mats, intr::ForceFreeStatesInternal, debug::Bool) - can_reintegrate = has_ua && ctrl !== nothing && equil !== nothing && ffit !== nothing + can_reintegrate = has_ua && ctrl !== nothing && equil !== nothing && mats !== nothing uShootR = Vector{Matrix{ComplexF64}}(undef, msing) uShootL = Vector{Matrix{ComplexF64}}(undef, msing) # uShootL[1] handled separately below @@ -261,7 +261,7 @@ function _build_S_axis_shooting_propagators( end if can_reintegrate && !isempty(shoot_range_R) uShootR[j] = integrate_fm_with_ua_ic(chunks, shoot_range_R, sing[j].ua_right, - ctrl, equil, ffit, intr; backward=false, psi_ua=sing[j].psi_ua_right) + ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right) else T_init = has_ua ? T_right_mats[j] : nothing uShootR[j] = assemble_fm_matrix(propagators, shoot_range_R; T_init=T_init) @@ -278,7 +278,7 @@ function _build_S_axis_shooting_propagators( end if can_reintegrate && !isempty(shoot_range_L) uShootL[j] = integrate_fm_with_ua_ic(chunks, shoot_range_L, sing[j].ua_left, - ctrl, equil, ffit, intr; backward=true, psi_ua=sing[j].psi_ua_left) + ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left) else T_init = has_ua ? T_left_mats[j] : nothing uShootL[j] = assemble_fm_matrix(propagators, shoot_range_L; T_init=T_init) @@ -288,7 +288,7 @@ function _build_S_axis_shooting_propagators( uAxis, i_axis_mid = _build_conditioned_axis_propagator(propagators, i_crossings, N) uShootL[1] = _build_uShootL_first(propagators, chunks, i_crossings, sing, T_left_mats, has_ua, can_reintegrate, i_axis_mid, - ctrl, equil, ffit, intr, N) + ctrl, equil, mats, intr, N) if debug shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1) @info " Axis propagator: $(i_axis_mid) chunks, cond=$(@sprintf("%.2e", cond(uAxis)))" @@ -358,11 +358,11 @@ function _build_uShootL_first(propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}}, has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int, - ctrl, equil, ffit, intr::ForceFreeStatesInternal, N::Int) + ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int) shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1) if can_reintegrate && !isempty(shoot_range_L1) return integrate_fm_with_ua_ic(chunks, shoot_range_L1, sing[1].ua_left, - ctrl, equil, ffit, intr; + ctrl, equil, mats, intr; backward=true, psi_ua=sing[1].psi_ua_left) elseif !isempty(shoot_range_L1) return assemble_fm_matrix(propagators, shoot_range_L1; diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl index 3defcdf24..18a23be2c 100644 --- a/src/ForceFreeStates/Riccati/Driver.jl +++ b/src/ForceFreeStates/Riccati/Driver.jl @@ -92,7 +92,7 @@ This is compatible with downstream code (which uses U₁/U₂ ratio): """ """ - riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) + riccati_eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, propagators, chunks, S_left) The Riccati/STRIDE integrator: a chunked fundamental matrix (propagator) driver for the EL integration. The trailing three return values feed the Δ' BVP in `compute_delta_prime_matrix!`; @@ -140,24 +140,24 @@ the DIIID-like example (N=26, n=1); bidirectional reduces this to within 2%. """ function riccati_eulerlagrange_integration( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, intr::ForceFreeStatesInternal + mats::MatrixSplines, intr::ForceFreeStatesInternal ) - odet = _initialize_parallel_odet(ctrl, equil, ffit, intr) + odet = _initialize_parallel_odet(ctrl, equil, mats, intr) chunks, propagators, odet_proxies = _setup_parallel_chunks_and_proxies(odet, ctrl, intr) _log_parallel_start(ctrl, odet, equil, chunks) - _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, ffit, intr, odet_proxies) + _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, mats, intr, odet_proxies) # Harvest solver-step counts accumulated thread-locally in each proxy during the BVP phase. # The outer re-integration below uses riccati_integrate_chunk!, which counts via its callback. odet.total_steps += sum(p.total_steps for p in odet_proxies) S_at_surface_left, last_crossing_step = - _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, ffit, intr) + _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, mats, intr) - _reintegrate_outer_plasma!(odet, last_crossing_step, ctrl, equil, ffit, intr) + _reintegrate_outer_plasma!(odet, last_crossing_step, ctrl, equil, mats, intr) - chunks, propagators = _handle_edge_dW_scan!(odet, chunks, propagators, ctrl, equil, ffit, intr) + chunks, propagators = _handle_edge_dW_scan!(odet, chunks, propagators, ctrl, equil, mats, intr) # 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, @@ -175,11 +175,11 @@ end # Build odet and initialize at the magnetic axis. Same path as serial eulerlagrange_integration. function _initialize_parallel_odet(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal) odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) if ctrl.sing_start <= 0 - initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr) + initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) elseif ctrl.sing_start <= intr.msing error("sing_start > 0 not implemented yet!") else @@ -223,11 +223,11 @@ end function _run_parallel_bvp_phase!(propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal, odet_proxies::Vector{OdeState}) Threads.@threads :static for i in eachindex(chunks) - integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, ffit, intr, + integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, mats, intr, odet_proxies[Threads.threadid()]) end end @@ -243,7 +243,7 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, intr::ForceFreeStatesInternal) + mats::MatrixSplines, intr::ForceFreeStatesInternal) N = intr.numpert_total S_at_surface_left = Matrix{ComplexF64}[] last_crossing_step = 1 @@ -265,7 +265,7 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu ctrl.kinetic_factor > 0 && error("kinetic_factor > 0 not implemented yet in Riccati!") # State is (S, I) from the renorm above — well-conditioned at the surface's left boundary. push!(S_at_surface_left, copy(odet.u[:, :, 1])) - riccati_cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr, chunk.ising) + riccati_cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, chunk.ising) last_crossing_step = odet.step - 1 else # Save non-crossing end-of-chunk state. These columns are FM/Riccati chunk @@ -291,7 +291,7 @@ end # before renormalization; we renorm here to (S_new, I) as the Riccati starting state. function _reintegrate_outer_plasma!(odet::OdeState, last_crossing_step::Int, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) N = intr.numpert_total odet.u .= odet.u_store[:, :, :, last_crossing_step] @@ -301,7 +301,7 @@ function _reintegrate_outer_plasma!(odet::OdeState, last_crossing_step::Int, renormalize_riccati_inplace!(odet.u, N) outer_chunk = IntegrationChunk(; psi_start=odet.psifac, psi_end=intr.psilim * (1 - eps), needs_crossing=false, ising=0) - riccati_integrate_chunk!(odet, ctrl, equil, ffit, intr, outer_chunk) + riccati_integrate_chunk!(odet, ctrl, equil, mats, intr, outer_chunk) # Post: odet.u is in (S, I) form; odet.step points to next empty slot. end @@ -316,7 +316,7 @@ end function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, propagators::Vector{ChunkPropagator}, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) N = intr.numpert_total odet.step -= 1 @@ -324,7 +324,7 @@ function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, ctrl.psiedge < intr.psilim || return chunks, propagators saved_psifac, saved_u = odet.psifac, copy(odet.u) - peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) + peak_step = findmax_dW_edge!(odet, ctrl, equil, mats, intr) if !ctrl.truncate_at_dW_peak odet.psifac = saved_psifac @@ -361,7 +361,7 @@ function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, chunks[last_chunk_idx] = new_chunk odet_proxy = OdeState(N, 1, 1, 0) integrate_propagator_chunk!(propagators[last_chunk_idx], new_chunk, - ctrl, equil, ffit, intr, odet_proxy) + ctrl, equil, mats, intr, odet_proxy) end n_dropped = 0 if last_chunk_idx < length(chunks) diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl index d6b43430a..ee638faed 100644 --- a/src/ForceFreeStates/Riccati/Propagators.jl +++ b/src/ForceFreeStates/Riccati/Propagators.jl @@ -140,11 +140,11 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form) du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, params::Tuple{ForceFreeStatesControl,Equilibrium.PlasmaEquilibrium, - FourFitVars,ForceFreeStatesInternal,OdeState,IntegrationChunk}, + MatrixSplines,ForceFreeStatesInternal,OdeState,IntegrationChunk}, psieval::Float64 ) - _, equil, ffit, intr, odet, _ = params + _, equil, mats, intr, odet, _ = params Npert = intr.numpert_total S = @view u[:, :, 1] @@ -167,9 +167,9 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form) tmp = similar!(pool, fmat_lower) # scratch # Evaluate F̄ (Cholesky factor), K̄, Ḡ splines at current ψ - ffit.ideal.fmats_lower(vec(fmat_lower), psieval; hint=ffit._hint) - ffit.ideal.kmats(vec(kmat), psieval; hint=ffit._hint) - ffit.ideal.gmats(vec(gmat), psieval; hint=ffit._hint) + mats.ideal.F_spline_lower(vec(fmat_lower), psieval; hint=mats._hint) + mats.ideal.K_spline(vec(kmat), psieval; hint=mats._hint) + mats.ideal.G_spline(vec(gmat), psieval; hint=mats._hint) # w = Q - K̄·S: w[i,j] = singfac_vec[i]·δ_ij - (K̄·S)[i,j] # Q is DIAGONAL (singfac_vec[i] only on i==j), so we cannot broadcast singfac_vec @@ -239,7 +239,7 @@ function riccati_integrator_callback!(integrator) end """ - riccati_integrate_chunk!(odet, ctrl, equil, ffit, intr, chunk) + riccati_integrate_chunk!(odet, ctrl, equil, mats, intr, chunk) Integrate the dual Riccati ODE from `chunk.psi_start` to `chunk.psi_end`. @@ -250,12 +250,12 @@ Ending state: u[:,:,1] = U₁, u[:,:,2] = U₂ (ratio S = U₁·U₂⁻¹ is the """ function riccati_integrate_chunk!( odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, intr::ForceFreeStatesInternal, chunk::IntegrationChunk + mats::MatrixSplines, intr::ForceFreeStatesInternal, chunk::IntegrationChunk ) cb = DiscreteCallback((u, t, integrator) -> true, riccati_integrator_callback!) rtol = ctrl.eulerlagrange_tolerance prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), - (ctrl, equil, ffit, intr, odet, chunk)) + (ctrl, equil, mats, intr, odet, chunk)) sol = solve(prob, Vern9(); reltol=rtol, callback=cb, save_everystep=false, save_end=true) odet.u .= sol.u[end] odet.psifac = sol.t[end] @@ -316,7 +316,7 @@ function renormalize_riccati_inplace!(u::Array{ComplexF64,3}, N::Int) end """ - integrate_propagator_chunk!(prop, chunk, ctrl, equil, ffit, intr, odet_proxy) + integrate_propagator_chunk!(prop, chunk, ctrl, equil, mats, intr, odet_proxy) Compute the fundamental matrix (propagator) for one integration chunk by solving the EL ODE twice from identity-block initial conditions. @@ -337,7 +337,7 @@ function integrate_propagator_chunk!( chunk::IntegrationChunk, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal, odet_proxy::OdeState ) @@ -349,7 +349,7 @@ function integrate_propagator_chunk!( (chunk.psi_start, chunk.psi_end) : (chunk.psi_end, chunk.psi_start) rtol = ctrl.eulerlagrange_tolerance - params = (ctrl, equil, ffit, intr, odet_proxy, chunk) + params = (ctrl, equil, mats, intr, odet_proxy, chunk) # Upper block IC: U₁ = I, U₂ = 0 u_upper = zeros(ComplexF64, N, N, 2) @@ -357,7 +357,7 @@ function integrate_propagator_chunk!( u_upper[i, i, 1] = 1 end odet_proxy.spline_hint[] = 1 - odet_proxy.ffit_hint[] = 1 + odet_proxy.mats_hint[] = 1 prob = ODEProblem(sing_der!, u_upper, tspan, params) sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) prop.block_upper_ic .= sol.u[end] @@ -369,7 +369,7 @@ function integrate_propagator_chunk!( u_lower[i, i, 2] = 1 end odet_proxy.spline_hint[] = 1 - odet_proxy.ffit_hint[] = 1 + odet_proxy.mats_hint[] = 1 prob = ODEProblem(sing_der!, u_lower, tspan, params) sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) prop.block_lower_ic .= sol.u[end] @@ -377,7 +377,7 @@ function integrate_propagator_chunk!( end """ - integrate_fm_with_ua_ic(chunks, chunk_range, ua, ctrl, equil, ffit, intr; + integrate_fm_with_ua_ic(chunks, chunk_range, ua, ctrl, equil, mats, intr; backward=false) -> Matrix{ComplexF64} Re-integrate a span of chunks using ua (asymptotic solution) as initial conditions, matching @@ -399,7 +399,7 @@ function integrate_fm_with_ua_ic( ua::Array{ComplexF64,3}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal; backward::Bool = false, psi_ua::Float64 = NaN @@ -422,14 +422,14 @@ function integrate_fm_with_ua_ic( result = zeros(ComplexF64, 2N, 2N) odet_proxy = OdeState(N, 1, 1, 0) dummy_chunk = IntegrationChunk(psi_start, psi_end, false, 0, backward ? -1 : 1) - params = (ctrl, equil, ffit, intr, odet_proxy, dummy_chunk) + params = (ctrl, equil, mats, intr, odet_proxy, dummy_chunk) # Batch 1: columns 1:N of T (big solutions) u0 = zeros(ComplexF64, N, N, 2) u0[:, :, 1] .= ua[:, 1:N, 1] u0[:, :, 2] .= ua[:, 1:N, 2] odet_proxy.spline_hint[] = 1 - odet_proxy.ffit_hint[] = 1 + odet_proxy.mats_hint[] = 1 prob = ODEProblem(sing_der!, u0, tspan, params) sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) result[1:N, 1:N] .= sol.u[end][:, :, 1] @@ -439,7 +439,7 @@ function integrate_fm_with_ua_ic( u0[:, :, 1] .= ua[:, N+1:2N, 1] u0[:, :, 2] .= ua[:, N+1:2N, 2] odet_proxy.spline_hint[] = 1 - odet_proxy.ffit_hint[] = 1 + odet_proxy.mats_hint[] = 1 prob = ODEProblem(sing_der!, u0, tspan, params) sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) result[1:N, N+1:2N] .= sol.u[end][:, :, 1] diff --git a/src/ForceFreeStates/Surfaces/Asymptotics.jl b/src/ForceFreeStates/Surfaces/Asymptotics.jl index cd93a026c..26746c0ca 100644 --- a/src/ForceFreeStates/Surfaces/Asymptotics.jl +++ b/src/ForceFreeStates/Surfaces/Asymptotics.jl @@ -1,7 +1,7 @@ # Frobenius asymptotics and resonant-basis evaluation at singular surfaces. """ - compute_sing_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) + compute_sing_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) Calculate asymptotic vmat and mmat matrices for a singular surface. Formerly `sing_vmat!`. Returns a `SingAsymptotics` struct with the computed data instead of @@ -24,7 +24,7 @@ function compute_sing_asymptotics( singp::SingType, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal; sig::Float64=1.0, alpha_override::Union{Nothing,Vector{ComplexF64}}=nothing, @@ -49,7 +49,7 @@ function compute_sing_asymptotics( # Compute mmat Taylor coefficients with direction parameter sig. # Fortran computes separate mmatl (sig=-1) and mmatr (sig=+1) — the sig flips # odd derivatives of all input quantities (q, F, G, K splines). - compute_sing_mmat!(mmat, singp, ctrl, equil.profiles, ffit, intr; sig=sig, sing_order=sing_order) + compute_sing_mmat!(mmat, singp, ctrl, equil.profiles, mats, intr; sig=sig, sing_order=sing_order) # Extract direction-specific m0mat from zeroth-order mmat m0mat = if length(r1) == 1 @@ -131,7 +131,7 @@ function compute_sing_asymptotics( end """ - compute_sing_mmat!(mmat::Array{ComplexF64,4}, singp::SingType, ctrl::ForceFreeStatesControl, profiles::Equilibrium.ProfileSplines, ffit::FourFitVars, intr::ForceFreeStatesInternal) + compute_sing_mmat!(mmat::Array{ComplexF64,4}, singp::SingType, ctrl::ForceFreeStatesControl, profiles::Equilibrium.ProfileSplines, mats::MatrixSplines, intr::ForceFreeStatesInternal) Calculate asymptotic mmat matrix for a singular surface. Formerly `sing_mmat!`. Performs the same function as `sing_mmat` in the Fortran code. Main differences are 1-indexing for @@ -166,7 +166,7 @@ Add a spline for F directly instead of the lower triangular factorization to avo singp::SingType, ctrl::ForceFreeStatesControl, profiles::Equilibrium.ProfileSplines, - ffit::FourFitVars, + mats::MatrixSplines, intr::ForceFreeStatesInternal; sig::Float64=1.0, sing_order::Int=ctrl.sing_order @@ -200,28 +200,28 @@ Add a spline for F directly instead of the lower triangular factorization to avo q_d2(singp.psifac), sig * q_d3(singp.psifac)) - # Evaluate fmats_lower and derivatives, applying sig to odd derivatives. + # Evaluate F_spline_lower and derivatives, applying sig to odd derivatives. # Fortran sing_mmat multiplies fmats_f1 and fmats_f3 by sig in the Taylor products. - ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.ideal.fmats_lower(vec(@view(f_lower_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + mats.ideal.F_spline_lower(vec(@view(f_lower_interp[:, :, 1])), singp.psifac; hint=mats._hint) + mats.ideal.F_spline_lower(vec(@view(f_lower_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + mats.ideal.F_spline_lower(vec(@view(f_lower_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + mats.ideal.F_spline_lower(vec(@view(f_lower_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views f_lower_interp[:, :, 2] .*= sig # 1st derivative @views f_lower_interp[:, :, 4] .*= sig # 3rd derivative - # Evaluate gmats and derivatives, applying sig to odd derivatives - ffit.ideal.gmats(vec(@view(g_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.ideal.gmats(vec(@view(g_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.ideal.gmats(vec(@view(g_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.ideal.gmats(vec(@view(g_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + # Evaluate G_spline and derivatives, applying sig to odd derivatives + mats.ideal.G_spline(vec(@view(g_interp[:, :, 1])), singp.psifac; hint=mats._hint) + mats.ideal.G_spline(vec(@view(g_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + mats.ideal.G_spline(vec(@view(g_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + mats.ideal.G_spline(vec(@view(g_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views g_interp[:, :, 2] .*= sig @views g_interp[:, :, 4] .*= sig - # Evaluate kmats and derivatives, applying sig to odd derivatives - ffit.ideal.kmats(vec(@view(k_interp[:, :, 1])), singp.psifac; hint=ffit._hint) - ffit.ideal.kmats(vec(@view(k_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) - ffit.ideal.kmats(vec(@view(k_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) - ffit.ideal.kmats(vec(@view(k_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) + # Evaluate K_spline and derivatives, applying sig to odd derivatives + mats.ideal.K_spline(vec(@view(k_interp[:, :, 1])), singp.psifac; hint=mats._hint) + mats.ideal.K_spline(vec(@view(k_interp[:, :, 2])), singp.psifac; deriv=DerivOp(1)) + mats.ideal.K_spline(vec(@view(k_interp[:, :, 3])), singp.psifac; deriv=DerivOp(2)) + mats.ideal.K_spline(vec(@view(k_interp[:, :, 4])), singp.psifac; deriv=DerivOp(3)) @views k_interp[:, :, 2] .*= sig @views k_interp[:, :, 4] .*= sig @@ -807,13 +807,13 @@ sing_get_dua_res(sing_asymp::SingAsymptotics, dpsi::Float64) = sing_get_dua_res!(Array{ComplexF64,3}(undef, size(sing_asymp.vmat, 1), 2, 2), sing_asymp, dpsi) """ - sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, ua, dua) -> matvec + sing_matvec(mats::MatrixSplines, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, ua, dua) -> matvec Apply the Euler-Lagrange residual operator `L u = -(F u' + K u)' + (K† u' + G u)` to the asymptotic solutions. Port of Fortran `sing_matvec` (sing.f). Returns `matvec`, shape `(numpert_total, size(ua,2))`. -Uses the reduced (Schur-complemented) `ffit.ideal.kmats` (= K̄) and `ffit.ideal.gmats` (= Ḡ) directly, with the +Uses the reduced (Schur-complemented) `mats.ideal.K_spline` (= K̄) and `mats.ideal.G_spline` (= Ḡ) directly, with the **direct** singular factor `singfac = m - n q` applied to `u' = dua[:,:,1]`. The second component `ua[:,:,2]` is the canonical momentum `F u' + K u`, so `-dua[:,:,2] = -(F u' + K u)'`. @@ -822,7 +822,7 @@ Uses the reduced (Schur-complemented) `ffit.ideal.kmats` (= K̄) and `ffit.ideal - `psi`: flux coordinate; `q`: safety factor at `psi` (passed in to avoid re-evaluating the spline) - `ua`, `dua`: asymptotic solution and its ψ-derivative from `sing_get_ua`/`sing_get_dua` """ -function sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, +function sing_matvec(mats::MatrixSplines, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, ua::Array{ComplexF64,3}, dua::Array{ComplexF64,3}) N = intr.numpert_total @@ -833,8 +833,8 @@ function sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Floa kmat = Matrix{ComplexF64}(undef, N, N) gmat = Matrix{ComplexF64}(undef, N, N) - ffit.ideal.kmats(vec(kmat), psi; hint=ffit._hint) - ffit.ideal.gmats(vec(gmat), psi; hint=ffit._hint) + mats.ideal.K_spline(vec(kmat), psi; hint=mats._hint) + mats.ideal.G_spline(vec(gmat), psi; hint=mats._hint) kdag = adjoint(kmat) matvec = zeros(ComplexF64, N, msol) @@ -849,7 +849,7 @@ function sing_matvec(ffit::FourFitVars, intr::ForceFreeStatesInternal, psi::Floa end """ - sing_matvec!(matvec, kmat, gmat, d1, tmp, sfvec, ffit, intr, psi, q, ua, dua) -> matvec + sing_matvec!(matvec, kmat, gmat, d1, tmp, sfvec, mats, intr, psi, q, ua, dua) -> matvec Allocation-free, in-place form of [`sing_matvec`](@ref) for the hot resonant-quadrature integrand. All scratch is caller-owned: `kmat`/`gmat` are `N×N`, `d1`/`tmp`/`sfvec` are length-`N`, `matvec` is @@ -859,7 +859,7 @@ All scratch is caller-owned: `kmat`/`gmat` are `N×N`, `d1`/`tmp`/`sfvec` are le 1e-16 change into ~1e-3 in Δ′. """ function sing_matvec!(matvec::AbstractMatrix{ComplexF64}, kmat::Matrix{ComplexF64}, gmat::Matrix{ComplexF64}, - d1::Vector{ComplexF64}, tmp::Vector{ComplexF64}, sfvec::Vector{Float64}, ffit::FourFitVars, + d1::Vector{ComplexF64}, tmp::Vector{ComplexF64}, sfvec::Vector{Float64}, mats::MatrixSplines, intr::ForceFreeStatesInternal, psi::Float64, q::Float64, ua::AbstractArray{ComplexF64,3}, dua::AbstractArray{ComplexF64,3}) N = intr.numpert_total @@ -872,8 +872,8 @@ function sing_matvec!(matvec::AbstractMatrix{ComplexF64}, kmat::Matrix{ComplexF6 sfvec[idx] = mm - q * nn end - ffit.ideal.kmats(vec(kmat), psi; hint=ffit._hint) - ffit.ideal.gmats(vec(gmat), psi; hint=ffit._hint) + mats.ideal.K_spline(vec(kmat), psi; hint=mats._hint) + mats.ideal.G_spline(vec(gmat), psi; hint=mats._hint) kdag = adjoint(kmat) for isol in 1:msol diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 6bb23e912..5154456ce 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -217,10 +217,10 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi p_vec = zeros(ComplexF64, np * np) pa_vec = zeros(ComplexF64, np * np) r1_vec = zeros(ComplexF64, np * np) - kin.f0mats(f0_vec, psi; hint=hint) - kin.pmats(p_vec, psi; hint=hint) - kin.paats(pa_vec, psi; hint=hint) - kin.r1mats(r1_vec, psi; hint=hint) + kin.F0_spline(f0_vec, psi; hint=hint) + kin.P_spline(p_vec, psi; hint=hint) + kin.P_spline_adj(pa_vec, psi; hint=hint) + kin.R1_spline(r1_vec, psi; hint=hint) f0mat = reshape(f0_vec, np, np) pmat = reshape(p_vec, np, np) paat = reshape(pa_vec, np, np) @@ -240,7 +240,7 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi end """ - find_kinetic_singular_surfaces!(ffit, equil, intr; ngrid=2000, cond_threshold=1e8) + find_kinetic_singular_surfaces!(mats, equil, intr; ngrid=2000, cond_threshold=1e8) Find kinetically-displaced singular surfaces — locations where cond(F̄) peaks, indicating near-singularity of the kinetic F̄ matrix in the ODE RHS. Populates @@ -257,8 +257,8 @@ Algorithm: 3. Refine each peak with golden-section minimization of -cond 4. Filter by threshold and resonance condition """ -function find_kinetic_singular_surfaces!(ffit::FourFitVars, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) - kin = ffit.kinetic +function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) + kin = mats.kinetic kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call make_kinetic_matrix first") psilow = equil.profiles.xs[1] psihigh = intr.psilim diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index 69c24616b..a73085c40 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -59,7 +59,7 @@ function store_ode_data!(odet::OdeState, psi::Float64, u) end """ - materialize_derivative_stores!(odet, equil, ffit, intr) -> Bool + materialize_derivative_stores!(odet, equil, mats, intr) -> Bool Fill `odet.du_store` (dΞ_ψ/dψ) and `odet.xi_s_store` from the stored solution, returning whether they hold valid data afterwards. Idempotent: a no-op when `du_store_populated` is @@ -71,20 +71,20 @@ fixup transforms and the free-boundary normalization is exact rather than merely Euler-Lagrange system is linear in `u`, and both operations right-multiply the solution by a mixing matrix `T`, so `du(ψ, u·T) = du(ψ, u)·T`. -Returns `false` without allocating when there is nothing to work from — no `ffit`, no stored +Returns `false` without allocating when there is nothing to work from — no `mats`, no stored steps, or a solution held in a basis the Euler-Lagrange kernel does not apply to (the sparse parallel path, which stores chunk-endpoint Riccati matrices). """ function materialize_derivative_stores!( odet::OdeState, equil::Equilibrium.PlasmaEquilibrium, - ffit::Union{FourFitVars,Nothing}, + mats::Union{MatrixSplines,Nothing}, intr::ModeSpace ) odet.du_store_populated && return true - (isnothing(ffit) || odet.step == 0 || isempty(odet.u_store) || !odet.u_store_el_basis) && return false + (isnothing(mats) || odet.step == 0 || isempty(odet.u_store) || !odet.u_store_el_basis) && return false - kinetic = ffit.kinetic !== nothing + kinetic = mats.kinetic !== nothing nstep = min(odet.step, size(odet.u_store, 4)) npert = odet.numpert_total odet.du_store = Array{ComplexF64}(undef, npert, npert, nstep) @@ -94,10 +94,10 @@ function materialize_derivative_stores!( u = zeros(ComplexF64, npert, npert, 2) @views for istep in 1:nstep u .= odet.u_store[:, :, :, istep] - odet.q = el_derivatives!(du, u, kinetic, equil, ffit, intr, odet.psi_store[istep], odet.spline_hint, odet.ffit_hint) + odet.q = el_derivatives!(du, u, kinetic, equil, mats, intr, odet.psi_store[istep], odet.spline_hint, odet.mats_hint) odet.du_store[:, :, istep] .= du[:, :, 1] - compute_node_xi_s!(odet.xi_s_store[:, :, istep], du[:, :, 1], u[:, :, 1], ffit, odet.psi_store[istep]; - kinetic=kinetic, hint=odet.ffit_hint) + compute_node_xi_s!(odet.xi_s_store[:, :, istep], du[:, :, 1], u[:, :, 1], mats, odet.psi_store[istep]; + kinetic=kinetic, hint=odet.mats_hint) end odet.du_store_populated = true diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 5d2dc3cec..1716454ae 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -241,8 +241,8 @@ function main_from_inputs( ffs_start = time() locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, ffit = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) - ffs_result = run_force_free_states(ctrl, equil, ffit, intr, metric) + metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) + ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric) if ctrl.write_outputs_to_HDF5 write_outputs_to_HDF5( @@ -470,7 +470,7 @@ function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.Pl end """ - prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, ffit) + prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, mats) Set up the force-free-states solve on `intr`: integration limits, the surviving singular surfaces and their GGJ coefficients, the poloidal mode range, and the metric plus @@ -558,8 +558,8 @@ function prepare_force_free_states!( @info "Computing F, G, and K matrices" end - # Compute matrices and populate FourFitVars struct - ffit = make_matrix(equil, intr, metric) + # Compute matrices and build the MatrixSplines container + mats = make_matrix(equil, intr, metric) if ctrl.kinetic_factor > 0 if ctrl.verbose @@ -572,22 +572,22 @@ function prepare_force_free_states!( KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles) - ffit = make_kinetic_matrix(ctrl, equil, ffit, intr, metric; + mats = make_kinetic_matrix(ctrl, equil, mats, intr, metric; calculated_source=calculated_cb) # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings. # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings; # singfac_min == 0 preserves single-chunk behavior. if ctrl.singfac_min > 0 - find_kinetic_singular_surfaces!(ffit, equil, intr) + find_kinetic_singular_surfaces!(mats, equil, intr) end end - return metric, ffit + return metric, mats end """ - run_force_free_states(ctrl, equil, ffit, intr, metric) -> ForceFreeStatesResult + run_force_free_states(ctrl, equil, mats, intr, metric) -> ForceFreeStatesResult Run the formalism selected by `ctrl.integrator` — the standalone Galerkin solve, or the Euler-Lagrange sweep with its free-boundary energies and Δ′ BVP — and publish its products as a @@ -596,7 +596,7 @@ Euler-Lagrange sweep with its free-boundary energies and Δ′ BVP — and publi function run_force_free_states( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit, + mats, intr::ForceFreeStatesInternal, metric ) @@ -614,14 +614,14 @@ function run_force_free_states( error("integrator = \"galerkin\" does not support kinetic runs (kinetic_factor > 0); use integrator = \"forward\".") gal_start = time() wv = ctrl.vac_flag ? first(ForceFreeStates.compute_scaled_wv(ctrl, equil, intr)) : nothing - gal_data, gal_dp = galerkin_solve(ctrl, equil, ffit, intr; wv=wv) + gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv) @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s" else # Integrate Euler-Lagrange Equation if ctrl.verbose @info "Integrating Euler-Lagrange equation" end - odet, fm_propagators, fm_chunks, fm_S_left = eulerlagrange_integration(ctrl, equil, ffit, intr) + odet, fm_propagators, fm_chunks, fm_S_left = eulerlagrange_integration(ctrl, equil, mats, intr) if odet.nzero > 0 && ctrl.verbose @warn "Fixed-boundary mode unstable for n = $nstring" end @@ -632,7 +632,7 @@ function run_force_free_states( wall_desc = intr.wall_settings.shape == "nowall" ? "no wall" : intr.wall_settings.shape @info "Computing free boundary energies ($wall_desc)" end - free_energies = free_run(odet, ctrl, equil, ffit, intr) + free_energies = free_run(odet, ctrl, equil, mats, intr) normalize_eigenfunctions!(odet, free_energies.wt, equil.psio) if real(free_energies.et[1]) < 0 if ctrl.verbose @@ -653,13 +653,13 @@ function run_force_free_states( ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; wv=free_energies.wv, psio=equil.psio, debug=ctrl.verbose, S_at_surface_left=fm_S_left, - ctrl=ctrl, equil=equil, ffit=ffit) + ctrl=ctrl, equil=equil, mats=mats) end end end # Publish the solve: from here on the downstream stages read the result, never `intr`. - return build_result(Symbol(ctrl.integrator), ctrl, equil, intr, metric, ffit, odet, free_energies, gal_data, gal_dp) + return build_result(Symbol(ctrl.integrator), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) end """ @@ -765,8 +765,8 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat end locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, ffit = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) - result = run_force_free_states(ctrl, equil, ffit, intr, metric) + metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) + result = run_force_free_states(ctrl, equil, mats, intr, metric) if ctrl.write_outputs_to_HDF5 write_outputs_to_HDF5(result; locstab=locstab, ballooning_boundary=ballooning_boundary) @@ -1017,7 +1017,7 @@ function write_outputs_to_HDF5( ctrl = result.control equil = result.equil - ffit = result.ffit + mats = result.mats free_energies = result.free_boundary gal_data = result.galerkin diag = result.diagnostics @@ -1287,27 +1287,27 @@ function write_outputs_to_HDF5( elm = "ForceFreeStates/EulerLagrangeMatrices" out_h5["$elm/psi"] = xs # Ideal primitive matrices (A, B, C, D, E, H) - out_h5["$elm/Ideal/A"] = _eval_mat_spline(ffit.ideal.amats) - out_h5["$elm/Ideal/B"] = _eval_mat_spline(ffit.ideal.bmats) - out_h5["$elm/Ideal/C"] = _eval_mat_spline(ffit.ideal.cmats) - out_h5["$elm/Ideal/D"] = _eval_mat_spline(ffit.ideal.dmats_prim) - out_h5["$elm/Ideal/E"] = _eval_mat_spline(ffit.ideal.emats_prim) - out_h5["$elm/Ideal/H"] = _eval_mat_spline(ffit.ideal.hmats) + out_h5["$elm/Ideal/A"] = _eval_mat_spline(mats.ideal.A_spline) + out_h5["$elm/Ideal/B"] = _eval_mat_spline(mats.ideal.B_spline) + out_h5["$elm/Ideal/C"] = _eval_mat_spline(mats.ideal.C_spline) + out_h5["$elm/Ideal/D"] = _eval_mat_spline(mats.ideal.D_spline_prim) + out_h5["$elm/Ideal/E"] = _eval_mat_spline(mats.ideal.E_spline_prim) + out_h5["$elm/Ideal/H"] = _eval_mat_spline(mats.ideal.H_spline) # Ideal derived matrices (F, K, G) - out_h5["$elm/Ideal/F"] = _eval_mat_spline(ffit.ideal.fmats_lower) - out_h5["$elm/Ideal/K"] = _eval_mat_spline(ffit.ideal.kmats) - out_h5["$elm/Ideal/G"] = _eval_mat_spline(ffit.ideal.gmats) + out_h5["$elm/Ideal/F"] = _eval_mat_spline(mats.ideal.F_spline_lower) + out_h5["$elm/Ideal/K"] = _eval_mat_spline(mats.ideal.K_spline) + out_h5["$elm/Ideal/G"] = _eval_mat_spline(mats.ideal.G_spline) # Kinetic-modified matrices - kin = ffit.kinetic + kin = mats.kinetic if kin !== nothing - out_h5["$elm/Kinetic/A"] = _eval_mat_spline(kin.amats) - out_h5["$elm/Kinetic/B"] = _eval_mat_spline(kin.bmats) - out_h5["$elm/Kinetic/C"] = _eval_mat_spline(kin.cmats) - out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(kin.f0mats) - out_h5["$elm/Kinetic/K"] = _eval_mat_spline(kin.kkmats) - out_h5["$elm/Kinetic/G"] = _eval_mat_spline(kin.gaats) + out_h5["$elm/Kinetic/A"] = _eval_mat_spline(kin.A_spline) + out_h5["$elm/Kinetic/B"] = _eval_mat_spline(kin.B_spline) + out_h5["$elm/Kinetic/C"] = _eval_mat_spline(kin.C_spline) + out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(kin.F0_spline) + out_h5["$elm/Kinetic/K"] = _eval_mat_spline(kin.Kk_spline) + out_h5["$elm/Kinetic/G"] = _eval_mat_spline(kin.G_spline_adj) end # Self-describing metadata pass (long_name/units/dims + dimension scales). diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index fe0f705af..455e515b2 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -8,7 +8,7 @@ injected from `GeneralizedPerturbedEquilibrium.main`). """ """ - compute_calculated_kinetic_matrices(ffs_ctrl, equil, ffs_intr, metric, ffit; + compute_calculated_kinetic_matrices(ffs_ctrl, equil, ffs_intr, metric, mats; kf_ctrl=KineticForcesControl(), kinetic_profiles) → (kw_flat, kt_flat) @@ -36,7 +36,7 @@ section. - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines - `ffs_intr`: ForceFreeStatesInternal (mode indexing) - `metric`: MetricData (provides ψ grid via `metric.xs`) - - `ffit`: FourFitVars (used only for `numpert_total` cross-check) + - `mats`: MatrixSplines (used only for `numpert_total` cross-check) # Keyword arguments @@ -56,7 +56,7 @@ function compute_calculated_kinetic_matrices( equil, ffs_intr, metric, - ffit; + mats; kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines ) @@ -66,7 +66,7 @@ function compute_calculated_kinetic_matrices( npert = ffs_intr.npert np = ffs_intr.numpert_total - @assert ffit.numpert_total == np "FourFitVars and ForceFreeStatesInternal disagree on numpert_total" + @assert mats.numpert_total == np "MatrixSplines and ForceFreeStatesInternal disagree on numpert_total" kw_flat = zeros(ComplexF64, mpsi, np^2, 6) kt_flat = zeros(ComplexF64, mpsi, np^2, 6) diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 74ec75937..8c349ae1d 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -37,7 +37,7 @@ Covariant components from metric tensor contraction (matches Fortran gpeq_cova): """ reconstruct_physical_fields( response_vector, flux_matrix, solution, - equil, ffs, intr, metric, ffit, ctrl + equil, ffs, intr, metric, mats, ctrl ) -> (xi_modes, b_modes) Reconstruct displacement and perturbed magnetic field from eigenmode response. @@ -74,7 +74,7 @@ function reconstruct_physical_fields( ffs::ForceFreeStatesResult, intr::PerturbedEquilibriumInternal, metric::MetricData, - ffit::FourFitVars, + mats::MatrixSplines, ctrl::PerturbedEquilibriumControl ) npsi = size(solution.u_store, 4) @@ -106,7 +106,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, ffit, ctrl + psi_grid, equil, ffs, mats, ctrl ) # Compute regularized (modified) b-field components (matches Fortran gpeq_sol bmt/bmz) @@ -321,7 +321,7 @@ end """ compute_clebsch_displacements( xi_psi_modes, xi_psi1_modes, xi_s_modes, - psi_grid, equil, ffs, ffit, ctrl + psi_grid, equil, ffs, mats, ctrl ) -> (clebsch_psi, clebsch_psi1, clebsch_alpha) Compute Clebsch displacement components for PENTRC output. @@ -335,7 +335,7 @@ Matches Fortran gpeq_sol regularization + gpout_xclebsch output convention: When reg_spot=0, clebsch_psi1 = xi_psi1 and clebsch_alpha = xi_s/χ₁ (no regularization). The regularized xms is computed as -A⁻¹(B·xmp1 + C·xsp) matching Fortran gpeq_sol, -where A, B, C are the stability matrices evaluated at each ψ via ffit interpolants. +where A, B, C are the stability matrices evaluated at each ψ from `mats`. """ function compute_clebsch_displacements( xi_psi_modes::Matrix{ComplexF64}, @@ -344,7 +344,7 @@ function compute_clebsch_displacements( psi_grid::Vector{Float64}, equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, - ffit::FourFitVars, + mats::MatrixSplines, ctrl::PerturbedEquilibriumControl ) npsi, mpert = size(xi_psi_modes) @@ -365,7 +365,7 @@ function compute_clebsch_displacements( end # A/B/C of the active model, matching what the ODE integrated. - active_mats = ffit.kinetic === nothing ? ffit.ideal : ffit.kinetic + active_mats = mats.kinetic === nothing ? mats.ideal : mats.kinetic # Per-thread workspaces: matrix ops and spline hints are not safe to share across threads. # Size by maxthreadid() and index by threadid() under :static scheduling (GPEC convention). @@ -400,16 +400,16 @@ function compute_clebsch_displacements( # Compute regularized xms = -A⁻¹(B·xmp1 + C·xsp) (matches Fortran gpeq_sol) # Evaluate stability matrices at this psi - active_mats.amats(view(amat, :), psi_norm; hint=hint) - active_mats.bmats(view(bmat, :), psi_norm; hint=hint) - active_mats.cmats(view(cmat_buf, :), psi_norm; hint=hint) + active_mats.A_spline(view(amat, :), psi_norm; hint=hint) + active_mats.B_spline(view(bmat, :), psi_norm; hint=hint) + active_mats.C_spline(view(cmat_buf, :), psi_norm; hint=hint) # xms = -(A\B)*xmp1 - (A\C)*xsp xsp_vec = view(xi_psi_modes, ipsi, :) mul!(xms_vec, bmat, xmp1_vec) # xms = B*xmp1 mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im) # xms += C*xsp # cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by - # active_mats.amats each surface), avoiding a fresh factorization allocation per surface. + # active_mats.A_spline each surface), avoiding a fresh factorization allocation per surface. # NOTE: this assumes the ideal A (positive-definite Newcomb kinetic-energy form). The # kinetic A is non-Hermitian and needs an LU, as compute_node_xi_s! does — see the # `active_mats` binding above. diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 245fa2949..66c1f3253 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: SolutionProfiles, ForceFreeStatesResult, FourFitVars, MetricData +import ..ForceFreeStates: SolutionProfiles, ForceFreeStatesResult, MatrixSplines, MetricData import ..Vacuum import ..ForcingTerms import ..ForcingTerms: ForcingMode, CoilSet, load_forcing_data!, convert_forcing_normalization! @@ -68,7 +68,7 @@ function compute_perturbed_equilibrium( state = PerturbedEquilibriumState() equil = ffs.equil - ffit = ffs.ffit + mats = ffs.mats mthvac = ffs.control.mthvac # Step 0: Initialize mode arrays for convenient indexing @@ -92,14 +92,14 @@ function compute_perturbed_equilibrium( if ctrl.compute_response && 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, mthvac, ffs, intr, ctrl, ffs.metric, ffit) + compute_plasma_response!(state, equil, solution, ffs.free_boundary.wt0, mthvac, ffs, intr, ctrl, ffs.metric, mats) end # Step 3: Compute 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, mthvac, ffs, intr, ctrl, ffit) + compute_singular_coupling_metrics!(state, equil, solution, mthvac, ffs, intr, ctrl, mats) end # Step 4: Output eigenmode fields (integrated into HDF5 output) diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index d305a079b..c82d86cd9 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -1,7 +1,7 @@ """ compute_plasma_response!( state, equil, solution, wt0, mthvac, ffs, - intr, ctrl, metric, ffit + intr, ctrl, metric, mats ) Compute plasma response to external forcing using ForceFreeStates eigenmode solutions. @@ -24,7 +24,7 @@ function compute_plasma_response!( intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl, metric::MetricData, - ffit::FourFitVars + mats::MatrixSplines ) if ctrl.verbose @info "Computing plasma response (wt0-based inductance)" @@ -99,7 +99,7 @@ function compute_plasma_response!( xi_modes, b_modes = reconstruct_physical_fields( response_flux, flux_matrix, solution, equil, ffs, intr, - metric, ffit, ctrl + metric, mats, ctrl ) npsi = size(solution.u_store, 4) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 04764f203..71d3faa95 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -156,12 +156,12 @@ function _solution_at( end """ - _el_solution_at(psi, resnum, odet, ffit, equil, ffs, nstep) -> (u, du) + _el_solution_at(psi, resnum, odet, mats, equil, ffs, nstep) -> (u, du) Evaluate the `resnum` row of Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution via the ideal Euler-Lagrange relation Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) [Glasser 2016 eqs. 22-24], with u₁, u₂ Hermite-interpolated to `psi`. Only valid for ideal runs where -`ffit.ideal.fmats_lower` and `kmats` generated the solution. +`mats.ideal.F_spline_lower` and `K_spline` generated the solution. The Hermite slopes need du₂ as well as du₁, and du₂ is not stored: both are evaluated here from the derivative kernel at the two bracketing nodes, which is where the handful of @@ -171,7 +171,7 @@ function _el_solution_at( psi::Float64, resnum::Int, odet::SolutionProfiles, - ffit::FourFitVars, + mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, nstep::Int @@ -190,8 +190,8 @@ function _el_solution_at( q_hint = Ref(1) du_a = zeros(ComplexF64, npert, npert, 2) du_b = zeros(ComplexF64, npert, npert, 2) - ForceFreeStates.el_derivatives!(du_a, odet.u_store[:, :, :, il], false, equil, ffit, ffs, psi_a, q_hint, hint) - ForceFreeStates.el_derivatives!(du_b, odet.u_store[:, :, :, ir], false, equil, ffit, ffs, psi_b, q_hint, hint) + ForceFreeStates.el_derivatives!(du_a, odet.u_store[:, :, :, il], false, equil, mats, ffs, psi_a, q_hint, hint) + ForceFreeStates.el_derivatives!(du_b, odet.u_store[:, :, :, ir], false, equil, mats, ffs, psi_b, q_hint, hint) du1_a = @view du_a[:, :, 1] du1_b = @view du_b[:, :, 1] du2_a = @view du_a[:, :, 2] @@ -206,8 +206,8 @@ function _el_solution_at( q_e = equil.profiles.q_spline(psi) singfac_inv = vec([1.0 / (m - q_e * n) for m in ffs.mlow:ffs.mhigh, n in ffs.nlow:ffs.nhigh]) fmat_lower = Matrix{ComplexF64}(undef, npert, npert) - ffit.ideal.fmats_lower(vec(fmat_lower), psi; hint=hint) - ffit.ideal.kmats(vec(kmat), psi; hint=hint) + mats.ideal.F_spline_lower(vec(fmat_lower), psi; hint=hint) + mats.ideal.K_spline(vec(kmat), psi; hint=hint) du1_e = u2_e .* singfac_inv du1_e .-= kmat * u1_e ldiv!(LowerTriangular(fmat_lower), du1_e) @@ -226,7 +226,7 @@ end ffs::ForceFreeStatesResult, intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl, - ffit::FourFitVars + mats::MatrixSplines ) Compute singular layer coupling matrices and applied resonant vectors. @@ -258,7 +258,7 @@ function compute_singular_coupling_metrics!( ffs::ForceFreeStatesResult, intr::PerturbedEquilibriumInternal, ctrl::PerturbedEquilibriumControl, - ffit::FourFitVars + mats::MatrixSplines ) ctrl.verbose && @info "Computing singular coupling metrics (GPEC method)" @@ -342,7 +342,7 @@ function compute_singular_coupling_metrics!( # @threads region. nstep = solution.step # ξ′ evaluation preference: the ideal EL relation, or the interpolated stored RHS for kinetic runs. - use_el = ffit.kinetic === nothing + use_el = mats.kinetic === nothing _blas_nthreads = BLAS.get_num_threads() BLAS.set_num_threads(1) try @@ -387,8 +387,8 @@ function compute_singular_coupling_metrics!( u_r, ud_r = _gal_solution_at(rpsi, resnum, solution, nstep) elseif use_el # interpolate u and evaluate dξ/dψ from the ideal EL relation - u_l, ud_l = _el_solution_at(lpsi, resnum, solution, ffit, equil, ffs, nstep) - u_r, ud_r = _el_solution_at(rpsi, resnum, solution, ffit, equil, ffs, nstep) + u_l, ud_l = _el_solution_at(lpsi, resnum, solution, mats, equil, ffs, nstep) + u_r, ud_r = _el_solution_at(rpsi, resnum, solution, mats, equil, ffs, nstep) else # interpolate u and the stored dξ/dψ, weighted to remove the resonant pole u_l, ud_l = _solution_at(lpsi, sing_surf.psifac, resnum, m_res, nn, solution, equil, nstep) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 6547352e6..ed12ef89c 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -450,19 +450,19 @@ end intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - ffit = FFS.make_matrix(equil, intr, metric) - odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) - return odet, ctrl, equil, ffit, intr + mats = FFS.make_matrix(equil, intr, metric) + odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, mats, intr) + return odet, ctrl, equil, mats, intr end - odet, ctrl, equil, ffit, intr = setup_solovev_run() + odet, ctrl, equil, mats, intr = setup_solovev_run() # Untouched copy of the solution, for the column-transform check further down. odet_pristine = deepcopy(odet) @testset "fills the stores once" begin @test isempty(odet.du_store) @test !odet.du_store_populated - @test FFS.materialize_derivative_stores!(odet, equil, ffit, intr) + @test FFS.materialize_derivative_stores!(odet, equil, mats, intr) @test odet.du_store_populated @test size(odet.du_store) == (intr.numpert_total, intr.numpert_total, odet.step) @test size(odet.xi_s_store) == (intr.numpert_total, intr.numpert_total, odet.step) @@ -471,7 +471,7 @@ end # Idempotent: a second call must not overwrite what is already there. du_first = copy(odet.du_store) - @test FFS.materialize_derivative_stores!(odet, equil, ffit, intr) + @test FFS.materialize_derivative_stores!(odet, equil, mats, intr) @test odet.du_store == du_first end @@ -482,8 +482,8 @@ end for istep in (1, odet.step ÷ 2, odet.step) psi = odet.psi_store[istep] u = odet.u_store[:, :, :, istep] - FFS.el_derivatives!(du, u, false, equil, ffit, intr, psi, Ref(1), Ref(1)) - FFS.compute_node_xi_s!(xi_s, @view(du[:, :, 1]), @view(u[:, :, 1]), ffit, psi) + FFS.el_derivatives!(du, u, false, equil, mats, intr, psi, Ref(1), Ref(1)) + FFS.compute_node_xi_s!(xi_s, @view(du[:, :, 1]), @view(u[:, :, 1]), mats, psi) @test odet.du_store[:, :, istep] == du[:, :, 1] @test odet.xi_s_store[:, :, istep] == xi_s end @@ -500,7 +500,7 @@ end odet_t.u_store[:, :, 1, istep] = odet_t.u_store[:, :, 1, istep] * T odet_t.u_store[:, :, 2, istep] = odet_t.u_store[:, :, 2, istep] * T end - @test FFS.materialize_derivative_stores!(odet_t, equil, ffit, intr) + @test FFS.materialize_derivative_stores!(odet_t, equil, mats, intr) for istep in (1, odet.step ÷ 2, odet.step) @test isapprox(odet_t.du_store[:, :, istep], odet.du_store[:, :, istep] * T; rtol=1e-10) @@ -512,7 +512,7 @@ end odet.du_store_populated = false odet.du_store = Array{ComplexF64}(undef, intr.numpert_total, intr.numpert_total, 0) odet.u_store_el_basis = false - @test !FFS.materialize_derivative_stores!(odet, equil, ffit, intr) + @test !FFS.materialize_derivative_stores!(odet, equil, mats, intr) @test isempty(odet.du_store) @test !odet.du_store_populated end diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 7e7737e06..9c39ee78a 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -126,10 +126,10 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) - GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr) + GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) @@ -213,10 +213,10 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) - GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr) + GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) # Default (bidirectional=false): all chunks should have direction=+1 chunks_fwd = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) @@ -275,9 +275,9 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert 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) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr end @@ -334,9 +334,9 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert 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) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr end @@ -384,11 +384,11 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) # Use the first chunk from chunk_el_integration_bounds: guaranteed rational-free interior odet_tmp = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, 10, 5, intr.msing) - GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet_tmp, ctrl, ffit, equil.profiles, intr) + GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet_tmp, ctrl, mats, equil.profiles, intr) chunks_tmp = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet_tmp, ctrl, intr) chunk1 = chunks_tmp[1] a = chunk1.psi_start @@ -436,10 +436,10 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert 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) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) # Derivatives are recomputed on demand; materialize so the stores can be compared. - GeneralizedPerturbedEquilibrium.ForceFreeStates.materialize_derivative_stores!(odet, equil, ffit, intr) + GeneralizedPerturbedEquilibrium.ForceFreeStates.materialize_derivative_stores!(odet, equil, mats, intr) return odet end @@ -498,14 +498,14 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = 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) + GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) + vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_delta_prime_matrix!( intr, fm_propagators, fm_chunks; wv=vac.wv, psio=equil.psio, - S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) + S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, mats=mats) msing = intr.msing dpm = intr.delta_prime_matrix diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index 5b2199df8..1f4dd7dc5 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -3,7 +3,7 @@ using LinearAlgebra, Random, TOML const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates # Configure a fresh ForceFreeStatesInternal from an already-built equilibrium. -# Cheap (sing_lim! + sing_find! + field assignment). Separate from equil/ffit +# Cheap (sing_lim! + sing_find! + field assignment). Separate from equil/mats # setup because intr is mutated by each integration (sing[s].delta_prime etc.). function make_solovev_intr(inputs, ctrl, equil, ex) intr = FFS.ForceFreeStatesInternal(; dir_path=ex) @@ -86,7 +86,7 @@ end # ── Shared Solovev setup ────────────────────────────────────────────────── # - # equil (Grad-Shafranov solve) and ffit (metric matrices) are expensive and + # equil (Grad-Shafranov solve) and mats (metric matrices) are expensive and # immutable after construction — built ONCE and shared across all tests below. # intr is cheap to (re)initialize but is mutated by each integration run # (sing[s].delta_prime etc.), so a fresh copy is made for each integration. @@ -107,25 +107,25 @@ end intr_tmp = make_solovev_intr(inputs, ctrl, equil, ex) metric = FFS.make_metric(equil, intr_tmp.mpert) - ffit = FFS.make_matrix(equil, intr_tmp, metric) + mats = FFS.make_matrix(equil, intr_tmp, metric) N = intr_tmp.numpert_total # Riccati integration. The driver returns (odet, propagators, chunks, S_at_surface_left); # only odet is used here. intr_ric = make_solovev_intr(inputs, ctrl, equil, ex) - odet_ric, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr_ric) + odet_ric, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr_ric) # Save inline Δ' values before any test that calls compute_delta_prime_from_ca! # (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, mats, intr_ric) et_ric = real(vac_ric.et[1]) # Forward integration (needed only for energy comparison). intr_fwd = make_solovev_intr(inputs, ctrl, equil, ex) - odet_fwd, _, _, _ = FFS.forward_eulerlagrange_integration(ctrl, equil, ffit, intr_fwd) - vac_fwd = FFS.free_run(odet_fwd, ctrl, equil, ffit, intr_fwd) + odet_fwd, _, _, _ = FFS.forward_eulerlagrange_integration(ctrl, equil, mats, intr_fwd) + vac_fwd = FFS.free_run(odet_fwd, ctrl, equil, mats, intr_fwd) et_fwd = real(vac_fwd.et[1]) # ───────────────────────────────────────────────────────────────────────── @@ -165,7 +165,7 @@ end # Use an initialized OdeState just for spline_hint and chunk bounds odet_tmp = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr_ric.msing) - FFS.initialize_el_at_axis!(odet_tmp, ctrl, ffit, equil.profiles, intr_ric) + FFS.initialize_el_at_axis!(odet_tmp, ctrl, mats, equil.profiles, intr_ric) chunks = FFS.chunk_el_integration_bounds(odet_tmp, ctrl, intr_ric) # 30% into each chunk: away from singularities at psi_end @@ -181,9 +181,9 @@ end L = zeros(ComplexF64, N, N) Kmat = zeros(ComplexF64, N, N) Gmat = zeros(ComplexF64, N, N) - ffit.ideal.fmats_lower(vec(L), psi; hint=ffit._hint) - ffit.ideal.kmats(vec(Kmat), psi; hint=ffit._hint) - ffit.ideal.gmats(vec(Gmat), psi; hint=ffit._hint) + mats.ideal.F_spline_lower(vec(L), psi; hint=mats._hint) + mats.ideal.K_spline(vec(Kmat), psi; hint=mats._hint) + mats.ideal.G_spline(vec(Gmat), psi; hint=mats._hint) q = equil.profiles.q_spline(psi) singfac = vec(1.0 ./ ((intr_ric.mlow:intr_ric.mhigh) .- q .* (intr_ric.nlow:intr_ric.nhigh)')) @@ -202,7 +202,7 @@ end u_ric[:, :, 1] .= S u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N) dummy = FFS.IntegrationChunk(psi, psi, false, 0, 1) - params = (ctrl, equil, ffit, intr_ric, odet_tmp, dummy) + params = (ctrl, equil, mats, intr_ric, odet_tmp, dummy) FFS.riccati_der!(du_ric, u_ric, params, psi) rel_err = norm(du_ric[:, :, 1] - dS_manual) / max(norm(dS_manual), 1e-10) diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index 9ae95d526..ddf42311f 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -108,19 +108,19 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex psifac_dummy = collect(range(0, 1, 10)); points = length(psifac_dummy) amat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/amat.dat")) - amats = copyForSplines(amat, psifac_dummy) + A_spline = copyForSplines(amat, psifac_dummy) bmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/bmat.dat")) - bmats = copyForSplines(bmat, psifac_dummy) + B_spline = copyForSplines(bmat, psifac_dummy) cmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/cmat.dat")) - cmats = copyForSplines(cmat, psifac_dummy) + C_spline = copyForSplines(cmat, psifac_dummy) fmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/fmat.dat")) fmat .= cholesky(Hermitian(fmat)).L; fmats = copyForSplines(fmat, psifac_dummy) kmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/kmat.dat")) - kmats = copyForSplines(kmat, psifac_dummy) + K_spline = copyForSplines(kmat, psifac_dummy) gmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/gmat.dat")) - gmats = copyForSplines(gmat, psifac_dummy) + G_spline = copyForSplines(gmat, psifac_dummy) umat_p1 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p1.dat")) umat_p2 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p2.dat")) @@ -132,20 +132,20 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex # Only the six matrices sing_der! reads are physical here; the rest are unused placeholders. unused = cubic_interp(psifac_dummy, Series(zeros(ComplexF64, points, intr.numpert_total^2)); itp_opts...) ideal = GeneralizedPerturbedEquilibrium.ForceFreeStates.IdealMatrices(; - amats=cubic_interp(psifac_dummy, Series(reshape(amats, points, :)); itp_opts...), - bmats=cubic_interp(psifac_dummy, Series(reshape(bmats, points, :)); itp_opts...), - cmats=cubic_interp(psifac_dummy, Series(reshape(cmats, points, :)); itp_opts...), - fmats_lower=cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); itp_opts...), - kmats=cubic_interp(psifac_dummy, Series(reshape(kmats, points, :)); itp_opts...), - gmats=cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); itp_opts...), - dmats_prim=unused, emats_prim=unused, hmats=unused, fmats_prim=unused, - fmats_gal=unused, jmats=unused) - ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.FourFitVars(; + A_spline=cubic_interp(psifac_dummy, Series(reshape(A_spline, points, :)); itp_opts...), + B_spline=cubic_interp(psifac_dummy, Series(reshape(B_spline, points, :)); itp_opts...), + C_spline=cubic_interp(psifac_dummy, Series(reshape(C_spline, points, :)); itp_opts...), + F_spline_lower=cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); itp_opts...), + K_spline=cubic_interp(psifac_dummy, Series(reshape(K_spline, points, :)); itp_opts...), + G_spline=cubic_interp(psifac_dummy, Series(reshape(G_spline, points, :)); itp_opts...), + D_spline_prim=unused, E_spline_prim=unused, H_spline=unused, F_spline_prim=unused, + F_spline_gal=unused, J_spline=unused) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.MatrixSplines(; numpert_total=intr.numpert_total, itp_opts, ideal) du = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) chunk = GeneralizedPerturbedEquilibrium.ForceFreeStates.IntegrationChunk(; psi_start=odet.psifac, psi_end=odet.psifac, needs_crossing=false) - params = (ctrl, equil, ffit, intr, odet, chunk) + params = (ctrl, equil, mats, intr, odet, chunk) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_der!(du, odet.u, params, odet.psifac) du_fortran = read_solutions_3d(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/sing_der_output_du.dat")) From 01a5c79e63d651375f2ca4a815a52cb7b1e902dc Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 18 Aug 2026 10:51:28 -0400 Subject: [PATCH 5/6] FFS - MINOR - removing unnecessary struct members from MatrixSplines --- src/ForceFreeStates/FixedKineticMatrices.jl | 7 +++--- src/ForceFreeStates/Fourfit.jl | 8 ++----- src/ForceFreeStates/Kinetic.jl | 22 +++++++----------- .../CalculatedKineticMatrices.jl | 2 -- test/runtests_sing.jl | 23 +++++++++---------- 5 files changed, 24 insertions(+), 38 deletions(-) diff --git a/src/ForceFreeStates/FixedKineticMatrices.jl b/src/ForceFreeStates/FixedKineticMatrices.jl index a854e2283..41f1f2245 100644 --- a/src/ForceFreeStates/FixedKineticMatrices.jl +++ b/src/ForceFreeStates/FixedKineticMatrices.jl @@ -43,7 +43,7 @@ function _build_x_matrix(mpert::Int, mlow::Int, sigma::Float64; hermitian::Bool= end """ - fixed_kinetic_matrices(mpert, mpsi, sigma, mlow, mats, xs) + fixed_kinetic_matrices(mpert, np, mpsi, sigma, mlow, mats, xs) Build X-shaped fixed kinetic energy matrices for testing all 6 components. @@ -69,10 +69,9 @@ Torque matrices (T) are all zero (torque requires finite rotation frequency). Returns `(kw_flat, kt_flat)` where each is `(mpsi, mpert^2, 6)`. """ function fixed_kinetic_matrices( - mpert::Int, mpsi::Int, sigma::Float64, mlow::Int, + mpert::Int, np::Int, mpsi::Int, sigma::Float64, mlow::Int, mats::MatrixSplines, xs::Vector{Float64} ) - np = mats.numpert_total kw_flat = zeros(ComplexF64, mpsi, np^2, 6) kt_flat = zeros(ComplexF64, mpsi, np^2, 6) @@ -100,7 +99,7 @@ function fixed_kinetic_matrices( # Scale: σ × ‖ideal(ψ)‖_F × unit X-pattern # For multi-n, tile the mpert×mpert X-pattern into the np×np block W = zeros(ComplexF64, np, np) - for jn in 0:(mats.numpert_total÷mpert-1) + for jn in 0:(np÷mpert-1) offset = jn * mpert W[(offset+1):(offset+mpert), (offset+1):(offset+mpert)] .= X end diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 86bf52752..98c677eaf 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -99,16 +99,12 @@ and for the Galerkin/vacuum paths that are defined only in the ideal basis. ## Fields - - `numpert_total::Int` - `mpert · npert`; each matrix carries `numpert_total^2` series. - - `itp_opts::NamedTuple` - interpolant options applied to every spline built for this fit. - `ideal::IdealMatrices` - always present. - `kinetic::Union{Nothing,KineticMatrices}` - present iff `ctrl.kinetic_factor > 0`. - `_hint::Base.RefValue{Int}` - shared bracket-search hint for sequential (single-threaded) evaluation. Not thread-safe: parallel paths must pass their own hint. """ -@kwdef struct MatrixSplines{S<:CubicSeriesInterpolant,Opts<:NamedTuple} - numpert_total::Int - itp_opts::Opts = (; extrap=ExtendExtrap()) +@kwdef struct MatrixSplines{S<:CubicSeriesInterpolant} ideal::IdealMatrices{S} kinetic::Union{Nothing,KineticMatrices{S}} = nothing _hint::Base.RefValue{Int} = Ref(1) @@ -633,5 +629,5 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates K_spline=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl J_spline=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) - return MatrixSplines(; numpert_total=intr.numpert_total, itp_opts, ideal) + return MatrixSplines(; ideal) end diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index f2d281741..c092fe664 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -34,7 +34,7 @@ function make_kinetic_matrix( # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" - kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs) + kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs) elseif ctrl.kinetic_source == "calculated" isnothing(calculated_source) && error( "kinetic_source=\"calculated\" requires the KineticForces callback. " * @@ -49,16 +49,12 @@ function make_kinetic_matrix( error("Unknown kinetic_source: $(ctrl.kinetic_source). Must be \"fixed\" or \"calculated\"") end - # Build splines for each of the 6 components - Kw_spline = [cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); mats.itp_opts...) for ic in 1:6] - Kt_spline = [cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); mats.itp_opts...) for ic in 1:6] - # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat, Kw_spline, Kt_spline) + return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat) end """ - _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat, Kw_spline, Kt_spline) -> MatrixSplines + _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat) -> MatrixSplines Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and return a new `MatrixSplines` whose `kinetic` field holds them alongside the kinetic-modified A/B/C; `mats.ideal` is carried over @@ -77,9 +73,7 @@ function _compute_fkg_matrices( intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, - kt_flat::Array{ComplexF64,3}, - Kw_spline::Vector, - Kt_spline::Vector + kt_flat::Array{ComplexF64,3} ) xs = metric.xs mpsi = length(xs) @@ -239,14 +233,14 @@ function _compute_fkg_matrices( end end - itp_opts = mats.itp_opts + itp_opts = (; extrap=ExtendExtrap()) kinetic = KineticMatrices(; # kinetic-modified A/B/C consumed by sing_der!; A is non-Hermitian here A_spline=cubic_interp(xs, Series(ak_flat); itp_opts...), B_spline=cubic_interp(xs, Series(bk_flat); itp_opts...), C_spline=cubic_interp(xs, Series(ck_flat); itp_opts...), - Kw_spline, - Kt_spline, + Kw_spline=[cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); itp_opts...) for ic in 1:6], + Kt_spline=[cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); itp_opts...) for ic in 1:6], F0_spline=cubic_interp(xs, Series(f0_flat); itp_opts...), P_spline=cubic_interp(xs, Series(p_flat); itp_opts...), P_spline_adj=cubic_interp(xs, Series(pa_flat); itp_opts...), @@ -257,5 +251,5 @@ function _compute_fkg_matrices( R3_spline=cubic_interp(xs, Series(r3_flat); itp_opts...), G_spline_adj=cubic_interp(xs, Series(ga_flat); itp_opts...)) - return MatrixSplines(mats.numpert_total, itp_opts, mats.ideal, kinetic, mats._hint) + return MatrixSplines(mats.ideal, kinetic, mats._hint) end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 455e515b2..50532341a 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -66,8 +66,6 @@ function compute_calculated_kinetic_matrices( npert = ffs_intr.npert np = ffs_intr.numpert_total - @assert mats.numpert_total == np "MatrixSplines and ForceFreeStatesInternal disagree on numpert_total" - kw_flat = zeros(ComplexF64, mpsi, np^2, 6) kt_flat = zeros(ComplexF64, mpsi, np^2, 6) diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index ddf42311f..d58d11c2b 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -108,19 +108,19 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex psifac_dummy = collect(range(0, 1, 10)); points = length(psifac_dummy) amat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/amat.dat")) - A_spline = copyForSplines(amat, psifac_dummy) + amats = copyForSplines(amat, psifac_dummy) bmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/bmat.dat")) - B_spline = copyForSplines(bmat, psifac_dummy) + bmats = copyForSplines(bmat, psifac_dummy) cmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/cmat.dat")) - C_spline = copyForSplines(cmat, psifac_dummy) + cmats = copyForSplines(cmat, psifac_dummy) fmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/fmat.dat")) fmat .= cholesky(Hermitian(fmat)).L; fmats = copyForSplines(fmat, psifac_dummy) kmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/kmat.dat")) - K_spline = copyForSplines(kmat, psifac_dummy) + kmats = copyForSplines(kmat, psifac_dummy) gmat = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/gmat.dat")) - G_spline = copyForSplines(gmat, psifac_dummy) + gmats = copyForSplines(gmat, psifac_dummy) umat_p1 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p1.dat")) umat_p2 = read_complex_fortran(joinpath(@__DIR__, "test_data/sing_der_testing/mat_dat/umat_p2.dat")) @@ -132,16 +132,15 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex # Only the six matrices sing_der! reads are physical here; the rest are unused placeholders. unused = cubic_interp(psifac_dummy, Series(zeros(ComplexF64, points, intr.numpert_total^2)); itp_opts...) ideal = GeneralizedPerturbedEquilibrium.ForceFreeStates.IdealMatrices(; - A_spline=cubic_interp(psifac_dummy, Series(reshape(A_spline, points, :)); itp_opts...), - B_spline=cubic_interp(psifac_dummy, Series(reshape(B_spline, points, :)); itp_opts...), - C_spline=cubic_interp(psifac_dummy, Series(reshape(C_spline, points, :)); itp_opts...), + A_spline=cubic_interp(psifac_dummy, Series(reshape(amats, points, :)); itp_opts...), + B_spline=cubic_interp(psifac_dummy, Series(reshape(bmats, points, :)); itp_opts...), + C_spline=cubic_interp(psifac_dummy, Series(reshape(cmats, points, :)); itp_opts...), F_spline_lower=cubic_interp(psifac_dummy, Series(reshape(fmats, points, :)); itp_opts...), - K_spline=cubic_interp(psifac_dummy, Series(reshape(K_spline, points, :)); itp_opts...), - G_spline=cubic_interp(psifac_dummy, Series(reshape(G_spline, points, :)); itp_opts...), + K_spline=cubic_interp(psifac_dummy, Series(reshape(kmats, points, :)); itp_opts...), + G_spline=cubic_interp(psifac_dummy, Series(reshape(gmats, points, :)); itp_opts...), D_spline_prim=unused, E_spline_prim=unused, H_spline=unused, F_spline_prim=unused, F_spline_gal=unused, J_spline=unused) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.MatrixSplines(; - numpert_total=intr.numpert_total, itp_opts, ideal) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.MatrixSplines(; ideal) du = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) chunk = GeneralizedPerturbedEquilibrium.ForceFreeStates.IntegrationChunk(; psi_start=odet.psifac, psi_end=odet.psifac, needs_crossing=false) From bc15b720e3e82b61eb6074fbd2e618b4c84ace71 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 18 Aug 2026 11:11:39 -0400 Subject: [PATCH 6/6] FFS - MINOR - renaming build_matrix functions to build_matrix_splines --- benchmarks/benchmark_delta_prime_methods.jl | 2 +- benchmarks/benchmark_riccati_der.jl | 2 +- benchmarks/benchmark_threads.jl | 2 +- docs/src/stability.md | 2 +- src/ForceFreeStates/FixedKineticMatrices.jl | 2 +- src/ForceFreeStates/Fourfit.jl | 16 ++++------------ src/ForceFreeStates/Kinetic.jl | 6 +++--- src/ForceFreeStates/Surfaces/Finding.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 6 +++--- src/KineticForces/CalculatedKineticMatrices.jl | 2 +- src/KineticForces/KineticForces.jl | 2 +- test/runtests_eulerlagrange.jl | 2 +- test/runtests_parallel_integration.jl | 14 +++++++------- test/runtests_riccati.jl | 2 +- 14 files changed, 27 insertions(+), 35 deletions(-) diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl index 04bb929eb..9390c1eb1 100644 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ b/benchmarks/benchmark_delta_prime_methods.jl @@ -41,7 +41,7 @@ function setup_and_run_solovev() intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - mats = FFS.make_matrix(equil, intr, metric) + mats = FFS.build_matrix_splines(equil, intr, metric) odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr) return ctrl, equil, mats, intr, odet end diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl index 4b9740a24..e9d88ef82 100644 --- a/benchmarks/benchmark_riccati_der.jl +++ b/benchmarks/benchmark_riccati_der.jl @@ -42,7 +42,7 @@ function setup_solovev() intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - mats = FFS.make_matrix(equil, intr, metric) + mats = FFS.build_matrix_splines(equil, intr, metric) return ctrl, equil, mats, intr end diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 12661706a..1b2912a70 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -30,7 +30,7 @@ function run_ffs(ex; integrator) intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr.numpert_total diff --git a/docs/src/stability.md b/docs/src/stability.md index 81379c4fd..84b7e177d 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -324,7 +324,7 @@ intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) -mats = FFS.make_matrix(equil, intr, metric) +mats = FFS.build_matrix_splines(equil, intr, metric) # Choose integration driver. The top-level `eulerlagrange_integration` dispatches # on ctrl.integrator and always returns a 4-tuple diff --git a/src/ForceFreeStates/FixedKineticMatrices.jl b/src/ForceFreeStates/FixedKineticMatrices.jl index 41f1f2245..92891961e 100644 --- a/src/ForceFreeStates/FixedKineticMatrices.jl +++ b/src/ForceFreeStates/FixedKineticMatrices.jl @@ -2,7 +2,7 @@ FixedKineticMatrices Test fixtures: synthetic X-shaped kinetic energy matrices keyed off the ideal -matrices' Frobenius norms. Used by `make_kinetic_matrix` when +matrices' Frobenius norms. Used by `build_kinetic_matrix_splines` when `ctrl.kinetic_source == "fixed"` to exercise the kinetic-MHD code path before the calculated NTV pipeline (see `KineticForces.compute_calculated_kinetic_matrices`) is wired in. diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 98c677eaf..f5d188abe 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -1,7 +1,7 @@ """ IdealMatrices -Ideal-MHD stability matrix ψ-splines assembled by [`make_matrix`](@ref). Each field flattens a +Ideal-MHD stability matrix ψ-splines assembled by [`build_matrix_splines`](@ref). Each field flattens a `numpert_total × numpert_total` matrix to `numpert_total^2` complex series (`J_spline` to `2·mpert−1`), following the appendix of Glasser Phys. Plasmas 2016 112506. @@ -47,7 +47,7 @@ end """ KineticMatrices -Kinetic-MHD matrix ψ-splines assembled by [`make_kinetic_matrix`](@ref). Present on a +Kinetic-MHD matrix ψ-splines assembled by [`build_kinetic_matrix_splines`](@ref). Present on a [`MatrixSplines`](@ref) only for kinetic runs; its presence means the solution obeys the FKG ODE rather than the ideal Euler-Lagrange relation. @@ -429,7 +429,7 @@ end """ - make_matrix(metric::MetricData, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) -> MatrixSplines + build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) -> MatrixSplines Constructs main ForceFreeStates matrices for a given toroidal mode number and returns them as a new `MatrixSplines` object. See the appendix of the Glasser Phys. Plasmas 2016 112506 @@ -447,13 +447,8 @@ later (i.e. `sing_der!`). ### Returns - `mats::MatrixSplines`: A struct holding cubic spline fits of the assembled matrices - -### TODOs - -Add kinetic metric tensor components for kinetic mode -Set powers if necessary """ -function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) +function build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) # --- Extract inputs --- profiles = equil.profiles @@ -604,10 +599,7 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates # Store factorized F matrix (lower triangular only) since we always will need F⁻¹ later # and this make computation more efficient via combined forward and back substitution - # TODO: does F stay Hermitian in the 3D case, allowing us to use the lower representation? fmat .= cholesky(Hermitian(fmat)).L - - # TODO: add kinetic matrices here end # --- Create Fourier coefficient splines (multi-quantity cubic interpolants) --- diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index c092fe664..0cd8e01dd 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,5 +1,5 @@ """ - make_kinetic_matrix(ctrl, equil, mats, intr, metric; + build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; calculated_source=nothing) Construct kinetic energy (W) and torque (T) matrices and pre-compute the FKG derived @@ -21,7 +21,7 @@ Dispatches on `ctrl.kinetic_source`: Both paths apply `ctrl.kinetic_factor` as a global scale before the FKG Schur reduction. """ -function make_kinetic_matrix( +function build_kinetic_matrix_splines( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, @@ -39,7 +39,7 @@ function make_kinetic_matrix( isnothing(calculated_source) && error( "kinetic_source=\"calculated\" requires the KineticForces callback. " * "Drive the run via `GeneralizedPerturbedEquilibrium.main` instead of " * - "calling make_kinetic_matrix directly, or pass " * + "calling build_kinetic_matrix_splines directly, or pass " * "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly." ) kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 5154456ce..2c7bdac14 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -259,7 +259,7 @@ Algorithm: """ function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) kin = mats.kinetic - kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call make_kinetic_matrix first") + kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call build_kinetic_matrix_splines first") psilow = equil.profiles.xs[1] psihigh = intr.psilim diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 1716454ae..26c5ef6e0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -73,7 +73,7 @@ include("Rerun.jl") using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings using .ForceFreeStates: ForceFreeStatesResult, build_result using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry -using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix +using .ForceFreeStates: make_metric, build_matrix_splines, build_kinetic_matrix_splines using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! using .ForceFreeStates: galerkin_solve, write_galerkin! @@ -559,7 +559,7 @@ function prepare_force_free_states!( end # Compute matrices and build the MatrixSplines container - mats = make_matrix(equil, intr, metric) + mats = build_matrix_splines(equil, intr, metric) if ctrl.kinetic_factor > 0 if ctrl.verbose @@ -572,7 +572,7 @@ function prepare_force_free_states!( KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles) - mats = make_kinetic_matrix(ctrl, equil, mats, intr, metric; + mats = build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; calculated_source=calculated_cb) # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings. diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 50532341a..51ebd2ef6 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -2,7 +2,7 @@ CalculatedKineticMatrices Bridge from KineticForces' bounce-averaged matrix kernels into the kinetic-MHD -stability path in ForceFreeStates. Used by `make_kinetic_matrix` when +stability path in ForceFreeStates. Used by `build_kinetic_matrix_splines` when `ctrl.kinetic_source == "calculated"` (via the `calculated_source` callback injected from `GeneralizedPerturbedEquilibrium.main`). """ diff --git a/src/KineticForces/KineticForces.jl b/src/KineticForces/KineticForces.jl index 1ac79ad22..cd27590ed 100644 --- a/src/KineticForces/KineticForces.jl +++ b/src/KineticForces/KineticForces.jl @@ -17,7 +17,7 @@ to MHD stability through torque and energy deposition calculations. - `Compute.jl`: Orchestration routines - compute_torque_all_methods!() - `CalculatedKineticMatrices.jl`: compute_calculated_kinetic_matrices() - ← Callback registered with ForceFreeStates.make_kinetic_matrix when + ← Callback registered with ForceFreeStates.build_kinetic_matrix_splines when kinetic_source="calculated" ### Supporting Functions diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index ed12ef89c..4da7402f8 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -450,7 +450,7 @@ end intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) - mats = FFS.make_matrix(equil, intr, metric) + mats = FFS.build_matrix_splines(equil, intr, metric) odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, mats, intr) return odet, ctrl, equil, mats, intr end diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 9c39ee78a..b427c3795 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -126,7 +126,7 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) @@ -213,7 +213,7 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) @@ -275,7 +275,7 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr @@ -334,7 +334,7 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr @@ -384,7 +384,7 @@ using TOML intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) # Use the first chunk from chunk_el_integration_bounds: guaranteed rational-free interior odet_tmp = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, 10, 5, intr.msing) @@ -436,7 +436,7 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) # Derivatives are recomputed on demand; materialize so the stores can be compared. GeneralizedPerturbedEquilibrium.ForceFreeStates.materialize_derivative_stores!(odet, equil, mats, intr) @@ -498,7 +498,7 @@ using TOML intr.mpert = intr.mhigh - intr.mlow + 1 intr.numpert_total = intr.mpert * intr.npert metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert) - mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric) + mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric) odet, fm_propagators, fm_chunks, fm_S_left = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index 1f4dd7dc5..e5989aed8 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -107,7 +107,7 @@ end intr_tmp = make_solovev_intr(inputs, ctrl, equil, ex) metric = FFS.make_metric(equil, intr_tmp.mpert) - mats = FFS.make_matrix(equil, intr_tmp, metric) + mats = FFS.build_matrix_splines(equil, intr_tmp, metric) N = intr_tmp.numpert_total # Riccati integration. The driver returns (odet, propagators, chunks, S_at_surface_left);