diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 858bb53a0..6cbde2a58 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -138,6 +138,42 @@ 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.) + +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. + +Measured on the DIII-D-like H-mode case (n = 1, C-coil drive), comparing the NTV torque against +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: + +| configuration | max ``|\xi^\alpha|`` | NTV torque [N·m] | EL dissipation [N·m] | +|---|---|---|---| +| ideal, `reg_spot = 0` | 643.8 | 6074.3 | — | +| ideal, `reg_spot = 0.05` | 0.058 | 0.554 | — | +| kinetic, `reg_spot = 0.05` | 0.055 | 0.1655 | 0.1322 | +| kinetic, `reg_spot = 0` | 0.059 | **0.1324** | **0.1322** | + +The ideal rows show why the knob exists; the kinetic rows show why it must be off there — the +two independent torques agree to 0.15 % with no regularization, and to 20 % with it. + ## 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 24b2b7f24..148d8ebf8 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -87,5 +87,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 where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) 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..e174c5ab9 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 where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index 47e2e603b..7854db0e1 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 where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..d9441bbe8 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -418,7 +418,7 @@ Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curv 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; +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 +428,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 +450,11 @@ 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 (e.g. kinetic-resonance surfaces): knot-at-node semantics via + # merge_mandatory_nodes — no cleared zone. Inserted before rational bracketing, so a + # pinned node inside a rational's bracket zone is cleared by it (the Δ′ clean-interval + # requirement wins locally; the rational's own dense floor resolves that neighbourhood). + 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 0f8095e3c..57cc28128 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -443,6 +443,7 @@ rational's resolution window, preserving the Δ′-stencil structure the equilib (`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) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 0cd8e01dd..e057276e1 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,3 +1,74 @@ + +# Knots laid across the envelope's [ψ_c, 2ψ_c] transition. The envelope is a quintic smoothstep, +# so a cubic spline needs several interior knots to follow it without overshoot; nine keeps the +# residual below the kernel's own tolerance 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 +98,19 @@ 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) + # The near-axis validity envelope (KineticForces) has structure on the scale of the + # suppression boundary; coarse equilibrium grids cannot represent env·(increment), and the + # spline overshoot can land on a rational surface. Pin the band ends (the smoothstep is + # only C² there) and resolve the transition with a fixed set of knots. + 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 +121,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 +136,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 +176,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 2c7bdac14..9cc7b6112 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 @@ -239,6 +241,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) @@ -257,7 +264,13 @@ Algorithm: 3. Refine each peak with golden-section minimization of -cond 4. Filter 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] @@ -282,11 +295,19 @@ 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 the threshold are not singular surfaces, but they mark where the kinetic F̄ comes + # closest to singular — the shifted/split resonances of Park & Logan Eq. (70). Report the + # strongest few so sharp kinetic structure is visible rather than silent (on a DIII-D-like case + # these track the NTV torque-density peaks at low collisionality/rotation). + 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 f2ed4ef65..aea43f817 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -85,8 +85,8 @@ using .Equilibrium: PlasmaEquilibrium using .ForcingTerms: RMPField const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim", - "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi", - "gal_flag") + "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi", + "gal_flag") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing @@ -210,7 +210,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; + pinned=KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr)) @info "Equilibrium construction completed in $(@sprintf("%.3f", time() - equil_start)) s" @@ -391,7 +392,8 @@ function maybe_reform_equilibrium( additional_input, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, - kinetic_profiles + kinetic_profiles; + pinned::Vector{Float64}=Float64[] ) Equilibrium.wants_two_pass(eq_config) || return equil @@ -399,8 +401,10 @@ function maybe_reform_equilibrium( # 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 @@ -579,12 +583,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; @@ -740,10 +750,10 @@ runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces] and errors here. ```julia -eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") +eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) -ffs = solve(prob, Riccati()) -ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # equivalent one-line form +ffs = solve(prob, Riccati()) +ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # equivalent one-line form ``` """ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrator) @@ -871,7 +881,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. @@ -948,6 +958,8 @@ function run_kinetic_forces( h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) + KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) end end else @@ -973,6 +985,8 @@ function run_kinetic_forces( end KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) + KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) end end end diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 274f54f14..3a88878e8 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -141,6 +141,25 @@ const MAIN_H5_ANNOTATIONS = [ "LocalStability/alpha_critical" => (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("psi_ballooning",), attach=(1 => "LocalStability/ballooning_psi",)), + # --- KineticForces/Validity/ --- + "KineticForces/Validity/psi" => (; long_name="normalized poloidal flux ψ_N of the kinetic validity profiles", scale="psi"), + "KineticForces/Validity/rho_i" => (; long_name="thermal ion gyroradius √(2mT)/(Z·e·B₀)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/rho_banana" => (; long_name="thermal banana orbit width q·ρ_i/√ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/rho_theta" => (; long_name="thermal poloidal gyroradius q·ρ_i/ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/w_potato" => (; long_name="potato orbit width (q²ρ_i²R₀)^(1/3)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/r_minor" => (; long_name="surface-average minor radius ⟨r⟩", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩", units="1"), + "KineticForces/Validity/envelope" => + (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/is_valid" => (; + long_name="1 where every zero-orbit-width ordering holds: max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, max orbit width < d_separatrix", + units="1", + attach=(1 => "KineticForces/Validity/psi",) + ), + # --- ForceFreeStates/Solutions/ForwardIntegration/ --- "ForceFreeStates/Solutions/ForwardIntegration/nstep" => (; long_name="number of saved solution snapshots"), "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" => (; long_name="total ODE solver steps taken"), @@ -245,7 +264,7 @@ const MAIN_H5_ANNOTATIONS = [ "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"), "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"), "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"), - "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"), + "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m") ] # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants. @@ -258,7 +277,7 @@ const _ELM_IDEAL_LETTERS = [ ("H", "Euler-Lagrange primitive coefficient matrix H"), ("F", "Euler-Lagrange derived coefficient matrix F"), ("K", "Euler-Lagrange derived coefficient matrix K"), - ("G", "Euler-Lagrange derived coefficient matrix G"), + ("G", "Euler-Lagrange derived coefficient matrix G") ] # The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared # unchanged from the ideal set and are not re-emitted. @@ -268,14 +287,19 @@ const _ELM_KINETIC_LETTERS = [ ("C", "Euler-Lagrange primitive coefficient matrix C"), ("K", "Euler-Lagrange derived coefficient matrix K"), ("G", "Euler-Lagrange derived coefficient matrix G"), - ("f0", "raw kinetic component matrix f0"), + ("f0", "raw kinetic component matrix f0") ] const ELM_H5_ANNOTATIONS = vcat( ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")], - ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => - (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS], - ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => - (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS] + [ + "ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => + (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS + ], + [ + "ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => + (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for + (l, d) in _ELM_KINETIC_LETTERS + ] ) """ diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 63a59aedc..f1bb1fc32 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -147,7 +147,7 @@ end """ compute_bounce_data(psi, n, l, q, bo, bmax, bmin, theta_bmax, - tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, + tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda=128, ntheta=128, smat=nothing, tmat=nothing, xmat=nothing, @@ -168,8 +168,8 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch). - `bmax, bmin`: Max/min of B(θ) at this ψ - `theta_bmax`: θ location of Bmax (nodal knot; the passing-transit start) - `tspl`: Periodic poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ] -- `B_extrap`: Endpoint-fit (non-periodic) cubic of B(θ) used for v_par and the - bounce-point roots (the Fortran `vspl` equivalent) +- `B_vpar`: Periodic cubic of B(θ) used for v_par and the bounce-point roots + (the Fortran `vspl` equivalent) - `mfac`: Poloidal mode numbers [mlow:mhigh] - `chi1`: 2π·ψ₀ flux normalization - `ro`: Major radius [m] @@ -190,7 +190,7 @@ function compute_bounce_data( psi::Float64, n::Int, l::Int, q::Float64, bo::Float64, bmax::Float64, bmin::Float64, theta_bmax::Float64, - tspl, B_extrap, mfac::Vector{Int}, chi1::Float64, ro::Float64, + tspl, B_vpar, mfac::Vector{Int}, chi1::Float64, ro::Float64, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, mass::Float64, chrg::Float64, @@ -241,12 +241,12 @@ function compute_bounce_data( # Find bounce points and build θ sub-grid _, _, tdt_pts, tdt_wts = _find_bounce_points_and_grid( - lmda, bo, sigma, B_extrap, theta_bmax, psi, ntheta) + lmda, bo, sigma, B_vpar, theta_bmax, psi, ntheta) # Bounce integrals over θ (Fortran lines 674-735) wbbar, wdbar, dJdJ_val, wmats_lmda = _bounce_integrate( tdt_pts, tdt_wts, lmda, lnq, sigma, n, q, bo, - tspl, B_extrap, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_vpar, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, do_matrices, mpert, smat, tmat, xmat, ymat, zmat, scr) # Physical frequencies (Fortran lines 744-745) @@ -409,12 +409,12 @@ end """ -Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the endpoint-fit cubic of B -(`B_extrap`, built where the surface interpolants are constructed), keeping v_par +Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the periodic cubic of B +(`B_vpar`, built where the surface interpolants are constructed), keeping v_par consistent with the bounce-point roots as in Fortran's `vspl`. """ -@inline _vpar_from_extrap(B_extrap, lmda::Float64, bo::Float64, θ::Float64) = - 1.0 - (lmda / bo) * B_extrap(mod(θ, 1.0)) +@inline _vpar_from_spline(B_vpar, lmda::Float64, bo::Float64, θ::Float64) = + 1.0 - (lmda / bo) * B_vpar(mod(θ, 1.0)) """ @@ -423,14 +423,14 @@ Returns (t1, t2, theta_points, theta_weights). """ function _find_bounce_points_and_grid( lmda::Float64, bo::Float64, sigma::Int, - B_extrap, theta_bmax::Float64, psi::Float64, + B_vpar, theta_bmax::Float64, psi::Float64, ntheta::Int ) if sigma == 0 # trapped - # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1), + # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_vpar(θ) in (0,1), # sorted descending — the same order as Fortran spline_roots, which the # marginally-trapped and deepest-well wrap logic below assume. - vpar_fn = θ -> _vpar_from_extrap(B_extrap, lmda, bo, θ) + vpar_fn = θ -> _vpar_from_spline(B_vpar, lmda, bo, θ) bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true) nbpts = length(bpts) @@ -443,7 +443,7 @@ function _find_bounce_points_and_grid( t1 = bpts[1] t2 = bpts[1] + 1.0 else - t1, t2 = _find_deepest_well(bpts, B_extrap, lmda, bo) + t1, t2 = _find_deepest_well(bpts, B_vpar, lmda, bo) end # Power-law grid refined near bounce points @@ -463,7 +463,7 @@ end Find the deepest potential well (largest midpoint v_par) among bounce-point pairs, handling pairs that wrap through θ = 0/1. """ -function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo::Float64) +function _find_deepest_well(bpts::Vector{Float64}, B_vpar, lmda::Float64, bo::Float64) nbpts = length(bpts) best_vpar = 0.0 best_t1 = 0.0 @@ -477,7 +477,7 @@ function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo:: else θmid = 0.5 * (bpts[i] + bpts[j]) end - vpar_mid = _vpar_from_extrap(B_extrap, lmda, bo, θmid) + vpar_mid = _vpar_from_spline(B_vpar, lmda, bo, θmid) if vpar_mid > best_vpar best_t1 = bpts[i] best_t2 = bpts[j] @@ -506,7 +506,7 @@ Ports Fortran torque.F90 lines 674-793. function _bounce_integrate( tdt_pts::Vector{Float64}, tdt_wts::Vector{Float64}, lmda::Float64, lnq::Float64, sigma::Int, n::Int, q::Float64, bo::Float64, - tspl, B_extrap, chi1::Float64, ro::Float64, + tspl, B_vpar, chi1::Float64, ro::Float64, mfac::Vector{Int}, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, do_matrices::Bool, mpert::Int, @@ -550,9 +550,9 @@ function _bounce_integrate( jac = tspl_f[4] djdpsi = tspl_f[5] - # v_par from the endpoint-fit cubic (consistent with the bounce points); + # v_par from the periodic cubic (consistent with the bounce points); # the periodic tspl B_val remains the numerator field in the integrands. - vpar = 1.0 - (lmda / bo) * B_extrap(θmod) + vpar = 1.0 - (lmda / bo) * B_vpar(θmod) if vpar <= 0 # Negative v_par near a bounce point: same fill rules as the Fortran bounce loop. @@ -632,9 +632,13 @@ function _bounce_integrate( return 0.0, 0.0, 0.0, nothing end - # Bounce-averaged frequencies + # Bounce-averaged frequencies. wbbar already carries one factor of ro that its own + # normalization bhat = sqrt(2T/m)/ro cancels; reusing it inside wdbar imports that ro + # a third time while dhat = (T/q)/(bo·ro²) removes only the two written explicitly, so + # the drift prefactor takes ro, not ro². (Otherwise ω_D = wdbar·dhat carries a surplus + # length: 4π·(I₂/I₁)·(T/q) is already V/Wb = 1/s, so the extra ro leaves m/s.) wbbar = ro * twopi / ((2 - sigma) * total_wb) - wdbar = ro^2 * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd + wdbar = ro * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd # Phase factor pl_i = exp(-2πi·lnq·fsi_wb(θ_i)/((2-σ)·total_wb)), using the # cumulative spline integral of the bounce action. diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 6c615a3c6..de916b0ac 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,19 @@ 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: suppress the drift-kinetic increments where the zero-orbit-width + # ordering fails, taking the widest-orbit species. Kernel evaluation is skipped where it is 0. + 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 +142,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) @@ -136,7 +155,9 @@ function compute_calculated_kinetic_matrices( z_s, m_s, kf_ctrl.wdfac, kf_ctrl.divxfac, el_s, equil, intr_t, prof_s; nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, - atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda + atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda, + atol_x=kf_ctrl.atol_x, rtol_x=kf_ctrl.rtol_x, + nested_tolerance_margin=kf_ctrl.nested_tolerance_margin ) full_w .+= block_w full_t .+= block_t @@ -146,8 +167,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 d0fdd591b..1ad347561 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -100,6 +100,18 @@ function integrate_psi_quadgk( panel_psis=Float64[], resonance_psis=Float64[]) end + # Near-axis validity suppression: same boundary and envelope as the EL kinetic + # matrices (one source of truth), applied to the torque density; the quadrature + # domain starts at the boundary since the integrand is identically zero below it. + psi_c = ctrl.axis_validity_suppression ? + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) : 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), @@ -137,7 +149,9 @@ function integrate_psi_quadgk( tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac, electron, method, equil, thread_intrs[tid], kinetic_profiles; op_wmats=w, - atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda, + atol_x=ctrl.atol_x, rtol_x=ctrl.rtol_x, + nested_tolerance_margin=ctrl.nested_tolerance_margin) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end @@ -147,6 +161,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) @@ -156,6 +171,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 @@ -164,7 +180,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) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index e4eca2581..8834a689f 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,11 +130,20 @@ builds a second control for its differing tolerance). nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances. - # *_xlmda: shared tolerances for inner λ (pitch) and x (energy) integrations - # *_psi: tolerances for outer ψ quadrature - atol_xlmda::Float64 = 1e-8 # Absolute tolerance for inner pitch + energy integrations - rtol_xlmda::Float64 = 1e-5 # Relative tolerance for inner pitch + energy integrations + # 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 + 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. + 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). @@ -149,6 +167,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" @@ -202,12 +221,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 @@ -297,17 +317,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 @@ -317,20 +337,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) # Copy mode numbers from FFS kf_intr.mlow = ffs.mlow kf_intr.mhigh = ffs.mhigh @@ -399,9 +421,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) @@ -421,8 +443,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 @@ -449,9 +471,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) @@ -521,8 +543,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..771873251 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/