diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 515186317..93df52d43 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -43,7 +43,7 @@ Top level (11 groups): | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` | | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` | -| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra | +| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/` (displacement and field spectra, plus the `dBoB_lagrangian` and `div_xi_perp` spectra that drive neoclassical transport), `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies///` with the same per-method layout, summing to the top-level total | | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) | | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_/` | diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index c04836766..14919aa46 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -1123,7 +1123,7 @@ function run_kinetic_forces( else # kf_ctrl and kinetic_profiles were loaded once before the equilibrium was re-formed. kf_intr = KineticForces.KineticForcesInternal(result.equil; verbose=kf_ctrl.verbose) - KineticForces.set_perturbation_data!(kf_intr, pe_state, result, result.equil, result.metric) + KineticForces.set_perturbation_data!(kf_intr, pe_state, result) if species === nothing kf_state = KineticForces.KineticForcesState() @@ -1376,7 +1376,7 @@ function efc_couplings( pe_intr.forcing_modes = modes pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(ffs, ForcingTerms.RMPField(ForcingTerms.ForcingTermsControl()), quiet, pe_intr) kf_intr = KineticForces.KineticForcesInternal(equil; verbose=verbose) - KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric) + KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs) return kf_intr end # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation. diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 9014dd9a4..2711ab9ef 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -319,27 +319,22 @@ function KineticForcesInternal(equil; verbose::Bool=false) end """ - set_perturbation_data!(kf_intr, pe_state, ffs, equil, metric) + set_perturbation_data!(kf_intr, pe_state, ffs) Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal. -Builds three interpolant sets from PE Clebsch displacements: +Builds three interpolant sets over ψ: - 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ - 2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) - 3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting + 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] from the PE Clebsch displacements + 2. `dbob_m` — Lagrangian δ|B|/B, taken from `pe_state.modb_modes` + 3. `divx_m` — ∇·ξ⊥, taken from `pe_state.modb_modes` -The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): - - 1. Apply geometric matrices S,T,X,Y,Z in m-space - 2. Inverse DFT to θ-space - 3. Divide by J·B² at each θ - 4. Forward DFT back to m-space +The perturbation spectra themselves are a property of the perturbed equilibrium and are +derived and written there (`PerturbedEquilibrium.compute_modb_modes`); this function only +interpolates them onto the torque quadrature. """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::PerturbedEquilibrium.PerturbedEquilibriumState, - ffs::ForceFreeStates.ForceFreeStatesResult, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + ffs::ForceFreeStates.ForceFreeStatesResult) ffs.npert == 1 || error("KineticForces supports a single toroidal mode; this run has n = $(ffs.nlow):$(ffs.nhigh).") # Copy mode numbers from FFS @@ -378,72 +373,19 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_mat); itp_opts...) kf_intr.xs_m = [xs_m_1, xs_m_2, xs_m_3] - # Build geometric matrices (S,T,X,Y,Z) for JBB deweighting - geom_mats = ForceFreeStates.build_kinetic_metric_matrices(equil, ffs, metric) - - # Build FourierTransform for the JBB deweighting DFT round-trip - mthsurf = kf_intr.mthsurf - ft = Utilities.FourierTransforms.FourierTransform(mthsurf, mpert, ffs.mlow) - - # JBB deweighting: convert Clebsch modes → physical δB/B and ∇·ξ⊥ modes - dbob_m_data = zeros(ComplexF64, npsi, mpert) - divx_m_data = zeros(ComplexF64, npsi, mpert) - - # Pre-allocate buffers - smat_flat = Vector{ComplexF64}(undef, mpert^2) - tmat_flat = Vector{ComplexF64}(undef, mpert^2) - xmat_flat = Vector{ComplexF64}(undef, mpert^2) - ymat_flat = Vector{ComplexF64}(undef, mpert^2) - zmat_flat = Vector{ComplexF64}(undef, mpert^2) - jbb_kapx = Vector{ComplexF64}(undef, mpert) - jbb_divx = Vector{ComplexF64}(undef, mpert) - jbb_dbob = Vector{ComplexF64}(undef, mpert) - theta_buf = Vector{ComplexF64}(undef, mthsurf) - - hint_s = Ref(1) - hint_t = Ref(1) - hint_x = Ref(1) - hint_y = Ref(1) - hint_z = Ref(1) - - for ipsi in 1:npsi - psi = psi_grid[ipsi] - - # Get Clebsch displacement vectors at this ψ - xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] - xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] - - # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert - geom_mats.smats(smat_flat, psi; hint=hint_s) - geom_mats.tmats(tmat_flat, psi; hint=hint_t) - geom_mats.xmats(xmat_flat, psi; hint=hint_x) - geom_mats.ymats(ymat_flat, psi; hint=hint_y) - geom_mats.zmats(zmat_flat, psi; hint=hint_z) - - smat = reshape(smat_flat, mpert, mpert) - tmat = reshape(tmat_flat, mpert, mpert) - xmat = reshape(xmat_flat, mpert, mpert) - ymat = reshape(ymat_flat, mpert, mpert) - zmat = reshape(zmat_flat, mpert, mpert) - - # Apply geometric matrices in m-space. xs_m ordering: (1)=∂ξ^ψ/∂ψ, (2)=ξ^ψ, (3)=ξ^α. - # jbb_kapx = smat·xsp + tmat·xms; jbb_divx = xmat·xmp1 + ymat·xsp + zmat·xms - mul!(jbb_kapx, smat, xsp) - mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += tmat * xms - mul!(jbb_divx, xmat, xmp1) - mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp - mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms - @. jbb_dbob = -(jbb_divx + jbb_kapx) - - # Inverse DFT to θ-space, divide by J·B², forward DFT back - _jbb_deweight!(view(dbob_m_data, ipsi, :), jbb_dbob, ft, psi, equil, mthsurf, theta_buf) - _jbb_deweight!(view(divx_m_data, ipsi, :), jbb_divx, ft, psi, equil, mthsurf, theta_buf) - end + # δB/B and ∇·ξ⊥ come from PerturbedEquilibrium, which derives them from the same + # displacement (`compute_modb_modes`, Fortran gpout_pmodb / set_peq) and writes them to + # `PerturbedEquilibrium/Response/`. Nothing kinetic enters them, so they are not rebuilt + # here; a run whose npert > 1 has already failed the single-mode check above. + modb = pe_state.modb_modes + modb === nothing && error("set_perturbation_data!: the perturbed equilibrium carries no " * + "δ|B|/B spectra; KineticForces cannot evaluate the torque without them.") + size(modb.dbob) == (npsi, mpert) || error("set_perturbation_data!: δ|B|/B spectra are " * + "$(size(modb.dbob)), expected ($npsi, $mpert).") # Build CubicSeriesInterpolants over ψ for dbob_m and divx_m - kf_intr.dbob_m = cubic_interp(psi_grid, Series(dbob_m_data); itp_opts...) - kf_intr.divx_m = cubic_interp(psi_grid, Series(divx_m_data); itp_opts...) + kf_intr.dbob_m = cubic_interp(psi_grid, Series(modb.dbob); itp_opts...) + kf_intr.divx_m = cubic_interp(psi_grid, Series(modb.divx); itp_opts...) if kf_intr.verbose @info "set_perturbation_data!: built dbob_m/divx_m/xs_m interpolants " * @@ -451,34 +393,7 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur end end -""" - _jbb_deweight!(out, jbb_modes, ft, psi, equil, mthsurf, theta_buf) - -JBB deweighting step: inverse DFT → divide by J(ψ,θ)·B(ψ,θ)² → forward DFT. -Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data -to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back. -""" -function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, - ft::Utilities.FourierTransforms.FourierTransform, - psi::Float64, equil::Equilibrium.PlasmaEquilibrium, - mthsurf::Int, theta_buf::Vector{ComplexF64}) - # Inverse DFT: m-space → θ-space - theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) - - # Divide by J·B² at each θ point - for j in 1:mthsurf - theta_norm = (j - 1) / mthsurf # θ ∈ [0, 1) - pt = (psi, theta_norm) - jac = equil.rzphi_jac(pt) - B = equil.eqfun_B(pt) - theta_buf[j] /= jac * B^2 - end - - # Forward DFT: θ-space → m-space - out .= ft(theta_buf) - return nothing -end # ============================================================================ diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index c4f073bc1..b27b4c927 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -1056,3 +1056,132 @@ function _apply_rzphi_transform( return R_modes, Z_modes, phi_modes end + +""" + compute_modb_modes(xi_modes, psi_grid, equil, ffs, metric) -> (; dbob, divx) + +Lagrangian δ|B|/B and ∇·ξ_⊥ Fourier modes, `[npsi, mpert]`, on `psi_grid`. + +These are the two spectra a neoclassical-transport calculation is driven by: NTV is a +functional of the Lagrangian perturbation of the field strength, not of the normal field. +They are a property of the perturbed equilibrium alone — no kinetic quantity enters — so they +are built here and merely consumed by `KineticForces`, which is what lets another NTV code be +run on GPEC's perturbation without running ours. The Fortran parent writes them from the same +stage (`gpec` `gpout_pmodb` under `pmodb_flag`, read back by `pentrc` `read_pmodb`); the +algorithm below is the one `pentrc` `set_peq` uses on the Clebsch displacements. + +The JBB deweighting (Fortran `pentrc/inputs.f90:828-868`): + + 1. apply the geometric matrices S,T,X,Y,Z in m-space, giving the J·B²-weighted κ·ξ and ∇·ξ_⊥ + 2. inverse DFT to θ-space + 3. divide by J·B² at each θ + 4. forward DFT back to m-space + +Only the two independent spectra are returned. The curvature term follows from +δB/B = −(∇·ξ_⊥ + κ·ξ), and the Eulerian δ|B|/B differs from the Lagrangian by ξ·∇|B|/|B|, +which `Response/xi_clebsch_psi`, `Response/xi_clebsch_alpha` and the equilibrium geometry +already determine — so neither is written. + +Single toroidal mode only: `build_kinetic_metric_matrices` assembles mpert×mpert matrices at +`n = ffs.nlow`. +""" +function compute_modb_modes(xi_modes::NamedTuple, psi_grid::Vector{Float64}, + equil::Equilibrium.PlasmaEquilibrium, + ffs::ForceFreeStatesResult, + metric::MetricData +) + npsi = length(psi_grid) + mpert = ffs.mpert + + # Geometric matrices (S,T,X,Y,Z) of the JBB-weighted Clebsch-to-physical transformation. + geom_mats = ForceFreeStates.build_kinetic_metric_matrices(equil, ffs, metric) + + # The deweighting round-trip runs on the equilibrium's own poloidal grid. + mthsurf = length(equil.rzphi_ys) - 1 + ft = Utilities.FourierTransforms.FourierTransform(mthsurf, mpert, ffs.mlow) + + dbob = zeros(ComplexF64, npsi, mpert) + divx = zeros(ComplexF64, npsi, mpert) + + smat_flat = Vector{ComplexF64}(undef, mpert^2) + tmat_flat = Vector{ComplexF64}(undef, mpert^2) + xmat_flat = Vector{ComplexF64}(undef, mpert^2) + ymat_flat = Vector{ComplexF64}(undef, mpert^2) + zmat_flat = Vector{ComplexF64}(undef, mpert^2) + jbb_kapx = Vector{ComplexF64}(undef, mpert) + jbb_divx = Vector{ComplexF64}(undef, mpert) + jbb_dbob = Vector{ComplexF64}(undef, mpert) + theta_buf = Vector{ComplexF64}(undef, mthsurf) + + # Serial over ψ: the rows are independent, but the sticky bracket hints below are shared + # mutable state, and a threaded reassociation here would make the spectra depend on the + # thread count. + hint_s = Ref(1) + hint_t = Ref(1) + hint_x = Ref(1) + hint_y = Ref(1) + hint_z = Ref(1) + + for ipsi in 1:npsi + psi = psi_grid[ipsi] + + xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] + xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] + xms = view(xi_modes.clebsch_alpha, ipsi, :) # ξ^α/χ₁ [mpert] + + geom_mats.smats(smat_flat, psi; hint=hint_s) + geom_mats.tmats(tmat_flat, psi; hint=hint_t) + geom_mats.xmats(xmat_flat, psi; hint=hint_x) + geom_mats.ymats(ymat_flat, psi; hint=hint_y) + geom_mats.zmats(zmat_flat, psi; hint=hint_z) + + smat = reshape(smat_flat, mpert, mpert) + tmat = reshape(tmat_flat, mpert, mpert) + xmat = reshape(xmat_flat, mpert, mpert) + ymat = reshape(ymat_flat, mpert, mpert) + zmat = reshape(zmat_flat, mpert, mpert) + + # jbb_kapx = smat·xsp + tmat·xms; jbb_divx = xmat·xmp1 + ymat·xsp + zmat·xms + mul!(jbb_kapx, smat, xsp) + mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) + mul!(jbb_divx, xmat, xmp1) + mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) + mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) + @. jbb_dbob = -(jbb_divx + jbb_kapx) + + _jbb_deweight!(view(dbob, ipsi, :), jbb_dbob, ft, psi, equil, mthsurf, theta_buf) + _jbb_deweight!(view(divx, ipsi, :), jbb_divx, ft, psi, equil, mthsurf, theta_buf) + end + + return (; dbob, divx) +end + +""" + _jbb_deweight!(out, jbb_modes, ft, psi, equil, mthsurf, theta_buf) + +JBB deweighting step: inverse DFT → divide by J(ψ,θ)·B(ψ,θ)² → forward DFT. + +Matches Fortran `set_peq` lines 859-868: transforms JBB-weighted m-space data to θ-space, +removes the J·B² weighting at each poloidal angle, and transforms back. +""" +function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, + ft::Utilities.FourierTransforms.FourierTransform, + psi::Float64, equil::Equilibrium.PlasmaEquilibrium, + mthsurf::Int, theta_buf::Vector{ComplexF64} +) + # Inverse DFT: m-space → θ-space + theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) + + # Divide by J·B² at each θ point + for j in 1:mthsurf + theta_norm = (j - 1) / mthsurf # θ ∈ [0, 1) + pt = (psi, theta_norm) + jac = equil.rzphi_jac(pt) + B = equil.eqfun_B(pt) + theta_buf[j] /= jac * B^2 + end + + # Forward DFT: θ-space → m-space + out .= ft(theta_buf) + return nothing +end diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 18cd783f3..9fb05b4eb 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -96,6 +96,8 @@ Response fields (mode space): - `b_modes::Union{Nothing, NamedTuple}` - Magnetic field; psi=b^ψ, b_psi_area_weighted=b^ψ/⟨J·|∇ψ|⟩_θ, theta/zeta=unregularized, theta_reg/zeta_reg=regularized [npsi, numpert_total] - `b_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal field b_n [npsi, numpert_total] - `xi_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal displacement xi_n [npsi, numpert_total] + - `modb_modes::Union{Nothing, NamedTuple}` - Perturbation spectra driving neoclassical transport; + dbob=Lagrangian δ|B|/B, divx=∇·ξ_⊥ [npsi, mpert]. Single-n runs only Coupling matrices [n_rational × numpert_total] — one row per resonant (surface, n) pair. Each row maps the full applied field to the resonant response at that surface. @@ -173,6 +175,8 @@ well-conditioned flux-space inductances L, Λ: b_modes::Union{Nothing,NamedTuple} = nothing b_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing # physical normal field b_n [npsi, numpert_total] xi_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing # physical normal displacement xi_n [npsi, numpert_total] + # Lagrangian δ|B|/B and ∇·ξ_⊥ spectra (see compute_modb_modes); the input an NTV code needs + modb_modes::Union{Nothing,NamedTuple} = nothing # Coupling matrices [n_rational × numpert_total] C_resonant_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index 2100e52c1..7a944e603 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -120,7 +120,21 @@ function compute_plasma_response!( state.b_n_modes = b_n_modes state.xi_n_modes = xi_n_modes + # Lagrangian δ|B|/B and ∇·ξ_⊥: the perturbation a neoclassical-transport calculation is + # driven by. Derived here, where the displacement is, rather than inside KineticForces -- + # it contains nothing kinetic, and another NTV code needs it without running ours. + # Restricted to a single toroidal mode by the geometric matrices it is built from. + modb_dt = 0.0 + if ffs.npert == 1 + modb_dt = @elapsed state.modb_modes = compute_modb_modes( + xi_modes, state.psi_grid, equil, ffs, metric + ) + end + if ctrl.verbose - @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" + @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = " * + "$(@sprintf("%.3e", maximum(abs.(response_flux))))" * + (state.modb_modes === nothing ? "; no δ|B|/B spectra (npert > 1)" : + "; δ|B|/B and ∇·ξ_⊥ spectra in $(@sprintf("%.2f", modb_dt)) s") end end diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index d2064aee4..7549517a8 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -70,7 +70,9 @@ PerturbedEquilibrium/ │ ├── b_n # Physical normal field b_n (ComplexF64 [npsi, numpert_total]) │ ├── xi_n # Physical normal displacement xi_n (ComplexF64 [npsi, numpert_total]) │ ├── Jb_theta -│ └── Jb_zeta +│ ├── Jb_zeta +│ ├── dBoB_lagrangian # Lagrangian δ|B|/B (ComplexF64 [npsi, numpert_total]) -- NTV input +│ └── div_xi_perp # ∇·ξ_⊥ (ComplexF64 [npsi, numpert_total]) -- NTV input ├── ResponseMatrices/ # [numpert_total × numpert_total], root-area-weighted field (b̃) space; R = S·A │ ├── plasma_inductance # Λ̃ = R⁻¹·Λ·R⁻† │ ├── surface_inductance # L̃ = R⁻¹·L·R⁻† @@ -166,6 +168,14 @@ function write_outputs_to_HDF5( response_group["xi_clebsch_alpha"] = state.xi_modes.clebsch_alpha end + # Lagrangian δ|B|/B and ∇·ξ_⊥ (gpout_pmodb): the perturbation a neoclassical-transport + # calculation is driven by, written whether or not this run also computes the torque. + # Empty for npert > 1, which the geometric matrices behind them do not support. + if !isnothing(state.modb_modes) + response_group["dBoB_lagrangian"] = state.modb_modes.dbob + response_group["div_xi_perp"] = state.modb_modes.divx + end + # Contravariant displacement (from gpeq_contra, all J-weighted) if have_xi response_group["Jxi_psi"] = state.xi_modes.psi_J @@ -291,6 +301,16 @@ const PE_H5_ANNOTATIONS = [ (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+regularization_width²))", dims=("psi", "mode"), attach=(1 => "Response/psi",)), "Response/xi_clebsch_alpha" => (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/dBoB_lagrangian" => + (; + long_name="Lagrangian perturbation of the field strength δ|B|/B (gpout_pmodb lagbpar/|B|), the NTV drive; " * + "δB/B(θ) = Σ_m value·exp(2πi·m·θ), θ ∈ [0,1). Eulerian form = this − ξ·∇|B|/|B| from xi_clebsch_*", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/div_xi_perp" => + (; + long_name="∇·ξ_⊥ (gpout_pmodb divxprp/|B|), the second NTV input; same angle and sign convention as " * + "dBoB_lagrangian, with κ·ξ = −(dBoB_lagrangian + div_xi_perp)", units="1/m", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), "Response/xi_Z" => (; long_name="cylindrical displacement component ξ_Z (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl index c743d6a7d..cecbc464a 100644 --- a/test/runtests_h5_schema.jl +++ b/test/runtests_h5_schema.jl @@ -105,3 +105,42 @@ end end end end + +@testset "PerturbedEquilibrium δ|B|/B output contract" begin + PE = GeneralizedPerturbedEquilibrium.PerturbedEquilibrium + npsi, mpert = 4, 3 + + # Written when the spectra exist. The NTV codes this serves read them as (m, ψ) — the + # transpose of the Julia shape — exactly like every other Response spectrum. + state = PE.PerturbedEquilibriumState(; + psi_grid=collect(range(0.1, 0.9; length=npsi)), + modb_modes=(dbob=ComplexF64[0.1i + 0.01im * j for i in 1:npsi, j in 1:mpert], + divx=ComplexF64[0.2i - 0.02im * j for i in 1:npsi, j in 1:mpert])) + f1 = tempname() * ".h5" + PE.write_outputs_to_HDF5(state, PE.PerturbedEquilibriumInternal(), f1) + HDF5.h5open(f1, "r") do f + g = f["PerturbedEquilibrium/Response"] + @test size(read(g["dBoB_lagrangian"])) == (npsi, mpert) + @test read(g["dBoB_lagrangian"]) ≈ state.modb_modes.dbob + @test read(g["div_xi_perp"]) ≈ state.modb_modes.divx + for d in ("dBoB_lagrangian", "div_xi_perp") + a = HDF5.attributes(g[d]) + @test haskey(a, "long_name") && haskey(a, "dims") + # Same ψ dimension scale as every other Response spectrum, which is what lets + # Analysis.PerturbedEquilibriumModes.modes_to_theta read these datasets unchanged. + # H5DS indexes file (row-major) dimensions, so Julia axis 1 of a 2-D dataset is + # C index 1 — see Utilities.HDF5Annotations.attach_scale!. + @test HDF5.API.h5ds_is_attached(g[d], g["psi"], 1) + end + end + rm(f1; force=true) + + # Absent, not empty, when the run could not produce them (npert > 1). + f2 = tempname() * ".h5" + PE.write_outputs_to_HDF5(PE.PerturbedEquilibriumState(), PE.PerturbedEquilibriumInternal(), f2) + HDF5.h5open(f2, "r") do f + @test !haskey(f["PerturbedEquilibrium/Response"], "dBoB_lagrangian") + @test !haskey(f["PerturbedEquilibrium/Response"], "div_xi_perp") + end + rm(f2; force=true) +end