PerturbedEquilibrium/KineticForces - FEATURE - Write the perturbation that drives neoclassical transport

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/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` 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_<k>/` |

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.

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,107 +373,27 @@ 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 " *
               "(npsi=$npsi, mpert=$mpert)"
     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
 
 
 # ============================================================================

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

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)

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