From 26fe4dea3d0ed6c5dcca76ce2b03ede99e159bfe Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 22 Aug 2026 11:31:49 -0400 Subject: [PATCH] PE/EQUIL - BUGFIX! - Withhold unvalidated kinetic tearing diagnostics; free the grid they constrained Delta' and its family (resonant current, resonant area-weighted field, island widths, Chirikov) all descend from the outer-solution jump across the rational surface, read through ideal outer-region asymptotic matching at an IDEAL rational surface. Kinetic terms displace the surfaces (kinsing) and change the layer response, so that reading has no validated basis -- withhold them in kinetic runs rather than deliver numbers with no established meaning. One opt-in (unvalidated_kinetic_tearing, default false) restores them for research use with a warning. Surface metadata and the inner-layer penetrated field are unaffected, as are all ideal runs. Rational-surface bracketing exists solely to keep that ideal Delta' stencil consistent across the surface, so when the Delta' is withheld the bracketing is skipped too and the fine near-rational knots survive to resolve the kinetic matrices' layer structure. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/Equilibrium/GridRefinement.jl | 11 +++++-- src/GeneralizedPerturbedEquilibrium.jl | 12 ++++++- .../PerturbedEquilibriumStructs.jl | 22 +++++++------ src/PerturbedEquilibrium/SingularCoupling.jl | 32 +++++++++++++++++-- 4 files changed, 63 insertions(+), 14 deletions(-) diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index d9441bbe8..bf8a695bd 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; `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); +whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket (`bracket_mandatory = false` skips the bracketing and keeps the fine knots inside instead); `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 @@ -429,6 +429,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium; kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], pinned::Vector{Float64}=Float64[], + bracket_mandatory::Bool=true, singfac_min::Float64=1e-4, n_min::Int=1, bracket_coef::Float64=BRACKET_COEF, @@ -455,7 +456,13 @@ function refined_psi_grid(equil::PlasmaEquilibrium; # 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 + # `bracket_mandatory = false` keeps the locally-uniform fine patch around each rational + # instead of clearing it: the cleared zone exists solely so the ideal Δ′ stencil sees a + # consistent 3rd derivative across the surface, so it is only worth its cost when that Δ′ is + # actually delivered (kinetic runs withhold it — see PerturbedEquilibriumControl + # `unvalidated_kinetic_tearing`). The knots inside are what resolve the kinetic matrices' + # near-rational structure. + (isempty(mandatory) || !bracket_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) end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index b95d2f263..7e0e2a14e 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -273,9 +273,19 @@ function main_from_inputs( isempty(pinned) || @info "Pinning $(length(pinned)) kinetic-resonance surfaces into the ψ grid: $(round.(sort(pinned); digits=3))" end + # Rational bracketing serves the ideal Δ′ stencil. Kinetic runs withhold that Δ′ (see + # PerturbedEquilibriumControl `unvalidated_kinetic_tearing`), so unless the user opts back + # in, skip the bracketing and let the fine near-rational knots resolve the kinetic + # matrices instead. Read from the raw inputs: the PE control struct is built later. + pe_raw = get(inputs, "PerturbedEquilibrium", Dict{String,Any}()) + wants_kinetic_tearing = Bool(get(pe_raw, "unvalidated_kinetic_tearing", false)) + bracket_mandatory = !(ctrl.kinetic_factor > 0) || wants_kinetic_tearing + bracket_mandatory || + @info "Kinetic run: rational-surface bracketing disabled (the ideal Δ′ it serves is withheld); " * + "near-rational knots retained for the kinetic matrices" psi_nodes = Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned, - singfac_min=ctrl.singfac_min, n_min=n_min) + bracket_mandatory=bracket_mandatory, 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 # path re-runs read_imas, which must resolve the same psihigh both passes; diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 1ce0895d6..329cb1c20 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -26,8 +26,10 @@ Medium Priority (defer for MWE): - `singular_point_method::String` - Method for singular point treatment (default: "standard") Regularization: - # High Priority (MWE) +# High Priority (MWE) + - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. + - `unvalidated_kinetic_tearing::Bool` - Emit the jump-derived tearing diagnostics (Δ′, resonant current, island widths, Chirikov) in **kinetic** runs. Default `false`: those quantities are withheld because their extraction assumes ideal outer-region asymptotic matching at an ideal rational surface, which kinetic terms invalidate — the surfaces are kinetically displaced (`kinsing`) and the layer response is not the ideal one, so the numbers would be delivered with no validated meaning. Setting `true` restores them for research use, at the user's risk, with a warning. Ideal runs are unaffected. Enabling it also restores rational-surface bracketing of the ψ grid (see `Equilibrium.refined_psi_grid`), which exists to keep the ideal Δ′ stencil consistent. """ @kwdef struct PerturbedEquilibriumControl # High Priority (MWE) @@ -48,6 +50,7 @@ Regularization: # Regularization width for singular surface smoothing (matches Fortran gpec.f reg_spot). # Set to 0 to disable regularization. Must be non-negative. reg_spot::Float64 = 5e-2 + unvalidated_kinetic_tearing::Bool = false end """ @@ -121,6 +124,7 @@ Metadata [n_rational] — identifies each (surface, n) row: Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field representations (all tesla; no flux/weber is stored): + - `forcing_b`/`response_b` - bare normal field b (Σ⁻¹·b̃) - `forcing_b_rootarea`/`response_b_rootarea` - root-area-weighted field b̃ (coordinate-invariant) - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄) @@ -185,18 +189,18 @@ well-conditioned flux-space inductances L, Λ: rational_surface_idx::Vector{Int} = Int[] # Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla - forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x) - forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant) - forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ - response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x) + forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x) + forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant) + forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ + response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x) response_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ - response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ + response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ # Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space - plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space) + plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space) surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # L̃ (field space) - permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R - reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R + permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R + reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # S = Σ/√A at psilim: b̃→b̄ recovery operator surface_area::Float64 = 0.0 # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index dd24dbfb1..ace45d58d 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -134,7 +134,7 @@ function _solution_at( # Same-side candidate nodes around the bracket, trimmed to the 4 nearest psi. side = sign(psi - psi_surf) - idxs = [j for j in max(1, il-3):min(nstep, ir+3) if sign(odet.psi_store[j] - psi_surf) == side] + idxs = [j for j in max(1, il - 3):min(nstep, ir + 3) if sign(odet.psi_store[j] - psi_surf) == side] while length(idxs) > 4 abs(odet.psi_store[idxs[1]] - psi) > abs(odet.psi_store[idxs[end]] - psi) ? popfirst!(idxs) : pop!(idxs) end @@ -307,7 +307,7 @@ function compute_singular_coupling_metrics!( state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) else state.C_penetrated_area_weighted_field = zeros(ComplexF64, 0, 0) - @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog=1 + @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog = 1 end state.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total) state.rational_psi = zeros(Float64, n_rational) @@ -493,6 +493,34 @@ function compute_singular_coupling_metrics!( # Phase 5: Island diagnostics from applied resonant vectors compute_island_diagnostics!(state, n_rational) + # Kinetic guard: every quantity above descends from the outer-solution jump across the + # rational surface, whose reading as a tearing index assumes ideal asymptotic matching at an + # ideal rational surface. Kinetic terms displace the surfaces (`kinsing`) and change the + # layer response, so that reading is unvalidated — withhold rather than deliver numbers with + # no established meaning. Surface metadata (ψ, q, m, n, area) and the inner-layer penetrated + # field are kept; they do not come from the jump. + if ffit.kinetic_populated && !ctrl.unvalidated_kinetic_tearing + state.C_delta_prime = zeros(ComplexF64, 0, 0) + state.delta_prime = ComplexF64[] + state.C_resonant_current = zeros(ComplexF64, 0, 0) + state.resonant_current = ComplexF64[] + state.C_resonant_area_weighted_field = zeros(ComplexF64, 0, 0) + state.resonant_area_weighted_field = ComplexF64[] + state.C_island_width_sq = zeros(ComplexF64, 0, 0) + state.island_width_sq = ComplexF64[] + state.island_half_width = Float64[] + state.chirikov_parameter = Float64[] + @warn "Kinetic run: withholding the jump-derived tearing diagnostics (Δ′, resonant current, " * + "resonant area-weighted field, island widths, Chirikov). Their extraction assumes ideal " * + "outer-region matching at an ideal rational surface, which kinetic terms invalidate; no " * + "validated kinetic formulation exists in this code. Set [PerturbedEquilibrium] " * + "unvalidated_kinetic_tearing = true to emit them anyway (research use, at your own risk)." + elseif ffit.kinetic_populated + @warn "Kinetic run with unvalidated_kinetic_tearing = true: Δ′ and the island/current diagnostics " * + "derived from the resonant-surface jump have NO validated kinetic formulation. Do not " * + "publish these numbers without independent verification." + end + if ctrl.verbose max_island = isempty(state.island_half_width) ? 0.0 : maximum(state.island_half_width) max_dp = isempty(state.delta_prime) ? 0.0 : maximum(abs.(real.(state.delta_prime)))