diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 575adf10e1..65071b9a21 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,106 @@ 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 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 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: 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 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 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 %. 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 (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 in its PEST coordinates +(`mpert = 8`): + +| `mpsi` | mode 1 effective / standard | mode 2 effective / standard | +|---|---|---| +| 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/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..358576ca48 --- /dev/null +++ b/src/Analysis/EnergyDecomposition.jl @@ -0,0 +1,351 @@ +""" + 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). 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 + +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 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 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 = 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 + +""" + 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 + +# 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_maxima(c_theta, c_zeta, jac, ft, n) -> (residual, scale) + +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_maxima(c_theta::Vector{ComplexF64}, c_zeta::Vector{ComplexF64}, jac::AbstractVector{Float64}, ft::FourierTransform, n::Int) + m = ft.mlow:(ft.mlow+ft.mpert-1) + 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 + +""" + 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 + - `relative_error_effective`, `relative_error_standard` - `(dW − real(ep[k])) / abs(ep[k])` of each total + - `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 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")) + 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 adds the equilibrium knots to 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)) + 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) + 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 + 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 + 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())) + 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)] + 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 + +""" + 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 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 diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl new file mode 100644 index 0000000000..5df8b66490 --- /dev/null +++ b/test/runtests_energy_decomposition.jl @@ -0,0 +1,130 @@ +using TOML +using Plots + +# 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 + ED = GPEC.Analysis.EnergyDecomposition + + # 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[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)] + 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[:, :, 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 + return standard, lagrangian + end + + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + 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 + mtheta = length(ed.theta) + knots = findall(in(equil.rzphi_xs), psi) + geo = ED.surface_geometry(ffs, psi) + + @testset "standard form equals the DCON Lagrangian on the knots" begin + 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 + # 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 7e-14), a check of the b_eff assembly. + @test all(ed.curl_residual .< 1e-10) + 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.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-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 + 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 + + @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=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) .< 5e-4) # measured 2e-5 + 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 + # 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 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",