diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 858bb53a0..95cf10315 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -138,6 +138,35 @@ terms respectively, and do not modify the stored kinetic profile splines. when computing the `toroidal_rotation_factor` back-solve. Julia uses a clean reimplementation with consistent pre-scaling derivatives throughout. +## Regularization: why kinetic runs set `reg_spot = 0` + +`[PerturbedEquilibrium] reg_spot` smooths the displacements before they drive the NTV +integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by +``Q^2/(Q^2 + \mathrm{reg\_spot}^2)`` with ``Q = m - nq``. It exists because **ideal** MHD is +singular at the rationals: ``\bar F_\mathrm{ideal} = Q F Q`` has ``\det \bar F = 0`` there, so +those two components diverge as ``1/Q`` and the torque integral does not converge. + +The **self-consistent kinetic** workflow has no such singularity. Park & Logan +([Phys. Plasmas 24, 032505 (2017)](https://doi.org/10.1063/1.4978562), §III D) decompose the +kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` where +``R_1 \neq 0`` at ``Q = 0``; with finite torque ``\det \bar F`` is complex, its zeros leave the +real ``\psi`` axis, and the singularity is removed from both the solution and the torque +integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to +``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case +for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmetries* (PhD thesis, +Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own +notation. + +GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override. +Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` +is never regularized, so damping the other two breaks their near-resonance cancellation in +``\delta B/B`` and leaves a spurious residue driving the NTV integrand. + +The check that this is right is that the NTV torque and the Euler–Lagrange solution's own +dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` — two independent calculations +of the same quantity — agree to a fraction of a percent with the regularization off, and disagree +by tens of percent with it on. + ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout The method level of `KineticForces//` reports the two physical scalars a user diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index a00db6260..757727d3f 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -94,5 +94,6 @@ f0type = "maxwellian" # Equilibrium distribution moment = "pressure" # Pressure-moment NTV torque atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = true # Enable verbose logging diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 39483c808..0e45a65fd 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -87,3 +87,4 @@ f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, [KineticForces] kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index 47e2e603b..7eabffae0 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -59,3 +59,4 @@ nutype = "harmonic" # Collision operator (zero, small, krook, harmoni f0type = "maxwellian" # Distribution function (maxwellian, jkp, cgl) atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 3ade63787..fa1be99c7 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -228,7 +228,51 @@ function direct_position!(raw_profile::DirectRunInput) end """ - direct_fieldline_int(psifac, raw_profile, ro, zo, rs2) + sfl_sample_time(sol, ts, sfl, sfl_index, x) -> (t, lo) + +Solver time `t` at which the normalised straight-fieldline angle reaches `x` ∈ [0, 1], and the index +`lo` of the saved step that starts its bracket. `sfl` is the SFL integral ∫jac·dl/Bp at the saved +times `ts`; it is monotone, so it brackets the root to one step, and Brent finishes on component +`sfl_index` of the dense solution. +""" +function sfl_sample_time(sol, ts::Vector{Float64}, sfl::AbstractVector{Float64}, sfl_index::Int, x::Float64) + n = length(ts) + x <= 0 && return ts[1], 1 + x >= 1 && return ts[n], n - 1 + target = x * sfl[n] + hi = clamp(searchsortedfirst(sfl, target), 2, n) + lo = hi - 1 + flo, fhi = sfl[lo] - target, sfl[hi] - target + flo == 0 && return ts[lo], lo + fhi == 0 && return ts[hi], lo + flo < 0 < fhi || error("sfl_sample_time: SFL integral is not monotone across a solver step (x = $x)") + return find_zero(t -> sol(t)[sfl_index] - target, (ts[lo], ts[hi]), Roots.Brent()), lo +end + +""" + sample_trace_at_sfl_angles(sol, y_out, sfl_index, theta_nodes, state_row) -> Matrix{Float64} + +Sample a traced flux surface at the straight-fieldline angles `theta_nodes` (0 to 1). Row `i` is +`state_row(sol(t), t, lo)` at the time `t` where the normalised SFL angle reaches `theta_nodes[i]` +(see `sfl_sample_time`), in the 5-column `y_out` layout. The last row is the saved closed endpoint, +so the surface totals `y_out[end, :]` are exact. Every surface sharing the same abscissae keeps the +resample error from varying between neighbouring surfaces. Every field-line tracer must return its +surface through this function. +""" +function sample_trace_at_sfl_angles(sol, y_out::Matrix{Float64}, sfl_index::Int, theta_nodes, state_row::F) where {F} + ts = sol.t::Vector{Float64} + sfl = view(y_out, :, 5) + samples = Matrix{Float64}(undef, length(theta_nodes), 5) + for i in 1:(length(theta_nodes)-1) + t, lo = sfl_sample_time(sol, ts, sfl, sfl_index, Float64(theta_nodes[i])) + samples[i, :] .= state_row(sol(t), t, lo) + end + samples[end, :] .= @view y_out[end, :] + return samples +end + +""" + direct_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield) Performs the field-line integration for a single flux surface. This is a Julia adaptation of the Fortran `direct_fl_int` subroutine. Note that the array `y_out` is now indexed @@ -240,10 +284,12 @@ from 1:5 rather than 0:4 as in Fortran. - `raw_profile`: `DirectRunInput` object containing splines and parameters. - `ro`, `zo`: Coordinates of the magnetic axis [m]. - `rs2`: R-coordinate of the outboard separatrix [m]. + - `theta_nodes`: straight-fieldline angles (0 to 1) at which to sample the surface. ## Returns: - - `y_out`: A matrix containing the integrated quantities vs. the geometric angle `η`. + - `y_out`: A matrix containing the integrated quantities at each of `theta_nodes` + (see `sample_trace_at_sfl_angles`). - `y_out[:, 1]`: η (geometric poloidal angle) - `y_out[:, 2]`: ∫(dl/Bp) - `y_out[:, 3]`: rfac (radial distance from magnetic axis) @@ -252,7 +298,8 @@ from 1:5 rather than 0:4 as in Fortran. - `bfield`: A `DirectBField` object with values at the integration start point. """ -function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)::Tuple{Matrix{Float64},DirectBField} +function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64, + theta_nodes)::Tuple{Matrix{Float64},DirectBField} # Find the starting point on the flux surface (outboard midplane) psi0_guess = raw_profile.psio * (1.0 - psifac) @@ -277,10 +324,11 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: u0 = zeros(Float64, 4) u0[2] = sqrt((r - ro)^2 + (z - zo)^2) - bfield = DirectBField() + # The ODE right-hand side overwrites its own field buffer, so keep `bfield` at the start point (F, P at ψ₀). + bfield_ode = DirectBField() equil_config = raw_profile.config params = FieldLineDerivParams(ro, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio, - equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield) + equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield_ode) # Use a callback to refine the solution at each step to stay on the flux surface function refine_affect!(integrator) @@ -291,7 +339,8 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false)) prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params) - sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=false) + # Dense output so the surface can be sampled at prescribed SFL angles between solver steps. + sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=true) # A failed solve returns a truncated solution instead of throwing; check both the retcode and # that the field line reached η = 2π, since a callback can end it early and still report Success. @@ -312,7 +361,9 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: end sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})' - return hcat(sol.t::Vector{Float64}, sol_matrix), bfield + y_out = hcat(sol.t::Vector{Float64}, sol_matrix) + # State is [∫dl/Bp, rfac, ∫dl/(R²Bp), ∫jac·dl/Bp] and the integration variable is η itself. + return sample_trace_at_sfl_angles(sol, y_out, 4, theta_nodes, (u, eta, _) -> (eta, u[1], u[2], u[3], u[4])), bfield end """ @@ -485,6 +536,9 @@ robustness. - `raw_profile`: A `DirectRunInput` object containing the initial splines (`psi_in`, `sq_in`) and run parameters (`equil_input`). + - `fieldline_int`: the flux-surface tracer, called as `fieldline_int(psifac, raw_profile, ro, zo, + rs2, theta_nodes) -> (y_out, bfield)` with `y_out` already sampled at `theta_nodes` through + `sample_trace_at_sfl_angles`. ## Returns: @@ -510,20 +564,18 @@ robustness. mpsi = length(psi_nodes) - 1 theta_nodes = range(0.0, 1.0; length=mtheta + 1) + ff_x_nodes = collect(theta_nodes) sq_nodes = zeros!(pool, Float64, mpsi + 1, 4) rzphi_nodes = zeros!(pool, Float64, mpsi + 1, mtheta + 1, 4) ff_val = zeros!(pool, Float64, 4) ff_deriv_val = zeros!(pool, Float64, 4) for ipsi in (mpsi+1):-1:1 # outermost to innermost - y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2) + # Every tracer returns the surface sampled at theta_nodes, so all surfaces share SFL abscissae. + y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2, theta_nodes) checkpoint!(pool, Float64) - # Fit data into temporary straight fieldline poloidal angle splines - ff_x_nodes = acquire!(pool, Float64, size(y_out, 1)) - @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5] - - ff_fs_nodes = acquire!(pool, Float64, size(y_out, 1), 4) + ff_fs_nodes = acquire!(pool, Float64, mtheta + 1, 4) @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2 @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4]) diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl index 32b93a0ac..e88559450 100644 --- a/src/Equilibrium/DirectEquilibriumArcLength.jl +++ b/src/Equilibrium/DirectEquilibriumArcLength.jl @@ -59,10 +59,11 @@ so the solver never restricts step size for them. end """ - arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2) + arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield) Arc-length-parameterized flux surface integration. Drop-in replacement for -`direct_fieldline_int` with identical return format: +`direct_fieldline_int` with identical return format, sampled at the straight-fieldline angles +`theta_nodes` (see `sample_trace_at_sfl_angles`): - `y_out[:, 1]`: geometric angle η ∈ 0 to 2π (CCW from outboard midplane) - `y_out[:, 2]`: accumulated ∫dl/Bp @@ -74,7 +75,7 @@ The ODE is terminated by a `ContinuousCallback` that detects the return to the outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. """ @with_pool pool function arclength_fieldline_int( - psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64 + psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64, theta_nodes )::Tuple{Matrix{Float64},DirectBField} psi0_guess = raw_profile.psio * (1.0 - psifac) @@ -114,13 +115,14 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. reltol_vec = [equil_config.etol, equil_config.etol, 1e20, 1e20, 1e20] abstol_vec = [1e-8, 1e-8, 1e20, 1e20, 1e20] sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec, - dt=2π / 200, adaptive=true, dense=false) + dt=2π / 200, adaptive=true, dense=true) n = length(sol.u) y_out = Matrix{Float64}(undef, n, 5) # Compute geometric angles and unwrap to monotone [0, 2π] - η_raw = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n] + η_atan = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n] + η_raw = copy(η_atan) for i in 2:n Δ = η_raw[i] - η_raw[i-1] if Δ > π @@ -141,7 +143,13 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. y_out[i, 5] = sol.u[i][5] end + # The integration variable is arc length, so η at a sample is unwrapped against its bracketing step. + function state_row(u, _s, lo) + η = y_out[lo, 1] + rem2pi(atan(u[2] - zo, u[1] - ro) - η_atan[lo], RoundNearest) + return (η, u[3], sqrt((u[1] - ro)^2 + (u[2] - zo)^2), u[4], u[5]) + end + # bfield at the starting point carries F and P for the surface-averaged quantities - return y_out, bfield + return sample_trace_at_sfl_angles(sol, y_out, 5, theta_nodes, state_row), bfield end diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..0b64337a7 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -419,6 +419,7 @@ density (`_knot_density`), equidistribution, a global minimum-spacing floor (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; +`pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the @@ -428,6 +429,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], + pinned::Vector{Float64}=Float64[], singfac_min::Float64=1e-4, n_min::Int=1, bracket_coef::Float64=BRACKET_COEF, @@ -449,6 +451,9 @@ function refined_psi_grid(equil::PlasmaEquilibrium; N == N_cap && M_total > N_cap && @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable" grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing) + # Pinned knots use knot-at-node semantics (no cleared zone); inserted before rational bracketing + # so a rational's bracket still wins locally where the two overlap. + grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned) isempty(mandatory) && return grid min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory] return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index f5d188abe..57cc28128 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -311,7 +311,7 @@ function build_kinetic_metric_matrices(equil::Equilibrium.PlasmaEquilibrium, v[3, 3] = 2π * rs # Raw g^ij (Fortran dcon_interface.f:933-937) - g12 = v[1, 1]*v[2, 1] + v[1, 2]*v[2, 2] + v[1, 3]*v[2, 3] + g12 = v[1, 1] * v[2, 1] + v[1, 2] * v[2, 2] + v[1, 3] * v[2, 3] g13 = v[3, 3] * v[1, 3] g22 = v[2, 1]^2 + v[2, 2]^2 + v[2, 3]^2 g23 = v[2, 3] * v[3, 3] @@ -428,6 +428,34 @@ function build_kinetic_metric_matrices(equil::Equilibrium.PlasmaEquilibrium, end +""" + core_capped_knots(xs, rationals) -> Vector{Int} + +Indices of `xs` the EL coefficient splines keep, decoupling their knots from the equilibrium grid +in the packed core. A cubic spline's third-derivative jumps scale as (node error)/Δψ³, so +equilibrium-grade core packing amplifies tolerance-level node error into jumps that slave the +Euler-Lagrange step size. Near the axis every component is a Frobenius power law in ψ, and power +laws are scale-free, so log-uniform sampling (Δψ ≥ 0.05·ψ) resolves them at constant relative +accuracy — cubic interpolation of ψ^p on that grid errs by ~(0.05p)⁴/384, resolving even the +steepest component (p = m_max/2) to ~2e-4, well below where the physics responds. The capped region +ends at the innermost rational (or ψ = 0.1, whichever is smaller) and never drops a knot inside a +rational's resolution window, preserving the Δ′-stencil structure the equilibrium grid encodes +(`Equilibrium.RATIONAL_RES_RADIUS`). +""" +function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int} + length(xs) < 3 && return collect(eachindex(xs)) + cap_edge = 0.1 + isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS)) + in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals) + keep = Int[1] + for i in 2:(length(xs)-1) + x = xs[i] + (x >= cap_edge || in_rational_window(x) || (x - xs[keep[end]]) >= 0.05 * x) && push!(keep, i) + end + push!(keep, length(xs)) + return keep +end + """ build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) -> MatrixSplines @@ -607,19 +635,23 @@ function build_matrix_splines(equil::Equilibrium.PlasmaEquilibrium, intr::ForceF # Create complex series interpolants with per-column extrap BC # TODO: set powers. Do we need this yet? Only called if power_flag = true itp_opts = (; extrap=ExtendExtrap()) + keep = core_capped_knots(metric.xs, [s.psifac for s in intr.sing]) + mxs = metric.xs[keep] + length(mxs) < length(metric.xs) && + @info "EL coefficient-spline grid: $(length(metric.xs)) -> $(length(mxs)) knots after core density cap" ideal = IdealMatrices(; - A_spline=cubic_interp(metric.xs, Series(amats_flat); itp_opts...), - B_spline=cubic_interp(metric.xs, Series(bmats_flat); itp_opts...), - C_spline=cubic_interp(metric.xs, Series(cmats_flat); itp_opts...), - D_spline_prim=cubic_interp(metric.xs, Series(dmats_flat); itp_opts...), - E_spline_prim=cubic_interp(metric.xs, Series(emats_flat); itp_opts...), - H_spline=cubic_interp(metric.xs, Series(hmats_flat); itp_opts...), - F_spline_lower=cubic_interp(metric.xs, Series(fmats_lower_flat); itp_opts...), - F_spline_prim=cubic_interp(metric.xs, Series(fmats_prim_flat); itp_opts...), - F_spline_gal=cubic_interp(metric.xs, Series(fmats_gal_flat); itp_opts...), - G_spline=cubic_interp(metric.xs, Series(gmats_flat); itp_opts...), - K_spline=cubic_interp(metric.xs, Series(kmats_flat); itp_opts...), + A_spline=cubic_interp(mxs, Series(amats_flat[keep, :]); itp_opts...), + B_spline=cubic_interp(mxs, Series(bmats_flat[keep, :]); itp_opts...), + C_spline=cubic_interp(mxs, Series(cmats_flat[keep, :]); itp_opts...), + D_spline_prim=cubic_interp(mxs, Series(dmats_flat[keep, :]); itp_opts...), + E_spline_prim=cubic_interp(mxs, Series(emats_flat[keep, :]); itp_opts...), + H_spline=cubic_interp(mxs, Series(hmats_flat[keep, :]); itp_opts...), + F_spline_lower=cubic_interp(mxs, Series(fmats_lower_flat[keep, :]); itp_opts...), + F_spline_prim=cubic_interp(mxs, Series(fmats_prim_flat[keep, :]); itp_opts...), + F_spline_gal=cubic_interp(mxs, Series(fmats_gal_flat[keep, :]); itp_opts...), + G_spline=cubic_interp(mxs, Series(gmats_flat[keep, :]); itp_opts...), + K_spline=cubic_interp(mxs, Series(kmats_flat[keep, :]); itp_opts...), # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl - J_spline=cubic_interp(metric.xs, Series(jmats_flat); itp_opts...)) + J_spline=cubic_interp(mxs, Series(jmats_flat[keep, :]); itp_opts...)) return MatrixSplines(; ideal) end diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 0cd8e01dd..d148b12e9 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,3 +1,73 @@ + +# Knots across the envelope's [ψ_c, 2ψ_c] quintic-smoothstep transition, enough for a cubic spline +# to follow it without overshoot even on coarse decks. +const BAND_KNOTS = 9 + +""" + refine_grid_at_fbar_peaks(xs, kw, kt, evaluate, kin, equil, intr, psi_c; + ngrid=1000, relaxed_frac=0.01, target=3, max_add=24) → (xs, kw, kt) + +Insert kinetic evaluation knots across near-singular structure of F̄ that the grid does not +resolve. Scans cond(F̄) (the same operator `find_kinetic_singular_surfaces!` searches — Park & +Logan Eq. 70, so shifted and split resonances are included), takes peaks between +`relaxed_frac`·threshold and the singular threshold, measures each peak's FWHM, and adds knots +only where fewer than `target` knots lie inside it. New points respect `MIN_KNOT_SPACING`, stay +above the near-axis validity band, and are capped at `max_add`; each costs one kernel evaluation +and the existing values are reused. A well-resolved grid inserts nothing. +""" +function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3}, + evaluate::Function, kin::KineticMatrices, equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal, psi_c::Float64; + ngrid::Int=1000, relaxed_frac::Float64=KINETIC_RELAXED_FRAC, target::Int=3, max_add::Int=24, + cond_threshold::Float64=KINETIC_SINGULAR_COND) + + lo, hi = xs[1], xs[end] + scan = collect(range(lo, hi; length=ngrid)) + hint = Ref(1) + cond_vals = [ + try + evaluate_fbar_condition(x, kin, equil, intr; hint=hint) + catch + Inf + end for x in scan + ] + + add = Float64[] + for i in 2:(ngrid-1) + c = cond_vals[i] + (c > cond_vals[i-1] && c > cond_vals[i+1] && relaxed_frac * cond_threshold < c <= cond_threshold) || continue + l = i + while l > 1 && cond_vals[l] > c / 2 + l -= 1 + end + r = i + while r < ngrid && cond_vals[r] > c / 2 + r += 1 + end + inside = count(x -> scan[l] <= x <= scan[r], xs) + inside >= target && continue + w = (scan[r] - scan[l]) / 3 + for x in (scan[i], scan[i] - w, scan[i] + w) + (lo < x < hi && x > 2 * psi_c) || continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, xs) && continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, add) && continue + push!(add, x) + end + end + isempty(add) && return xs, kw, kt + length(add) > max_add && (add = sort(add)[1:max_add]) + + sort!(add) + @info "Kinetic grid: $(length(add)) knot(s) added across unresolved near-singular F̄ structure at " * + "ψ=$(round.(add; digits=4)) (cond peaks below the singular threshold)" + kw_new, kt_new = evaluate(add) + allxs = vcat(xs, add) + perm = sortperm(allxs) + kw_all = cat(kw, kw_new; dims=1)[perm, :, :] + kt_all = cat(kt, kt_new; dims=1)[perm, :, :] + return allxs[perm], kw_all, kt_all +end + """ build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; calculated_source=nothing) @@ -27,11 +97,17 @@ function build_kinetic_matrix_splines( mats::MatrixSplines, intr::ForceFreeStatesInternal, metric::MetricData; - calculated_source::Union{Nothing,Function}=nothing + calculated_source::Union{Nothing,Function}=nothing, + axis_validity_psi_c::Float64=0.0 ) xs = metric.xs mpsi = length(xs) + # Coarse equilibrium grids cannot represent env·(increment) across the validity band, and the + # spline overshoot can land on a rational; pin the band ends (C² there) and resolve the middle. + band_knots(lo, hi) = axis_validity_psi_c > 0 ? + [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=BAND_KNOTS) if lo < x < hi] : Float64[] + # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs) @@ -42,7 +118,14 @@ function build_kinetic_matrix_splines( "calling build_kinetic_matrix_splines directly, or pass " * "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly." ) - kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) + band = band_knots(xs[1], xs[end]) + if isempty(band) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) + else + xs = sort!(unique!(vcat(collect(xs), band))) + mpsi = length(xs) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats; psis=xs) + end kw_flat .*= ctrl.kinetic_factor kt_flat .*= ctrl.kinetic_factor else @@ -50,7 +133,24 @@ function build_kinetic_matrix_splines( end # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat) + mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs) + + # The FKG splines now exist, so F̄ can be scanned: add knots only where near-singular structure + # (shifted/split kinetic resonances) falls in an interval that does not resolve it. + if ctrl.kinetic_source == "calculated" && calculated_source !== nothing && mats.kinetic !== nothing + xs2, kw_flat, kt_flat = refine_grid_at_fbar_peaks( + collect(xs), kw_flat, kt_flat, + psis -> begin + kwn, ktn = calculated_source(ctrl, equil, intr, metric, mats; psis=psis) + (kwn .* ctrl.kinetic_factor, ktn .* ctrl.kinetic_factor) + end, + mats.kinetic, equil, intr, axis_validity_psi_c) + if length(xs2) != length(xs) + xs = xs2 + mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs) + end + end + return mats end """ @@ -73,9 +173,9 @@ function _compute_fkg_matrices( intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, - kt_flat::Array{ComplexF64,3} + kt_flat::Array{ComplexF64,3}; + xs::Vector{Float64}=metric.xs ) - xs = metric.xs mpsi = length(xs) np = intr.numpert_total mpert = intr.mpert diff --git a/src/ForceFreeStates/Surfaces/Asymptotics.jl b/src/ForceFreeStates/Surfaces/Asymptotics.jl index 26746c0ca..87cac444a 100644 --- a/src/ForceFreeStates/Surfaces/Asymptotics.jl +++ b/src/ForceFreeStates/Surfaces/Asymptotics.jl @@ -74,7 +74,7 @@ function compute_sing_asymptotics( # This is the parameter α but for all modes - α = 0 for non-resonant modes power[ipert_res] .= -alpha - power[ipert_res .+ intr.numpert_total] .= alpha + power[ipert_res.+intr.numpert_total] .= alpha # Zeroth-order non-resonant solutions for ipert in 1:intr.numpert_total @@ -111,7 +111,7 @@ function compute_sing_asymptotics( msg *= @sprintf(" m0mat(1,2)= %+.12e %+.12ei\n", real(m0mat[1, 2]), imag(m0mat[1, 2])) msg *= @sprintf(" m0mat(2,1)= %+.12e %+.12ei\n", real(m0mat[2, 1]), imag(m0mat[2, 1])) msg *= @sprintf(" m0mat(2,2)= %+.12e %+.12ei\n", real(m0mat[2, 2]), imag(m0mat[2, 2])) - di = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2] + di = m0mat[1, 1] * m0mat[2, 2] - m0mat[2, 1] * m0mat[1, 2] msg *= @sprintf(" di= %+.12e, alpha= %+.12e %+.12ei\n", real(di), real(alpha[1]), imag(alpha[1])) msg *= @sprintf(" psifac= %+.12e, r1=%d, ipert0=%d\n", singp.psifac, r1[1], ipert0) msg *= @sprintf(" vmat(ip,ip,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0, 2, 1]), imag(vmat[ipert0, ipert0, 2, 1])) @@ -579,7 +579,7 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64) # Restore powers (unshear v→u) — matches Fortran STRIDE sing_get_ua for i in eachindex(r1) - pfac = pfac_base ^ sing_asymp.alpha[i] # dpsi^α + pfac = pfac_base^sing_asymp.alpha[i] # dpsi^α ua[:, r2[2*i-1], :] ./= pfac # big solution column: /dpsi^α ua[:, r2[2*i], :] .*= pfac # small solution column: *dpsi^α ua[r1[i], :, 1] ./= sqrtfac # resonant row ξ: /√dpsi diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 18c45def7..9a5ec118e 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -125,9 +125,11 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends # on which n is used — and fall back to qhigh / psihigh truncation with a warning. if ctrl.set_psilim_via_dmlim && intr.nlow <= 0 - error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * - "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * - "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.") + error( + "sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * + "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * + "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead." + ) elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation." elseif ctrl.set_psilim_via_dmlim @@ -256,6 +258,11 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi return cond(fbar) end +# Kinetic F̄ counts as singular above this condition number; structure within +# KINETIC_RELAXED_FRAC of it is reported as near-singular (the shifted/split resonances). +const KINETIC_SINGULAR_COND = 1.0e8 +const KINETIC_RELAXED_FRAC = 0.01 + """ find_kinetic_singular_surfaces!(mats, equil, intr; ngrid=2000, cond_threshold=1e8) @@ -267,14 +274,16 @@ Mirrors the intent of Fortran `ksing_find` (`sing.f:1486-1616`) which finds zero det(F̄) via adaptive bisection. Here we use condition number peaks instead of determinant zeros for better numerical robustness and scale invariance. -Algorithm: - - 1. Evaluate cond(F̄) on a dense ψ grid - 2. Find local maxima (peaks where gradient changes from + to -) - 3. Refine each peak with golden-section minimization of -cond - 4. Filter by threshold and resonance condition +Evaluates cond(F̄) on a dense ψ grid, finds its local maxima, refines each with golden-section +minimization of -cond, then filters by threshold and resonance condition. """ -function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) +function find_kinetic_singular_surfaces!( + mats::MatrixSplines, + equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal; + ngrid::Int=2000, + cond_threshold::Float64=KINETIC_SINGULAR_COND +) kin = mats.kinetic kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call build_kinetic_matrix_splines first") psilow = equil.profiles.xs[1] @@ -299,11 +308,17 @@ function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium intr.kinsing_scan_threshold = cond_threshold # Find local maxima of cond(F̄): points where cond increases then decreases - peak_indices = Int[] - for i in 2:(ngrid-1) - if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && cond_vals[i] > cond_threshold - push!(peak_indices, i) - end + local_maxima = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1]] + peak_indices = filter(i -> cond_vals[i] > cond_threshold, local_maxima) + + # Peaks below threshold mark near-singular structure (the shifted/split resonances of Park & + # Logan Eq. 70); report the strongest few so sharp kinetic structure is visible rather than silent. + subthreshold = filter(i -> KINETIC_RELAXED_FRAC * cond_threshold < cond_vals[i] <= cond_threshold, local_maxima) + if !isempty(subthreshold) + top = sort(subthreshold; by=i -> cond_vals[i], rev=true)[1:min(3, length(subthreshold))] + @info "Kinetic F̄ near-singular structure below the singular threshold at " * + join(["ψ=$(round(psi_grid[i]; digits=4)) (cond=$(round(cond_vals[i]; sigdigits=3)))" for i in top], ", ") * + " — full scan in SingularSurfaces/Kinetic/scan_cond; check the ψ grid resolves these if results look grid-sensitive" end # Refine each peak to find the precise ψ location diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 27e29ed54..fade8405f 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -220,7 +220,8 @@ function main_from_inputs( equil = Equilibrium.setup_equilibrium(eq_config, additional_input) kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) - equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) + equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; + kf_ctrl=kf_ctrl, kf_species=kf_species) equil_dt = time() - equil_start push!(runtimes, "equilibrium" => equil_dt) @@ -269,6 +270,15 @@ function main_from_inputs( @info "Results written to $(ctrl.HDF5_filename)" end + # Validity describes the kinetic model, not the torque contraction, so it is written wherever + # kinetic profiles were loaded — including FFS-only runs that never reach run_kinetic_forces. + if ctrl.write_outputs_to_HDF5 && kinetic_profiles !== nothing + h5open(joinpath(ffs_result.dir_path, ctrl.HDF5_filename), "cw") do h5file + KineticForces.write_validity!(h5file, kf_ctrl, kf_species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) + end + end + ffs_dt = time() - ffs_start push!(runtimes, "force_free_states" => ffs_dt) @info "Force-Free States completed in $(@sprintf("%.3f", ffs_dt)) s" @@ -441,7 +451,8 @@ function load_kinetic_context( end """ - maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) -> equil + maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; + kf_ctrl=nothing, kf_species=nothing) -> equil Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory @@ -453,16 +464,24 @@ function maybe_reform_equilibrium( additional_input, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, - kinetic_profiles + kinetic_profiles; + kf_ctrl::Union{Nothing,KineticForces.KineticForcesControl}=nothing, + kf_species=nothing ) Equilibrium.wants_two_pass(eq_config) || return equil mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) + # Kinetic-resonance nodes are located here rather than at the call site so the validity + # boundary that filters them is cleared against this pass's rationals (intr.sing is empty yet). + pinned = kf_ctrl === nothing ? Float64[] : + KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr; rationals=mandatory) # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), # so the rational-surface brackets clear a zone large enough for every mode. n_min = minimum(abs(n) for n in intr.nlow:intr.nhigh if n != 0) + isempty(pinned) || + @info "Pinning $(length(pinned)) kinetic-resonance surfaces into the ψ grid: $(round.(sort(pinned); digits=3))" psi_nodes = Equilibrium.refined_psi_grid(equil; - tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, + tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned, singfac_min=ctrl.singfac_min, n_min=n_min) rerun_input = if additional_input !== nothing # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS @@ -628,12 +647,18 @@ function prepare_force_free_states!( # Inject the KineticForces callback so the "calculated" source can # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). - calculated_cb = (c, e, i, m, f) -> - KineticForces.compute_calculated_kinetic_matrices( - c, e, i, m, f; - kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species) + # One boundary for the whole run: the widest-orbit species, moved clear of any rational + # surface whose window the transition band would otherwise cut through. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) + calculated_cb = + (c, e, i, m, f; psis::Vector{Float64}=Float64[]) -> + KineticForces.compute_calculated_kinetic_matrices( + c, e, i, m, f; + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species, psis=psis, + axis_psi_c=axis_psi_c) mats = build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; - calculated_source=calculated_cb) + calculated_source=calculated_cb, axis_validity_psi_c=axis_psi_c) # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings. # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings; @@ -915,7 +940,7 @@ function perturbed_equilibrium( kwargs... ) ctrl = ffs.control - pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; kwargs...) + pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; PerturbedEquilibrium.kinetic_regularization_kwargs(ffs, kwargs)...) pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=ffs.dir_path) # Inner-layer penetrated resonant field; zeros under ideal closure. @@ -976,6 +1001,10 @@ function run_kinetic_forces( @info "\n KineticForces\n$_SECTION" kf_start = time() + # One boundary for the whole run, shared with the EL kinetic matrices and the Validity output. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) + # Standalone NTV torque diagnostics need a PE state (they contract kinetic operators # against ξ). The self-consistent kinetic_source="calculated" path produces none — skip. if pe_state === nothing @@ -987,7 +1016,8 @@ function run_kinetic_forces( if species === nothing kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles) + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles; + axis_psi_c=axis_psi_c) if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; @@ -1004,7 +1034,8 @@ function run_kinetic_forces( (f => getfield(kf_ctrl, f) for f in fieldnames(KineticForces.KineticForcesControl))..., zi=sp.z, mi=sp.m, electron=sp.electron, ion_species=KineticForces.IonSpecies[]) st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles) + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles; + axis_psi_c=axis_psi_c) push!(states, st) push!(labels, sp.label) end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index bb61135e4..d6436b710 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -37,6 +37,7 @@ section. - `ffs_intr`: ForceFreeStatesInternal (mode indexing) - `metric`: MetricData (provides ψ grid via `metric.xs`) - `mats`: MatrixSplines (used only for `numpert_total` cross-check) + - `psis`: ψ grid for the kernel; empty (default) evaluates on `metric.xs` # Keyword arguments @@ -59,9 +60,13 @@ function compute_calculated_kinetic_matrices( mats; kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, - species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing + species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing, + psis::Vector{Float64}=Float64[], + axis_psi_c::Float64=0.0 ) - xs = metric.xs + # The kernel is a pure function of psi (it evaluates equilibrium splines), so it can be + # driven over any knot list; default is the full equilibrium grid. + xs = isempty(psis) ? metric.xs : psis mpsi = length(xs) mpert = ffs_intr.mpert npert = ffs_intr.npert @@ -110,6 +115,20 @@ function compute_calculated_kinetic_matrices( splist = species === nothing ? [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species + + # Near-axis validity envelope: suppresses the drift-kinetic increments where the zero-orbit-width + # ordering fails (kernel evaluation skipped where it is 0). The same real scalar must scale kw and + # kt at the same ψ — that is what keeps the blocks Hermitian and the energy/torque relation intact. + env = ones(Float64, mpsi) + if kf_ctrl.axis_validity_suppression + psi_c = axis_psi_c + if psi_c > 0 + env .= kinetic_axis_validity_envelope.(xs, psi_c) + @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " * + "psi=$(round(2 * psi_c; sigdigits=3)) (kernel evaluation skipped below psi_c)" maxlog = 1 + end + end + for sp in splist # Hoisted per-species scalars: the closure then captures concrete values, not a # union-typed `sp` (the ternary above mixes the resolved vector with the fallback). @@ -124,6 +143,7 @@ function compute_calculated_kinetic_matrices( block_w = thread_block_w[tid] block_t = thread_block_t[tid] psi = xs[ipsi] + env[ipsi] == 0.0 && continue # inside the near-axis validity band: kernel skipped for in_idx in 1:npert n = ffs_intr.nlow + in_idx - 1 fill!(full_w, 0) @@ -148,8 +168,8 @@ function compute_calculated_kinetic_matrices( row_offset = (in_idx - 1) * mpert for k in 1:6, j in 1:mpert, i in 1:mpert idx = (row_offset + j - 1) * np + (row_offset + i) - kw_flat[ipsi, idx, k] += full_w[i, j, k] - kt_flat[ipsi, idx, k] += full_t[i, j, k] + kw_flat[ipsi, idx, k] += env[ipsi] * full_w[i, j, k] + kt_flat[ipsi, idx, k] += env[ipsi] * full_t[i, j, k] end end end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 62c55d7f1..3ed4fdd9d 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -83,7 +83,7 @@ function integrate_psi_quadgk( wdfac::Float64, divxfac::Float64, electron::Bool, method::String, equil, intr::KineticForcesInternal, ctrl::KineticForcesControl, kinetic_profiles::Equilibrium.KineticProfileSplines; - psi_min::Float64=0.0, psi_max::Float64=1.0 + psi_min::Float64=0.0, psi_max::Float64=1.0, axis_psi_c::Float64=0.0 ) is_matrix_method = occursin("mm", method) mpert = intr.mpert @@ -100,6 +100,17 @@ function integrate_psi_quadgk( panel_psis=Float64[], resonance_psis=Float64[]) end + # Near-axis validity suppression on the torque density. `axis_psi_c` is the run's single + # boundary (widest-orbit species, cleared of rationals), so the quadrature domain matches the + # EL kinetic matrices and the Validity output; below it the integrand is identically zero. + psi_c = ctrl.axis_validity_suppression ? axis_psi_c : 0.0 + x0 = max(x0, psi_c) + if x0 >= xout + return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, + panel_psis=Float64[], resonance_psis=Float64[]) + end + psi_c > 0 && @info "Kinetic axis-validity suppression in ψ torque quadrature: domain starts at psi_c=$(round(psi_c; sigdigits=3))" + # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine # loop invokes the callback with small batches (~15 nodes), so threading it is # fork-join-bound. Instead thread the inner bounce-harmonic loop (2·nl+1 harmonics), @@ -149,6 +160,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm total += harm_vals[ell_idx] end + total *= kinetic_axis_validity_envelope(psi, psi_c) y[k] = total push!(logged_psi, psi) @@ -158,6 +170,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm elems_accum .+= harm_elems[ell_idx] end + elems_accum .*= kinetic_axis_validity_envelope(psi, psi_c) push!(logged_elems, elems_accum) end end @@ -166,7 +179,7 @@ function integrate_psi_quadgk( # Panels at the rational surfaces the run resolved plus the kinetic-resonance # surfaces (thermal-energy Ω_ℓ = 0 for ℓ ∈ -nl:nl) — both are torque-density peaks. resonance_psis = kinetic_resonance_psi_nodes(kinetic_profiles, equil; n, nl, zi, mi, electron, wdfac) - pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis), x0, xout) + pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis, [2 * psi_c]), x0, xout) bi = QuadGK.BatchIntegrand(psi_batch!, ComplexF64[], Float64[]) total, quad_err = quadgk(bi, pts...; atol=ctrl.atol_psi, rtol=ctrl.rtol_psi, maxevals=ctrl.maxevals_psi) @@ -291,10 +304,13 @@ block-diagonal kinetic matrices. - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run - `equil`: PlasmaEquilibrium with 2D interpolants - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + - `axis_psi_c`: the run's near-axis validity boundary, suppressing the torque density below it + (`KineticForces.axis_validity_boundary`); 0 disables the suppression """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + kinetic_profiles::Equilibrium.KineticProfileSplines; + axis_psi_c::Float64=0.0) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -337,7 +353,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor n, ctrl.nl, ctrl.zi, ctrl.mi, ctrl.wdfac, ctrl.divxfac, ctrl.electron, method, equil, intr, ctrl, kinetic_profiles; - psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2]) + psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2], axis_psi_c=axis_psi_c) total_torque += result.total psi_nsteps_total += result.psi_nsteps diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index ed9ef14a6..30725ed9d 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -3,12 +3,13 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTuple `(name, flag, kind, doc)`: -- `name` — short method identifier used as the HDF5 group key and in `intr.method` -- `flag` — the `KineticForcesControl` field symbol that enables the method -- `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` - (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the - three special-cased methods) -- `doc` — one-line description printed in verbose output + + - `name` — short method identifier used as the HDF5 group key and in `intr.method` + - `flag` — the `KineticForcesControl` field symbol that enables the method + - `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` + (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the + three special-cased methods) + - `doc` — one-line description printed in verbose output The method names/docs and the `Compute.jl` enable list are all derived from this tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. @@ -77,6 +78,7 @@ User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -84,6 +86,13 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce Immutable: vary a field by building a new control rather than assigning to one (the multi-species loop does this per species, and `check_psi_quadrature_convergence`'s test builds a second control for its differing tolerance). + + - `axis_validity_suppression::Bool` - Suppress the calculated kinetic terms where the + zero-orbit-width ordering fails near the axis: below ψ_c — the outermost ψ at which a thermal + orbit width of the widest-orbit species reaches ⟨r⟩, moved clear of any rational surface the + transition would otherwise cut through — the increments are zeroed and the kernel is skipped, + rising to full strength at 2ψ_c through a C² envelope. Boundary and envelope come from the + equilibrium and kinetic profiles; there are no tuning parameters. Default `true`. """ @kwdef struct KineticForcesControl # Moment type @@ -121,32 +130,19 @@ builds a second control for its differing tolerance). nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy). - # Each level must be resolved more tightly than the one enclosing it, or the outer - # integrator chases its integrand's own quadrature noise instead of converging. - # *_xlmda: tolerances for the λ (pitch) integration - # *_x: tolerances for the x (energy) integration nested inside it; NaN ⇒ derive as - # nested_tolerance_margin × the pitch tolerances - # *_psi: tolerances for the outer ψ quadrature + # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy); each level must + # resolve tighter than the one enclosing it, or the outer integrator chases quadrature noise. atol_xlmda::Float64 = 1e-8 # Absolute tolerance for the inner pitch integration rtol_xlmda::Float64 = 1e-5 # Relative tolerance for the inner pitch integration atol_x::Float64 = NaN # Absolute tolerance for the energy integration (NaN ⇒ derived) rtol_x::Float64 = NaN # Relative tolerance for the energy integration (NaN ⇒ derived) - # The pitch integrand IS the energy integral, so the energy level is resolved this much - # tighter than the pitch level by default. + # Energy tolerances default to this fraction of the pitch ones (the pitch integrand is itself an energy integral). nested_tolerance_margin::Float64 = 1e-2 # Factor relating derived energy tolerances to the pitch ones - # rtol_psi is the primary convergence knob: ~2 significant figures matches the validity - # of the NTV model approximations. Do not set it tighter than the noise floor of the - # inner integrals (keep rtol_psi ≳ 10 × rtol_xlmda). + # Primary convergence knob (~2 sig figs matches NTV model validity); keep rtol_psi ≳ 10 × rtol_xlmda or it chases inner-integral noise. rtol_psi::Float64 = 1e-2 # Relative tolerance for outer ψ quadrature - # atol_psi is in N·m and therefore amplitude-sensitive: NTV scales as δB², so a 10× - # weaker applied field gives a 100× smaller torque and any fixed absolute tolerance can - # silently dominate termination with O(1) relative error. Default 0 (rtol-only); a - # nonzero value is an expert opt-out for near-net-zero-torque cases and triggers a - # warning when it dominates. + # Amplitude-sensitive (NTV ∝ δB²); default 0 (rtol-only). Nonzero is an expert opt-out for near-zero-torque cases and warns when it dominates. atol_psi::Float64 = 0.0 # Absolute tolerance for outer ψ quadrature [N·m] - # Runaway guard for the outer quadrature (e.g. sign-cancelling torque density with tiny - # net torque under rtol-only control); a convergence warning fires when hit. + # Runaway guard for rtol-only control (e.g. sign-cancelling torque density); warns when hit. maxevals_psi::Int = 2000 # Max integrand evaluations for outer ψ quadrature # Scaling factors @@ -158,6 +154,7 @@ builds a second control for its differing tolerance). nufac::Float64 = 1.0 # Collisionality scaling divxfac::Float64 = 1.0 # div(xi_perp) scaling + axis_validity_suppression::Bool = true # documented in the `## Fields` docstring # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" @@ -211,12 +208,13 @@ Internal working state for KineticForces calculations. Holds equilibrium-derived quantities, profile interpolants, and integration results. Fields replacing former module-level globals: -- `ro`, `bo`, `chi1`: Equilibrium geometry parameters -- `mthsurf`, `mfac`: Poloidal grid info -- `dbob_m`, `divx_m`: Perturbation mode interpolants -- `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as - panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on - Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection + + - `ro`, `bo`, `chi1`: Equilibrium geometry parameters + - `mthsurf`, `mfac`: Poloidal grid info + - `dbob_m`, `divx_m`: Perturbation mode interpolants + - `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as + panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on + Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection Equilibrium and kinetic profile data are read directly from the `PlasmaEquilibrium` (`equil.profiles`, `equil.geometry`) and the @@ -306,17 +304,17 @@ function KineticForcesInternal(equil; verbose::Bool=false) # Axis toroidal field F(0)/ro that normalizes λ = μ·bo/E; F_spline stores 2πF. bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro) KineticForcesInternal(; - ro = equil.ro, - bo = bo_axis, - chi1 = 2π * equil.psio, + ro=equil.ro, + bo=bo_axis, + chi1=2π * equil.psio, mthsurf, - tpsi_xs = collect(range(0.0, 1.0, length=nth)), - tpsi_B = Vector{Float64}(undef, nth), - tpsi_dBdpsi = Vector{Float64}(undef, nth), - tpsi_dBdtheta = Vector{Float64}(undef, nth), - tpsi_jac = Vector{Float64}(undef, nth), - tpsi_djdpsi = Vector{Float64}(undef, nth), - verbose, + tpsi_xs=collect(range(0.0, 1.0; length=nth)), + tpsi_B=Vector{Float64}(undef, nth), + tpsi_dBdpsi=Vector{Float64}(undef, nth), + tpsi_dBdtheta=Vector{Float64}(undef, nth), + tpsi_jac=Vector{Float64}(undef, nth), + tpsi_djdpsi=Vector{Float64}(undef, nth), + verbose ) end @@ -326,20 +324,22 @@ end Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal. Builds three interpolant sets from PE Clebsch displacements: -1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ -2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) -3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting + + 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ + 2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) + 3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): -1. Apply geometric matrices S,T,X,Y,Z in m-space -2. Inverse DFT to θ-space -3. Divide by J·B² at each θ -4. Forward DFT back to m-space + + 1. Apply geometric matrices S,T,X,Y,Z in m-space + 2. Inverse DFT to θ-space + 3. Divide by J·B² at each θ + 4. Forward DFT back to m-space """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::PerturbedEquilibrium.PerturbedEquilibriumState, - ffs::ForceFreeStates.ForceFreeStatesResult, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + ffs::ForceFreeStates.ForceFreeStatesResult, + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) ffs.npert == 1 || error("KineticForces supports a single toroidal mode; this run has n = $(ffs.nlow):$(ffs.nhigh).") # Copy mode numbers from FFS @@ -410,9 +410,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur psi = psi_grid[ipsi] # Get Clebsch displacement vectors at this ψ - xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] + xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] + xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert geom_mats.smats(smat_flat, psi; hint=hint_s) @@ -432,8 +432,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur mul!(jbb_kapx, smat, xsp) mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += tmat * xms mul!(jbb_divx, xmat, xmp1) - mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp - mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms + mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp + mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms @. jbb_dbob = -(jbb_divx + jbb_kapx) # Inverse DFT to θ-space, divide by J·B², forward DFT back @@ -460,9 +460,9 @@ Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back. """ function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, - ft::Utilities.FourierTransforms.FourierTransform, - psi::Float64, equil::Equilibrium.PlasmaEquilibrium, - mthsurf::Int, theta_buf::Vector{ComplexF64}) + ft::Utilities.FourierTransforms.FourierTransform, + psi::Float64, equil::Equilibrium.PlasmaEquilibrium, + mthsurf::Int, theta_buf::Vector{ComplexF64}) # Inverse DFT: m-space → θ-space theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) @@ -532,8 +532,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "KineticForces" group. """ @kwdef mutable struct KineticForcesState - method_results::Dict{String, MethodResult} = Dict{String, MethodResult}() + method_results::Dict{String,MethodResult} = Dict{String,MethodResult}() # Block-diagonal kinetic matrices: key=method, value=(numpert_total, numpert_total, 6) - kinetic_matrices::Dict{String, Array{ComplexF64,3}} = Dict{String, Array{ComplexF64,3}}() + kinetic_matrices::Dict{String,Array{ComplexF64,3}} = Dict{String,Array{ComplexF64,3}}() completed::Bool = false end diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 75f527ae2..867fafff4 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -13,15 +13,16 @@ then write to gpec.h5 in a single pass. Write KineticForces results to the "KineticForces" group in gpec.h5. # Arguments -- `h5file::HDF5.File`: Open HDF5 file handle -- `state::KineticForcesState`: Accumulated computation results -- `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is - written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) - is directly available -- `species_label`: `nothing` writes the run total to `KineticForces//`; a label - (e.g. `"ion_z1_m2"`, `"electron"`) writes one species' contribution to - `KineticForces/PerSpecies/