From dfd26eca983bfd2f1dbf0da52f2b6e6d7d2ce108 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 14:51:07 -0400 Subject: [PATCH 1/8] =?UTF-8?q?ForceFreeStates=20-=20FEATURE=20-=20Integra?= =?UTF-8?q?te=20the=20Euler-Lagrange=20state=20to=20any=20=CF=88=20from=20?= =?UTF-8?q?the=20nearest=20saved=20step?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add ForceFreeStates.solution_state_at(ffs, psi), which returns the stored basis state of a forward run at an arbitrary ψ. A saved step comes back as stored; otherwise the linear Euler-Lagrange system is integrated to psi from the nearest saved step that no singular surface separates from it, so the state does not depend on how densely the sweep saved. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/Result.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index 4c03700f90..9782be13e3 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -249,3 +249,29 @@ function build_result( solution, odet, wp, free_energies, delta_prime, gal_data ) end + +""" + solution_state_at(ffs, psi; rtol=ffs.control.eulerlagrange_tolerance) -> Array{ComplexF64,3} + +State `(u₁, u₂)` of a forward run's stored Euler-Lagrange solution at `psi`, in the basis of +`ffs.solution.u_store`. A saved step is returned as stored; otherwise the system is integrated to +`psi` from the nearest saved step that no singular surface separates from it. The stored basis is a +constant right-multiple of the integrated one and the system is linear in `u`, so this is the state +the forward sweep would have saved at `psi`, without interpolating between saved steps. +""" +function solution_state_at(ffs::ForceFreeStatesResult, psi::Float64; rtol::Float64=ffs.control.eulerlagrange_tolerance) + sol = ffs.solution + psis = @view sol.psi_store[1:sol.step] + i = searchsortedlast(psis, psi) + (i >= 1 && psis[i] == psi) && return sol.u_store[:, :, :, i] + sing = [s.psifac for s in ffs.surfaces] + reachable(j) = 1 <= j <= sol.step && !any(s -> min(psis[j], psi) < s < max(psis[j], psi), sing) + from = filter(reachable, [i, i + 1]) + isempty(from) && throw(ArgumentError("ψ = $psi is cut off from both neighbouring saved steps by a singular surface")) + j = argmin(j -> abs(psis[j] - psi), from) + kinetic = ffs.mats.kinetic !== nothing + hints = (Ref(1), Ref(1)) + rhs!(du, u, _, t) = (el_derivatives!(du, u, kinetic, ffs.equil, ffs.mats, ffs, t, hints...); nothing) + prob = ODEProblem(rhs!, sol.u_store[:, :, :, j], (psis[j], psi)) + return solve_el(prob, rtol).u[end] +end From a969f779f9ecdfd84e96f29be4ea2209facd8066 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 14:51:11 -0400 Subject: [PATCH 2/8] =?UTF-8?q?Analysis=20-=20FEATURE=20-=20Decompose=20th?= =?UTF-8?q?e=20plasma=20energy=20of=20free-boundary=20eigenmodes=20on=20th?= =?UTF-8?q?e=20(=CF=88,=CE=B8)=20grid?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add Analysis.EnergyDecomposition, a scripted diagnostic that evaluates δW_p of the free-boundary DCON eigenmodes of a forward, ideal, single-n solve in the effective-field form ½∫[|b_eff|²/μ₀ − (K₁+K₂+K₃)|ξ_n|²] and the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)], scaled to compare with ep[k]. decompose_plasma_energy returns a NamedTuple of the seven pointwise integrands, their surface integrals and both totals; plot_energy_decomposition draws the radial profiles and cumulative δW_p. All geometry comes from ffs.metric: the contravariant products follow from the covariant metric by cofactors and B² is the metric field strength. ξ at each knot comes from ForceFreeStates.solution_state_at. It refuses kinetic, multi-n, non-forward and fixed-boundary results. On the same solve it reproduces the former PerturbedEquilibrium stage to 1e-13 in both totals. Co-Authored-By: Claude Opus 5.5 --- src/Analysis/Analysis.jl | 2 + src/Analysis/EnergyDecomposition.jl | 291 ++++++++++++++++++++++++++++ 2 files changed, 293 insertions(+) create mode 100644 src/Analysis/EnergyDecomposition.jl diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl index 59309953da..06f38936ce 100644 --- a/src/Analysis/Analysis.jl +++ b/src/Analysis/Analysis.jl @@ -11,6 +11,7 @@ Post-processing and visualization utilities for GPEC simulation outputs. - `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium and singular coupling results - `PerturbedEquilibriumModes`: Data helpers to convert modal GPEC output to (ψ, θ) and (ψ, θ, φ) grids - `ErrorFields`: Plotting functions for the error-field assessment (coil sensitivities, tolerance Monte Carlo, locking risk, phasing maps) + - `EnergyDecomposition`: Decomposition of the plasma energy of free-boundary eigenmodes on the (ψ, θ) grid """ module Analysis @@ -20,5 +21,6 @@ include("CoilForcing.jl") include("PerturbedEquilibrium.jl") include("PerturbedEquilibriumModes.jl") include("ErrorFields.jl") +include("EnergyDecomposition.jl") end # module Analysis diff --git a/src/Analysis/EnergyDecomposition.jl b/src/Analysis/EnergyDecomposition.jl new file mode 100644 index 0000000000..2e713bbdc3 --- /dev/null +++ b/src/Analysis/EnergyDecomposition.jl @@ -0,0 +1,291 @@ +""" + EnergyDecomposition + +Decomposition of the plasma potential energy δW_p of the free-boundary DCON eigenmodes of a forward +ForceFreeStates (FFS) solve on the (ψ, θ) grid, in two forms: + + effective-field form 2δW_p = ∫ [ |b_eff|²/μ₀ − K|ξ_n|² ] dV, K = K₁ + K₂ + K₃ + standard form 2δW_p = ∫ [ |b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p) ] dV + +with b_eff = b + ξ_n(μ₀ j × n̂) and K₁ = |∇ψ|²σS (shear current), K₂ = μ₀B²σ² (parallel current), +K₃ = 2p'κ_ψ (pressure curvature) (Bernstein et al. 1958; Bernstein 1983; Chance et al. 1992; +Boozer 2004 Sec. V.C). The Fortran counterpart is `gpout_recon` with `effective_field_flag` and +`standard_form_flag`. Both totals are scaled by 2μ₀/ψ₀² so they read in the normalization of +`ffs.free_boundary.ep`. + +Call [`decompose_plasma_energy`](@ref) on a `ForceFreeStatesResult` and plot the result with +[`plot_energy_decomposition`](@ref). +""" +module EnergyDecomposition + +using FastInterpolations +using LinearAlgebra +using Plots +import ...ForceFreeStates as FFS +import ...PerturbedEquilibrium as PE +import ...Utilities.FourierTransforms: FourierTransform, inverse +import ...Utilities.PhysicalConstants: MU_0 + +export decompose_plasma_energy, plot_energy_decomposition + +""" +Names of the seven density terms, in the order of the third axis of `density` and the second of `profile`. +""" +const TERMS = (:b_eff_sq, :K1, :K2, :K3, :b_sq, :j_b_xi, :divxi_gradp) + +# ψ_N gap below which the last equilibrium knot counts as the solution edge. +const EDGE_KNOT_TOL = 1e-10 + +""" + radial_nodes(ffs) -> Vector{Float64} + +The equilibrium ψ knots inside the ξ solution domain, plus the solution edge when the last knot +falls short of it. On the knots the metric, the Euler-Lagrange matrices and the geometry are exact. +""" +function radial_nodes(ffs::FFS.ForceFreeStatesResult) + sol = ffs.solution + lo, hi = sol.psi_store[1], sol.psi_store[sol.step] + psi = [x for x in ffs.equil.rzphi_xs if lo <= x <= hi] + (isempty(psi) || psi[end] < hi - EDGE_KNOT_TOL) && push!(psi, hi) + return psi +end + +""" + surface_geometry(ffs, psi) -> NamedTuple of (length(psi), mtheta) matrices + +Geometry, equilibrium current and curvature-kernel coefficients at the nodes `psi` on the θ knots, +all from `ffs.metric` (columns e_i·e_j/J, J, ∂J/∂ψ of `make_metric`). The contravariant products +follow from the covariant metric by cofactors, |∇ψ|² = (g₂₂g₃₃ − g₂₃²)/J², ∇ψ·∇θ = (g₁₃g₂₃ − g₁₂g₃₃)/J², +∇ψ·∇ζ = (g₁₂g₂₃ − g₁₃g₂₂)/J², and B² = (χ'/J)²(g₂₂ + 2q g₂₃ + q² g₃₃). The kernel follows Fortran +`gpeq_shear`, `gpeq_curvature` and `gpeq_K`: + + S = (χ'²/J) [ q' + ∂_θ( (q ∇ψ·∇θ − ∇ψ·∇ζ)/|∇ψ|² ) ] + κ_ψ = (|∇ψ|²/B²) [ (μ₀p)' + ½∂_ψB² + ½∂_θB² (∇ψ·∇θ)/|∇ψ|² ] + σ = j·B/B², μ₀j^θ = −F'/J, μ₀j^ζ = −(μ₀p)'/χ' − qF'/J +""" +function surface_geometry(ffs::FFS.ForceFreeStatesResult, psi::Vector{Float64}) + equil = ffs.equil + profiles = equil.profiles + fs = ffs.metric.fs + xs = equil.rzphi_xs + theta_ext = collect(equil.rzphi_ys) + mtheta = length(theta_ext) - 1 + chi1 = 2π * equil.psio + r2 = equil.rzphi_rsquared.nodal_derivs.partials + deta = equil.rzphi_offset.nodal_derivs.partials + + # Knot rows of the eight metric columns, B², R and Z, splined in ψ to reach the edge node off the knots. + knot_vals = zeros(length(xs), 11 * mtheta) + for i in eachindex(xs), k in 1:mtheta + J = fs[i, k, 7] + q = profiles.q_spline(xs[i]) + rfac = sqrt(abs(r2[1, i, k])) + eta = 2π * (theta_ext[k] + deta[1, i, k]) + # Deviation from Fortran gpeq_K/gpeq_curvature: B² of the metric, not eqfun B², so K₂|ξ_n|² cancels |b_eff|² exactly. + bsq = (chi1 / J)^2 * J * (fs[i, k, 2] + 2q * fs[i, k, 4] + q^2 * fs[i, k, 3]) + for (c, v) in enumerate((fs[i, k, 1:8]..., bsq, equil.ro + rfac * cos(eta), equil.zo + rfac * sin(eta))) + knot_vals[i, (c-1)*mtheta+k] = v + end + end + itp = cubic_interp(xs, Series(knot_vals); bc=CubicFit(), extrap=ExtendExtrap()) + + names = (:R, :Z, :jac, :jac_psi, :g11, :g12, :g13, :g22, :g23, :g33, :delpsi2, :dpdt, :dpdz, :mu0_j_theta, :mu0_j_zeta, :shear, :K1, :K2, :K3) + geo = NamedTuple{names}(ntuple(_ -> zeros(length(psi), mtheta), length(names))) + ratio = zeros(mtheta + 1) + for (i, p) in enumerate(psi) + v = itp(p) + dv = itp(p; deriv=DerivOp(1)) + col(c, k) = v[(c-1)*mtheta+k] + q, q1 = profiles.q_spline(p), profiles.q_deriv(p) + F1, mu0p1 = profiles.F_deriv(p), profiles.P_deriv(p) + bsq = [col(9, k) for k in 1:mtheta] + dbsq_theta = deriv1(cubic_interp(theta_ext, [bsq; bsq[1]]; bc=PeriodicBC())) + for k in 1:mtheta + J = col(7, k) + g11, g22, g33, g23, g13, g12 = (col(c, k) * J for c in 1:6) + d2 = (g22 * g33 - g23^2) / J^2 + dpdt = (g13 * g23 - g12 * g33) / J^2 + dpdz = (g12 * g23 - g13 * g22) / J^2 + geo.R[i, k], geo.Z[i, k] = col(10, k), col(11, k) + geo.jac[i, k], geo.jac_psi[i, k] = J, col(8, k) + geo.g11[i, k], geo.g12[i, k], geo.g13[i, k], geo.g22[i, k], geo.g23[i, k], geo.g33[i, k] = g11, g12, g13, g22, g23, g33 + geo.delpsi2[i, k], geo.dpdt[i, k], geo.dpdz[i, k] = d2, dpdt, dpdz + geo.mu0_j_theta[i, k] = -F1 / J + geo.mu0_j_zeta[i, k] = -mu0p1 / chi1 - F1 * q / J + ratio[k] = (q * dpdt - dpdz) / d2 + end + ratio[end] = ratio[1] + dratio = deriv1(cubic_interp(theta_ext, ratio; bc=PeriodicBC())) + for k in 1:mtheta + J, d2, b2 = geo.jac[i, k], geo.delpsi2[i, k], bsq[k] + S = chi1^2 / J * (q1 + dratio(theta_ext[k])) + kappa = d2 / b2 * (mu0p1 + 0.5 * dv[8*mtheta+k] + 0.5 * dbsq_theta(theta_ext[k]) * geo.dpdt[i, k] / d2) + jt, jz = geo.mu0_j_theta[i, k] / MU_0, geo.mu0_j_zeta[i, k] / MU_0 + bt, bz = chi1 / J, q * chi1 / J + sigma = (geo.g22[i, k] * bt * jt + geo.g33[i, k] * bz * jz + geo.g23[i, k] * (bt * jz + bz * jt)) / b2 + geo.shear[i, k] = S + geo.K1[i, k] = d2 * sigma * S + geo.K2[i, k] = MU_0 * b2 * sigma^2 + geo.K3[i, k] = 2 * mu0p1 / MU_0 * kappa + end + end + return geo +end + +""" + node_states(ffs, psi) -> (u1, du1, xi_s) + +Euler-Lagrange basis Ξ_ψ, dΞ_ψ/dψ and Ξ_s = −A⁻¹(BΞ'_ψ + CΞ_ψ) at the nodes `psi`, each +`(numpert_total, numpert_total, length(psi))`. Each node state is integrated from the nearest saved +step (`ForceFreeStates.solution_state_at`), so it does not depend on how densely the sweep saved. +""" +function node_states(ffs::FFS.ForceFreeStatesResult, psi::Vector{Float64}) + N = ffs.numpert_total + u1, du1, xi_s = (zeros(ComplexF64, N, N, length(psi)) for _ in 1:3) + du = zeros(ComplexF64, N, N, 2) + for (i, p) in enumerate(psi) + u = FFS.solution_state_at(ffs, p) + FFS.el_derivatives!(du, u, false, ffs.equil, ffs.mats, ffs, p, Ref(1), Ref(1)) + u1[:, :, i] .= @view u[:, :, 1] + du1[:, :, i] .= @view du[:, :, 1] + FFS.compute_node_xi_s!(view(xi_s, :, :, i), view(du1, :, :, i), view(u1, :, :, i), ffs.mats, p) + end + return u1, du1, xi_s +end + +""" + eigenmode_fields(ffs, k, psi, states, ft) -> NamedTuple of (length(psi), mtheta) matrices + +θ-space fields of free-boundary eigenmode `k`: ξ^ψ (`xi`), ∂ξ^ψ/∂ψ (`xi1`), Jξ^θ, Jξ^ζ, ∂_θ(Jξ^θ) +and Jb^ψ, Jb^θ, Jb^ζ. The edge condition is the unit-power eigenvector `wt[:, k]`; the mode-space +chain is the one `PerturbedEquilibrium` runs for a forced response. +""" +function eigenmode_fields(ffs::FFS.ForceFreeStatesResult, k::Int, psi::Vector{Float64}, states, ft::FourierTransform) + u1, du1, xs = states + sol = ffs.solution + alpha = sol.u_store[:, :, 1, sol.step] \ Vector{ComplexF64}(ffs.free_boundary.wt[:, k]) + rows(A) = permutedims(reduce(hcat, (A[:, :, i] * alpha for i in eachindex(psi)))) + xi, xi1, xi_s = rows(u1), rows(du1), rows(xs) + chi1 = 2π * ffs.equil.psio + Jbp, Jbt, Jbz = PE.compute_perturbed_field_modes(xi, xi1, xi_s, psi, ffs.equil, ffs) + # Deviation from Fortran gpeq_contra: always unregularized (regularization width 0). + _, Jxt, Jxz, _, _ = PE.compute_contra_displacements(xi, xi1, xi_s ./ chi1, psi, ffs.equil, ffs, ffs.metric, 0.0) + # Deviation from Fortran gpeq_firstform: ∂_θ(Jξ^θ) is spectral, not a periodic-spline derivative. + dJxt = Jxt .* transpose(2π * im .* (ffs.mlow:ffs.mhigh)) + theta(X) = permutedims(inverse(ft, permutedims(X))) + return (; xi=theta(xi), xi1=theta(xi1), Jxt=theta(Jxt), Jxz=theta(Jxz), dJxt=theta(dJxt), Jbp=theta(Jbp), Jbt=theta(Jbt), Jbz=theta(Jbz)) +end + +""" + effective_field(geo, f, i, k) -> (Jc_psi, Jc_theta, Jc_zeta, c_psi, c_theta, c_zeta, b_psi, b_theta, b_zeta, xi_n) + +Contravariant J b_eff^i and covariant b_eff,i of b_eff = b + ξ_n(μ₀ j × n̂) at node `i`, θ point `k` +(Fortran `gpeq_c`), with the covariant field b_i = g_ij Jb^j/J and ξ_n = ξ^ψ/|∇ψ|: + + J b_eff^θ = J b^θ + Jξ^ψ/(J|∇ψ|²) (μ₀j^θ g_θζ + μ₀j^ζ g_ζζ) b_eff,θ = b_θ + Jξ^ψ μ₀j^ζ + J b_eff^ζ = J b^ζ − Jξ^ψ/(J|∇ψ|²) (μ₀j^θ g_θθ + μ₀j^ζ g_θζ) b_eff,ζ = b_ζ − Jξ^ψ μ₀j^θ + b_eff,ψ = b_ψ + Jξ^ψ/|∇ψ|² (μ₀j^θ ∇ψ·∇ζ − μ₀j^ζ ∇ψ·∇θ), J b_eff^ψ = J b^ψ +""" +function effective_field(geo::NamedTuple, f::NamedTuple, i::Int, k::Int) + J, d2 = geo.jac[i, k], geo.delpsi2[i, k] + g11, g12, g13, g22, g23, g33 = geo.g11[i, k], geo.g12[i, k], geo.g13[i, k], geo.g22[i, k], geo.g23[i, k], geo.g33[i, k] + jt, jz = geo.mu0_j_theta[i, k], geo.mu0_j_zeta[i, k] + Jbp, Jbt, Jbz = f.Jbp[i, k], f.Jbt[i, k], f.Jbz[i, k] + # Deviation from Fortran gpeq_cova/gpeq_normal/gpeq_c: b_i, Jξ^ψ, ξ_n and b_eff are pointwise in θ, not band-truncated. + xw = J * f.xi[i, k] + bp = (g11 * Jbp + g12 * Jbt + g13 * Jbz) / J + bt = (g12 * Jbp + g22 * Jbt + g23 * Jbz) / J + bz = (g13 * Jbp + g23 * Jbt + g33 * Jbz) / J + Jct = Jbt + xw / (d2 * J) * (jt * g23 + jz * g33) + Jcz = Jbz - xw / (d2 * J) * (jt * g22 + jz * g23) + cp = bp + xw / d2 * (jt * geo.dpdz[i, k] - jz * geo.dpdt[i, k]) + return (Jbp, Jct, Jcz, cp, bt + xw * jz, bz - xw * jt, bp, bt, bz, xw / (J * sqrt(d2))) +end + +""" + decompose_plasma_energy(ffs; eigenmodes=[1]) -> NamedTuple + +Decompose δW_p of the free-boundary eigenmodes `eigenmodes` of a forward, ideal, single-n FFS solve +with the vacuum response (`vac_flag = true`). Fields of the result: + + - `psi` - radial nodes (`radial_nodes`); `theta` - θ knots in [0, 1); `R`, `Z` - node positions in m + - `terms` - `TERMS`, the names along the term axis + - `eigenmodes` - the decomposed eigenmode indices + - `density` - `(npsi, mtheta, 7, neig)` pointwise integrands, scaled by 2μ₀/ψ₀² + - `profile` - `(npsi, 7, neig)` surface integrals ∮J·density dθ + - `dW_effective`, `dW_standard` - ½∫ of the effective-field and standard-form combinations + +The relative error against the DCON energy is `(dW − real(ep[k])) / abs(ep[k])` with +`ep = ffs.free_boundary.ep`. +""" +function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::AbstractVector{<:Integer}=[1]) + ffs.integrator === :forward && ffs.solution !== nothing || throw(ArgumentError("energy decomposition needs a forward-integrator solve")) + ffs.free_boundary === nothing && throw(ArgumentError("energy decomposition needs the free-boundary eigenmodes (vac_flag = true)")) + ffs.nlow == ffs.nhigh || throw(ArgumentError("energy decomposition handles a single toroidal mode number")) + ffs.mats.kinetic === nothing || throw(ArgumentError("energy decomposition handles ideal solves only, not kinetic ones")) + nwt = size(ffs.free_boundary.wt, 2) + all(k -> 1 <= k <= nwt, eigenmodes) || throw(ArgumentError("eigenmodes must lie in 1:$nwt")) + + equil = ffs.equil + # Deviation from Fortran gpout_recon: the radial integral runs over the equilibrium knots, not the ODE psifac grid. + psi = radial_nodes(ffs) + mtheta = length(equil.rzphi_ys) - 1 + geo = surface_geometry(ffs, psi) + states = node_states(ffs, psi) + ft = FourierTransform(mtheta, ffs.mpert, ffs.mlow) + n = ffs.nlow + scale = 2 * MU_0 / equil.psio^2 + + density = zeros(ComplexF64, length(psi), mtheta, length(TERMS), length(eigenmodes)) + for (m, k) in enumerate(eigenmodes) + f = eigenmode_fields(ffs, k, psi, states, ft) + for i in eachindex(psi) + p1 = equil.profiles.P_deriv(psi[i]) / MU_0 + for j in 1:mtheta + Jcp, Jct, Jcz, cp, ct, cz, bp, bt, bz, xi_n = effective_field(geo, f, i, j) + J, xw = geo.jac[i, j], geo.jac[i, j] * f.xi[i, j] + jt, jz = geo.mu0_j_theta[i, j], geo.mu0_j_zeta[i, j] + divxi = f.xi1[i, j] + geo.jac_psi[i, j] / J * f.xi[i, j] + f.dJxt[i, j] / J - 2π * im * n * f.Jxz[i, j] / J + terms = (real(conj(Jcp) * cp + conj(Jct) * ct + conj(Jcz) * cz) / (MU_0 * J), + geo.K1[i, j] * abs2(xi_n), geo.K2[i, j] * abs2(xi_n), geo.K3[i, j] * abs2(xi_n), + real(conj(f.Jbp[i, j]) * bp + conj(f.Jbt[i, j]) * bt + conj(f.Jbz[i, j]) * bz) / (MU_0 * J), + (-jt * (f.Jbp[i, j] * conj(f.Jxz[i, j]) - f.Jbz[i, j] * conj(xw)) + jz * (f.Jbp[i, j] * conj(f.Jxt[i, j]) - f.Jbt[i, j] * conj(xw))) / (MU_0 * J), + conj(divxi) * p1 * f.xi[i, j]) + for t in eachindex(terms) + density[i, j, t, m] = scale * terms[t] + end + end + end + end + profile = dropdims(sum(density .* geo.jac; dims=2); dims=2) ./ mtheta + # The imaginary parts of the surface integrals are roundoff, so the ψ integrals take the real part. + total(t, m) = FastInterpolations.integrate(cubic_interp(psi, real.(profile[:, t, m]); bc=CubicFit())) + dW_effective = [0.5 * (total(1, m) - total(2, m) - total(3, m) - total(4, m)) for m in eachindex(eigenmodes)] + dW_standard = [0.5 * (total(5, m) - total(6, m) + total(7, m)) for m in eachindex(eigenmodes)] + return (; psi, theta=equil.rzphi_ys[1:mtheta], R=geo.R, Z=geo.Z, terms=TERMS, eigenmodes=collect(eigenmodes), density, profile, dW_effective, dW_standard) +end + +""" + plot_energy_decomposition(ed; eigenmode=first(ed.eigenmodes), save_path=nothing) -> Plots.Plot + +Radial profiles of the seven terms of `ed` (from `decompose_plasma_energy`) for one eigenmode, over +the cumulative effective-field and standard-form δW_p. +""" +function plot_energy_decomposition(ed::NamedTuple; eigenmode::Integer=first(ed.eigenmodes), save_path=nothing) + m = findfirst(==(eigenmode), ed.eigenmodes) + m === nothing && throw(ArgumentError("eigenmode $eigenmode was not decomposed")) + prof = real.(ed.profile[:, :, m]) + cum = reduce(hcat, (vec(cumulative_integrate(cubic_interp(ed.psi, prof[:, t]; bc=CubicFit()))) for t in eachindex(ed.terms))) + eff = 0.5 .* (cum[:, 1] .- cum[:, 2] .- cum[:, 3] .- cum[:, 4]) + std = 0.5 .* (cum[:, 5] .- cum[:, 6] .+ cum[:, 7]) + p1 = plot(ed.psi, prof; label=permutedims(collect(string.(ed.terms))), xlabel="ψ_N", ylabel="∮ J·density dθ", title="Eigenmode $eigenmode") + p2 = plot(ed.psi, [eff std]; label=["effective field" "standard form"], xlabel="ψ_N", ylabel="cumulative δW_p") + p = plot(p1, p2; layout=(2, 1), size=(800, 800), left_margin=12Plots.mm, bottom_margin=4Plots.mm) + if save_path !== nothing + savefig(p, save_path) + println("Saved energy decomposition plot to ", abspath(save_path)) + end + return p +end + +end # module EnergyDecomposition From 34199838d8d2b5f31b40b17339a316d769fdf9e3 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 14:51:13 -0400 Subject: [PATCH 3/8] Analysis - TEST - Check the energy decomposition against the DCON Lagrangian and the plasma energy MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit On the Solovev fixture: the standard-form surface integral equals 2/ψ₀² times the Euler-Lagrange quadratic form on every knot; both totals match each other and ep[k] within the measured quadrature error; the (∇×b_eff)·∇ψ identity, the shear surface average, the off-knot metric and solution_state_at are checked; the plot is drawn. Co-Authored-By: Claude Opus 5.5 --- test/runtests_energy_decomposition.jl | 108 ++++++++++++++++++++++++++ test/suites.jl | 1 + 2 files changed, 109 insertions(+) create mode 100644 test/runtests_energy_decomposition.jl diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl new file mode 100644 index 0000000000..bfe2760401 --- /dev/null +++ b/test/runtests_energy_decomposition.jl @@ -0,0 +1,108 @@ +using TOML +using Plots + +# Analysis.EnergyDecomposition on the coarse Solovev fixture (mpsi = 16, forward integrator, vac_flag = true): +# one solve and one decomposition, shared by every testset. +@testset "Analysis.EnergyDecomposition" begin + ENV["GKSwstype"] = "100" + GPEC = GeneralizedPerturbedEquilibrium + ED = GPEC.Analysis.EnergyDecomposition + FT = GPEC.Utilities.FourierTransforms + + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), + GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])) + wall = GPEC.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in deck["Wall"])...) + ffs_kwargs = Dict(Symbol(k) => v for (k, v) in deck["ForceFreeStates"] if !(k in ("integrator", "nn_low", "nn_high"))) + ffs = solve(equil, Forward(); nn=1, wall=wall, dir_path=mktempdir(), ffs_kwargs...) + + ed = ED.decompose_plasma_energy(ffs; eigenmodes=[1, 2]) + psi = ed.psi + mtheta = length(ed.theta) + knots = findall(in(equil.rzphi_xs), psi) + geo = ED.surface_geometry(ffs, psi) + states = ED.node_states(ffs, psi) + ep = real.(ffs.free_boundary.ep) + + @testset "standard form equals the DCON Lagrangian on the knots" begin + # On a knot the standard-form surface integral is 2/ψ₀² times the Euler-Lagrange quadratic form ξ'†Fξ' + 2Re(ξ'†Kξ) + ξ†Gξ. + N = ffs.numpert_total + F, K, G = (zeros(ComplexF64, N, N) for _ in 1:3) + u1, du1, _ = states + sol = ffs.solution + for (m, k) in enumerate(ed.eigenmodes) + alpha = sol.u_store[:, :, 1, sol.step] \ Vector{ComplexF64}(ffs.free_boundary.wt[:, k]) + for i in knots + sf = (ffs.mlow:ffs.mhigh) .- ffs.nlow * equil.profiles.q_spline(psi[i]) + ffs.mats.ideal.F_spline_lower(vec(G), psi[i]) + F .= LowerTriangular(G) * LowerTriangular(G)' .* sf .* sf' + ffs.mats.ideal.K_spline(vec(K), psi[i]) + K .*= sf + ffs.mats.ideal.G_spline(vec(G), psi[i]) + xi, dxi = u1[:, :, i] * alpha, du1[:, :, i] * alpha + lagrangian = real(dot(dxi, F, dxi) + 2 * dot(dxi, K, xi) + dot(xi, G, xi)) + standard = real(ed.profile[i, 5, m] - ed.profile[i, 6, m] + ed.profile[i, 7, m]) + @test standard ≈ 2 / equil.psio^2 * lagrangian rtol = 1e-9 + end + end + end + + @testset "both forms reproduce the DCON plasma energy" begin + for (m, k) in enumerate(ed.eigenmodes) + # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.3e-5). + @test ed.dW_effective[m] ≈ ed.dW_standard[m] rtol = 2e-4 + # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 2.0 %); it converges with mpsi. + @test abs(ed.dW_standard[m] - ep[k]) / abs(ep[k]) < 0.025 + end + end + + @testset "(∇×b_eff)·∇ψ vanishes on the knots" begin + # An algebraic identity of the b_eff construction on the knots, so this checks the effective-field formulas only. + ft = FT.FourierTransform(mtheta, ffs.mpert, ffs.mlow) + mvals = ffs.mlow:ffs.mhigh + f = ED.eigenmode_fields(ffs, 1, psi, states, ft) + for i in knots + c = [ED.effective_field(geo, f, i, k) for k in 1:mtheta] + c_theta, c_zeta = ft(getindex.(c, 5)), ft(getindex.(c, 6)) + curl = FT.inverse(ft, 2π * im .* (mvals .* c_zeta .+ ffs.nlow .* c_theta)) ./ geo.jac[i, :] + dtheta_c_zeta = FT.inverse(ft, 2π * im .* mvals .* c_zeta) ./ geo.jac[i, :] + @test maximum(abs, curl) < 1e-10 * maximum(abs, dtheta_c_zeta) + end + end + + @testset "flux-surface average of the shear" begin + # The θ-derivative term of S averages out, so ⟨J S⟩ = χ'² q'. + i = knots[end÷2] + @test sum(geo.shear[i, :] .* geo.jac[i, :]) / mtheta ≈ (2π * equil.psio)^2 * equil.profiles.q_deriv(psi[i]) rtol = 1e-12 + end + + @testset "metric between the knots" begin + # Off the knots the metric is splined in ψ; it must match the rzphi bicubic to interpolation accuracy (measured 1.2e-5). + xs = equil.rzphi_xs + off = [(xs[end÷2] + xs[end÷2+1]) / 2, psi[end]] + g = ED.surface_geometry(ffs, off) + ref = [GPEC.Equilibrium.flux_surface_metric(equil, p, t) for p in off, t in ed.theta] + @test maximum(abs, g.jac ./ getfield.(ref, :jac) .- 1) < 1e-4 + @test maximum(abs, sqrt.(g.delpsi2) ./ getfield.(ref, :delpsi) .- 1) < 1e-4 + @test maximum(abs, g.R ./ getfield.(ref, :r) .- 1) < 1e-4 + end + + @testset "solution state between saved steps" begin + sol = ffs.solution + @test ForceFreeStates.solution_state_at(ffs, sol.psi_store[5]) == sol.u_store[:, :, :, 5] + # Just left of a mid-interval point is integrated from the left step, just right of it from the right one; they must meet. + sing = [s.psifac for s in ffs.surfaces] + i = 9 + findfirst(j -> !any(s -> sol.psi_store[j] < s < sol.psi_store[j+1], sing), 10:(sol.step-1)) + a, b = sol.psi_store[i], sol.psi_store[i+1] + from_left = ForceFreeStates.solution_state_at(ffs, (a + b) / 2 - 1e-9 * (b - a)) + from_right = ForceFreeStates.solution_state_at(ffs, (a + b) / 2 + 1e-9 * (b - a)) + @test maximum(abs, from_left - from_right) / maximum(abs, from_left) < 1e-6 + end + + @testset "plot" begin + path = joinpath(mktempdir(), "energy_decomposition.png") + @test ED.plot_energy_decomposition(ed; eigenmode=2, save_path=path) isa Plots.Plot + @test isfile(path) + end +end diff --git a/test/suites.jl b/test/suites.jl index 225eb2fd8a..70c356cad8 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -53,6 +53,7 @@ const TEST_SUITES = [ "runtests_multiion.jl" => "core", "runtests_fullruns.jl" => "api", "runtests_torque_response.jl" => "api", + "runtests_energy_decomposition.jl" => "api", "runtests_coils.jl" => "api", "runtests_imas.jl" => "api", "runtests_rerun_from_h5.jl" => "api", From 5226d82c40c2b38329f050116099ac912bd52640 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 14:51:15 -0400 Subject: [PATCH 4/8] Docs - DOCS - Document the energy decomposition and cite its energy-principle references Add the EnergyDecomposition section to the Analysis page (both forms, the K split, usage, grid and accuracy, Fortran deviations, the mpsi convergence table) with its @autodocs block, and cite Bernstein et al. 1958, Bernstein 1983, Chance et al. 1992 and Boozer 2004 in citations.md. Co-Authored-By: Claude Opus 5.5 --- docs/src/analysis.md | 83 +++++++++++++++++++++++++++++++++++++++++++ docs/src/citations.md | 33 +++++++++++++++++ 2 files changed, 116 insertions(+) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 575adf10e1..f1d135ad1c 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -9,6 +9,7 @@ The Analysis module provides post-processing and visualization utilities for GPE - `CoilForcing`: Plotting functions for coil geometry and normal field spectra - `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium results - `PerturbedEquilibriumModes`: Data helpers to convert modal GPEC output to (ψ, θ) and (ψ, θ, φ) grids +- `EnergyDecomposition`: Decomposition of the plasma energy of free-boundary eigenmodes on the (ψ, θ) grid ```@docs GeneralizedPerturbedEquilibrium.Analysis @@ -49,3 +50,85 @@ Modules = [GeneralizedPerturbedEquilibrium.Analysis.PerturbedEquilibriumModes] ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Analysis.ErrorFields] ``` + +## EnergyDecomposition + +`Analysis.EnergyDecomposition` evaluates the ideal-MHD plasma potential energy ``\delta W_p`` of the +free-boundary eigenmodes of a forward ForceFreeStates solve on the ``(\psi,\theta)`` grid and splits it +into its physical terms. It is a scripted diagnostic, the Julia counterpart of the Fortran GPEC +`gpout_recon` with `effective_field_flag` and `standard_form_flag`: + +```julia +using GeneralizedPerturbedEquilibrium +const ED = GeneralizedPerturbedEquilibrium.Analysis.EnergyDecomposition +ffs = GeneralizedPerturbedEquilibrium.main(["examples/Solovev_ideal_example"]).ffs +ed = ED.decompose_plasma_energy(ffs; eigenmodes=[1, 2]) +ED.plot_energy_decomposition(ed; eigenmode=1, save_path="energy_decomposition.png") +``` + +It needs `integrator = "forward"`, `vac_flag = true`, a single toroidal mode number and an ideal solve. + +### The two forms + +The energy principle (Bernstein et al. 1958; Boozer 2004, Sec. V.C) gives the *standard form* + +```math +2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B|^2}{\mu_0} +- \mathbf j\cdot(\delta\mathbf B\times\boldsymbol\xi^*) + (\nabla\cdot\boldsymbol\xi)^*\,(\boldsymbol\xi\cdot\nabla p) \right], +``` + +and completing the square with the equilibrium current (Bernstein 1983; Chance et al. 1992, Eq. 12c; +Boozer 2004, Eq. 100) gives the *effective-field form* + +```math +2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B_{\rm eff}|^2}{\mu_0} - K\,|\xi_n|^2 \right], +\qquad \delta\mathbf B_{\rm eff} = \delta\mathbf B + \xi_n\,(\mu_0 \mathbf j\times\hat n), +``` + +```math +K = K_1 + K_2 + K_3,\qquad K_1 = |\nabla\psi|^2\,\sigma S,\quad K_2 = \mu_0 B^2 \sigma^2,\quad K_3 = 2 p'\,\kappa_\psi, +\qquad \sigma = \frac{\mathbf j\cdot\mathbf B}{B^2}, +``` + +with ``S`` the local shear and ``\kappa_\psi`` the curvature projection (Chance et al. 1992, Eqs. 15–17): +the shear–current, parallel-current and pressure–curvature drives. Both forms omit the compressional +``\gamma p|\nabla\cdot\boldsymbol\xi|^2`` term, as the Fortran does. The code keeps the shorthand `b` for +``\delta\mathbf B``. + +The result holds the seven integrands `(b_eff_sq, K1, K2, K3, b_sq, j_b_xi, divxi_gradp)` per +``(\psi,\theta)`` node, their surface integrals and the two totals. Everything is scaled by +``2\mu_0/\psi_0^2``, so the totals compare directly with the DCON plasma energy +``e_{p,k}`` = `ffs.free_boundary.ep[k]` of the unit-power eigenvector ``\mathbf w_k``. + +### Grid and accuracy + +The radial nodes are the equilibrium ``\psi`` knots inside the solution domain, plus the solution edge. +On the knots all geometry is read from `ffs.metric`, and the standard-form surface integral equals +``2/\psi_0^2`` times the DCON Lagrangian ``\xi'^\dagger F\xi' + 2\,\mathrm{Re}(\xi'^\dagger K\xi) + \xi^\dagger G\xi`` +to roundoff, so the standard-form total differs from ``e_{p,k}`` only by the ``\psi`` quadrature of the knot +grid. ``\xi`` at each node is integrated from the nearest saved ODE step that no singular surface +separates from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. + +``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi = 0`` (Boozer 2004, Eq. 99) is an algebraic identity of +the construction on the knots. It checks how ``\delta\mathbf B_{\rm eff}`` is assembled from ``\delta\mathbf B``, +``\xi_n`` and the metric, not the radial ``\xi`` solution, the kernel ``K`` or the radial integral. + +Choices that differ from the Fortran: ``B^2`` is the metric field strength +``(\chi'/\mathcal J)^2(g_{\theta\theta} + 2q g_{\theta\zeta} + q^2 g_{\zeta\zeta})``, which makes ``K_2|\xi_n|^2`` +cancel against ``|\delta\mathbf B_{\rm eff}|^2``; the covariant field, ``\xi_n`` and ``\delta\mathbf B_{\rm eff}`` +are formed pointwise in ``\theta`` rather than band-truncated; ``\partial_\theta(\mathcal J\xi^\theta)`` is +spectral; the reconstruction is unregularized; and the totals integrate over the knots rather than the ODE +`psifac` grid. + +Relative error of both forms against ``e_{p,k}`` on the Solovev test fixture (`mpert = 8`): + +| `mpsi` | mode 1 effective / standard | mode 2 effective / standard | +|---|---|---| +| 16 | +1.39e-2 / +1.39e-2 | −1.98e-2 / −1.97e-2 | +| 32 | −1.42e-3 / −1.42e-3 | +1.07e-2 / +1.07e-2 | +| 64 | +2.10e-4 / +2.10e-4 | +6.43e-4 / +6.43e-4 | +| 128 | +1.62e-4 / +1.62e-4 | −4.18e-5 / −4.18e-5 | + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Analysis.EnergyDecomposition] +``` diff --git a/docs/src/citations.md b/docs/src/citations.md index 9bf7f5a03e..6cb92461e4 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -110,6 +110,39 @@ Establishes the area weightings that make flux-surface Fourier spectra of pertur --- +### Energy decomposition of free-boundary eigenmodes + +> I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M. Kulsrud, "An energy principle for hydromagnetic stability problems," +> *Proceedings of the Royal Society A* **244**, 17 (1958). +> DOI: [10.1098/rspa.1958.0023](https://doi.org/10.1098/rspa.1958.0023) + +The ideal-MHD energy principle: the standard form of the plasma potential energy δW_p that `Analysis.EnergyDecomposition` evaluates on the (ψ,θ) grid — see the [Analysis Module](analysis.md). + +--- + +> I. B. Bernstein, in *Basic Plasma Physics I*, Handbook of Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan, +> North-Holland, Amsterdam (1983), p. 421. + +Derives the effective-field (Bernstein) form of δW_p, in which the current and pressure terms are completed into |δB_eff|²/μ₀ − K|ξ_n|². + +--- + +> M. S. Chance, Y.-C. Sun, S. C. Jardin, C. E. Kessel, and M. Okabayashi, "MHD stability of tokamak plasmas," +> PPPL-CFP-2687, 2nd Symposium on Plasma Dynamics, Trieste (1992). +> Link: [OSTI 7234275](https://www.osti.gov/biblio/7234275) + +Writes both forms of δW_p side by side (Eqs. 12a–12c) and defines the destabilizing kernel K = K₁ + K₂ + K₃ with its shear–current, parallel-current and pressure–curvature pieces (Eqs. 15–17), the `K1`, `K2` and `K3` terms of the decomposition. + +--- + +> A. H. Boozer, "Physics of magnetically confined plasmas," +> *Reviews of Modern Physics* **76**, 1071 (2004). +> DOI: [10.1103/RevModPhys.76.1071](https://doi.org/10.1103/RevModPhys.76.1071) + +Sec. V.C derives the energy principle in the form used here, introduces the effective field C = δB + ξ_n(μ₀ j × n̂) (Eq. 98), the constraint (∇×C)·∇ψ = 0 that the perturbed current lie in the perturbed flux surfaces (Eq. 99) and the split w = C²/2μ₀ − w_d (Eq. 100). + +--- + ## InnerLayer Module > A. H. Glasser, Z. R. Wang, and J.-K. Park, "Computation of resistive instabilities by matched asymptotic expansions," From 3d53fb30c0bdb20e82a9e0ee5ad38673f5683df3 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 15:25:11 -0400 Subject: [PATCH 5/8] Analysis - DOCS - Cite the energy-principle references in the module docstring, not the citations page The citations page lists the papers to cite for GPEC itself. Move Bernstein et al. 1958, Bernstein 1983, Chance et al. 1992 and Boozer 2004 out of citations.md and the Analysis page text into a References section of the EnergyDecomposition module docstring, next to the code that uses them. Co-Authored-By: Claude Opus 5.5 --- docs/src/analysis.md | 11 +++++----- docs/src/citations.md | 33 ----------------------------- src/Analysis/EnergyDecomposition.jl | 14 +++++++++--- 3 files changed, 16 insertions(+), 42 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index f1d135ad1c..5575c726eb 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -70,15 +70,14 @@ It needs `integrator = "forward"`, `vac_flag = true`, a single toroidal mode num ### The two forms -The energy principle (Bernstein et al. 1958; Boozer 2004, Sec. V.C) gives the *standard form* +The energy principle gives the *standard form* ```math 2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B|^2}{\mu_0} - \mathbf j\cdot(\delta\mathbf B\times\boldsymbol\xi^*) + (\nabla\cdot\boldsymbol\xi)^*\,(\boldsymbol\xi\cdot\nabla p) \right], ``` -and completing the square with the equilibrium current (Bernstein 1983; Chance et al. 1992, Eq. 12c; -Boozer 2004, Eq. 100) gives the *effective-field form* +and completing the square with the equilibrium current gives the *effective-field form* ```math 2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B_{\rm eff}|^2}{\mu_0} - K\,|\xi_n|^2 \right], @@ -90,8 +89,8 @@ K = K_1 + K_2 + K_3,\qquad K_1 = |\nabla\psi|^2\,\sigma S,\quad K_2 = \mu_0 B^2 \qquad \sigma = \frac{\mathbf j\cdot\mathbf B}{B^2}, ``` -with ``S`` the local shear and ``\kappa_\psi`` the curvature projection (Chance et al. 1992, Eqs. 15–17): -the shear–current, parallel-current and pressure–curvature drives. Both forms omit the compressional +with ``S`` the local shear and ``\kappa_\psi`` the curvature projection: the shear–current, +parallel-current and pressure–curvature drives. Both forms omit the compressional ``\gamma p|\nabla\cdot\boldsymbol\xi|^2`` term, as the Fortran does. The code keeps the shorthand `b` for ``\delta\mathbf B``. @@ -109,7 +108,7 @@ to roundoff, so the standard-form total differs from ``e_{p,k}`` only by the ``\ grid. ``\xi`` at each node is integrated from the nearest saved ODE step that no singular surface separates from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. -``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi = 0`` (Boozer 2004, Eq. 99) is an algebraic identity of +``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi = 0`` is an algebraic identity of the construction on the knots. It checks how ``\delta\mathbf B_{\rm eff}`` is assembled from ``\delta\mathbf B``, ``\xi_n`` and the metric, not the radial ``\xi`` solution, the kernel ``K`` or the radial integral. diff --git a/docs/src/citations.md b/docs/src/citations.md index 6cb92461e4..9bf7f5a03e 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -110,39 +110,6 @@ Establishes the area weightings that make flux-surface Fourier spectra of pertur --- -### Energy decomposition of free-boundary eigenmodes - -> I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M. Kulsrud, "An energy principle for hydromagnetic stability problems," -> *Proceedings of the Royal Society A* **244**, 17 (1958). -> DOI: [10.1098/rspa.1958.0023](https://doi.org/10.1098/rspa.1958.0023) - -The ideal-MHD energy principle: the standard form of the plasma potential energy δW_p that `Analysis.EnergyDecomposition` evaluates on the (ψ,θ) grid — see the [Analysis Module](analysis.md). - ---- - -> I. B. Bernstein, in *Basic Plasma Physics I*, Handbook of Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan, -> North-Holland, Amsterdam (1983), p. 421. - -Derives the effective-field (Bernstein) form of δW_p, in which the current and pressure terms are completed into |δB_eff|²/μ₀ − K|ξ_n|². - ---- - -> M. S. Chance, Y.-C. Sun, S. C. Jardin, C. E. Kessel, and M. Okabayashi, "MHD stability of tokamak plasmas," -> PPPL-CFP-2687, 2nd Symposium on Plasma Dynamics, Trieste (1992). -> Link: [OSTI 7234275](https://www.osti.gov/biblio/7234275) - -Writes both forms of δW_p side by side (Eqs. 12a–12c) and defines the destabilizing kernel K = K₁ + K₂ + K₃ with its shear–current, parallel-current and pressure–curvature pieces (Eqs. 15–17), the `K1`, `K2` and `K3` terms of the decomposition. - ---- - -> A. H. Boozer, "Physics of magnetically confined plasmas," -> *Reviews of Modern Physics* **76**, 1071 (2004). -> DOI: [10.1103/RevModPhys.76.1071](https://doi.org/10.1103/RevModPhys.76.1071) - -Sec. V.C derives the energy principle in the form used here, introduces the effective field C = δB + ξ_n(μ₀ j × n̂) (Eq. 98), the constraint (∇×C)·∇ψ = 0 that the perturbed current lie in the perturbed flux surfaces (Eq. 99) and the split w = C²/2μ₀ − w_d (Eq. 100). - ---- - ## InnerLayer Module > A. H. Glasser, Z. R. Wang, and J.-K. Park, "Computation of resistive instabilities by matched asymptotic expansions," diff --git a/src/Analysis/EnergyDecomposition.jl b/src/Analysis/EnergyDecomposition.jl index 2e713bbdc3..c22127b1eb 100644 --- a/src/Analysis/EnergyDecomposition.jl +++ b/src/Analysis/EnergyDecomposition.jl @@ -8,13 +8,21 @@ ForceFreeStates (FFS) solve on the (ψ, θ) grid, in two forms: standard form 2δW_p = ∫ [ |b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p) ] dV with b_eff = b + ξ_n(μ₀ j × n̂) and K₁ = |∇ψ|²σS (shear current), K₂ = μ₀B²σ² (parallel current), -K₃ = 2p'κ_ψ (pressure curvature) (Bernstein et al. 1958; Bernstein 1983; Chance et al. 1992; -Boozer 2004 Sec. V.C). The Fortran counterpart is `gpout_recon` with `effective_field_flag` and -`standard_form_flag`. Both totals are scaled by 2μ₀/ψ₀² so they read in the normalization of +K₃ = 2p'κ_ψ (pressure curvature). The Fortran counterpart is `gpout_recon` with `effective_field_flag` +and `standard_form_flag`. Both totals are scaled by 2μ₀/ψ₀² so they read in the normalization of `ffs.free_boundary.ep`. Call [`decompose_plasma_energy`](@ref) on a `ForceFreeStatesResult` and plot the result with [`plot_energy_decomposition`](@ref). + +## References + + - I. B. Bernstein, E. A. Frieman, M. D. Kruskal and R. M. Kulsrud, Proc. R. Soc. A 244, 17 (1958): + the energy principle and its standard form. + - I. B. Bernstein, Handbook of Plasma Physics Vol. 1, p. 421 (North-Holland, 1983): the effective-field form. + - M. S. Chance et al., PPPL-CFP-2687 (1992), Eqs. 12a–12c and 15–17: both forms and the split of K. + - A. H. Boozer, Rev. Mod. Phys. 76, 1071 (2004), Sec. V.C, Eqs. 98–100: the effective field and + (∇×b_eff)·∇ψ = 0. """ module EnergyDecomposition From fb2ab23b2bc15f4323d6271bd85f6ce18c02fed1 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 15:25:12 -0400 Subject: [PATCH 6/8] Analysis - FEATURE - Report the curl residual and the plasma-energy mismatch of each decomposed eigenmode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit decompose_plasma_energy now returns and checks two success measures per eigenmode: relative_error_effective/relative_error_standard, the mismatch of each δW_p total against ep[k] (warned above 5 %; it is the ψ quadrature of the knot grid and flags an mpsi that is too coarse), and curl_residual, the largest relative (∇×b_eff)·∇ψ over the knots (warned above 1e-8; an algebraic identity there that checks the b_eff assembly). The tests bound the returned values instead of recomputing the curl, and the Analysis page lists the checks. Co-Authored-By: Claude Opus 5.5 --- docs/src/analysis.md | 16 +++++++--- src/Analysis/EnergyDecomposition.jl | 45 +++++++++++++++++++++++++-- test/runtests_energy_decomposition.jl | 31 +++++------------- 3 files changed, 62 insertions(+), 30 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 5575c726eb..40c5535e0a 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -108,10 +108,6 @@ to roundoff, so the standard-form total differs from ``e_{p,k}`` only by the ``\ grid. ``\xi`` at each node is integrated from the nearest saved ODE step that no singular surface separates from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. -``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi = 0`` is an algebraic identity of -the construction on the knots. It checks how ``\delta\mathbf B_{\rm eff}`` is assembled from ``\delta\mathbf B``, -``\xi_n`` and the metric, not the radial ``\xi`` solution, the kernel ``K`` or the radial integral. - Choices that differ from the Fortran: ``B^2`` is the metric field strength ``(\chi'/\mathcal J)^2(g_{\theta\theta} + 2q g_{\theta\zeta} + q^2 g_{\zeta\zeta})``, which makes ``K_2|\xi_n|^2`` cancel against ``|\delta\mathbf B_{\rm eff}|^2``; the covariant field, ``\xi_n`` and ``\delta\mathbf B_{\rm eff}`` @@ -119,6 +115,18 @@ are formed pointwise in ``\theta`` rather than band-truncated; ``\partial_\theta spectral; the reconstruction is unregularized; and the totals integrate over the knots rather than the ODE `psifac` grid. +### Success checks + +Every call reports two checks per eigenmode and warns when either fails: + +- `relative_error_effective`, `relative_error_standard`: the mismatch of each total against + ``e_{p,k}``, warned above 5 %. It is the ``\psi`` quadrature error of the knot grid, so a large value + means `mpsi` is too coarse; it converges as in the table below. +- `curl_residual`: the largest relative ``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi`` over the + knots, warned above 1e-8. It is an algebraic identity there, so it checks how ``\delta\mathbf B_{\rm eff}`` + is assembled from ``\delta\mathbf B``, ``\xi_n`` and the metric, not the radial ``\xi`` solution, the + kernel ``K`` or the radial integral. + Relative error of both forms against ``e_{p,k}`` on the Solovev test fixture (`mpert = 8`): | `mpsi` | mode 1 effective / standard | mode 2 effective / standard | diff --git a/src/Analysis/EnergyDecomposition.jl b/src/Analysis/EnergyDecomposition.jl index c22127b1eb..8b33e955a1 100644 --- a/src/Analysis/EnergyDecomposition.jl +++ b/src/Analysis/EnergyDecomposition.jl @@ -210,6 +210,27 @@ function effective_field(geo::NamedTuple, f::NamedTuple, i::Int, k::Int) return (Jbp, Jct, Jcz, cp, bt + xw * jz, bz - xw * jt, bp, bt, bz, xw / (J * sqrt(d2))) end +# Relative (∇×b_eff)·∇ψ residual above which the effective field is reported as inconsistent. +const CURL_RESIDUAL_TOL = 1e-8 +# Relative mismatch of a δW_p total against ep[k] above which the radial grid is reported as too coarse. +const ENERGY_MISMATCH_TOL = 0.05 + +""" + curl_residual(c_theta, c_zeta, jac, ft, n) -> Float64 + +Relative (∇×b_eff)·∇ψ = (∂_θ b_eff,ζ − ∂_ζ b_eff,θ)/J on one surface (Fortran `gpeq_cveri`), from the +covariant components `c_theta`, `c_zeta` on the θ grid of `ft`, normalized by max|∂_θ b_eff,ζ/J|. On +the knots it is an algebraic identity, so it checks how b_eff is assembled from b, ξ_n and the metric, +not the radial ξ solution, the kernel K or the radial integral. +""" +function curl_residual(c_theta::Vector{ComplexF64}, c_zeta::Vector{ComplexF64}, jac::AbstractVector{Float64}, ft::FourierTransform, n::Int) + m = ft.mlow:(ft.mlow+ft.mpert-1) + theta_modes, zeta_modes = ft(c_theta), ft(c_zeta) + curl = inverse(ft, 2π * im .* (m .* zeta_modes .+ n .* theta_modes)) ./ jac + dtheta_c_zeta = inverse(ft, 2π * im .* m .* zeta_modes) ./ jac + return maximum(abs, curl) / maximum(abs, dtheta_c_zeta) +end + """ decompose_plasma_energy(ffs; eigenmodes=[1]) -> NamedTuple @@ -222,9 +243,12 @@ with the vacuum response (`vac_flag = true`). Fields of the result: - `density` - `(npsi, mtheta, 7, neig)` pointwise integrands, scaled by 2μ₀/ψ₀² - `profile` - `(npsi, 7, neig)` surface integrals ∮J·density dθ - `dW_effective`, `dW_standard` - ½∫ of the effective-field and standard-form combinations + - `relative_error_effective`, `relative_error_standard` - `(dW − real(ep[k])) / abs(ep[k])` of each total + - `curl_residual` - largest relative (∇×b_eff)·∇ψ over the knots (`curl_residual`) -The relative error against the DCON energy is `(dW − real(ep[k])) / abs(ep[k])` with -`ep = ffs.free_boundary.ep`. +These last two are the success checks of the decomposition: a curl residual above `CURL_RESIDUAL_TOL` +means b_eff was assembled inconsistently, and a mismatch against `ep[k]` above `ENERGY_MISMATCH_TOL` +means the radial grid is too coarse for the ψ quadrature (it converges with `mpsi`). Both are warned. """ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::AbstractVector{<:Integer}=[1]) ffs.integrator === :forward && ffs.solution !== nothing || throw(ArgumentError("energy decomposition needs a forward-integrator solve")) @@ -245,12 +269,15 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs scale = 2 * MU_0 / equil.psio^2 density = zeros(ComplexF64, length(psi), mtheta, length(TERMS), length(eigenmodes)) + curl = zeros(length(eigenmodes)) + c_theta, c_zeta = zeros(ComplexF64, mtheta), zeros(ComplexF64, mtheta) for (m, k) in enumerate(eigenmodes) f = eigenmode_fields(ffs, k, psi, states, ft) for i in eachindex(psi) p1 = equil.profiles.P_deriv(psi[i]) / MU_0 for j in 1:mtheta Jcp, Jct, Jcz, cp, ct, cz, bp, bt, bz, xi_n = effective_field(geo, f, i, j) + c_theta[j], c_zeta[j] = ct, cz J, xw = geo.jac[i, j], geo.jac[i, j] * f.xi[i, j] jt, jz = geo.mu0_j_theta[i, j], geo.mu0_j_zeta[i, j] divxi = f.xi1[i, j] + geo.jac_psi[i, j] / J * f.xi[i, j] + f.dJxt[i, j] / J - 2π * im * n * f.Jxz[i, j] / J @@ -263,6 +290,7 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs density[i, j, t, m] = scale * terms[t] end end + psi[i] in equil.rzphi_xs && (curl[m] = max(curl[m], curl_residual(c_theta, c_zeta, view(geo.jac, i, :), ft, n))) end end profile = dropdims(sum(density .* geo.jac; dims=2); dims=2) ./ mtheta @@ -270,7 +298,18 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs total(t, m) = FastInterpolations.integrate(cubic_interp(psi, real.(profile[:, t, m]); bc=CubicFit())) dW_effective = [0.5 * (total(1, m) - total(2, m) - total(3, m) - total(4, m)) for m in eachindex(eigenmodes)] dW_standard = [0.5 * (total(5, m) - total(6, m) + total(7, m)) for m in eachindex(eigenmodes)] - return (; psi, theta=equil.rzphi_ys[1:mtheta], R=geo.R, Z=geo.Z, terms=TERMS, eigenmodes=collect(eigenmodes), density, profile, dW_effective, dW_standard) + ep = real.(ffs.free_boundary.ep[eigenmodes]) + relative_error_effective = (dW_effective .- ep) ./ abs.(ep) + relative_error_standard = (dW_standard .- ep) ./ abs.(ep) + for (m, k) in enumerate(eigenmodes) + curl[m] > CURL_RESIDUAL_TOL && @warn "Energy decomposition: (∇×b_eff)·∇ψ of eigenmode $k does not vanish (relative residual $(curl[m]))" + for (form, err) in (("effective-field", relative_error_effective[m]), ("standard-form", relative_error_standard[m])) + abs(err) > ENERGY_MISMATCH_TOL && + @warn "Energy decomposition: $form δW_p of eigenmode $k differs from ep[$k] by $(round(100 * err; sigdigits=3)) %; refine mpsi" + end + end + return (; psi, theta=equil.rzphi_ys[1:mtheta], R=geo.R, Z=geo.Z, terms=TERMS, eigenmodes=collect(eigenmodes), density, profile, + dW_effective, dW_standard, relative_error_effective, relative_error_standard, curl_residual=curl) end """ diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl index bfe2760401..e4ffe63758 100644 --- a/test/runtests_energy_decomposition.jl +++ b/test/runtests_energy_decomposition.jl @@ -7,7 +7,6 @@ using Plots ENV["GKSwstype"] = "100" GPEC = GeneralizedPerturbedEquilibrium ED = GPEC.Analysis.EnergyDecomposition - FT = GPEC.Utilities.FourierTransforms template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") deck = TOML.parsefile(joinpath(template, "gpec.toml")) @@ -23,7 +22,6 @@ using Plots knots = findall(in(equil.rzphi_xs), psi) geo = ED.surface_geometry(ffs, psi) states = ED.node_states(ffs, psi) - ep = real.(ffs.free_boundary.ep) @testset "standard form equals the DCON Lagrangian on the knots" begin # On a knot the standard-form surface integral is 2/ψ₀² times the Euler-Lagrange quadratic form ξ'†Fξ' + 2Re(ξ'†Kξ) + ξ†Gξ. @@ -48,27 +46,14 @@ using Plots end end - @testset "both forms reproduce the DCON plasma energy" begin - for (m, k) in enumerate(ed.eigenmodes) - # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.3e-5). - @test ed.dW_effective[m] ≈ ed.dW_standard[m] rtol = 2e-4 - # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 2.0 %); it converges with mpsi. - @test abs(ed.dW_standard[m] - ep[k]) / abs(ep[k]) < 0.025 - end - end - - @testset "(∇×b_eff)·∇ψ vanishes on the knots" begin - # An algebraic identity of the b_eff construction on the knots, so this checks the effective-field formulas only. - ft = FT.FourierTransform(mtheta, ffs.mpert, ffs.mlow) - mvals = ffs.mlow:ffs.mhigh - f = ED.eigenmode_fields(ffs, 1, psi, states, ft) - for i in knots - c = [ED.effective_field(geo, f, i, k) for k in 1:mtheta] - c_theta, c_zeta = ft(getindex.(c, 5)), ft(getindex.(c, 6)) - curl = FT.inverse(ft, 2π * im .* (mvals .* c_zeta .+ ffs.nlow .* c_theta)) ./ geo.jac[i, :] - dtheta_c_zeta = FT.inverse(ft, 2π * im .* mvals .* c_zeta) ./ geo.jac[i, :] - @test maximum(abs, curl) < 1e-10 * maximum(abs, dtheta_c_zeta) - end + @testset "success checks: δW_p against ep and the curl residual" begin + # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 2.0 %); it converges with mpsi. + @test all(abs.(ed.relative_error_effective) .< 0.025) + @test all(abs.(ed.relative_error_standard) .< 0.025) + # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.3e-5). + @test ed.dW_effective ≈ ed.dW_standard rtol = 2e-4 + # (∇×b_eff)·∇ψ is an algebraic identity on the knots (measured 8e-13), a check of the b_eff assembly. + @test all(ed.curl_residual .< 1e-10) end @testset "flux-surface average of the shear" begin From dd146be108c23329517464933b69265ea8682878 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 19:47:50 -0400 Subject: [PATCH 7/8] Analysis - BUGFIX - Integrate the energy decomposition from the ODE start and normalize its curl check as Fortran does MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The ψ integrals started at the first saved ODE step. The forward sweep stores that step one ODE step after its start, so the integrals dropped [start, first knot]: 1.25 % of ep for a qlow = 2.5 start (Solovev fixture, mode 2, mpsi 16 to 64). They now start at max(psilow, first profile knot), where initialize_el_at_axis! starts. curl_residual was normalized surface by surface by max|∂θ b_eff,ζ/J|, which inflates it where the mode is weak. It is now the largest residual on the knots over the largest (|∂θ b_eff,ζ| + |∂ζ b_eff,θ|)/J, as in Fortran gpeq_cveri. On the DIII-D-like example with mpsi = 512 that is 2.2e-5 instead of 3.4e-3. The tests could not see the F′ terms or the g_θζ and g_ψζ terms: the Solovev fixture has F′ = 0 and uses PEST coordinates, and 5 of 33 injected errors survived. Changes to the tests: - The fixture now runs in Boozer coordinates. - A DIII-D-like case checks the standard form against the DCON Lagrangian and the curl residual. - A qlow case checks where the integration starts. 4 of the 5 are now caught. The fifth scales (∇·ξ)*(ξ·∇p), which vanishes for DCON displacements. The docs now say what each success check can and cannot detect. Co-Authored-By: Claude Opus 5.5 --- docs/src/analysis.md | 31 +++++---- src/Analysis/EnergyDecomposition.jl | 56 ++++++++------- test/runtests_energy_decomposition.jl | 99 ++++++++++++++++++--------- 3 files changed, 119 insertions(+), 67 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 40c5535e0a..2abc5d23b1 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -101,9 +101,10 @@ The result holds the seven integrands `(b_eff_sq, K1, K2, K3, b_sq, j_b_xi, divx ### Grid and accuracy -The radial nodes are the equilibrium ``\psi`` knots inside the solution domain, plus the solution edge. -On the knots all geometry is read from `ffs.metric`, and the standard-form surface integral equals -``2/\psi_0^2`` times the DCON Lagrangian ``\xi'^\dagger F\xi' + 2\,\mathrm{Re}(\xi'^\dagger K\xi) + \xi^\dagger G\xi`` +The radial nodes are the equilibrium ``\psi`` knots inside the solution domain, plus each end of the +domain that falls between knots: the start of the ODE sweep (`max(psilow, first knot)`, or a `qlow` +crossing) and the solution edge. On the knots all geometry is read from `ffs.metric`, and the +standard-form surface integral equals ``2/\psi_0^2`` times the DCON Lagrangian ``\xi'^\dagger F\xi' + 2\,\mathrm{Re}(\xi'^\dagger K\xi) + \xi^\dagger G\xi`` to roundoff, so the standard-form total differs from ``e_{p,k}`` only by the ``\psi`` quadrature of the knot grid. ``\xi`` at each node is integrated from the nearest saved ODE step that no singular surface separates from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. @@ -121,20 +122,26 @@ Every call reports two checks per eigenmode and warns when either fails: - `relative_error_effective`, `relative_error_standard`: the mismatch of each total against ``e_{p,k}``, warned above 5 %. It is the ``\psi`` quadrature error of the knot grid, so a large value - means `mpsi` is too coarse; it converges as in the table below. -- `curl_residual`: the largest relative ``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi`` over the - knots, warned above 1e-8. It is an algebraic identity there, so it checks how ``\delta\mathbf B_{\rm eff}`` - is assembled from ``\delta\mathbf B``, ``\xi_n`` and the metric, not the radial ``\xi`` solution, the - kernel ``K`` or the radial integral. + means `mpsi` is too coarse; it converges as in the table below. It does not check the algebra: an + error that moves a total by a few percent, such as a factor 2 lost in ``K_3`` (4.8 % on the test + fixture), stays below 5 %. +- `curl_residual`: the largest ``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi`` over the knots, + relative to the largest ``(|\partial_\theta \delta B_{{\rm eff},\zeta}| + |\partial_\zeta \delta B_{{\rm eff},\theta}|)/\mathcal J`` + it cancels, warned above 1e-8; normalization and threshold are those of the Fortran `gpeq_cveri`. It is + an algebraic identity on the knots, so it checks how ``\delta\mathbf B_{\rm eff}`` is assembled from + ``\delta\mathbf B``, ``\xi_n``, the metric and the equilibrium current, not the radial ``\xi`` solution, + the kernel ``K`` or the radial integral. It is roundoff (4e-14) on the Solovev test fixture and 4e-12 + on the DIII-D-like example, but that example reaches 2e-5 with a uniform `mpsi = 512` grid and 5e-5 in + PEST coordinates, all from knots inside ``\psi = 0.1``, without an error in ``\delta\mathbf B_{\rm eff}``. Relative error of both forms against ``e_{p,k}`` on the Solovev test fixture (`mpert = 8`): | `mpsi` | mode 1 effective / standard | mode 2 effective / standard | |---|---|---| -| 16 | +1.39e-2 / +1.39e-2 | −1.98e-2 / −1.97e-2 | -| 32 | −1.42e-3 / −1.42e-3 | +1.07e-2 / +1.07e-2 | -| 64 | +2.10e-4 / +2.10e-4 | +6.43e-4 / +6.43e-4 | -| 128 | +1.62e-4 / +1.62e-4 | −4.18e-5 / −4.18e-5 | +| 16 | +1.38e-2 / +1.38e-2 | −1.99e-2 / −1.98e-2 | +| 32 | −1.53e-3 / −1.53e-3 | +1.01e-2 / +1.01e-2 | +| 64 | +2.40e-4 / +2.40e-4 | +8.36e-4 / +8.36e-4 | +| 128 | +1.66e-4 / +1.66e-4 | −1.63e-5 / −1.63e-5 | ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Analysis.EnergyDecomposition] diff --git a/src/Analysis/EnergyDecomposition.jl b/src/Analysis/EnergyDecomposition.jl index 8b33e955a1..57d9e00e31 100644 --- a/src/Analysis/EnergyDecomposition.jl +++ b/src/Analysis/EnergyDecomposition.jl @@ -41,20 +41,23 @@ Names of the seven density terms, in the order of the third axis of `density` an """ const TERMS = (:b_eff_sq, :K1, :K2, :K3, :b_sq, :j_b_xi, :divxi_gradp) -# ψ_N gap below which the last equilibrium knot counts as the solution edge. +# ψ_N gap below which an equilibrium knot counts as an end of the solution domain. const EDGE_KNOT_TOL = 1e-10 """ radial_nodes(ffs) -> Vector{Float64} -The equilibrium ψ knots inside the ξ solution domain, plus the solution edge when the last knot -falls short of it. On the knots the metric, the Euler-Lagrange matrices and the geometry are exact. +The equilibrium ψ knots inside the ξ solution domain, plus each end of the domain (the ODE start, +e.g. a `qlow` crossing, and the solution edge) that falls between knots. On the knots the metric, +the Euler-Lagrange matrices and the geometry are exact. """ function radial_nodes(ffs::FFS.ForceFreeStatesResult) sol = ffs.solution - lo, hi = sol.psi_store[1], sol.psi_store[sol.step] + # The sweep starts at max(psilow, first profile knot), as initialize_el_at_axis! does, but saves its first state one step later. + lo, hi = max(ffs.psilow, ffs.equil.profiles.xs[1]), sol.psi_store[sol.step] psi = [x for x in ffs.equil.rzphi_xs if lo <= x <= hi] - (isempty(psi) || psi[end] < hi - EDGE_KNOT_TOL) && push!(psi, hi) + (isempty(psi) || psi[1] > lo + EDGE_KNOT_TOL) && pushfirst!(psi, lo) + psi[end] < hi - EDGE_KNOT_TOL && push!(psi, hi) return psi end @@ -210,25 +213,25 @@ function effective_field(geo::NamedTuple, f::NamedTuple, i::Int, k::Int) return (Jbp, Jct, Jcz, cp, bt + xw * jz, bz - xw * jt, bp, bt, bz, xw / (J * sqrt(d2))) end -# Relative (∇×b_eff)·∇ψ residual above which the effective field is reported as inconsistent. +# Relative (∇×b_eff)·∇ψ residual above which the effective field is reported as inconsistent (Fortran cveri_tol). const CURL_RESIDUAL_TOL = 1e-8 # Relative mismatch of a δW_p total against ep[k] above which the radial grid is reported as too coarse. const ENERGY_MISMATCH_TOL = 0.05 """ - curl_residual(c_theta, c_zeta, jac, ft, n) -> Float64 + curl_maxima(c_theta, c_zeta, jac, ft, n) -> (residual, scale) -Relative (∇×b_eff)·∇ψ = (∂_θ b_eff,ζ − ∂_ζ b_eff,θ)/J on one surface (Fortran `gpeq_cveri`), from the -covariant components `c_theta`, `c_zeta` on the θ grid of `ft`, normalized by max|∂_θ b_eff,ζ/J|. On -the knots it is an algebraic identity, so it checks how b_eff is assembled from b, ξ_n and the metric, -not the radial ξ solution, the kernel K or the radial integral. +Largest |(∇×b_eff)·∇ψ| = |∂_θ b_eff,ζ − ∂_ζ b_eff,θ|/J on one surface and the largest +(|∂_θ b_eff,ζ| + |∂_ζ b_eff,θ|)/J that it cancels, from the covariant components `c_theta`, `c_zeta` +on the θ grid of `ft`. Perturbations vary as exp(2πi(mθ − nζ)), so ∂_ζ → −2πin. As in Fortran +`gpeq_cveri`, the reported residual is the largest `residual` over all surfaces divided by the +largest `scale`. """ -function curl_residual(c_theta::Vector{ComplexF64}, c_zeta::Vector{ComplexF64}, jac::AbstractVector{Float64}, ft::FourierTransform, n::Int) +function curl_maxima(c_theta::Vector{ComplexF64}, c_zeta::Vector{ComplexF64}, jac::AbstractVector{Float64}, ft::FourierTransform, n::Int) m = ft.mlow:(ft.mlow+ft.mpert-1) - theta_modes, zeta_modes = ft(c_theta), ft(c_zeta) - curl = inverse(ft, 2π * im .* (m .* zeta_modes .+ n .* theta_modes)) ./ jac - dtheta_c_zeta = inverse(ft, 2π * im .* m .* zeta_modes) ./ jac - return maximum(abs, curl) / maximum(abs, dtheta_c_zeta) + dtheta_c_zeta = inverse(ft, 2π * im .* m .* ft(c_zeta)) ./ jac + dzeta_c_theta = inverse(ft, -2π * im * n .* ft(c_theta)) ./ jac + return maximum(abs, dtheta_c_zeta .- dzeta_c_theta), maximum(abs.(dtheta_c_zeta) .+ abs.(dzeta_c_theta)) end """ @@ -244,11 +247,14 @@ with the vacuum response (`vac_flag = true`). Fields of the result: - `profile` - `(npsi, 7, neig)` surface integrals ∮J·density dθ - `dW_effective`, `dW_standard` - ½∫ of the effective-field and standard-form combinations - `relative_error_effective`, `relative_error_standard` - `(dW − real(ep[k])) / abs(ep[k])` of each total - - `curl_residual` - largest relative (∇×b_eff)·∇ψ over the knots (`curl_residual`) - -These last two are the success checks of the decomposition: a curl residual above `CURL_RESIDUAL_TOL` -means b_eff was assembled inconsistently, and a mismatch against `ep[k]` above `ENERGY_MISMATCH_TOL` -means the radial grid is too coarse for the ψ quadrature (it converges with `mpsi`). Both are warned. + - `curl_residual` - largest (∇×b_eff)·∇ψ on the knots relative to the largest term it cancels (`curl_maxima`) + +These last two are the success checks of the decomposition, and both are warned. On the knots +(∇×b_eff)·∇ψ vanishes identically, so a curl residual above `CURL_RESIDUAL_TOL` flags b_eff inconsistent +with the metric and the equilibrium current; near the axis of a reconstructed equilibrium the radial grid +and the coordinates can raise it too. A mismatch against `ep[k]` above `ENERGY_MISMATCH_TOL` means the +radial grid is too coarse for the ψ quadrature (it converges with `mpsi`); it is not a check of the +algebra, since an error that moves a total by a few percent stays below it. """ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::AbstractVector{<:Integer}=[1]) ffs.integrator === :forward && ffs.solution !== nothing || throw(ArgumentError("energy decomposition needs a forward-integrator solve")) @@ -269,7 +275,7 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs scale = 2 * MU_0 / equil.psio^2 density = zeros(ComplexF64, length(psi), mtheta, length(TERMS), length(eigenmodes)) - curl = zeros(length(eigenmodes)) + curl_max, curl_scale = zeros(length(eigenmodes)), zeros(length(eigenmodes)) c_theta, c_zeta = zeros(ComplexF64, mtheta), zeros(ComplexF64, mtheta) for (m, k) in enumerate(eigenmodes) f = eigenmode_fields(ffs, k, psi, states, ft) @@ -290,9 +296,13 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs density[i, j, t, m] = scale * terms[t] end end - psi[i] in equil.rzphi_xs && (curl[m] = max(curl[m], curl_residual(c_theta, c_zeta, view(geo.jac, i, :), ft, n))) + if psi[i] in equil.rzphi_xs + residual, cancelled = curl_maxima(c_theta, c_zeta, view(geo.jac, i, :), ft, n) + curl_max[m], curl_scale[m] = max(curl_max[m], residual), max(curl_scale[m], cancelled) + end end end + curl = curl_max ./ curl_scale profile = dropdims(sum(density .* geo.jac; dims=2); dims=2) ./ mtheta # The imaginary parts of the surface integrals are roundoff, so the ψ integrals take the real part. total(t, m) = FastInterpolations.integrate(cubic_interp(psi, real.(profile[:, t, m]); bc=CubicFit())) diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl index e4ffe63758..41f1573f8e 100644 --- a/test/runtests_energy_decomposition.jl +++ b/test/runtests_energy_decomposition.jl @@ -1,58 +1,75 @@ using TOML using Plots -# Analysis.EnergyDecomposition on the coarse Solovev fixture (mpsi = 16, forward integrator, vac_flag = true): -# one solve and one decomposition, shared by every testset. +# Analysis.EnergyDecomposition on the coarse Solovev fixture (mpsi = 16, forward integrator, vac_flag = true) in Boozer +# coordinates, which unlike the deck's PEST coordinates have nonzero g_θζ and g_ψζ: one solve and one decomposition, +# shared by every testset but the last two. @testset "Analysis.EnergyDecomposition" begin ENV["GKSwstype"] = "100" GPEC = GeneralizedPerturbedEquilibrium ED = GPEC.Analysis.EnergyDecomposition - template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - deck = TOML.parsefile(joinpath(template, "gpec.toml")) - equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), - GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])) - wall = GPEC.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in deck["Wall"])...) - ffs_kwargs = Dict(Symbol(k) => v for (k, v) in deck["ForceFreeStates"] if !(k in ("integrator", "nn_low", "nn_high"))) - ffs = solve(equil, Forward(); nn=1, wall=wall, dir_path=mktempdir(), ffs_kwargs...) + # Forward solve of a deck for toroidal mode number 1, with `[Equilibrium]` and `[ForceFreeStates]` overrides. + function forward_solve(dir; eq_overrides=Dict{String,Any}(), ffs_overrides=Dict{Symbol,Any}()) + deck = TOML.parsefile(joinpath(dir, "gpec.toml")) + merge!(deck["Equilibrium"], eq_overrides) + config = GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], dir) + equil = + haskey(deck, "SOL_INPUT") ? GPEC.Equilibrium.PlasmaEquilibrium(config, GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])) : + GPEC.Equilibrium.PlasmaEquilibrium(config) + wall = GPEC.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in deck["Wall"])...) + ffs_kwargs = Dict{Symbol,Any}(Symbol(k) => v for (k, v) in deck["ForceFreeStates"] if !(k in ("integrator", "nn_low", "nn_high"))) + return solve(equil, Forward(); nn=1, wall=wall, dir_path=mktempdir(), merge(ffs_kwargs, ffs_overrides)...) + end + + # Standard-form surface integrals and 2/ψ₀² times the Euler-Lagrange quadratic form ξ'†Fξ' + 2Re(ξ'†Kξ) + ξ†Gξ, + # (knot, eigenmode) matrices over the knots among the radial nodes. + function standard_and_lagrangian(ffs, ed) + knots = findall(in(ffs.equil.rzphi_xs), ed.psi) + u1, du1, _ = ED.node_states(ffs, ed.psi) + N = ffs.numpert_total + F, F_lower, K, G = (zeros(ComplexF64, N, N) for _ in 1:4) + standard = [real(ed.profile[i, 5, m] - ed.profile[i, 6, m] + ed.profile[i, 7, m]) for i in knots, m in eachindex(ed.eigenmodes)] + lagrangian = similar(standard) + for (m, k) in enumerate(ed.eigenmodes) + alpha = ffs.solution.u_store[:, :, 1, ffs.solution.step] \ Vector{ComplexF64}(ffs.free_boundary.wt[:, k]) + for (j, i) in enumerate(knots) + m_nq = (ffs.mlow:ffs.mhigh) .- ffs.nlow * ffs.equil.profiles.q_spline(ed.psi[i]) + ffs.mats.ideal.F_spline_lower(vec(F_lower), ed.psi[i]) + F .= LowerTriangular(F_lower) * LowerTriangular(F_lower)' .* m_nq .* m_nq' + ffs.mats.ideal.K_spline(vec(K), ed.psi[i]) + K .*= m_nq + ffs.mats.ideal.G_spline(vec(G), ed.psi[i]) + xi, dxi = u1[:, :, i] * alpha, du1[:, :, i] * alpha + lagrangian[j, m] = 2 / ffs.equil.psio^2 * real(dot(dxi, F, dxi) + 2 * dot(dxi, K, xi) + dot(xi, G, xi)) + end + end + return standard, lagrangian + end + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + boozer = Dict{String,Any}("jac_type" => "boozer") + ffs = forward_solve(template; eq_overrides=boozer) + equil = ffs.equil ed = ED.decompose_plasma_energy(ffs; eigenmodes=[1, 2]) psi = ed.psi mtheta = length(ed.theta) knots = findall(in(equil.rzphi_xs), psi) geo = ED.surface_geometry(ffs, psi) - states = ED.node_states(ffs, psi) @testset "standard form equals the DCON Lagrangian on the knots" begin - # On a knot the standard-form surface integral is 2/ψ₀² times the Euler-Lagrange quadratic form ξ'†Fξ' + 2Re(ξ'†Kξ) + ξ†Gξ. - N = ffs.numpert_total - F, K, G = (zeros(ComplexF64, N, N) for _ in 1:3) - u1, du1, _ = states - sol = ffs.solution - for (m, k) in enumerate(ed.eigenmodes) - alpha = sol.u_store[:, :, 1, sol.step] \ Vector{ComplexF64}(ffs.free_boundary.wt[:, k]) - for i in knots - sf = (ffs.mlow:ffs.mhigh) .- ffs.nlow * equil.profiles.q_spline(psi[i]) - ffs.mats.ideal.F_spline_lower(vec(G), psi[i]) - F .= LowerTriangular(G) * LowerTriangular(G)' .* sf .* sf' - ffs.mats.ideal.K_spline(vec(K), psi[i]) - K .*= sf - ffs.mats.ideal.G_spline(vec(G), psi[i]) - xi, dxi = u1[:, :, i] * alpha, du1[:, :, i] * alpha - lagrangian = real(dot(dxi, F, dxi) + 2 * dot(dxi, K, xi) + dot(xi, G, xi)) - standard = real(ed.profile[i, 5, m] - ed.profile[i, 6, m] + ed.profile[i, 7, m]) - @test standard ≈ 2 / equil.psio^2 * lagrangian rtol = 1e-9 - end + for (standard, lagrangian) in zip(standard_and_lagrangian(ffs, ed)...) + @test standard ≈ lagrangian rtol = 1e-9 end end @testset "success checks: δW_p against ep and the curl residual" begin - # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 2.0 %); it converges with mpsi. + # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 1.9 %); it converges with mpsi. @test all(abs.(ed.relative_error_effective) .< 0.025) @test all(abs.(ed.relative_error_standard) .< 0.025) - # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.3e-5). + # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.2e-5). @test ed.dW_effective ≈ ed.dW_standard rtol = 2e-4 - # (∇×b_eff)·∇ψ is an algebraic identity on the knots (measured 8e-13), a check of the b_eff assembly. + # (∇×b_eff)·∇ψ is an algebraic identity on the knots (measured 4e-14), a check of the b_eff assembly. @test all(ed.curl_residual .< 1e-10) end @@ -90,4 +107,22 @@ using Plots @test ED.plot_energy_decomposition(ed; eigenmode=2, save_path=path) isa Plots.Plot @test isfile(path) end + + @testset "ψ integrals start at the ODE start" begin + # A qlow start falls between knots, and the sweep saves its first state one ODE step after it. + qffs = forward_solve(template; eq_overrides=boozer, ffs_overrides=Dict{Symbol,Any}(:qlow => 2.5)) + qed = ED.decompose_plasma_energy(qffs; eigenmodes=[1, 2]) + @test qed.psi[1] == qffs.psilow < qffs.solution.psi_store[1] + @test all(abs.(qed.relative_error_standard) .< 0.025) + end + + @testset "F′ ≠ 0: the DIII-D-like reconstruction" begin + # The Solovev fixture has F′ = 0, so the poloidal current and the qF′ part of j^ζ are exercised only here. + dffs = forward_solve(joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example")) + ded = ED.decompose_plasma_energy(dffs) + standard, lagrangian = standard_and_lagrangian(dffs, ded) + # Judged against the largest density (measured 7e-11), since single knots where the mode is weak agree less closely. + @test maximum(abs, standard .- lagrangian) < 1e-8 * maximum(abs, lagrangian) + @test all(ded.curl_residual .< 1e-10) # measured 4e-12 + end end From 4daa3b1daa4cb466690e8da87a61af773dc78e68 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Fri, 30 Oct 2026 22:51:53 -0400 Subject: [PATCH 8/8] Analysis - BUGFIX - Integrate the energy decomposition over the saved ODE steps as well as the equilibrium knots MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The ψ integrals ran over the equilibrium knots only. The knots follow the equilibrium profiles, not the eigenmode, so near the rational surfaces a spline through the knot densities misses the structure of ξ. Against ep, the standard form was off by: - 0.3–0.5 % on the DIII-D-like example, at every knot grid; - 35 % for eigenmode 1 of the LAR reversed-shear example at mpsi = 128; - 5 % on the a10 example at mpsi = 64. On DIII-D the whole error was knot quadrature: Re(ξ†u₂) at the edge equals ep, the jumps at the five rational surfaces sum to 5e-6, and ∫L over the saved ODE steps equals ep to 1e-5. The nodes are now the ODE start and every saved step, as in Fortran gpout_recon, plus the knots, where the geometry is exact and the curl check runs. The standard form now matches ep to: - 1e-5 on DIII-D; - 6e-5 on a10 (mpsi 64); - 1e-4 on TJ (mpsi 256); - 7e-4 on LAR (mpsi 128). A saved step needs no re-integration, so the extra nodes are cheap. The density arrays grow with them; DIII-D goes from 281 to 1444 nodes. Test changes: - The fixture runs at mpsi = 32, since between its own 16 knots the splined metric moves the totals by 3 %. - The success-check tolerances tighten to the new measured values. - A DIII-D check guards the fix. The effective-field total still differs from the standard-form one by a grid-independent amount: −0.4 % on DIII-D and −9.5 % for LAR eigenmode 1. The docs now say so. Co-Authored-By: Claude Opus 5.5 --- docs/src/analysis.md | 45 ++++++++++++++++----------- src/Analysis/EnergyDecomposition.jl | 29 +++++++++-------- test/runtests_energy_decomposition.jl | 38 +++++++++++----------- 3 files changed, 62 insertions(+), 50 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 2abc5d23b1..65071b9a21 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -101,47 +101,54 @@ The result holds the seven integrands `(b_eff_sq, K1, K2, K3, b_sq, j_b_xi, divx ### Grid and accuracy -The radial nodes are the equilibrium ``\psi`` knots inside the solution domain, plus each end of the -domain that falls between knots: the start of the ODE sweep (`max(psilow, first knot)`, or a `qlow` -crossing) and the solution edge. On the knots all geometry is read from `ffs.metric`, and the -standard-form surface integral equals ``2/\psi_0^2`` times the DCON Lagrangian ``\xi'^\dagger F\xi' + 2\,\mathrm{Re}(\xi'^\dagger K\xi) + \xi^\dagger G\xi`` -to roundoff, so the standard-form total differs from ``e_{p,k}`` only by the ``\psi`` quadrature of the knot -grid. ``\xi`` at each node is integrated from the nearest saved ODE step that no singular surface -separates from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. +The radial nodes are the start of the forward ODE sweep (`max(psilow, first knot)`, or a `qlow` crossing) +and every step it saved, plus the equilibrium ``\psi`` knots inside that domain. The saved steps follow +``\xi`` into the rational surfaces, where the densities vary too fast for the knots alone: integrated over +the knots only, the standard form misses ``e_{p,k}`` by 0.3–0.5 % on the DIII-D-like example whatever the +knot grid, and by 35 % for eigenmode 1 of the LAR reversed-shear example at `mpsi = 128`. On the knots +all geometry is read from `ffs.metric`, and the standard-form surface integral equals ``2/\psi_0^2`` times +the DCON Lagrangian +``\xi'^\dagger F\xi' + 2\,\mathrm{Re}(\xi'^\dagger K\xi) + \xi^\dagger G\xi`` to roundoff; between the knots +the geometry is splined in ``\psi`` from the knot values. ``\xi`` at a saved step is the stored state, and +at a knot or the ODE start it is integrated from the nearest saved step that no singular surface separates +from it (`ForceFreeStates.solution_state_at`), not splined between saved steps. Choices that differ from the Fortran: ``B^2`` is the metric field strength ``(\chi'/\mathcal J)^2(g_{\theta\theta} + 2q g_{\theta\zeta} + q^2 g_{\zeta\zeta})``, which makes ``K_2|\xi_n|^2`` cancel against ``|\delta\mathbf B_{\rm eff}|^2``; the covariant field, ``\xi_n`` and ``\delta\mathbf B_{\rm eff}`` are formed pointwise in ``\theta`` rather than band-truncated; ``\partial_\theta(\mathcal J\xi^\theta)`` is -spectral; the reconstruction is unregularized; and the totals integrate over the knots rather than the ODE -`psifac` grid. +spectral; the reconstruction is unregularized; and the totals add the equilibrium knots to the ODE `psifac` +grid. ### Success checks Every call reports two checks per eigenmode and warns when either fails: - `relative_error_effective`, `relative_error_standard`: the mismatch of each total against - ``e_{p,k}``, warned above 5 %. It is the ``\psi`` quadrature error of the knot grid, so a large value - means `mpsi` is too coarse; it converges as in the table below. It does not check the algebra: an - error that moves a total by a few percent, such as a factor 2 lost in ``K_3`` (4.8 % on the test - fixture), stays below 5 %. + ``e_{p,k}``, warned above 5 %. For the standard form it is 1e-5 on the DIII-D-like example and grows + when `mpsi` is too coarse for the metric splined between the knots; it converges as in the table + below. The effective-field total can also differ from the standard-form one by a grid-independent + amount whose cause is still open: −0.4 % on the DIII-D-like example and −9.5 % for eigenmode 1 of the + LAR reversed-shear example. The mismatch does not check the algebra: an error that moves a total by a + few percent, such as a factor 2 lost in ``K_3`` (3.9 % on the test fixture), stays below 5 %. - `curl_residual`: the largest ``(\nabla\times\delta\mathbf B_{\rm eff})\cdot\nabla\psi`` over the knots, relative to the largest ``(|\partial_\theta \delta B_{{\rm eff},\zeta}| + |\partial_\zeta \delta B_{{\rm eff},\theta}|)/\mathcal J`` it cancels, warned above 1e-8; normalization and threshold are those of the Fortran `gpeq_cveri`. It is an algebraic identity on the knots, so it checks how ``\delta\mathbf B_{\rm eff}`` is assembled from ``\delta\mathbf B``, ``\xi_n``, the metric and the equilibrium current, not the radial ``\xi`` solution, - the kernel ``K`` or the radial integral. It is roundoff (4e-14) on the Solovev test fixture and 4e-12 + the kernel ``K`` or the radial integral. It is roundoff (7e-14) on the Solovev test fixture and 4e-12 on the DIII-D-like example, but that example reaches 2e-5 with a uniform `mpsi = 512` grid and 5e-5 in PEST coordinates, all from knots inside ``\psi = 0.1``, without an error in ``\delta\mathbf B_{\rm eff}``. -Relative error of both forms against ``e_{p,k}`` on the Solovev test fixture (`mpert = 8`): +Relative error of both forms against ``e_{p,k}`` on the Solovev test fixture in its PEST coordinates +(`mpert = 8`): | `mpsi` | mode 1 effective / standard | mode 2 effective / standard | |---|---|---| -| 16 | +1.38e-2 / +1.38e-2 | −1.99e-2 / −1.98e-2 | -| 32 | −1.53e-3 / −1.53e-3 | +1.01e-2 / +1.01e-2 | -| 64 | +2.40e-4 / +2.40e-4 | +8.36e-4 / +8.36e-4 | -| 128 | +1.66e-4 / +1.66e-4 | −1.63e-5 / −1.63e-5 | +| 16 | −4.44e-3 / −4.50e-3 | −3.19e-2 / −3.23e-2 | +| 32 | −2.17e-6 / −8.19e-6 | −3.35e-5 / −6.86e-5 | +| 64 | −1.07e-6 / −1.37e-6 | −9.14e-6 / −1.11e-5 | +| 128 | −2.01e-6 / −1.98e-6 | −1.72e-5 / −1.72e-5 | ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Analysis.EnergyDecomposition] diff --git a/src/Analysis/EnergyDecomposition.jl b/src/Analysis/EnergyDecomposition.jl index 57d9e00e31..358576ca48 100644 --- a/src/Analysis/EnergyDecomposition.jl +++ b/src/Analysis/EnergyDecomposition.jl @@ -41,23 +41,25 @@ Names of the seven density terms, in the order of the third axis of `density` an """ const TERMS = (:b_eff_sq, :K1, :K2, :K3, :b_sq, :j_b_xi, :divxi_gradp) -# ψ_N gap below which an equilibrium knot counts as an end of the solution domain. -const EDGE_KNOT_TOL = 1e-10 +# ψ_N gap below which a saved ODE step or the ODE start counts as the same radial node as a knot. +const SAME_NODE_TOL = 1e-10 """ radial_nodes(ffs) -> Vector{Float64} -The equilibrium ψ knots inside the ξ solution domain, plus each end of the domain (the ODE start, -e.g. a `qlow` crossing, and the solution edge) that falls between knots. On the knots the metric, -the Euler-Lagrange matrices and the geometry are exact. +The start and every saved step of the forward ODE sweep (e.g. from a `qlow` crossing to the edge), +plus the equilibrium ψ knots inside that domain. The saved steps follow ξ into the rational surfaces, +where the densities vary too fast for the knots alone; on the knots the metric, the Euler-Lagrange +matrices and the geometry are exact. """ function radial_nodes(ffs::FFS.ForceFreeStatesResult) sol = ffs.solution + steps = sol.psi_store[1:sol.step] # The sweep starts at max(psilow, first profile knot), as initialize_el_at_axis! does, but saves its first state one step later. - lo, hi = max(ffs.psilow, ffs.equil.profiles.xs[1]), sol.psi_store[sol.step] - psi = [x for x in ffs.equil.rzphi_xs if lo <= x <= hi] - (isempty(psi) || psi[1] > lo + EDGE_KNOT_TOL) && pushfirst!(psi, lo) - psi[end] < hi - EDGE_KNOT_TOL && push!(psi, hi) + lo = max(ffs.psilow, ffs.equil.profiles.xs[1]) + knots = [x for x in ffs.equil.rzphi_xs if lo <= x <= steps[end]] + psi = sort!([knots; [p for p in steps if all(x -> abs(p - x) > SAME_NODE_TOL, knots)]]) + psi[1] > lo + SAME_NODE_TOL && pushfirst!(psi, lo) return psi end @@ -252,9 +254,10 @@ with the vacuum response (`vac_flag = true`). Fields of the result: These last two are the success checks of the decomposition, and both are warned. On the knots (∇×b_eff)·∇ψ vanishes identically, so a curl residual above `CURL_RESIDUAL_TOL` flags b_eff inconsistent with the metric and the equilibrium current; near the axis of a reconstructed equilibrium the radial grid -and the coordinates can raise it too. A mismatch against `ep[k]` above `ENERGY_MISMATCH_TOL` means the -radial grid is too coarse for the ψ quadrature (it converges with `mpsi`); it is not a check of the -algebra, since an error that moves a total by a few percent stays below it. +and the coordinates can raise it too. A standard-form mismatch against `ep[k]` above `ENERGY_MISMATCH_TOL` +means the knots are too coarse for the metric splined between them (it converges with `mpsi`); the +effective-field total can also differ by a grid-independent amount on some equilibria. The mismatch is +not a check of the algebra, since an error that moves a total by a few percent stays below it. """ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::AbstractVector{<:Integer}=[1]) ffs.integrator === :forward && ffs.solution !== nothing || throw(ArgumentError("energy decomposition needs a forward-integrator solve")) @@ -265,7 +268,7 @@ function decompose_plasma_energy(ffs::FFS.ForceFreeStatesResult; eigenmodes::Abs all(k -> 1 <= k <= nwt, eigenmodes) || throw(ArgumentError("eigenmodes must lie in 1:$nwt")) equil = ffs.equil - # Deviation from Fortran gpout_recon: the radial integral runs over the equilibrium knots, not the ODE psifac grid. + # Deviation from Fortran gpout_recon: the radial integral adds the equilibrium knots to the ODE psifac grid. psi = radial_nodes(ffs) mtheta = length(equil.rzphi_ys) - 1 geo = surface_geometry(ffs, psi) diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl index 41f1573f8e..5df8b66490 100644 --- a/test/runtests_energy_decomposition.jl +++ b/test/runtests_energy_decomposition.jl @@ -1,9 +1,9 @@ using TOML using Plots -# Analysis.EnergyDecomposition on the coarse Solovev fixture (mpsi = 16, forward integrator, vac_flag = true) in Boozer -# coordinates, which unlike the deck's PEST coordinates have nonzero g_θζ and g_ψζ: one solve and one decomposition, -# shared by every testset but the last two. +# Analysis.EnergyDecomposition on the Solovev fixture (forward integrator, vac_flag = true) at mpsi = 32, since its own 16 knots +# leave the metric between them too coarse, and in Boozer coordinates, which unlike the deck's PEST coordinates have nonzero +# g_θζ and g_ψζ: one solve and one decomposition, shared by every testset but the last two. @testset "Analysis.EnergyDecomposition" begin ENV["GKSwstype"] = "100" GPEC = GeneralizedPerturbedEquilibrium @@ -26,7 +26,7 @@ using Plots # (knot, eigenmode) matrices over the knots among the radial nodes. function standard_and_lagrangian(ffs, ed) knots = findall(in(ffs.equil.rzphi_xs), ed.psi) - u1, du1, _ = ED.node_states(ffs, ed.psi) + u1, du1, _ = ED.node_states(ffs, ed.psi[knots]) N = ffs.numpert_total F, F_lower, K, G = (zeros(ComplexF64, N, N) for _ in 1:4) standard = [real(ed.profile[i, 5, m] - ed.profile[i, 6, m] + ed.profile[i, 7, m]) for i in knots, m in eachindex(ed.eigenmodes)] @@ -40,7 +40,7 @@ using Plots ffs.mats.ideal.K_spline(vec(K), ed.psi[i]) K .*= m_nq ffs.mats.ideal.G_spline(vec(G), ed.psi[i]) - xi, dxi = u1[:, :, i] * alpha, du1[:, :, i] * alpha + xi, dxi = u1[:, :, j] * alpha, du1[:, :, j] * alpha lagrangian[j, m] = 2 / ffs.equil.psio^2 * real(dot(dxi, F, dxi) + 2 * dot(dxi, K, xi) + dot(xi, G, xi)) end end @@ -48,8 +48,8 @@ using Plots end template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - boozer = Dict{String,Any}("jac_type" => "boozer") - ffs = forward_solve(template; eq_overrides=boozer) + overrides = Dict{String,Any}("jac_type" => "boozer", "mpsi" => 32) + ffs = forward_solve(template; eq_overrides=overrides) equil = ffs.equil ed = ED.decompose_plasma_energy(ffs; eigenmodes=[1, 2]) psi = ed.psi @@ -64,12 +64,12 @@ using Plots end @testset "success checks: δW_p against ep and the curl residual" begin - # ψ-quadrature error of the mpsi = 16 knot grid (measured 1.4 % and 1.9 %); it converges with mpsi. - @test all(abs.(ed.relative_error_effective) .< 0.025) - @test all(abs.(ed.relative_error_standard) .< 0.025) - # The forms agree only after integration, so they differ by the knot-grid quadrature (measured 1.3e-5 and 8.2e-5). + # Measured at most 4e-5 (effective) and 7e-5 (standard); at mpsi = 16 the metric splined between the knots raises it to 3 %. + @test all(abs.(ed.relative_error_effective) .< 5e-4) + @test all(abs.(ed.relative_error_standard) .< 5e-4) + # The forms agree only after integration, so they differ by the ψ quadrature (measured 6e-6 and 3e-5). @test ed.dW_effective ≈ ed.dW_standard rtol = 2e-4 - # (∇×b_eff)·∇ψ is an algebraic identity on the knots (measured 4e-14), a check of the b_eff assembly. + # (∇×b_eff)·∇ψ is an algebraic identity on the knots (measured 7e-14), a check of the b_eff assembly. @test all(ed.curl_residual .< 1e-10) end @@ -80,14 +80,14 @@ using Plots end @testset "metric between the knots" begin - # Off the knots the metric is splined in ψ; it must match the rzphi bicubic to interpolation accuracy (measured 1.2e-5). + # Off the knots the metric is splined in ψ; it must match the rzphi bicubic to interpolation accuracy (measured 1.5e-7). xs = equil.rzphi_xs off = [(xs[end÷2] + xs[end÷2+1]) / 2, psi[end]] g = ED.surface_geometry(ffs, off) ref = [GPEC.Equilibrium.flux_surface_metric(equil, p, t) for p in off, t in ed.theta] - @test maximum(abs, g.jac ./ getfield.(ref, :jac) .- 1) < 1e-4 - @test maximum(abs, sqrt.(g.delpsi2) ./ getfield.(ref, :delpsi) .- 1) < 1e-4 - @test maximum(abs, g.R ./ getfield.(ref, :r) .- 1) < 1e-4 + @test maximum(abs, g.jac ./ getfield.(ref, :jac) .- 1) < 1e-6 + @test maximum(abs, sqrt.(g.delpsi2) ./ getfield.(ref, :delpsi) .- 1) < 1e-6 + @test maximum(abs, g.R ./ getfield.(ref, :r) .- 1) < 1e-6 end @testset "solution state between saved steps" begin @@ -110,10 +110,10 @@ using Plots @testset "ψ integrals start at the ODE start" begin # A qlow start falls between knots, and the sweep saves its first state one ODE step after it. - qffs = forward_solve(template; eq_overrides=boozer, ffs_overrides=Dict{Symbol,Any}(:qlow => 2.5)) + qffs = forward_solve(template; eq_overrides=overrides, ffs_overrides=Dict{Symbol,Any}(:qlow => 2.5)) qed = ED.decompose_plasma_energy(qffs; eigenmodes=[1, 2]) @test qed.psi[1] == qffs.psilow < qffs.solution.psi_store[1] - @test all(abs.(qed.relative_error_standard) .< 0.025) + @test all(abs.(qed.relative_error_standard) .< 5e-4) # measured 2e-5 end @testset "F′ ≠ 0: the DIII-D-like reconstruction" begin @@ -124,5 +124,7 @@ using Plots # Judged against the largest density (measured 7e-11), since single knots where the mode is weak agree less closely. @test maximum(abs, standard .- lagrangian) < 1e-8 * maximum(abs, lagrangian) @test all(ded.curl_residual .< 1e-10) # measured 4e-12 + # The knots alone miss the structure at the pedestal rational surfaces (0.5 % here); with the saved steps it is 1e-5. + @test all(abs.(ded.relative_error_standard) .< 1e-4) end end