From a47f409f0e896317e6fdd3802746afb0d2f94f62 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 28 Aug 2026 17:07:47 -0400 Subject: [PATCH 1/2] ForceFreeStates - REFACTOR - Single-source the auto chunk-count target The nchunks=0 target and the crossing floor were mirrored by hand in two test files. They are now min_crossing_chunks(msing) and auto_chunk_target(msing), used by balance_integration_chunks and both tests; one deliberate assertion pins the formula so a src change still fails loudly instead of propagating silently. Co-Authored-By: Claude Fable 5 --- src/ForceFreeStates/EulerLagrange.jl | 33 +++++++++++++++++------ test/runtests_decomposition_invariance.jl | 9 +++---- test/runtests_parallel_integration.jl | 12 ++++----- 3 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 88c1d9e4a..de8d8a743 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -61,6 +61,29 @@ function ode_itime_cost(psi1::Float64, psi2::Float64, intr::ForceFreeStatesInter return cost end +""" + min_crossing_chunks(msing) -> Int + +Floor on the chunk count set by the `msing` rational-surface crossings; explicit `nchunks` +requests below it are clamped up (with a warning) by `balance_integration_chunks`. +""" +min_crossing_chunks(msing::Integer) = 2 * msing + 3 + +""" + auto_chunk_target(msing) -> Int + +The chunk-count target `balance_integration_chunks` resolves `nchunks = 0` to. Beyond the +crossing floor, BVP propagator conditioning needs at least 8 non-crossing sub-chunks per +segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge) plus the crossing chunks — without them, +`assemble_fm_matrix(condition=true)` cannot keep accumulated products well-conditioned +because single long-span propagators may already have cond ~ 10²⁴ (STRIDE uses 33 intervals +for comparable problems). Derived from `msing` alone — never from the thread count — which +is what keeps the chunk list, and hence every Riccati output, identical whatever `julia -t` +provides. Tests that steer decompositions relative to the auto target call this instead of +mirroring the formula. +""" +auto_chunk_target(msing::Integer) = max(min_crossing_chunks(msing), 8 * (msing + 1) + msing) + """ balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk} @@ -79,20 +102,14 @@ each original chunk retains `needs_crossing=true` and the original `ising`, so t rational surface crossing still fires at the correct ψ in the serial assembly phase. """ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal) - min_chunks = 2 * intr.msing + 3 - # Ensure enough sub-chunks for BVP propagator conditioning: at least 5 non-crossing - # sub-chunks per segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge), plus crossing - # chunks. STRIDE uses 33 intervals for comparable problems. Without enough sub-chunks, - # assemble_fm_matrix(condition=true) can't keep accumulated products well-conditioned - # because single long-span propagators may already have cond ~ 10²⁴. - min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing + min_chunks = min_crossing_chunks(intr.msing) if ctrl.nchunks > 0 if ctrl.nchunks < min_chunks @warn "nchunks = $(ctrl.nchunks) is below the $min_chunks chunks required by $(intr.msing) singular surfaces; clamping up." end target_n = max(ctrl.nchunks, min_chunks) else - target_n = max(min_chunks, min_bvp_intervals) + target_n = auto_chunk_target(intr.msing) end result = collect(chunks) diff --git a/test/runtests_decomposition_invariance.jl b/test/runtests_decomposition_invariance.jl index 097af281d..84bcd9052 100644 --- a/test/runtests_decomposition_invariance.jl +++ b/test/runtests_decomposition_invariance.jl @@ -49,11 +49,10 @@ end @test auto.delta_prime !== nothing msing = size(auto.delta_prime.matrix, 1) - # The nchunks=0 target, mirroring balance_integration_chunks' internal formula (as - # runtests_parallel_integration.jl does). Invariance is asserted ABOVE this floor only: - # fewer chunks than the msing-derived minimum is not a different decomposition but a - # structurally deficient one (the floor gives the rational-surface crossings room). - auto_target = max(2 * msing + 3, 8 * (msing + 1) + msing) + # The nchunks=0 target, from its single source. Invariance is asserted ABOVE this floor + # only: fewer chunks than the msing-derived minimum is not a different decomposition but + # a structurally deficient one (the floor gives the rational-surface crossings room). + auto_target = GP_TI.ForceFreeStates.auto_chunk_target(msing) finer = _solve_at_nchunks(dir, auto_target + 11) @test finer.delta_prime !== nothing diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 7e7737e06..0b3f4fe28 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -134,11 +134,11 @@ using TOML base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) - # Must mirror balance_integration_chunks' internal target_n formula for nchunks = 0 - # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. The formula reads only - # intr.msing — no thread count enters it, which is what makes Riccati outputs - # independent of how many threads `julia -t` provides. - target_n = max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) + # The nchunks = 0 target, from its single source. It reads only intr.msing — no + # thread count enters it, which is what makes Riccati outputs independent of how + # many threads `julia -t` provides. + target_n = GeneralizedPerturbedEquilibrium.ForceFreeStates.auto_chunk_target(intr.msing) + @test target_n == max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) # After balancing, chunk count equals target_n: the while-loop adds exactly one # chunk per iteration (a bisection split) and exits when length(result) >= target_n, @@ -183,7 +183,7 @@ using TOML @test length(balanced_more) == target_n + 7 # An nchunks below the singular-surface floor is clamped up, with a warning. - min_chunks = 2 * intr.msing + 3 + min_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.min_crossing_chunks(intr.msing) ctrl_few = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., nchunks=1) balanced_few = @test_logs (:warn,) match_mode=:any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr) From b7ace78939bc05d1dca40700689a308b57d58a62 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 23 Sep 2026 16:39:45 -0400 Subject: [PATCH 2/2] ForceFreeStates - MINOR - Finish single-sourcing the chunk target and pin it by value MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The struct docstring, docs/src/stability.md and the test comment now point at auto_chunk_target / min_crossing_chunks instead of restating the formula. - The test pins auto_chunk_target and min_crossing_chunks by value (a spec), not by a copy of the expression. - auto_chunk_target drops a max whose floor could never win (9m+8 > 2m+3 for m ≥ 0); value-identical for every msing. Co-Authored-By: Claude Opus 5.5 --- docs/src/stability.md | 6 +- src/ForceFreeStates/EulerLagrange.jl | 63 ++++++------- src/ForceFreeStates/ForceFreeStatesStructs.jl | 30 +----- test/runtests_parallel_integration.jl | 93 ++++++++++++------- 4 files changed, 90 insertions(+), 102 deletions(-) diff --git a/docs/src/stability.md b/docs/src/stability.md index 529304a90..6b5589ddb 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -108,9 +108,9 @@ nchunks = 0 # 0 = auto: derived from the singular-surface count alone #### Chunking and thread independence `balance_integration_chunks` splits the base chunks until the count reaches a target derived -from the number of singular surfaces, `max(2 m_s + 3, 8(m_s + 1) + m_s)`, where ``m_s`` is -`msing`. Setting `nchunks` overrides that target; a value below the ``2 m_s + 3`` floor is -clamped up with a warning. The target never consults `Threads.nthreads()`, and every chunk +from the number of singular surfaces alone ([`auto_chunk_target`](@ref)). Setting `nchunks` +overrides that target; a value below the [`min_crossing_chunks`](@ref) floor is clamped up with a +warning. The target never consults `Threads.nthreads()`, and every chunk integrates independently from identity initial conditions, so the results do not depend on how many threads `julia -t` provides — threads change wall-clock only. diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index de8d8a743..ec33e5025 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -37,9 +37,9 @@ end # at the interval endpoints. Coefficients are ported from STRIDE's ode_itime cost model # (Fortran reference) and unchanged here. Tune only after re-fitting against a per-chunk # step-count sweep; touching these affects parallel-chunk load balancing. -const ODE_COST_AXIS = (a = 39695.0, b = 212830.0) -const ODE_COST_RAT = (a = 17147.0, b = 470710.0) -const ODE_COST_EDGE = (a = 1646.0, b = 4683.0) +const ODE_COST_AXIS = (a=39695.0, b=212830.0) +const ODE_COST_RAT = (a=17147.0, b=470710.0) +const ODE_COST_EDGE = (a=1646.0, b=4683.0) """ ode_itime_cost(psi1, psi2, intr) -> Float64 @@ -72,17 +72,12 @@ min_crossing_chunks(msing::Integer) = 2 * msing + 3 """ auto_chunk_target(msing) -> Int -The chunk-count target `balance_integration_chunks` resolves `nchunks = 0` to. Beyond the -crossing floor, BVP propagator conditioning needs at least 8 non-crossing sub-chunks per -segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge) plus the crossing chunks — without them, -`assemble_fm_matrix(condition=true)` cannot keep accumulated products well-conditioned -because single long-span propagators may already have cond ~ 10²⁴ (STRIDE uses 33 intervals -for comparable problems). Derived from `msing` alone — never from the thread count — which -is what keeps the chunk list, and hence every Riccati output, identical whatever `julia -t` -provides. Tests that steer decompositions relative to the auto target call this instead of -mirroring the formula. +The chunk-count target `balance_integration_chunks` resolves `nchunks = 0` to: 8 non-crossing +sub-chunks per segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge) plus the crossing chunks, so that +`assemble_fm_matrix(condition=true)` can keep the accumulated propagator products well-conditioned +(STRIDE uses 33 intervals for comparable problems). Always above [`min_crossing_chunks`](@ref). """ -auto_chunk_target(msing::Integer) = max(min_crossing_chunks(msing), 8 * (msing + 1) + msing) +auto_chunk_target(msing::Integer) = 8 * (msing + 1) + msing """ balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk} @@ -95,7 +90,7 @@ never from `Threads.nthreads()` — so the chunk list, and hence every Riccati o identical whatever thread count `julia -t` provides. Each split finds the equal-cost midpoint ψ_mid via bisection: - ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2 +ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2 Sub-chunks inherit `needs_crossing=false` and `ising=0`. Only the LAST sub-chunk of each original chunk retains `needs_crossing=true` and the original `ising`, so the @@ -148,10 +143,10 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc psi_mid = (lo + hi) / 2.0 left = IntegrationChunk(; psi_start=chunk.psi_start, psi_end=psi_mid, - needs_crossing=false, ising=0, direction=1) + needs_crossing=false, ising=0, direction=1) right = IntegrationChunk(; psi_start=psi_mid, psi_end=chunk.psi_end, - needs_crossing=chunk.needs_crossing, ising=chunk.ising, - direction=chunk.direction) + needs_crossing=chunk.needs_crossing, ising=chunk.ising, + direction=chunk.direction) splice!(result, best_idx, [left, right]) end @@ -301,20 +296,20 @@ is identified by dominant |U₁| component, giving the physically correct consta Frobenius solution and avoiding the spurious logarithmic irregularity. """ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplines, - intr::ForceFreeStatesInternal, psi_low::Float64) - N = intr.numpert_total + intr::ForceFreeStatesInternal, psi_low::Float64) + N = intr.numpert_total hint = Ref(1) # Evaluate stability matrices at psi_low F_lower = zeros(ComplexF64, N, N) - kmat = zeros(ComplexF64, N, N) - gmat = zeros(ComplexF64, N, N) + kmat = zeros(ComplexF64, N, N) + gmat = zeros(ComplexF64, N, N) ffit.fmats_lower(vec(F_lower), psi_low; hint=hint) - ffit.kmats(vec(kmat), psi_low; hint=hint) - ffit.gmats(vec(gmat), psi_low; hint=hint) + ffit.kmats(vec(kmat), psi_low; hint=hint) + ffit.gmats(vec(gmat), psi_low; hint=hint) # singfac[j] = 1 / (m_j − n_j · q) for each mode j - q0 = profiles.q_spline(psi_low; hint=hint) + q0 = profiles.q_spline(psi_low; hint=hint) singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q0 .* (intr.nlow:intr.nhigh)')) # F̄⁻¹ = (F_lower · F_lower')⁻¹ via the Cholesky factor @@ -328,15 +323,15 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin for j in 1:N sf = singfac[j] fi = Finv[j, j] - k = kmat[j, j] + k = kmat[j, j] kd = conj(k) # K̄†[j,j] - g = gmat[j, j] + g = gmat[j, j] # 2×2 ODE matrix block for mode j [Glasser 2016 Eq. 22-24, diagonal approximation] m11 = -sf * fi * k - m12 = sf^2 * fi - m21 = g - kd * fi * k - m22 = sf * kd * fi + m12 = sf^2 * fi + m21 = g - kd * fi * k + m22 = sf * kd * fi # Frobenius matrix A₀_j = ψ_low · M_j [Glasser 2016 Eq. 51] #! format: off @@ -363,7 +358,7 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin if abs(v2) > Base.sqrt(Base.eps(Float64)) * abs(v1) U1_init[j, j] = v1 / v2 else - U1_init[j, j] = one(ComplexF64) + U1_init[j, j] = one(ComplexF64) U2_init[j, j] = zero(ComplexF64) end end @@ -382,7 +377,7 @@ Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, an Move ising_start logic to chunk_el_integration_bounds? """ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ffit::FourFitVars, - profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) + profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) # Default psifac to minimum equilibrium psi value odet.psifac = profiles.xs[1] @@ -559,7 +554,7 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro psi_end=psi_end, needs_crossing=true, ising=ising_current, - direction = bidirectional ? -1 : 1 + direction=bidirectional ? -1 : 1 )) # After crossing, we jump to the other side of the singular surface @@ -992,7 +987,7 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal) temp[ksol, jsol] = odet.fixfac[ksol, jsol, ifix] end end - mul!(gauss_buffer, view(gauss,:,:,ifix), temp) + mul!(gauss_buffer, view(gauss, :, :, ifix), temp) gauss[:, :, ifix] .= gauss_buffer end # Account for zeroed indices at singular surfaces in `ode_ideal_cross` @@ -1009,7 +1004,7 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal) # and mfix + 1 is the for the region after the last fixup and before the edge transforms[:, :, end] .= identity for ifix in odet.ifix:-1:1 - mul!(view(transforms,:,:,ifix), view(gauss,:,:,ifix), view(transforms,:,:,(ifix+1))) + mul!(view(transforms, :, :, ifix), view(gauss, :, :, ifix), view(transforms, :, :, (ifix + 1))) end # Now that we have the transform matrices, we can apply them to the solution vectors diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 79735b895..50d99dfd7 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -262,7 +262,7 @@ gpec.toml. - `save_interval::Int` - Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. (Same as `euler_step` in the Fortran) - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations) - `integrator::String` - Which formalism integrates the Euler-Lagrange system. `"forward"` sweeps the plasma serially with Gaussian reduction and returns `u_store` / `du_store` / `xi_s_store` dense in the axis (EL) basis — the only convention PerturbedEquilibrium and FieldReconstruction consume correctly, and the only path that supports `kinetic_factor > 0`. `"riccati"` (default) runs the chunked fundamental-matrix propagator driver (Glasser 2018 Phys. Plasmas 25, 032507): chunks are integrated independently from identity initial conditions and assembled serially with Riccati-style crossings, which is the only way to obtain the singular-surface Δ' matrix for the tearing-mode solvers downstream, but leaves `u_store` as sparse chunk-endpoint Riccati states, so dense ξ profiles are unavailable. `"galerkin"` solves the same Euler-Lagrange system variationally instead of by radial ODE integration — the RDCON outer-region singular Galerkin method (Glasser, Wang & Park 2016 Phys. Plasmas 23, 112506), which discretizes the displacement on packed Hermite-cubic elements and solves one global banded system — producing the resistive Δ′ matrix and, when `gal_match_flag` is set, the RPEC inner-layer-matched ξ; it computes its own vacuum response and returns no free-boundary energies, and does not support `kinetic_factor > 0`. Requires `singfac_min != 0` for `"riccati"`. - - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone: `max(2·msing + 3, 8·(msing + 1) + msing)`, enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below `2·msing + 3` is clamped up with a warning. Chunk sizing never consults `Threads.nthreads()`, so Riccati outputs are identical whatever thread count `julia -t` provides; threads only change wall-clock. + - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone ([`auto_chunk_target`](@ref)), enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below [`min_crossing_chunks`](@ref) is clamped up with a warning. Chunk sizing never consults `Threads.nthreads()`, so Riccati outputs are identical whatever thread count `julia -t` provides; threads only change wall-clock. - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64. """ @kwdef struct ForceFreeStatesControl @@ -493,71 +493,46 @@ and a small set of temporary matrices and factors used to compute singular-layer - `numpert_total::Int` - Total number of Fourier mode combinations (m × n) used in the calculation. - `numunorms_init::Int` - Initial allocation size for the number of normalization operations recorded. - - `msing::Int` - Number of singular surfaces in the equilibrium (used to size asymptotic coefficient arrays). - - `numsteps_init::Int` - Initial allocation size for the number of integration steps to store. - - `step::Int` - Current integration step index (1-based, like `istep` in the original Fortran). - - `psi_store::Vector{Float64}` - Stored psi values at each saved integration step (length `numsteps_init`). - - `q_store::Vector{Float64}` - Stored q values at each saved integration step (length `numsteps_init`). - - `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape `(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver). - - `du_store::Array{ComplexF64,3}` - dΞ_ψ/dψ (the u₁ block only) at each saved step, shape `(numpert_total, numpert_total, step)`. Empty until `materialize_derivative_stores!` fills it, except on the galerkin-matched path which supplies the analytic derivative at construction. du₂/dψ is never stored densely — its only consumer evaluates it on demand at bracket nodes. - - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s at each saved step, eq. 18 of Glasser 2016, shape `(numpert_total, numpert_total, step)`. Empty until materialized, same as `du_store`. - - `u_store_el_basis::Bool` - True when `u_store` holds the Euler-Lagrange state `(u₁, u₂)`, so the derivative kernel can be re-applied to it. False on the sparse parallel path, whose stored columns are chunk-endpoint Riccati matrices; `materialize_derivative_stores!` refuses to run there. - - `du_store_populated::Bool` - True once `du_store`/`xi_s_store` hold valid data in the final (post-transform, post-normalization) basis. Set by `materialize_derivative_stores!` or by the galerkin-matched constructor; stays false where the stores cannot be materialized, e.g. the sparse parallel path whose solution is in the Riccati basis. - - `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`). - - `ca_r::Array{ComplexF64,4}` - Asymptotic coefficients just to the right of each singular surface with shape `(numpert_total, numpert_total, 2, msing)`. - - `ca_l::Array{ComplexF64,4}` - Asymptotic coefficients just to the left of each singular surface with shape `(numpert_total, numpert_total, 2, msing)`. - - `edge_scan::EdgeScanState` - Edge dW scan state and results. Initialized as a disabled sentinel (N_edge=0) and replaced by `findmax_dW_edge!` when a scan runs. - - `psifac::Float64` - Current normalized flux coordinate for the integrator. - - `q::Float64` - Safety factor value at `psifac` (current q during integration). - - `u::Array{ComplexF64,3}` - Current working solution arrays with shape `(numpert_total, numpert_total, 2)`. - - `ising_start::Int` - Index of the starting singular surface to be crossed during integration. - - `psimax::Float64` - Maximum psi value for which the integrator is allowed to run in next integration region. - - `needs_crossing::Bool` - Flag indicating whether a rational surface needs to be crossed after the current integration region. - - `nzero::Int` - Count of detected zero crossings (used for diagnostics). - - `new::Bool` - Flag indicating whether a new `unorm0` should be computed after a fixup. # Initialization parameters - - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`). - - `unorm0::Vector{Float64}` - Reference/initial norms of the solution vectors (length `numpert_total`). # Saved data throughout integration - - `ifix::Int` - Number of normalization operations performed (index into normalization arrays). # Total ODE solver steps taken (all steps, not just saved ones) @@ -566,11 +541,8 @@ and a small set of temporary matrices and factors used to compute singular-layer - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) (length `numunorms_init`). - - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed. # Data for integrator - - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`. - - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`). """ @kwdef mutable struct OdeState diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 0b3f4fe28..804ff1859 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -28,8 +28,8 @@ using TOML GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator!(odet, prop) - @test odet.u[:, :, 1] ≈ u1_in rtol=1e-12 - @test odet.u[:, :, 2] ≈ u2_in rtol=1e-12 + @test odet.u[:, :, 1] ≈ u1_in rtol = 1e-12 + @test odet.u[:, :, 2] ≈ u2_in rtol = 1e-12 end @testset "apply_propagator! linearity" begin @@ -65,8 +65,8 @@ using TOML u1_expected = U1_upper * u1_in + U1_lower * u2_in u2_expected = U2_upper * u1_in + U2_lower * u2_in - @test odet.u[:, :, 1] ≈ u1_expected rtol=1e-12 - @test odet.u[:, :, 2] ≈ u2_expected rtol=1e-12 + @test odet.u[:, :, 1] ≈ u1_expected rtol = 1e-12 + @test odet.u[:, :, 2] ≈ u2_expected rtol = 1e-12 end @testset "apply_propagator_inverse! is inverse of apply_propagator!" begin @@ -100,13 +100,13 @@ using TOML GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator_inverse!(odet, prop) GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator!(odet, prop) - @test odet.u[:, :, 1] ≈ u1_in rtol=1e-12 - @test odet.u[:, :, 2] ≈ u2_in rtol=1e-12 + @test odet.u[:, :, 1] ≈ u1_in rtol = 1e-12 + @test odet.u[:, :, 2] ≈ u2_in rtol = 1e-12 end @testset "balance_integration_chunks produces target count" begin - # Verify that balance_integration_chunks creates at least - # max(2*msing+3, 8*(msing+1)+msing) chunks from a small set of base chunks. + # Verify that balance_integration_chunks reaches auto_chunk_target(msing) chunks from a + # small set of base chunks. ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false @@ -114,9 +114,12 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -134,11 +137,11 @@ using TOML base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) - # The nchunks = 0 target, from its single source. It reads only intr.msing — no - # thread count enters it, which is what makes Riccati outputs independent of how - # many threads `julia -t` provides. - target_n = GeneralizedPerturbedEquilibrium.ForceFreeStates.auto_chunk_target(intr.msing) - @test target_n == max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) + # nchunks = 0 target, pinned by value so a formula change is a visible spec change. + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + @test FFS.auto_chunk_target.((0, 1, 5)) == (8, 17, 53) + @test FFS.min_crossing_chunks.((0, 1, 5)) == (3, 5, 13) + target_n = FFS.auto_chunk_target(intr.msing) # After balancing, chunk count equals target_n: the while-loop adds exactly one # chunk per iteration (a bisection split) and exits when length(result) >= target_n, @@ -156,7 +159,7 @@ using TOML # gap of ≈2·singfac_min/|n·q1| between the pre-crossing and post-crossing intervals. for i in eachindex(balanced)[2:end] if !balanced[i-1].needs_crossing - @test balanced[i].psi_start ≈ balanced[i-1].psi_end rtol=1e-10 + @test balanced[i].psi_start ≈ balanced[i-1].psi_end rtol = 1e-10 else # Inner-layer gap: post-crossing chunk starts AFTER the rational surface @test balanced[i].psi_start > balanced[i-1].psi_end @@ -186,7 +189,7 @@ using TOML min_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.min_crossing_chunks(intr.msing) ctrl_few = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., nchunks=1) - balanced_few = @test_logs (:warn,) match_mode=:any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr) + balanced_few = @test_logs (:warn,) match_mode = :any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr) @test length(balanced_few) == max(min_chunks, length(base_chunks)) end @@ -201,9 +204,12 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -262,11 +268,14 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -311,7 +320,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) # Apply the two-pass auto grid exactly as the main driver does (the example ships # grid_type="auto", mpsi=0): measured-curvature refinement with rational surfaces # pinned as mandatory knots, re-formed from the captured ingest. The pinned values @@ -324,8 +336,8 @@ using TOML end intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -372,9 +384,12 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -423,11 +438,14 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], example_dir) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil) @@ -476,7 +494,10 @@ using TOML ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex) - equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config, haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing) + equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium( + eq_config, + haskey(inputs, "SOL_INPUT") ? GeneralizedPerturbedEquilibrium.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + ) # Apply the two-pass auto grid (measured-curvature refinement, rational surfaces pinned # as mandatory knots) exactly as the main driver does (see the FM testset above); the # pinned values below are for this grid, the production default. @@ -488,8 +509,8 @@ using TOML end intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - intr.nlow = ctrl.nn_low; - intr.nhigh = ctrl.nn_high; + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high intr.npert = 1 GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)