From 1d07e9e7e93b7d8355d0cb4140c381586a0c874e Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Thu, 8 Oct 2026 16:16:18 -0400 Subject: [PATCH] ForceFreeStates - BUGFIX! - Default the axis start to the fixed initialization The diagonal Frobenius start does not recover the axis limit it documented, and its m = 0 branch is an imaginary pair whose sign follows roundoff. frobenius_psi_max now defaults to 0, so that start is opt-in. --- .../gpec.toml | 3 - src/ForceFreeStates/CoreTypes.jl | 4 +- src/ForceFreeStates/EulerLagrange.jl | 90 ++++++++++++------- src/ForceFreeStates/Riccati/Driver.jl | 8 +- test/runtests_eulerlagrange.jl | 30 ++++--- 5 files changed, 81 insertions(+), 54 deletions(-) diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index 581ca3d47..002596c21 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml @@ -4,9 +4,6 @@ # small edge singular factor. The q = 2 and q = 3 surfaces resonate at both n, so SingularCoupling # repeats them across rows. psihigh and delta_mhigh keep the case small for the regression harness; # every n block spans m = -8..20, which a single-n n = 1 run matches only with delta_mhigh = 8. -# frobenius_psi_max = 0 uses the fixed axis start: the Frobenius start's m = 0 branch is picked by -# roundoff, which makes the n = 1 anti-Hermitian part of W (and so the torque) differ between -# otherwise identical runs. [Equilibrium] eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index c982268ba..5544b0e60 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -123,7 +123,7 @@ gpec.toml. - `mthvac::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput) - `nzvac::Int` - Number of vacuum toroidal grid points (corresponds to `nzeta` in VacuumInput3D) - `sing_start::Int` - Start integration at the `sing_start`-th singular surface - - `frobenius_psi_max::Float64` - Largest starting ψ_N at which the Frobenius start is used. The Frobenius start [Glasser 2016 Eq. 51] is a power series about the magnetic axis, so an integration beginning above this value (an interior start from a raised `psilow`, or `qlow` above q₀) instead starts every solution with zero displacement and unit "momentum" (U₁ = 0, U₂ = I), the initialization Fortran DCON always uses, and warns that it did. `0` selects the fixed start everywhere, silently. Default 0.01, which keeps every axis-started case on the Frobenius start. + - `frobenius_psi_max::Float64` - Largest starting ψ_N at which the opt-in Frobenius start is used. Default `0` selects the fixed start everywhere and silently: U₁ = 0, U₂ = I, the initialization Fortran DCON uses. A positive value uses the diagonal Frobenius 2×2 when the integration starts at or below it, and warns that this start is currently unreliable. A start above the threshold uses the fixed start with no warning. That 2×2 does not recover a displacement that vanishes as a positive power of the starting ψ, and its m = 0 branch is an imaginary pair whose sign follows roundoff. Leave this at 0 unless the Frobenius start itself is under test. - `nn_low::Int` - Lower bound for toroidal modes - `nn_high::Int` - Upper bound for toroidal modes - `delta_mlow::Int` - Expands lower bound of Fourier harmonics by delta_mlow @@ -221,5 +221,5 @@ gpec.toml. gal_rho::Vector{Float64} = Float64[] # per-surface mass density ρ [kg/m³] (length msing, core→edge); Fortran rmatch `massden` gal_rotation::Vector{Float64} = Float64[] # per-surface rotation frequency f [Hz] (length msing, core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation` gal_gamma::Float64 = 5 / 3 # ratio of specific heats Γ for the resistive-layer coefficients (resist_eval G term) - frobenius_psi_max::Float64 = 0.01 + frobenius_psi_max::Float64 = 0.0 end diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 48e08b062..2224c53e2 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -45,49 +45,81 @@ 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)`. + - `ca_populated::Bool` - True once an ideal singular-surface crossing has filled `ca_l`/`ca_r`; kinetic and galerkin-matched runs never populate them and leave this false, and the HDF5 writer then emits zero-extent `ca_left`/`ca_right` datasets instead of unpopulated arrays. + - `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. + + # Initialization parameters + - `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 + + # Saved data throughout integration + - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`). + +# Total ODE solver steps taken (all steps, not just saved ones) + - `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) @@ -96,8 +128,15 @@ 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`). + + # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) + + # Initialization parameters + - `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)`. + + - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`. # Data for integrator + - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`). """ @kwdef mutable struct OdeState @@ -404,20 +443,15 @@ end """ compute_axis_init(mats, profiles, intr, psi_low) -> (U1_init, U2_init) -Compute axis initial conditions for the Euler-Lagrange ODE via the Frobenius -leading-coefficient eigenvalue problem [Glasser Phys. Plasmas 2016 112506 Eq. 51]: - - lim_{ψ→0} [ψ M(ψ) − a I] v = 0 - -For each mode j, solves the 2×2 Frobenius eigenvalue problem for the diagonal block of -A₀ = ψ_low · M(ψ_low). The eigenvector with Re(a) ≥ 0 is the regular (non-singular) -Frobenius solution. Returns U₁_init and U₂_init normalized so that U₂_init = I (consistent -with the N independent solutions convention). +Opt-in diagonal approximation to the axis Frobenius problem [Glasser Phys. Plasmas 2016 112506 Eq. 51], +evaluated at finite `psi_low` rather than in the limit ψ → 0. Used only when `frobenius_psi_max` is +positive and the integration starts at or below it. -For m≠0 the regular eigenvector has a negligible U₁ component (~ψ_low^(|m|/2)), recovering -the Glasser [0, I] limit as ψ_low → 0. For m=0 (degenerate a≈0), the regular eigenvector -is identified by dominant |U₁| component, giving the physically correct constant-displacement -Frobenius solution and avoiding the spurious logarithmic irregularity. +For each mode, forms the 2×2 diagonal block of ψ_low · M from the diagonal entries of F̄⁻¹, K̄, and Ḡ, +discards mode coupling, and returns that block's eigenvector normalized so that U₂ = I (or U₁ = 1, U₂ = 0 +when the momentum component underflows). The component ratio does not scale as a positive power of +`psi_low`. A negative Ḡ diagonal, including the m = 0 block, gives an imaginary conjugate pair, so the +returned U₁ is non-real and which branch is stored follows the eigen-solver's column order. """ function compute_axis_init(mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal, psi_low::Float64) @@ -495,12 +529,9 @@ end Initialize the OdeState struct for the case of sing_start = 0 (axis initialization). Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, and `u`. -The Frobenius start is used only when the starting surface lies at or below `ctrl.frobenius_psi_max`; -an interior start (or `ctrl.frobenius_psi_max = 0`) uses U₁ = 0, U₂ = I. - -### TODOs - -Move ising_start logic to chunk_el_integration_bounds? +The default `ctrl.frobenius_psi_max = 0` uses the fixed start, U₁ = 0, U₂ = I, with no warning. The +Frobenius start is used only when `frobenius_psi_max` is positive and the starting surface lies at or +below it. """ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, mats::MatrixSplines, profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal) @@ -533,13 +564,8 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ma odet.ising_start = searchsortedfirst(getfield.(intr.sing, :psifac), odet.psifac) - 1 end - # The Frobenius start is a series about the axis; an interior start uses the fixed start Fortran DCON always uses. + # Default (frobenius_psi_max = 0) is the fixed start Fortran DCON uses. A positive threshold opts into the diagonal Frobenius start. fixed_start = odet.psifac > ctrl.frobenius_psi_max - if fixed_start && ctrl.frobenius_psi_max > 0 - @warn "Integration starts at ψ_N = $(round(odet.psifac; sigdigits=4)), above frobenius_psi_max = $(ctrl.frobenius_psi_max): " * - "the axis Frobenius start is only valid near the axis, so the fixed start (U₁ = 0, U₂ = I) is used, as in Fortran DCON. " * - "Set frobenius_psi_max = 0 to select it silently, or raise it to keep the Frobenius start." - end if fixed_start # Original Glasser initialization: U₁=0, U₂=I [Glasser 2016 §VI]. # Constrains the displacement ξ^ψ=0 for all modes at the starting surface (fixed axis or rigid core). @@ -547,9 +573,10 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ma odet.u[ipert, ipert, 2] = 1 end else - # Frobenius initialization [Glasser 2016 §VI Eq. 51]: selects the regular - # (non-logarithmic) solution for each mode, including the correct constant - # displacement solution for the degenerate m=0 case (free magnetic axis). + # Opt-in diagonal Frobenius start. See compute_axis_init. + @warn "Frobenius axis start at ψ_N = $(round(odet.psifac; sigdigits=4)). " * + "This start is currently unreliable - use only if explicitly requested." * + "Set frobenius_psi_max = 0 for the fixed start (U₁ = 0, U₂ = I)." U1_init, U2_init = compute_axis_init(mats, profiles, intr, odet.psifac) odet.u[:, :, 1] .= U1_init odet.u[:, :, 2] .= U2_init @@ -897,7 +924,6 @@ making it clear what region is being integrated. - `mats::MatrixSplines` - Fourier fit variables - `intr::ForceFreeStatesInternal` - Internal data - `chunk::IntegrationChunk` - Integration chunk containing start and end ψ for integration - """ function integrate_el_region!( odet::OdeState, @@ -1177,7 +1203,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` @@ -1194,7 +1220,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/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl index 8111b92e9..1fb29d2ea 100644 --- a/src/ForceFreeStates/Riccati/Driver.jl +++ b/src/ForceFreeStates/Riccati/Driver.jl @@ -76,10 +76,10 @@ After renormalization (at crossing or when norms exceed ucrit): This is compatible with downstream code (which uses U₁/U₂ ratio): - Free.jl: wp = u[:,:,2] / u[:,:,1] = I · S⁻¹ = P ✓ (post-renorm) - FixedBoundaryStability.jl: crit = min_eigval(u[:,:,1] / u[:,:,2]) = min_eigval(S) ✓ - - Axis init: determined by `ctrl.frobenius_psi_max`. When the start lies above it (or it is 0), U₁=0, U₂=I → S(ψ₀)=0 (original - Glasser fixed-axis BC). Otherwise (default near the axis), Frobenius eigenvalue init [Glasser 2016 Eq. 51] - sets U₂=I and U₁ to the regular Frobenius eigenvector per mode → S(ψ₀) = U₁_Frobenius is - nonzero in general. Riccati S-evolution remains well-defined either way. + - Axis init: determined by `ctrl.frobenius_psi_max`. The default 0, and any start above a positive threshold, + uses U₁=0, U₂=I → S(ψ₀)=0 (Fortran DCON's fixed-axis BC). A positive threshold at or above the start + opts into the diagonal Frobenius init, which sets U₂=I and a per-mode U₁, so S(ψ₀) is nonzero. + Riccati S-evolution remains well-defined either way. ## Key Differences from Standard Integration diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 3d45922df..98fca1a63 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -136,8 +136,10 @@ end @testset "truncate_chunks! cuts the Riccati chunks back to the edge" begin F = GeneralizedPerturbedEquilibrium.ForceFreeStates # Crossing chunks end just short of surfaces 1-3; the edge drops surface 3 (msing = 2). - chunks() = [F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0) - for (i, (a, b)) in enumerate(((0.10, 0.29), (0.31, 0.59), (0.61, 0.89), (0.91, 0.99)))] + chunks() = [ + F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0) + for (i, (a, b)) in enumerate(((0.10, 0.29), (0.31, 0.59), (0.61, 0.89), (0.91, 0.99))) + ] cs = chunks() @test F.truncate_chunks!(cs, 0.75, 2) == 3 # straddling chunk shortened, so re-integrate it @test length(cs) == 3 && cs[3].psi_end == 0.75 && count(c -> c.needs_crossing, cs) == 2 @@ -405,7 +407,7 @@ end @test length(odet.unorm0) == numpert_total end - @testset "interior start falls back to the fixed initialization" begin + @testset "axis start defaults to the fixed initialization" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) @@ -430,19 +432,21 @@ end odet = FFS.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) return odet, ctrl, mats, equil, intr end - # Near the axis the Frobenius start gives the regular solution: U₂ = I with a nonzero U₁. + # The default threshold is 0: near the axis the fixed start is used, with no warning. odet, ctrl, mats, equil, intr = axis_state(1e-4) - FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) - @test odet.u[:, :, 2] ≈ I - @test any(!iszero, odet.u[:, :, 1]) - # An interior start switches to the fixed start (U₁ = 0, U₂ = I) and says so. - odet, ctrl, mats, equil, intr = axis_state(0.3) - @test_logs (:warn, r"fixed start") FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + @test ctrl.frobenius_psi_max == 0 + @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I - # The threshold is a control: raising it keeps the Frobenius start, and zero selects the fixed start silently. - odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.5) + # A start above a positive threshold uses the fixed start, still with no warning. + odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.01) @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I + # Raising the threshold to cover the start opts into the Frobenius start and warns. + odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.5) + @test_logs (:warn, r"Frobenius axis start") FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + @test odet.u[:, :, 2] ≈ I @test any(!iszero, odet.u[:, :, 1]) + # An explicit zero selects the fixed start silently even near the axis. odet, ctrl, mats, equil, intr = axis_state(1e-4; frobenius_psi_max=0.0) @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I @@ -613,7 +617,7 @@ end # materialize after the Gaussian fixups and free-boundary normalization rather than # transforming stored derivatives alongside u_store. npert = intr.numpert_total - T = Matrix{ComplexF64}(I, npert, npert) .+ 0.25 .* ComplexF64.(reshape(sin.(1:npert^2), npert, npert)) + T = Matrix{ComplexF64}(I, npert, npert) .+ 0.25 .* ComplexF64.(reshape(sin.(1:(npert^2)), npert, npert)) odet_t = deepcopy(odet_pristine) for istep in 1:odet_t.step odet_t.u_store[:, :, 1, istep] = odet_t.u_store[:, :, 1, istep] * T