From 780de2f68037ba76dd228b4e1afc4cc7b1dd661e Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 1 Oct 2026 16:47:35 -0400 Subject: [PATCH 1/3] KineticForces - FEATURE - Per-bounce-harmonic torque density output The torque density is already computed harmonic by harmonic -- the l loop in integrate_psi_quadgk fills harm_vals and then reduces it -- but only the sum was ever kept. Multi-code NTV benchmarking wants the breakdown: the ITPA TC-24 activity compares contributions from l = 0, +/-1 ... +/-5 across codes, and the Fortran parent already emits dTdpsi_fgar(psi, ell), so Julia could not answer that request at all. KineticForces// gains dTdpsi_ell on a new ell dimension scale, shape (npsi, nell) in Julia and (nell, npsi) to a C-ordered reader. Summing over l reproduces dTdpsi. Three things needed care: - The per-l log is permuted by the same sortperm as psi. QuadGK visits nodes in refinement order, not in psi order, so an unsorted log would misalign with the grid silently rather than erroring. - harm_vals is reused at every psi node, so the log stores a copy. - combine_species_states interpolates onto a union grid, and the benchmark runs multi-species, so the summed total needs per-l as well as the per-species groups. The interpolation is now a shared helper used by both the total and each l column, and the combined dtdpsi is taken as the row sum of the per-l matrix, which makes the sum rule exact by construction rather than true to interpolation roundoff. A species without per-l data leaves the combined matrix empty instead of erroring. Quadrature adapts on the summed integrand, so a sharply resonant individual harmonic can be less well resolved than the total; the annotation says so. Profiles still come from the first toroidal mode only, so for npert > 1 the sum over l matches dtdpsi but neither integrates to total_torque. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FabuzB6rkkv638ZUCzaAA2 --- docs/development/hdf5-conventions.md | 2 +- src/KineticForces/Compute.jl | 75 +++++++++++++++++++---- src/KineticForces/KineticForcesStructs.jl | 4 ++ src/KineticForces/Output.jl | 6 ++ test/runtests_kinetic.jl | 1 + test/runtests_multiion.jl | 34 +++++++++- 6 files changed, 107 insertions(+), 15 deletions(-) diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c4631c984..858e99671 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -44,7 +44,7 @@ Top level (11 groups): | `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), `Kinetic/` | | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `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 | +| `KineticForces/` | `/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `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_/` | | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds | diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 62c55d7f1..ef33be14a 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -121,6 +121,7 @@ function integrate_psi_quadgk( logged_psi = Float64[] logged_dtdpsi = ComplexF64[] + logged_dtdpsi_ell = Vector{Vector{ComplexF64}}() logged_elems = is_matrix_method ? Vector{Array{ComplexF64,3}}() : nothing function psi_batch!(y::AbstractVector{ComplexF64}, x::AbstractVector) @@ -153,6 +154,8 @@ function integrate_psi_quadgk( push!(logged_psi, psi) push!(logged_dtdpsi, total) + # `harm_vals` is overwritten at every ψ node, so the per-ℓ log must own a copy. + push!(logged_dtdpsi_ell, copy(harm_vals)) if is_matrix_method elems_accum = zeros(ComplexF64, mpert, mpert, 6) for ell_idx in 1:nharm @@ -177,6 +180,11 @@ function integrate_psi_quadgk( perm = sortperm(logged_psi) sorted_psi = logged_psi[perm] sorted_dtdpsi = logged_dtdpsi[perm] + # Per-ℓ profile, permuted by the SAME `perm`: QuadGK visits nodes in refinement order, not in + # ψ order, so a per-ℓ log left unsorted would misalign with `sorted_psi` without any error. + # `reduce(hcat, ...)` builds (nharm, npsi); transpose so ψ is the first axis, as `dtdpsi` is. + sorted_dtdpsi_ell = isempty(logged_dtdpsi_ell) ? zeros(ComplexF64, 0, 0) : + permutedims(reduce(hcat, logged_dtdpsi_ell[perm])) # Cumulative trapezoidal integration for T(ψ) profile npts = length(sorted_psi) @@ -189,7 +197,9 @@ function integrate_psi_quadgk( end end - torque_profile = npts > 1 ? (psi=sorted_psi, dtdpsi=sorted_dtdpsi, t_cumulative=t_cumulative) : nothing + torque_profile = npts > 1 ? + (psi=sorted_psi, dtdpsi=sorted_dtdpsi, t_cumulative=t_cumulative, + dtdpsi_ell=sorted_dtdpsi_ell, ell=collect(-nl:nl)) : nothing # Trapezoidal integration of kinetic matrices over ψ (matrix methods only) matrix_integrated = nothing @@ -217,6 +227,28 @@ end # High-level orchestration # ============================================================================ +""" + _accumulate_interp!(acc, grid, xs, ys) + +Add the linear interpolant of `(xs, ys)` onto `grid` into `acc`, contributing zero outside +`[xs[1], xs[end]]`. `xs` must be sorted ascending. A duplicate node or the right endpoint takes +`ys[j]` directly rather than dividing by a zero interval. Shared by the total and per-ℓ profiles in +[`combine_species_states`](@ref) so both use identical interpolation. +""" +function _accumulate_interp!(acc, grid, xs::AbstractVector{Float64}, ys::AbstractVector{ComplexF64}) + for (k, ψ) in enumerate(grid) + (ψ < xs[1] || ψ > xs[end]) && continue + j = searchsortedlast(xs, ψ) + if j == length(xs) || xs[j+1] == xs[j] # endpoint or duplicate node: no interpolation + acc[k] += ys[j] + else + t = (ψ - xs[j]) / (xs[j+1] - xs[j]) + acc[k] += (1 - t) * ys[j] + t * ys[j+1] + end + end + return acc +end + """ combine_species_states(states) -> KineticForcesState @@ -246,19 +278,30 @@ function combine_species_states(states::AbstractVector{KineticForcesState}) for r in results length(r.psi_grid) >= 2 || continue o = sortperm(r.psi_grid) - xs = r.psi_grid[o] - ys = r.dtdpsi[o] - for (k, ψ) in enumerate(grid) - (ψ < xs[1] || ψ > xs[end]) && continue - j = searchsortedlast(xs, ψ) - if j == length(xs) || xs[j+1] == xs[j] # endpoint or duplicate node: no interpolation - dtdpsi[k] += ys[j] - else - t = (ψ - xs[j]) / (xs[j+1] - xs[j]) - dtdpsi[k] += (1 - t) * ys[j] + t * ys[j+1] + _accumulate_interp!(dtdpsi, grid, r.psi_grid[o], r.dtdpsi[o]) + end + + # Per-ℓ profile, interpolated onto the same union grid column by column. Only combined + # when every species carries it on an identical ℓ range; otherwise left empty. + ell = first(results).ell + has_ell = !isempty(ell) && all(r -> r.ell == ell && size(r.dtdpsi_ell, 2) == length(ell) && + size(r.dtdpsi_ell, 1) == length(r.psi_grid), results) + dtdpsi_ell = zeros(ComplexF64, has_ell ? length(grid) : 0, has_ell ? length(ell) : 0) + if has_ell + for r in results + length(r.psi_grid) >= 2 || continue + o = sortperm(r.psi_grid) + xs = r.psi_grid[o] + for j in eachindex(ell) + _accumulate_interp!(view(dtdpsi_ell, :, j), grid, xs, r.dtdpsi_ell[o, j]) end end + # Take the total from the per-ℓ row sum so Σ_ℓ ≡ dtdpsi holds exactly rather than to + # interpolation roundoff. Linear interpolation commutes with the sum, so this is the + # same quantity the loop above produced. + dtdpsi = vec(sum(dtdpsi_ell; dims=2)) end + tcum = zeros(ComplexF64, length(grid)) for j in 2:length(grid) tcum[j] = tcum[j-1] + 0.5 * (dtdpsi[j] + dtdpsi[j-1]) * (grid[j] - grid[j-1]) @@ -268,6 +311,7 @@ function combine_species_states(states::AbstractVector{KineticForcesState}) total_torque=sum(r.total_torque for r in results), total_energy=sum(r.total_energy for r in results), psi_grid=grid, dtdpsi=dtdpsi, t_cumulative=tcum, + dtdpsi_ell=dtdpsi_ell, ell=(has_ell ? ell : Int[]), psi_nsteps=sum(r.psi_nsteps for r in results)) end combined.completed = true @@ -322,6 +366,8 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor psi_grid_out = Float64[] dtdpsi_out = ComplexF64[] t_cum_out = ComplexF64[] + dtdpsi_ell_out = zeros(ComplexF64, 0, 0) + ell_out = Int[] psi_nsteps_total = 0 panel_psis_out = Float64[] resonance_psis_out = Float64[] @@ -349,6 +395,11 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor psi_grid_out = result.torque_profile.psi dtdpsi_out = result.torque_profile.dtdpsi t_cum_out = result.torque_profile.t_cumulative + # Profiles come from the first n only, while `total_torque` sums over all n. + # So sum(dtdpsi_ell; dims=2) == dtdpsi always, but neither integrates to + # `total_torque` when npert > 1. + dtdpsi_ell_out = result.torque_profile.dtdpsi_ell + ell_out = result.torque_profile.ell end end @@ -371,6 +422,8 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor psi_grid=psi_grid_out, dtdpsi=dtdpsi_out, t_cumulative=t_cum_out, + dtdpsi_ell=dtdpsi_ell_out, + ell=ell_out, psi_nsteps=psi_nsteps_total, panel_psis=panel_psis_out, resonance_psis=resonance_psis_out diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index ed9ef14a6..450bfb9d5 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -519,6 +519,10 @@ Results for one NTV computation method across all flux surfaces. psi_grid::Vector{Float64} = Float64[] dtdpsi::Vector{ComplexF64} = ComplexF64[] t_cumulative::Vector{ComplexF64} = ComplexF64[] + # dT/dψ resolved by bounce harmonic ℓ, shape (length(psi_grid), length(ell)) with ell = -nl:nl. + # Summing over ℓ reproduces `dtdpsi`. Empty when no ψ profile was produced. + dtdpsi_ell::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + ell::Vector{Int} = Int[] psi_nsteps::Int = 0 # ψ-quadrature panel boundaries and the located kinetic-resonance surfaces (first n) panel_psis::Vector{Float64} = Float64[] diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 75f527ae2..9ed89e0f3 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -66,6 +66,10 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spl mg["psi"] = result.psi_grid mg["dTdpsi"] = result.dtdpsi mg["T"] = result.t_cumulative + if !isempty(result.dtdpsi_ell) + mg["ell"] = result.ell + mg["dTdpsi_ell"] = result.dtdpsi_ell + end if dVdpsi_spline !== nothing mg["dVdpsi"] = [dVdpsi_spline(p) for p in result.psi_grid] end @@ -114,6 +118,8 @@ const KF_METHOD_H5_ANNOTATIONS = [ "psi" => (; long_name="normalized poloidal flux ψ_N at quadrature evaluation points", scale="psi"), "dTdpsi" => (; long_name="complex torque density dT/dψ_N = dT_φ/dψ_N + 2i·n·dδW_k/dψ_N; divide by dVdpsi for dT/dV", units="N*m", dims=("psi",), attach=(1 => "psi",)), "T" => (; long_name="cumulative complex torque T(ψ_N) = T_φ + 2i·n·δW_k (trapezoidal)", units="N*m", dims=("psi",), attach=(1 => "psi",)), + "ell" => (; long_name="bounce harmonic ℓ ∈ -nl:nl of the per-harmonic torque density", scale="ell"), + "dTdpsi_ell" => (; long_name="complex torque density dT/dψ_N resolved by bounce harmonic ℓ; summing over ℓ reproduces dTdpsi", units="N*m", dims=("psi", "ell"), attach=(1 => "psi", 2 => "ell")), "dVdpsi" => (; long_name="flux-surface volume derivative dV/dψ_N at quadrature points", units="m^3", dims=("psi",), attach=(1 => "psi",)), "EnergyIntegrals/psi" => (; long_name="ψ_N of each energy-integration record"), "EnergyIntegrals/lambda" => (; long_name="pitch λ = μB0/E of each record"), diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index d502385f7..fa3655fac 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -590,6 +590,7 @@ @test mr.total_torque == 0.0 + 0.0im @test mr.total_energy == 0.0 @test isempty(mr.records) + @test isempty(mr.dtdpsi_ell) && isempty(mr.ell) end # ========================================================================= diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index 773174c9e..e951380e8 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -182,7 +182,8 @@ using HDF5 st = KF.KineticForcesState(); st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, total_torque=ComplexF64(sum(dt)), psi_grid=collect(Float64, psi), - dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt))); + dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt)), + dtdpsi_ell=ComplexF64.(hcat(0.2 .* dt, 0.5 .* dt, 0.3 .* dt)), ell=collect(-1:1)); st ) s1 = mk(0.0:0.1:1.0, fill(1.0, 11)) @@ -197,6 +198,23 @@ using HDF5 s3 = mk([0.0, 0.5, 0.5, 1.0], fill(1.0, 4)) c2 = KF.combine_species_states([s1, s3]) @test all(isfinite, real.(c2.method_results["fgar"].dtdpsi)) + @test all(isfinite, real.(c2.method_results["fgar"].dtdpsi_ell)) + + # Per-ℓ survives the combine on the same union grid, and summing over ℓ reproduces the + # total exactly (it is derived as the row sum, so this is not merely interpolation noise). + @test r.ell == collect(-1:1) + @test size(r.dtdpsi_ell) == (length(r.psi_grid), 3) + @test vec(sum(r.dtdpsi_ell; dims=2)) ≈ r.dtdpsi + ov = [real(d) for (p, d) in zip(r.psi_grid, r.dtdpsi_ell[:, 2]) if 0.1 <= p <= 0.9] + @test all(x -> isapprox(x, 1.5; atol=1e-9), ov) # 0.5·(1 + 2) on the overlap + # A species with no per-ℓ data must leave the combined matrix empty, not error. + s4 = KF.KineticForcesState() + s4.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, + psi_grid=[0.0, 0.5, 1.0], dtdpsi=fill(ComplexF64(1.0), 3), + t_cumulative=zeros(ComplexF64, 3)) + c3 = KF.combine_species_states([s1, s4]) + @test isempty(c3.method_results["fgar"].dtdpsi_ell) + @test isempty(c3.method_results["fgar"].ell) end @testset "per-species HDF5 layout" begin @@ -206,7 +224,8 @@ using HDF5 st = KF.KineticForcesState(); st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, total_torque=ComplexF64(tt), psi_grid=collect(0.0:0.25:1.0), - dtdpsi=fill(ComplexF64(tt), 5), t_cumulative=zeros(ComplexF64, 5)); + dtdpsi=fill(ComplexF64(tt), 5), t_cumulative=zeros(ComplexF64, 5), + dtdpsi_ell=fill(ComplexF64(tt / 3), 5, 3), ell=collect(-1:1)); st ) h5 = tempname() * ".h5" @@ -221,10 +240,19 @@ using HDF5 @test read(f["KineticForces"]["fgar"]["total_torque"]) == 3.0 @test read(f["KineticForces"]["PerSpecies"]["ion_z1_m2"]["fgar"]["total_torque"]) == 1.0 # The metadata pass must reach the per-species method level, not stop at PerSpecies. - for d in ("total_torque", "dTdpsi") + for d in ("total_torque", "dTdpsi", "dTdpsi_ell") a = HDF5.attributes(f["KineticForces"]["PerSpecies"]["electron"]["fgar"][d]) @test haskey(a, "long_name") && haskey(a, "units") end + # HDF5.jl round-trips the Julia shape, so a read here gives (npsi, nell); a + # C-ordered reader (h5py, MATLAB) sees the transpose, (nell, npsi), which is the + # orientation the Fortran dTdpsi_fgar(psi, ell) implies. + g = f["KineticForces"]["fgar"] + @test size(read(g["dTdpsi_ell"])) == (5, 3) + @test read(g["ell"]) == [-1, 0, 1] + # Summing over ℓ must reproduce dTdpsi. + @test vec(sum(read(g["dTdpsi_ell"]); dims=2)) ≈ read(g["dTdpsi"]) + @test haskey(HDF5.attributes(g["dTdpsi_ell"]), "dims") end rm(h5) end From 3696f809183c2dfdc15f272fbaef7893ec8ec6bf Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Tue, 20 Oct 2026 13:19:21 +0200 Subject: [PATCH 2/3] KineticForces - TEST - Correct the per-harmonic orientation comment about the Fortran layout PENTRC writes dTdpsi_ as (i, ell, psi) in Fortran order, so the h5py view of dTdpsi_ell, (nell, npsi), is not the Fortran orientation. MATLAB reverses HDF5 dimensions too, so it is not a C-ordered reader. Co-Authored-By: Claude Opus 5.5 --- test/runtests_multiion.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index e951380e8..95b0d58ee 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -245,8 +245,7 @@ using HDF5 @test haskey(a, "long_name") && haskey(a, "units") end # HDF5.jl round-trips the Julia shape, so a read here gives (npsi, nell); a - # C-ordered reader (h5py, MATLAB) sees the transpose, (nell, npsi), which is the - # orientation the Fortran dTdpsi_fgar(psi, ell) implies. + # C-ordered reader (h5py, xarray) sees the transpose, (nell, npsi). g = f["KineticForces"]["fgar"] @test size(read(g["dTdpsi_ell"])) == (5, 3) @test read(g["ell"]) == [-1, 0, 1] From 9aeee4d673d54ab963d12ad8d7a3e2a1fca3e673 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Tue, 20 Oct 2026 13:43:41 +0200 Subject: [PATCH 3/3] KineticForces - BUGFIX - Apply the axis-validity envelope to the per-harmonic torque log The near-axis validity suppression scales each node's summed torque density by an envelope, but the per-harmonic log kept the unscaled harmonics. Between psi_c and 2 psi_c, sum over ell of dTdpsi_ell then differed from dTdpsi. The multi-ion combined total takes its dtdpsi from that row sum, so it lost the suppression as well. Scale the per-harmonic log by the same envelope. Co-Authored-By: Claude Opus 5.5 --- src/KineticForces/Compute.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 3ec1b1039..a4f8a4ec8 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -161,19 +161,20 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm total += harm_vals[ell_idx] end - total *= kinetic_axis_validity_envelope(psi, psi_c) + envelope = kinetic_axis_validity_envelope(psi, psi_c) + total *= envelope y[k] = total push!(logged_psi, psi) push!(logged_dtdpsi, total) - # `harm_vals` is overwritten at every ψ node, so the per-ℓ log must own a copy. - push!(logged_dtdpsi_ell, copy(harm_vals)) + # Same envelope as `total` so Σ_ℓ still reproduces it; the broadcast also copies the reused `harm_vals`. + push!(logged_dtdpsi_ell, harm_vals .* envelope) if is_matrix_method elems_accum = zeros(ComplexF64, mpert, mpert, 6) for ell_idx in 1:nharm elems_accum .+= harm_elems[ell_idx] end - elems_accum .*= kinetic_axis_validity_envelope(psi, psi_c) + elems_accum .*= envelope push!(logged_elems, elems_accum) end end