diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index 04c4aba22..344d643e0 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -20,6 +20,12 @@ control surface (when `vac_flag`) and produces no free-boundary energies or ODE Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution that `PerturbedEquilibrium` consumes in place of a forward solution. +`gal_solver` picks the banded solver. `"auto"` (the default) tries a Hermitian Cholesky +factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved +to the right-hand side; the diagonal's imaginary part is ignored) and falls back to Jacobi-scaled +LU when the matrix is not positive definite. `"LU"` skips the Cholesky attempt; +`"cholesky"` stops with an error instead of falling back. + The implementation lives in `src/ForceFreeStates/Galerkin/`: | File | Content | diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index 71d16afcb..890137170 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -57,7 +57,7 @@ dmlim = 0.2 # Truncate integration at (last_rational_q + dmli # Outer-region singular Galerkin Δ′ solver with rpec coil columns. # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking. -gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_solver = "auto" # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky" gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces gal_nq = 6 # Gauss-Lobatto quadrature order per element gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index 71e584905..4249bbe12 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -57,7 +57,7 @@ dmlim = 0.2 # Truncate integration at (last_rational_q + dmli # Outer-region singular Galerkin Δ′ solver with rpec coil columns. # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking. -gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_solver = "auto" # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky" gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces gal_nq = 6 # Gauss-Lobatto quadrature order per element gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index 0706863b7..6b4a808cf 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -57,7 +57,7 @@ force_termination = true # Terminate after the stability stage (no perturb write_outputs_to_HDF5 = true # Write stability outputs to HDF5 # Outer-region singular Galerkin Delta-prime solver with rpec coil columns. -gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_solver = "LU" # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky" gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces gal_nq = 6 # Gauss-Lobatto quadrature order per element gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 0571301c7..94106ee0d 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -58,7 +58,7 @@ force_termination = true # Terminate after the stability stage (no perturb write_outputs_to_HDF5 = true # Write stability outputs to HDF5 # Outer-region singular Galerkin Delta-prime solver with rpec coil columns. -gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_solver = "LU" # Banded solver: "auto" (Cholesky if Hermitian positive definite, else LU), "LU" or "cholesky" gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces gal_nq = 6 # Gauss-Lobatto quadrature order per element gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 4dde78771..9c31d151c 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -153,7 +153,7 @@ gpec.toml. - `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. - - `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. + - `extended_precision_bvp::Bool` - Solve the Δ' BVP and its PEST3 combination in `Complex{Double64}` (~31 digits). Default `false` (Float64, 1.5–5.4× faster on large matrices); set `true` to check Float64 at small `singfac_min`, where the four-term combination loses digits. """ @kwdef struct ForceFreeStatesControl verbose::Bool = true @@ -191,10 +191,10 @@ gpec.toml. force_termination::Bool = false integrator::String = "riccati" # Default: unlocks SingularSurfaces/Delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. Use "forward" for dense ξ (PerturbedEquilibrium) or kinetic runs. nchunks::Int = 0 # Riccati chunk-count target; 0 = auto (derived from msing alone, never from Threads.nthreads()). - extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). + extended_precision_bvp::Bool = false # Δ' BVP in Complex{Double64}; opt-in check of Float64 at small singfac_min. # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port); selected by integrator = "galerkin" --- - gal_solver::String = "LU" # "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs) + gal_solver::String = "auto" # "auto" (banded Cholesky if the matrix is positive definite, else LU), "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs, error if not positive definite) gal_nx::Int = 256 # elements per interval between singular surfaces gal_nq::Int = 6 # Gauss-Lobatto quadrature order per element gal_pfac::Float64 = 0.001 # grid packing ratio near singular surfaces diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 2224c53e2..339d5b919 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -769,7 +769,7 @@ function cross_ideal_singular_surf!( # Compute direction-specific asymptotic power series for this singular surface singp = intr.sing[ising] sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) - sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha) + sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0) dpsi = singp.psifac - odet.psifac # ψ_res - ψ (positive) # Get asymptotic coefficients before crossing (left side) diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl index 1ad84cb8b..593a15ef2 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl @@ -6,7 +6,7 @@ # gal_gauss_quad (gal.f) nonresonant stiffness via Gauss-Lobatto quadrature # gal_extension (gal.f) extension-cell matrix/rhs (big-solution driving + surface terms) # gal_resonant (gal.f) resonant-cell integral (QuadGK in place of Fortran LSODE) -# gal_assemble_mat (gal.f) scatter cell blocks into the banded global matrix (LU/Cholesky) +# gal_assemble_mat (gal.f) scatter cell blocks into the full-band global matrix # gal_assemble_rhs (gal.f) scatter cell rhs into the global RHS # gal_set_boundary (gal.f) axis identity + edge (vacuum or fixed) BCs # @@ -359,15 +359,13 @@ end """ gal_assemble_mat!(ws, mpert) -Scatter every cell's `mat`/`emat`/`ediag` into the banded global matrix `ws.mat`. Port of -`gal_assemble_mat` (gal.f). Supports both `"LU"` (full band, `offset = kl+ku+1`) and -`"cholesky"` (lower band only, `offset = 1`). +Scatter every cell's `mat`/`emat`/`ediag` into the full-band global matrix `ws.mat` +(`offset = kl+ku+1`). Port of `gal_assemble_mat` (gal.f). """ function gal_assemble_mat!(ws::GalWorkspace, mpert::Int) np = ws.np msing = length(ws.intvl) - 1 - lu = ws.solver == "LU" - offset = lu ? ws.kl + ws.ku + 1 : 1 + offset = ws.kl + ws.ku + 1 for ising in 0:msing for cell in ws.intvl[ising+1].cells @@ -377,24 +375,14 @@ function gal_assemble_mat!(ws::GalWorkspace, mpert::Int) # nonresonant stiffness for jp in 0:np, jpert in 1:mpert j = cell.map[jpert, jp+1] - if lu - ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1] - elseif j <= i - ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1] - end + ws.mat[offset+i-j, j] += cell.mat[ipert, jpert, ip+1, jp+1] end # resonant coupling (res/ext only) if res j = cell.emap e = cell.emat[ipert, ip+1] - if lu - ws.mat[offset+i-j, j] += e - ws.mat[offset+j-i, i] += conj(e) - elseif j < i - ws.mat[offset+i-j, j] += e - elseif j > i - ws.mat[offset+j-i, i] += conj(e) - end + ws.mat[offset+i-j, j] += e + ws.mat[offset+j-i, i] += conj(e) end end if res @@ -452,14 +440,12 @@ including its full three-way edge branch (rpec count `ncoil = ws.nsol - 2*msing` `wv_edge` is `nothing` for the rpec and fixed cases. """ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing,Matrix{ComplexF64}}) - chol = ws.solver == "cholesky" msing = length(ws.intvl) - 1 ncoil = ws.nsol - 2 * msing # axis u(0): first cell of interval 0, Hermite value DOF ip=0 (Julia index 1) cell = ws.intvl[1].cells[1] cell.mat[:, :, 1, :] .= 0 - chol && (cell.mat[:, :, :, 1] .= 0) for idx in 1:mpert cell.mat[idx, idx, 1, 1] = 1 end @@ -469,7 +455,6 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing, if ncoil > 0 # rpec: identity edge + unit coil sources (gal.f) cell.mat[:, :, 4, :] .= 0 - chol && (cell.mat[:, :, :, 4] .= 0) for idx in 1:mpert cell.mat[idx, idx, 4, 4] = 1 end @@ -482,7 +467,6 @@ function gal_set_boundary!(ws::GalWorkspace, mpert::Int, wv_edge::Union{Nothing, else # fixed boundary: identity edge (gal.f) cell.mat[:, :, 4, :] .= 0 - chol && (cell.mat[:, :, :, 4] .= 0) for idx in 1:mpert cell.mat[idx, idx, 4, 4] = 1 end diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index a06e37e82..e25765a54 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -39,7 +39,9 @@ function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil, return ws end -"""Empty `GalerkinResult` for a domain with no resonant surfaces.""" +""" +Empty `GalerkinResult` for a domain with no resonant surfaces. +""" function empty_galerkin_result() return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing, nothing) end @@ -60,7 +62,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines intr::ForceFreeStatesInternal; wv=nothing) intr.npert == 1 || error("galerkin_solve: only single-n (npert == 1) is supported") - ctrl.gal_solver in ("LU", "cholesky") || error("galerkin_solve: gal_solver must be \"LU\" or \"cholesky\"") + ctrl.gal_solver in ("auto", "LU", "cholesky") || + error("galerkin_solve: gal_solver must be \"auto\", \"LU\" or \"cholesky\"") nn = intr.nlow mpert = intr.numpert_total np = GAL_NP @@ -74,22 +77,20 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines ctrl.verbose && @info "Starting outer-region Galerkin Δ′ solve (msing=$msing, solver=$(ctrl.gal_solver))" - # Per-surface two-sided asymptotic series matching Fortran sing.f vmatr/vmatl: - # right = sig=+1, left = sig=-1, no √det normalization. The Mercier exponent α is a - # property of the surface, so the left series reuses the right's α (alpha_override). The order is - # raised to gal_sing_order + ceil(2·Re(α)) for high-Mercier-index surfaces (Fortran sing1_vmat). + # Per-surface two-sided asymptotic series (sig = ±1). Each side uses its own resonant block. asymps = GalSingAsymp[] for s in sings sing_order = ctrl.gal_sing_order ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order) if ctrl.gal_sing_order_ceiling + # The order is raised for high-Mercier-index surfaces order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1])) if order > ctrl.gal_sing_order sing_order = order ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order) end end - al = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order) + al = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=-1.0, sing_order=sing_order) push!(asymps, GalSingAsymp(ar, al)) end @@ -97,14 +98,10 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines nx = ctrl.gal_nx kl = mpert * (np + 1) ku = kl - ldab = ctrl.gal_solver == "LU" ? 2kl + ku + 1 : kl + 1 + ldab = 2kl + ku + 1 # rpec_flag (RDCON, gal.f): append mpert coil-response columns. Each is a unit source at # the edge value DOF in one poloidal mode; the recorded plasma response is the coil block of Δ_gw. - # Cholesky's lower-only edge zeroing can't represent the rpec identity edge, so require LU. ncoil = ctrl.gal_rpec_flag ? mpert : 0 - if ncoil > 0 && ctrl.gal_solver != "LU" - error("galerkin_solve: gal_rpec_flag=true requires gal_solver=\"LU\" (coil edge BC needs the full-band path)") - end nsol = 2 * msing + ncoil intvl = [GalInterval(zeros(Float64, nx + 1), zeros(Float64, nx + 1), [GalCell(mpert) for _ in 1:nx]) for _ in 0:msing] @@ -131,7 +128,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) if ctrl.verbose - offdbg = ws.solver == "LU" ? ws.kl + ws.ku + 1 : 1 + offdbg = ws.kl + ws.ku + 1 @info "Galerkin assembly: ndim=$(ws.ndim) nsol=$nsol kl=$(ws.kl) ldab=$(ws.ldab) |rhs|=$(norm(ws.rhs))" for ising in 0:msing for cell in ws.intvl[ising+1].cells @@ -144,14 +141,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines # --- banded solve (gal.f) --- ws.sol .= ws.rhs - if ws.solver == "LU" - ctrl.verbose && @info "Galerkin LU banded factorization + solve" - gal_scaled_lu_solve!(ws.mat, ws.sol, ws.kl, ws.ku) - else - ctrl.verbose && @info "Galerkin Cholesky banded factorization + solve" - LinearAlgebra.LAPACK.pbtrf!('L', ws.kl, ws.mat) - LinearAlgebra.LAPACK.pbtrs!('L', ws.kl, ws.mat, ws.sol) - end + ws.solver = gal_banded_solve!(ws.mat, ws.sol, ws.kl, ws.ku, ctrl.gal_solver) + ctrl.verbose && @info "Galerkin banded solve: $(ws.solver)" # --- extract Δ′ from the small resonant coefficients (gal.f) --- delta = zeros(ComplexF64, nsol, 2 * msing) @@ -167,8 +158,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines Ap, Bp, Gammap, Deltap = gal_pest3_blocks(delta, msing) - # Mercier index D_I = -Re(α²); α is the small/large solution exponent (a surface property, taken - # from the right series — the left series reuses the same α via alpha_override above). + # Mercier index D_I = -Re(α²). α is the half-separation of the resonant exponents. di = [real(-asymps[i].right.alpha[1]^2) for i in 1:msing] alpha = [asymps[i].right.alpha[1] for i in 1:msing] @@ -215,7 +205,7 @@ function gal_scaled_lu_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{Com n = size(mat, 2) off = kl + ku + 1 # band row holding the diagonal d = [iszero(mat[off, j]) ? 1.0 : 1 / sqrt(abs(mat[off, j])) for j in 1:n] - for j in 1:n, i in max(1, j - ku):min(n, j + kl) + for j in 1:n, i in max(1, j-ku):min(n, j+kl) mat[off+i-j, j] *= d[i] * d[j] end sol .*= d @@ -320,11 +310,19 @@ const GALERKIN_H5_ANNOTATIONS = [ "ForceFreeStates/Solutions/GalerkinIntegration/q" => (; long_name="safety factor q on the Galerkin solution grid", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)), "ForceFreeStates/Solutions/GalerkinIntegration/xi_psi" => - (; long_name="closed axis-to-edge ξ^ψ profiles (ideal or inner-layer closure; identity-at-edge basis)", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)), + (; + long_name="closed axis-to-edge ξ^ψ profiles (ideal or inner-layer closure; identity-at-edge basis)", + dims=("mode", "solution", "psi"), + attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",) + ), "ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi" => (; long_name="analytic ψ_N derivative of the closed ξ^ψ profiles", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)), "ForceFreeStates/Solutions/GalerkinIntegration/xi_s" => - (; long_name="surface displacement ξ_s from the outer ideal-MHD relation", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)), + (; + long_name="surface displacement ξ_s from the outer ideal-MHD relation", + dims=("mode", "solution", "psi"), + attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",) + ), "ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi" => (; long_name="normalized poloidal flux ψ_N grid of the raw Galerkin basis", scale="psi_gal_basis"), "ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational" => (; long_name="flag: grid node lies on a rational surface", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi",)), @@ -361,3 +359,71 @@ function annotate_galerkin!(out_h5) Utilities.HDF5Annotations.annotate!(out_h5, GALERKIN_H5_ANNOTATIONS) return nothing end + +""" + gal_banded_solve!(mat, sol, kl, ku, solver) -> String + +Solve the full-band system `mat` (LAPACK `gbtrf` storage) in place of `sol` and return the +solver used. `"auto"` tries `gal_cholesky_solve!` and falls back to `gal_scaled_lu_solve!` when +the matrix is not positive definite; `"cholesky"` errors instead of falling back. +""" +function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String) + if solver != "LU" + gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky" + solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not positive definite; use \"auto\" or \"LU\"") + end + gal_scaled_lu_solve!(mat, sol, kl, ku) + return "LU" +end + +""" + gal_cholesky_solve!(mat, sol, kl, ku) -> Bool + +Banded Cholesky solve of the full-band system `mat`, leaving `mat` untouched. Rows that are +identity rows (Dirichlet DOFs: axis, fixed or coil edge) have their columns moved to the +right-hand side, which keeps the system Hermitian. Only the lower band and the real part of the +diagonal are used. Returns `false`, with `sol` unchanged, when the matrix is not positive definite. +""" +function gal_cholesky_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int) + kl == ku || return false + n = size(mat, 2) + off = 2kl + 1 + a(i, j) = mat[off+i-j, j] + band(i) = max(1, i - kl):min(n, i + kl) + fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n] + ab = zeros(ComplexF64, kl + 1, n) + for j in 1:n, i in j:min(n, j + kl) + ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j) + end + gal_zpbtrf!(ab, kl) || return false + for c in 1:n + fixed[c] || continue + for r in band(c) + r == c || fixed[r] || (@views sol[r, :] .-= a(r, c) .* sol[c, :]) + end + end + gal_zpbtrs!(ab, kl, sol) + return true +end + +# Lower-band Hermitian Cholesky as rdcon gal.f: zpbtrf('L') then zpbtrs('L'). LinearAlgebra.LAPACK +# has no pbtrf!/pbtrs! wrapper, so call LAPACK directly. Returns false if not positive definite. +function gal_zpbtrf!(AB::Matrix{ComplexF64}, kd::Int) + info = Ref{LinearAlgebra.BlasInt}() + ccall((LinearAlgebra.BLAS.@blasfunc(zpbtrf_), LinearAlgebra.libblastrampoline), Cvoid, + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, + Ref{LinearAlgebra.BlasInt}, Clong), + 'L', size(AB, 2), kd, AB, max(1, stride(AB, 2)), info, 1) + LAPACK.chklapackerror(min(info[], 0)) + return info[] == 0 +end + +function gal_zpbtrs!(AB::Matrix{ComplexF64}, kd::Int, B::AbstractVecOrMat{ComplexF64}) + info = Ref{LinearAlgebra.BlasInt}() + ccall((LinearAlgebra.BLAS.@blasfunc(zpbtrs_), LinearAlgebra.libblastrampoline), Cvoid, + (Ref{UInt8}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, + Ref{LinearAlgebra.BlasInt}, Ptr{ComplexF64}, Ref{LinearAlgebra.BlasInt}, Ref{LinearAlgebra.BlasInt}, Clong), + 'L', size(AB, 2), kd, size(B, 2), AB, max(1, stride(AB, 2)), B, max(1, stride(B, 2)), info, 1) + LAPACK.chklapackerror(info[]) + return B +end diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl index f4e79deb0..a1de3187a 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl @@ -84,10 +84,8 @@ end Global workspace for the Galerkin solve. Port of Fortran `gal_type` (gal.f). `intvl` has length `msing+1` (Fortran `intvl(0:msing)`, accessed as `intvl[ising+1]`). -`kl = ku = mpert*(np+1)`. The banded `mat` storage layout depends on `solver` (gal.f): - - - `"LU"` → `ldab = 2*kl + ku + 1`, full band (LAPACK `gbtrf!`/`gbtrs!`) - - `"cholesky"` → `ldab = kl + 1`, lower band only (LAPACK `pbtrf!`/`pbtrs!`, `'L'`) +`kl = ku = mpert*(np+1)`. `mat` is stored full band, `ldab = 2*kl + ku + 1` (LAPACK `gbtrf!` +layout); `solver` records the factorization `gal_banded_solve!` used (`"cholesky"` or `"LU"`). """ mutable struct GalWorkspace solver::String diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl index 8b123b2d0..05a79f81d 100644 --- a/src/ForceFreeStates/Integrators.jl +++ b/src/ForceFreeStates/Integrators.jl @@ -37,7 +37,7 @@ and never depends on the thread count. Maps onto `integrator = "riccati"`. end """ - Galerkin(; solver="LU", nx=256, ...) + Galerkin(; solver="auto", nx=256, ...) RDCON outer-region singular Galerkin solver: solves the same Euler-Lagrange system variationally on a finite-element grid packed around the rational surfaces, producing Δ′ @@ -46,7 +46,7 @@ matching `gal_*` control key without the prefix. ## Fields - - `solver::String` - Banded linear solver, `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (zpbtrf/zpbtrs). Matching requires `"LU"`. + - `solver::String` - Banded linear solver: `"auto"` (Cholesky, zpbtrf/zpbtrs, when the matrix is positive definite, else LU), `"LU"` (zgbtrf/zgbtrs) or `"cholesky"` (error if not positive definite). - `nx::Int` - Elements per interval between singular surfaces. - `nq::Int` - Gauss-Lobatto quadrature order per element. - `pfac::Float64` - Grid packing ratio near singular surfaces. @@ -63,7 +63,7 @@ matching `gal_*` control key without the prefix. - `edge_onesided::Bool` - Pack the two end intervals one-sided toward their single rational end instead of the Fortran symmetric pack. """ @kwdef struct Galerkin <: AbstractIntegrator - solver::String = "LU" + solver::String = "auto" nx::Int = 256 nq::Int = 6 pfac::Float64 = 0.001 diff --git a/src/ForceFreeStates/Riccati/Crossings.jl b/src/ForceFreeStates/Riccati/Crossings.jl index 4962ee240..e4b80fd56 100644 --- a/src/ForceFreeStates/Riccati/Crossings.jl +++ b/src/ForceFreeStates/Riccati/Crossings.jl @@ -6,14 +6,14 @@ Cross a singular surface for the Riccati formulation. Replaces `cross_ideal_singular_surf!` for the Riccati integration path with two key differences: -1. **No Gaussian reduction**: `cross_ideal_singular_surf!` calls `compute_solution_norms!` - which applies Gaussian reduction to (S, I). This divides by pivot elements of S, which - can be near-zero (S = 0 at axis and grows slowly), producing NaN/Inf in U₂. For Riccati, - S is bounded so Gaussian reduction is unnecessary. + 1. **No Gaussian reduction**: `cross_ideal_singular_surf!` calls `compute_solution_norms!` + which applies Gaussian reduction to (S, I). This divides by pivot elements of S, which + can be near-zero (S = 0 at axis and grows slowly), producing NaN/Inf in U₂. For Riccati, + S is bounded so Gaussian reduction is unnecessary. -2. **Direct column zeroing**: Instead of using the GR-sorted `odet.index` to identify the - column to zero, we use `ipert_res` directly (the resonant mode index). This is valid since - without GR there is no permutation applied to the columns of S. + 2. **Direct column zeroing**: Instead of using the GR-sorted `odet.index` to identify the + column to zero, we use `ipert_res` directly (the resonant mode index). This is valid since + without GR there is no permutation applied to the columns of S. **Δ' normalization**: This function expects `odet.u` in the bounded (U₁, U₂) form produced by `riccati_integrate_chunk!` with `needs_crossing=true` (final renorm skipped). ca_l is computed @@ -54,16 +54,13 @@ end """ _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) -> (left, right) -Compute left- (`sig=-1`) and right- (`sig=+1`) side singular asymptotics matching -Fortran STRIDE's separate vmatl/vmatr (sing_vmat). Alpha is taken from the right -side and shared with the left. +Left (`sig=-1`) and right (`sig=+1`) singular asymptotics, each from its own resonant block. """ function _two_sided_singular_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal) + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, + intr::ForceFreeStatesInternal) sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) - sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, - alpha_override=sing_asymp_right.alpha) + sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0) return sing_asymp_left, sing_asymp_right end @@ -83,7 +80,7 @@ end # Capture left-side asymptotic data into odet.ca_l and singp.ua_left/psi_ua_left. function _capture_left_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_left, - dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int) + dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int) ua = sing_get_ua(sing_asymp_left, dpsi) singp.ua_left = copy(ua) singp.psi_ua_left = odet.psifac @@ -95,9 +92,9 @@ end # odet.psifac to the right side. The zeroed columns stay zero through the predictor # since du[:, ipert_res, :] = 0 when u[:, ipert_res, :] = 0. function _predict_across_singular_surface!(odet::OdeState, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal, ising::Int, - ipert_res, dpsi::Float64, sing_asymp_right) + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, + intr::ForceFreeStatesInternal, ising::Int, + ipert_res, dpsi::Float64, sing_asymp_right) if ctrl.kinetic_factor == 0 for i in eachindex(sing_asymp_right.r1) odet.u[:, ipert_res[i], :] .= 0 @@ -117,8 +114,8 @@ end # Column ipert_res of [U₁_new; U₂_new] = ua[:, ipert_res+N, :] (the introduced small asymptotic), # so ca_r[ipert_res, ipert_res, 2] = 1 regardless of other columns' normalization. function _capture_right_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_right, - dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int, - ipert_res, ctrl::ForceFreeStatesControl) + dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int, + ipert_res, ctrl::ForceFreeStatesControl) ua = sing_get_ua(sing_asymp_right, dpsi) singp.ua_right = copy(ua) singp.psi_ua_right = odet.psifac @@ -135,9 +132,9 @@ end # STUB: per-surface ca-based Δ' (not physically valid; see SingType.delta_prime docstring). # The canonical Δ' is intr.delta_prime_matrix from compute_delta_prime_matrix!. function _stash_per_surface_delta_prime_stub!(odet::OdeState, intr::ForceFreeStatesInternal, - ising::Int, ipert_res, sing_asymp_right, - equil::Equilibrium.PlasmaEquilibrium, - ctrl::ForceFreeStatesControl) + ising::Int, ipert_res, sing_asymp_right, + equil::Equilibrium.PlasmaEquilibrium, + ctrl::ForceFreeStatesControl) ctrl.kinetic_factor == 0 || return denom = (2π)^2 * equil.psio n_res = length(sing_asymp_right.r1) diff --git a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl index 0a955fde7..c878628dd 100644 --- a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl +++ b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl @@ -7,6 +7,7 @@ Compute the inter-surface tearing stability matrix (msing × msing) using the STRIDE global BVP formulation [Glasser 2018 Phys. Plasmas 25, 032501, Sec. III.B]. The BVP encodes the full plasma response with unknowns at each surface boundary: + ``` x_axis (N): free IC parameters at the axis (U₁ = 0 regular solutions) x_left[j] (2N): state at left inner-layer boundary of surface j @@ -19,9 +20,11 @@ The BVP encodes the full plasma response with unknowns at each surface boundary: When `wv` is provided (the vacuum response matrix, singfac-scaled), the edge BC follows the Fortran STRIDE convention: + ``` U₁ = c, U₂ = -wv·ψ₀²·c ``` + which is the free-boundary condition `wp + wv = 0` at the edge. When `wv` is `nothing`, a conducting wall BC (`U₁ = 0`) is used. @@ -38,9 +41,11 @@ keeps the BVP matrix full-rank and well-conditioned. The raw BVP solution is a 2·msing × 2·msing matrix `dp` with left/right sub-indices at each surface. The PEST3-convention Δ' matrix is the linear combination [Chance, PPPL-2527]: + ``` deltap(i,j) = dp(2i,2j) - dp(2i,2j-1) - dp(2i-1,2j) + dp(2i-1,2j-1) ``` + stored in `intr.delta_prime_matrix` (msing × msing). ## Limitations @@ -65,13 +70,13 @@ function compute_delta_prime_matrix!( intr::ForceFreeStatesInternal, propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}; - wv::Union{Nothing,Matrix{ComplexF64}} = nothing, - psio::Float64 = 0.0, - debug::Bool = false, - S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}} = nothing, - ctrl::Union{Nothing,ForceFreeStatesControl} = nothing, - equil::Union{Nothing,Equilibrium.PlasmaEquilibrium} = nothing, - mats::Union{Nothing,MatrixSplines} = nothing + wv::Union{Nothing,Matrix{ComplexF64}}=nothing, + psio::Float64=0.0, + debug::Bool=false, + S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}}=nothing, + ctrl::Union{Nothing,ForceFreeStatesControl}=nothing, + equil::Union{Nothing,Equilibrium.PlasmaEquilibrium}=nothing, + mats::Union{Nothing,MatrixSplines}=nothing ) intr.msing == 0 && return _has_unsupported_multi_resonance(intr) && return @@ -91,8 +96,10 @@ function compute_delta_prime_matrix!( if !use_S_axis for ic in i_crossings chunks[ic].direction == 1 || - error("compute_delta_prime_matrix!: FM-axis fallback (use_S_axis=false) requires forward crossing chunks; " * - "chunk $ic has direction=$(chunks[ic].direction). Either provide S_at_surface_left or use bidirectional=false.") + error( + "compute_delta_prime_matrix!: FM-axis fallback (use_S_axis=false) requires forward crossing chunks; " * + "chunk $ic has direction=$(chunks[ic].direction). Either provide S_at_surface_left or use bidirectional=false." + ) end end @@ -105,15 +112,15 @@ function compute_delta_prime_matrix!( _build_asymptotic_basis_matrices(sing, has_ua, N, msing) debug && _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, - Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) + Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) if use_S_axis uShootR, uShootL, uAxis = _build_S_axis_shooting_propagators( propagators, chunks, i_crossings, sing, msing, N, T_left_mats, T_right_mats, has_ua, ctrl, equil, mats, intr, debug) debug && _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, - S_at_surface_left, T_left_mats, - ipert_all, has_ua, msing, N) + S_at_surface_left, T_left_mats, + ipert_all, has_ua, msing, N) M, nMat, col_edge = _assemble_bvp_S_axis( uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio) else @@ -140,13 +147,13 @@ function compute_delta_prime_matrix!( # The raw matrix is consumed by `pest3_decompose` to recover (A', B', Γ', Δ') for the full # det(D' − D(γ)) = 0 eigenvalue problem; see the `delta_prime_raw` docstring in CoreTypes.jl. intr.delta_prime_matrix = deltap - intr.delta_prime_raw = dp_raw_persisted + intr.delta_prime_raw = dp_raw_persisted end # Column index helpers for the BVP matrix. j is the 1-based singular-surface index, # N is numpert_total. Layout: c_axis(N), c_left[1](2N), c_right[1](2N), ..., c_edge(N). -_col_left(j::Int, N::Int) = (N + 4N*(j-1) + 1):(N + 4N*(j-1) + 2N) -_col_right(j::Int, N::Int) = (N + 4N*(j-1) + 2N + 1):(N + 4N*j) +_col_left(j::Int, N::Int) = (N+4N*(j-1)+1):(N+4N*(j-1)+2N) +_col_right(j::Int, N::Int) = (N+4N*(j-1)+2N+1):(N+4N*j) # Multi-resonance surfaces (one q value satisfying multiple (m,n) tuples in a multi-n run) # are not yet handled by the inter-surface BVP. Returns true if any surface has >1 modes; @@ -185,28 +192,28 @@ end # Midpoint splitting halves each inter-surface span's condition number — STRIDE's trick: # cond(full) = 10¹⁵ → cond(half) ≈ 10⁷·⁵, an 8-digit accuracy gain. function _assemble_segment_propagators(propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, - i_crossings::Vector{Int}, msing::Int, N::Int, - use_S_axis::Bool) + chunks::Vector{IntegrationChunk}, + i_crossings::Vector{Int}, msing::Int, N::Int, + use_S_axis::Bool) Phi_L_mats = [assemble_fm_matrix(propagators, i_crossings[j]:i_crossings[j]) for j in 1:msing] Phi_R_mats = Vector{Matrix{ComplexF64}}(undef, msing + 1) if !use_S_axis - Phi_R_mats[1] = assemble_fm_matrix(propagators, 1:i_crossings[1]-1; condition=true) + Phi_R_mats[1] = assemble_fm_matrix(propagators, 1:(i_crossings[1]-1); condition=true) end for j in 2:msing - Phi_R_mats[j] = assemble_fm_matrix(propagators, i_crossings[j-1]+1:i_crossings[j]-1) + Phi_R_mats[j] = assemble_fm_matrix(propagators, (i_crossings[j-1]+1):(i_crossings[j]-1)) end - Phi_R_mats[msing+1] = assemble_fm_matrix(propagators, i_crossings[msing]+1:length(chunks)) + Phi_R_mats[msing+1] = assemble_fm_matrix(propagators, (i_crossings[msing]+1):length(chunks)) Phi_R_halves = Vector{Tuple{Matrix{ComplexF64},Matrix{ComplexF64}}}(undef, msing - 1) - for j in 1:msing-1 + for j in 1:(msing-1) chunk_start = i_crossings[j] + 1 - chunk_end = i_crossings[j+1] - 1 - n_chunks = chunk_end - chunk_start + 1 + chunk_end = i_crossings[j+1] - 1 + n_chunks = chunk_end - chunk_start + 1 if n_chunks >= 2 i_mid = chunk_start + div(n_chunks, 2) - 1 - Phi_left_half = assemble_fm_matrix(propagators, chunk_start:i_mid) - Phi_right_half = assemble_fm_matrix(propagators, i_mid+1:chunk_end) + Phi_left_half = assemble_fm_matrix(propagators, chunk_start:i_mid) + Phi_right_half = assemble_fm_matrix(propagators, (i_mid+1):chunk_end) Phi_R_halves[j] = (Phi_left_half, Phi_right_half) else Phi_R_halves[j] = (Matrix{ComplexF64}(I, 2N, 2N), Phi_R_mats[j+1]) @@ -220,17 +227,17 @@ end # N+1:2N = small solutions (z^{+α}, bounded). Fortran STRIDE bakes T into the shooting # propagators (uFM_sing_init); we multiply T into the BVP propagator blocks at each surface. function _build_asymptotic_basis_matrices(sing::Vector{SingType}, has_ua::Bool, N::Int, msing::Int) - T_left_mats = Vector{Matrix{ComplexF64}}(undef, msing) + T_left_mats = Vector{Matrix{ComplexF64}}(undef, msing) T_right_mats = Vector{Matrix{ComplexF64}}(undef, msing) - T_left_inv = Vector{Matrix{ComplexF64}}(undef, msing) - T_right_inv = Vector{Matrix{ComplexF64}}(undef, msing) + T_left_inv = Vector{Matrix{ComplexF64}}(undef, msing) + T_right_inv = Vector{Matrix{ComplexF64}}(undef, msing) if has_ua for j in 1:msing sp = sing[j] - T_left_mats[j] = [sp.ua_left[:,:,1]; sp.ua_left[:,:,2]] - T_right_mats[j] = [sp.ua_right[:,:,1]; sp.ua_right[:,:,2]] - T_left_inv[j] = inv(T_left_mats[j]) - T_right_inv[j] = inv(T_right_mats[j]) + T_left_mats[j] = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]] + T_right_mats[j] = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]] + T_left_inv[j] = inv(T_left_mats[j]) + T_right_inv[j] = inv(T_right_mats[j]) end end return T_left_mats, T_right_mats, T_left_inv, T_right_inv @@ -261,7 +268,7 @@ function _build_S_axis_shooting_propagators( end if can_reintegrate && !isempty(shoot_range_R) uShootR[j] = integrate_fm_with_ua_ic(chunks, shoot_range_R, sing[j].ua_right, - ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right) + ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right) else T_init = has_ua ? T_right_mats[j] : nothing uShootR[j] = assemble_fm_matrix(propagators, shoot_range_R; T_init=T_init) @@ -278,7 +285,7 @@ function _build_S_axis_shooting_propagators( end if can_reintegrate && !isempty(shoot_range_L) uShootL[j] = integrate_fm_with_ua_ic(chunks, shoot_range_L, sing[j].ua_left, - ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left) + ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left) else T_init = has_ua ? T_left_mats[j] : nothing uShootL[j] = assemble_fm_matrix(propagators, shoot_range_L; T_init=T_init) @@ -287,10 +294,10 @@ function _build_S_axis_shooting_propagators( uAxis, i_axis_mid = _build_conditioned_axis_propagator(propagators, i_crossings, N) uShootL[1] = _build_uShootL_first(propagators, chunks, i_crossings, sing, - T_left_mats, has_ua, can_reintegrate, i_axis_mid, - ctrl, equil, mats, intr, N) + T_left_mats, has_ua, can_reintegrate, i_axis_mid, + ctrl, equil, mats, intr, N) if debug - shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1) + shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1) @info " Axis propagator: $(i_axis_mid) chunks, cond=$(@sprintf("%.2e", cond(uAxis)))" @info " uShootL[1]: range=$(shoot_range_L1), cond=$(@sprintf("%.2e", cond(uShootL[1])))" end @@ -303,25 +310,25 @@ end # chunk+1 to chunk(i_crossings[j]-1). The ψ midpoint is used (not the chunk-index midpoint) # because chunks near singularities are packed tighter in ψ — Fortran convention. function _midpoint_shoot_range(chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, - j::Int, msing::Int; side::Symbol) + j::Int, msing::Int; side::Symbol) if side === :right - j == msing && return (i_crossings[msing] + 1):length(chunks) + j == msing && return (i_crossings[msing]+1):length(chunks) chunk_start = i_crossings[j] + 1 - chunk_end = i_crossings[j+1] - 1 + chunk_end = i_crossings[j+1] - 1 else # :left, j >= 2 chunk_start = i_crossings[j-1] + 1 - chunk_end = i_crossings[j] - 1 + chunk_end = i_crossings[j] - 1 end psi_mid_target = (chunks[chunk_start].psi_start + chunks[chunk_end].psi_end) / 2 i_mid_inter = chunk_start - for ic in chunk_start:chunk_end-1 + for ic in chunk_start:(chunk_end-1) if chunks[ic].psi_end >= psi_mid_target i_mid_inter = ic break end i_mid_inter = ic end - return side === :right ? (chunk_start:i_mid_inter) : ((i_mid_inter + 1):chunk_end) + return side === :right ? (chunk_start:i_mid_inter) : ((i_mid_inter+1):chunk_end) end # Build a well-conditioned axis propagator by forward-propagating [0; I] through the @@ -329,7 +336,7 @@ end # midpoint is placed one chunk before the first surface so that uShootL[1] covers only the # last chunk, keeping it well-conditioned. function _build_conditioned_axis_propagator(propagators::Vector{ChunkPropagator}, - i_crossings::Vector{Int}, N::Int) + i_crossings::Vector{Int}, N::Int) n_pre_cross = i_crossings[1] - 1 i_axis_mid = max(1, n_pre_cross - 1) uAxis = zeros(ComplexF64, 2N, N) @@ -339,9 +346,9 @@ function _build_conditioned_axis_propagator(propagators::Vector{ChunkPropagator} for ic in 1:i_axis_mid prop = propagators[ic] upper_old = uAxis[1:N, :] - lower_old = uAxis[N+1:2N, :] - uAxis[1:N, :] .= prop.block_upper_ic[:,:,1] * upper_old .+ prop.block_lower_ic[:,:,1] * lower_old - uAxis[N+1:2N, :] .= prop.block_upper_ic[:,:,2] * upper_old .+ prop.block_lower_ic[:,:,2] * lower_old + lower_old = uAxis[(N+1):2N, :] + uAxis[1:N, :] .= prop.block_upper_ic[:, :, 1] * upper_old .+ prop.block_lower_ic[:, :, 1] * lower_old + uAxis[(N+1):2N, :] .= prop.block_upper_ic[:, :, 2] * upper_old .+ prop.block_lower_ic[:, :, 2] * lower_old Q, _ = qr(uAxis) uAxis .= Matrix(Q)[:, 1:N] end @@ -355,18 +362,18 @@ end # Falls back to T_left_mats[1] (or identity if no ua) when there's only 1 chunk before the # first crossing. function _build_uShootL_first(propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, - sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}}, - has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int, - ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int) - shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1) + chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, + sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}}, + has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int, + ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int) + shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1) if can_reintegrate && !isempty(shoot_range_L1) return integrate_fm_with_ua_ic(chunks, shoot_range_L1, sing[1].ua_left, - ctrl, equil, mats, intr; - backward=true, psi_ua=sing[1].psi_ua_left) + ctrl, equil, mats, intr; + backward=true, psi_ua=sing[1].psi_ua_left) elseif !isempty(shoot_range_L1) return assemble_fm_matrix(propagators, shoot_range_L1; - T_init=has_ua ? T_left_mats[1] : nothing) + T_init=has_ua ? T_left_mats[1] : nothing) else return has_ua ? T_left_mats[1] : Matrix{ComplexF64}(I, 2N, 2N) end @@ -377,19 +384,19 @@ end # the catastrophically ill-conditioned axis FM. Fortran-matched structure with # nMat = (2 + 4·msing)·N. Returns (M, nMat, col_edge). function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, - uShootL::Vector{Matrix{ComplexF64}}, - uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int}, - msing::Int, N::Int, - wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) + uShootL::Vector{Matrix{ComplexF64}}, + uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int}, + msing::Int, N::Int, + wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) # STRIDE global BVP block structure [Glasser-Kolemen 2018 PoP 25, 032501 Eq. 37]. nMat = (2 + 4 * msing) * N col_axis = 1:N - col_edge = (nMat - N + 1):nMat + col_edge = (nMat-N+1):nMat M = zeros(ComplexF64, nMat, nMat) # Axis matching: uShootL[1] · c_left[1] = uAxis · c_axis (2N equations) M[1:2N, _col_left(1, N)] .= uShootL[1] - M[1:2N, col_axis] .= -uAxis + M[1:2N, col_axis] .= -uAxis row_offset = 2N for j in 1:msing @@ -398,24 +405,24 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, for i in 1:2N if i != ipert_j && i != ipert_j + N row_offset += 1 - M[row_offset, _col_left(j, N)[i]] = 1 + M[row_offset, _col_left(j, N)[i]] = 1 M[row_offset, _col_right(j, N)[i]] = -1 end end - junc_rows = (row_offset + 1):(row_offset + 2N) + junc_rows = (row_offset+1):(row_offset+2N) if j < msing # Midpoint matching between consecutive surfaces - M[junc_rows, _col_right(j, N)] .= -uShootR[j] - M[junc_rows, _col_left(j+1, N)] .= uShootL[j+1] + M[junc_rows, _col_right(j, N)] .= -uShootR[j] + M[junc_rows, _col_left(j+1, N)] .= uShootL[j+1] else # Edge junction M[junc_rows, _col_right(msing, N)] .= uShootR[msing] if wv !== nothing - M[junc_rows[1:N], col_edge] .= -I(N) - M[junc_rows[N+1:end], col_edge] .= wv .* psio^2 + M[junc_rows[1:N], col_edge] .= -I(N) + M[junc_rows[(N+1):end], col_edge] .= wv .* psio^2 else - M[junc_rows[N+1:end], col_edge] .= -I(N) + M[junc_rows[(N+1):end], col_edge] .= -I(N) end end row_offset = last(junc_rows) @@ -425,7 +432,7 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, for j in 1:msing ipert_j = ipert_all[j] row_offset += 1 - M[row_offset, _col_left(j, N)[ipert_j]] = 1 + M[row_offset, _col_left(j, N)[ipert_j]] = 1 row_offset += 1 M[row_offset, _col_right(j, N)[ipert_j]] = 1 end @@ -460,18 +467,18 @@ end # Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available). # Uses the conditioned axis propagator Phi_R[1][:,N+1:2N] in place of S-axis matching. function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}}, - Phi_R_mats::Vector{Matrix{ComplexF64}}, ipert_all::Vector{Int}, - msing::Int, N::Int, - T_left_inv::Vector{Matrix{ComplexF64}}, - T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool, - wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) + Phi_R_mats::Vector{Matrix{ComplexF64}}, ipert_all::Vector{Int}, + msing::Int, N::Int, + T_left_inv::Vector{Matrix{ComplexF64}}, + T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool, + wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) nMat = (2 + 4 * msing) * N col_axis = 1:N - col_edge = (N + 4N*msing + 1):nMat + col_edge = (N+4N*msing+1):nMat M = zeros(ComplexF64, nMat, nMat) M[1:2N, (N+1):(N+2N)] .= Phi_L_mats[1] - M[1:2N, col_axis] .= -view(Phi_R_mats[1], :, N+1:2N) + M[1:2N, col_axis] .= -view(Phi_R_mats[1], :, (N+1):2N) row_drive_base = 2N + (4N-2)*msing for j in 1:msing @@ -482,45 +489,41 @@ function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}}, for i in 1:2N if i != ipert_j && i != ipert_j + N row_cont += 1 - M[row_cont, cl[i]] = 1 + M[row_cont, cl[i]] = 1 M[row_cont, cr[i]] = -1 end end - junc_rows = (row_cont + 1):(2N + (4N-2)*j) + junc_rows = (row_cont+1):(2N+(4N-2)*j) if j < msing - M[junc_rows, cr] .= Phi_R_mats[j+1] + M[junc_rows, cr] .= Phi_R_mats[j+1] M[junc_rows, _col_left(j+1, N)] .= -Phi_L_mats[j+1] else M[junc_rows, cr] .= Phi_R_mats[msing+1] if wv !== nothing - M[junc_rows[1:N], col_edge] .= -I(N) - M[junc_rows[N+1:end], col_edge] .= wv .* psio^2 + M[junc_rows[1:N], col_edge] .= -I(N) + M[junc_rows[(N+1):end], col_edge] .= wv .* psio^2 else - M[junc_rows[N+1:end], col_edge] .= -I(N) + M[junc_rows[(N+1):end], col_edge] .= -I(N) end end if has_ua - M[row_drive_base + 2j-1, cl] .= T_left_inv[j][ipert_j, :] - M[row_drive_base + 2j, cr] .= T_right_inv[j][ipert_j, :] + M[row_drive_base+2j-1, cl] .= T_left_inv[j][ipert_j, :] + M[row_drive_base+2j, cr] .= T_right_inv[j][ipert_j, :] else - M[row_drive_base + 2j-1, cl[ipert_j]] = 1 - M[row_drive_base + 2j, cr[ipert_j]] = 1 + M[row_drive_base+2j-1, cl[ipert_j]] = 1 + M[row_drive_base+2j, cr[ipert_j]] = 1 end end return M, nMat, col_edge end -# Solve the BVP for each driving configuration and apply the PEST3 four-term combination. -# Promotes to Complex{Double64} if ctrl.extended_precision_bvp (default true) — the PEST3 -# combination subtracts dp_raw entries up to ~3×10⁴ larger than the result, and Float64 -# precision lets the imaginary part drift 2–5× on DIIID-class equilibria. +# Solve the BVP and apply the PEST3 four-term combination. function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int, - use_S_axis::Bool, ipert_all::Vector{Int}, col_edge, - ctrl, debug::Bool) + use_S_axis::Bool, ipert_all::Vector{Int}, col_edge, + ctrl, debug::Bool) s2 = 2 * msing - Tc = (ctrl === nothing || ctrl.extended_precision_bvp) ? Complex{Double64} : ComplexF64 + Tc = (ctrl !== nothing && ctrl.extended_precision_bvp) ? Complex{Double64} : ComplexF64 M_solve = Tc.(M) - M_lu = lu(M_solve; check=false) use_lu = issuccess(M_lu) M_pinv = use_lu ? nothing : pinv(M_solve) @@ -538,32 +541,30 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, x = use_lu ? (M_lu \ b) : (M_pinv * b) debug && _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, - ipert_all, col_edge, use_S_axis) + ipert_all, col_edge, use_S_axis) for ksing in 1:msing ipert_k = ipert_all[ksing] dp_raw[dRow, 2ksing-1] = x[_col_left(ksing, N)[ipert_k+N]] - dp_raw[dRow, 2ksing] = x[_col_right(ksing, N)[ipert_k+N]] + dp_raw[dRow, 2ksing] = x[_col_right(ksing, N)[ipert_k+N]] end end # PEST3 four-term combination [Chance PPPL-2527; Glasser-Kolemen 2018 PoP 25, 032501 Eq. 31]. - # Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw, in extended precision. - deltap_ext = zeros(Tc, msing, msing) + # Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw. + deltap = zeros(Tc, msing, msing) for i in 1:msing, j in 1:msing - deltap_ext[i, j] = dp_raw[2i, 2j] - dp_raw[2i, 2j-1] - dp_raw[2i-1, 2j] + dp_raw[2i-1, 2j-1] + deltap[i, j] = dp_raw[2i, 2j] - dp_raw[2i, 2j-1] - dp_raw[2i-1, 2j] + dp_raw[2i-1, 2j-1] end - deltap = ComplexF64.(deltap_ext) - debug && _log_bvp_pest3(dp_raw, deltap, s2, msing, Tc) - # Return the PEST3-combined matrix AND the raw 2msing×2msing D' matrix (ComplexF64 - # for compatibility with downstream pest3_decompose / HDF5 writer). - return deltap, ComplexF64.(dp_raw) + debug && _log_bvp_pest3(dp_raw, deltap, s2, msing) + # Return the PEST3-combined matrix AND the raw 2msing×2msing D' matrix. + return ComplexF64.(deltap), ComplexF64.(dp_raw) end # Logging helpers for `compute_delta_prime_matrix!`. Called only when debug=true. function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, - Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) + Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) @info "Δ' BVP: $(length(chunks)) chunks, $msing surfaces, N=$N" @info "Δ' BVP: Axis BC: $(use_S_axis ? "S-based (Riccati)" : "FM-based (conditioned)")" @info "Δ' BVP: Asymptotic basis: $(has_ua ? "available" : "NOT available (raw basis driving)")" @@ -575,8 +576,8 @@ function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, if has_ua for j in 1:msing sp = sing[j] - T_l = [sp.ua_left[:,:,1]; sp.ua_left[:,:,2]] - T_r = [sp.ua_right[:,:,1]; sp.ua_right[:,:,2]] + T_l = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]] + T_r = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]] @info " Surface $j: cond(T_left)=$(@sprintf("%.2e", cond(T_l))), cond(T_right)=$(@sprintf("%.2e", cond(T_r)))" ipert_j = ipert_all[j] @info " Surface $j ua_left (ipert=$ipert_j, psi_ua_left=$(@sprintf("%.8f", sp.psi_ua_left))):" @@ -586,7 +587,7 @@ function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, @info " small: ua(1,$(ipert_j+N),1)=$(@sprintf("%16.8e %16.8e", real(sp.ua_left[1,ipert_j+N,1]), imag(sp.ua_left[1,ipert_j+N,1])))" end end - for j in 1:msing-1 + for j in 1:(msing-1) Phi_L_h, Phi_R_h = Phi_R_halves[j] @info " Inter-surface $j→$(j+1): half_L cond=$(@sprintf("%.2e",cond(Phi_L_h))), half_R cond=$(@sprintf("%.2e",cond(Phi_R_h))), full cond=$(@sprintf("%.2e",cond(Phi_R_mats[j+1])))" end @@ -603,7 +604,7 @@ function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, end function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_left, - T_left_mats, ipert_all, has_ua, msing, N) + T_left_mats, ipert_all, has_ua, msing, N) @info " Shooting propagators (S-based axis BC, no axis unknowns):" for j in 1:msing shoot_R_str = @sprintf("%.2e", cond(uShootR[j])) @@ -613,7 +614,7 @@ function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_ S1 = S_at_surface_left[1] if has_ua T1 = T_left_mats[1] - axis_BC = T1[1:N, :] - S1 * T1[N+1:2N, :] + axis_BC = T1[1:N, :] - S1 * T1[(N+1):2N, :] @info " S-axis BC matrix: cond=$(@sprintf("%.2e", cond(axis_BC)))" end for j in 1:msing @@ -627,7 +628,7 @@ function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_ @info " uShootL[$j] col ipert=$ipert_j norm=$(@sprintf("%.2e", col_norms_L[ipert_j])), col ipert+N=$(ipert_j+N) norm=$(@sprintf("%.2e", col_norms_L[ipert_j+N]))" end end - for j in 1:msing-1 + for j in 1:(msing-1) mid_block = hcat(uShootR[j], -uShootL[j+1]) @info " Midpoint $j→$(j+1): cond([uShootR[$j] | -uShootL[$(j+1)]]) = $(@sprintf("%.2e", cond(mid_block)))" col_norms_Ljp1 = [norm(view(uShootL[j+1], :, k)) for k in 1:2N] @@ -636,7 +637,7 @@ function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_ end function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, - ipert_all, col_edge, use_S_axis) + ipert_all, col_edge, use_S_axis) residual = norm(ComplexF64.(M_solve * x - b)) side_str = side == 1 ? "left" : "right" @info " BVP solve: jsing=$jsing side=$side_str (dRow=$dRow): ||Mx-b||=$(@sprintf("%.2e", residual)), ||x||=$(@sprintf("%.2e", Float64(norm(x))))" @@ -644,9 +645,9 @@ function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, ipert_ks = ipert_all[ks] cl = _col_left(ks, N) cr = _col_right(ks, N) - xl_big = ComplexF64(x[cl[ipert_ks]]) + xl_big = ComplexF64(x[cl[ipert_ks]]) xl_small = ComplexF64(x[cl[ipert_ks+N]]) - xr_big = ComplexF64(x[cr[ipert_ks]]) + xr_big = ComplexF64(x[cr[ipert_ks]]) xr_small = ComplexF64(x[cr[ipert_ks+N]]) @info " surf $ks: x_left[big]=$(@sprintf("%+.4e%+.4ei", real(xl_big), imag(xl_big))), x_left[small]=$(@sprintf("%+.4e%+.4ei", real(xl_small), imag(xl_small)))" @info " surf $ks: x_right[big]=$(@sprintf("%+.4e%+.4ei", real(xr_big), imag(xr_big))), x_right[small]=$(@sprintf("%+.4e%+.4ei", real(xr_small), imag(xr_small)))" @@ -657,10 +658,10 @@ function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, end end -function _log_bvp_pest3(dp_raw, deltap, s2, msing, Tc) - @info "Δ' BVP: Full dp_raw matrix ($(s2)×$(s2)) [$(Tc)]:" +function _log_bvp_pest3(dp_raw, deltap, s2, msing) + @info "Δ' BVP: Full dp_raw matrix ($(s2)×$(s2)):" for i in 1:s2 - row_str = join([@sprintf("%+.6e", Float64(real(dp_raw[i,j]))) for j in 1:s2], " ") + row_str = join([@sprintf("%+.6e", Float64(real(dp_raw[i, j]))) for j in 1:s2], " ") @info " dp_raw[$i,:] = $row_str" end @info "Δ' BVP: Raw dp diagonal = $([@sprintf("%.4f%+.4fi", Float64(real(dp_raw[i,i])), Float64(imag(dp_raw[i,i]))) for i in 1:s2])" @@ -720,8 +721,8 @@ function pest3_decompose(dp_raw::AbstractMatrix) for i in 1:m, j in 1:m LL = dp_raw[2i-1, 2j-1] LR = dp_raw[2i-1, 2j] - RL = dp_raw[2i, 2j-1] - RR = dp_raw[2i, 2j] + RL = dp_raw[2i, 2j-1] + RR = dp_raw[2i, 2j] Ap[i, j] = RR + RL + LR + LL Bp[i, j] = RR - RL + LR - LL Gp[i, j] = RR + RL - LR - LL diff --git a/src/ForceFreeStates/Surfaces/Asymptotics.jl b/src/ForceFreeStates/Surfaces/Asymptotics.jl index 26746c0ca..eb15bd6fc 100644 --- a/src/ForceFreeStates/Surfaces/Asymptotics.jl +++ b/src/ForceFreeStates/Surfaces/Asymptotics.jl @@ -1,5 +1,35 @@ # Frobenius asymptotics and resonant-basis evaluation at singular surfaces. +""" + resonant_block_exponents(m0, sig) -> (α, p_big, p_small) + +Eigenvalues λ = tr/2 ∓ √(tr²/4 − det) of one resonant 2×2 M₀ block. +α = √(tr²/4 − det) is the Mercier half-separation. Powers are sig·λ, big solution +(more negative real part) first. A traceless block returns α = √(−det) and (−α, +α). +When tr²/4 − det is negative (Mercier-unstable), α is taken with Im α > 0, so both sides +of the surface get the same α and the same big root whatever the round-off sign. +""" +function resonant_block_exponents(m0::AbstractMatrix, sig::Real) + tr = m0[1, 1] + m0[2, 2] + det = m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2] + α = sqrt(ComplexF64(tr * tr / 4 - det)) + abs(real(α)) <= 64 * Base.eps(Float64) * abs(α) && (α = ComplexF64(0, abs(imag(α)))) + p_minus = sig * (tr / 2 - α) + p_plus = sig * (tr / 2 + α) + # Equal real parts: sig ≥ 0 keeps big = sig·(tr/2 − α). + if real(p_minus) < real(p_plus) || (real(p_minus) == real(p_plus) && sig >= 0) + return α, p_minus, p_plus + else + return α, p_plus, p_minus + end +end + +# Scale the big and small resonant columns by dpsi^power. +function _apply_column_powers!(big_cols, small_cols, p_big, p_small, dpsi) + big_cols .*= dpsi^p_big + small_cols .*= dpsi^p_small +end + """ compute_sing_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) @@ -27,7 +57,6 @@ function compute_sing_asymptotics( mats::MatrixSplines, intr::ForceFreeStatesInternal; sig::Float64=1.0, - alpha_override::Union{Nothing,Vector{ComplexF64}}=nothing, sing_order::Int=ctrl.sing_order ) @@ -58,41 +87,31 @@ function compute_sing_asymptotics( Matrix(vcat([transpose(mmat[r1[i], r2, :, 1]) for i in eachindex(r1)]...)) end - # Alpha (Mercier index) — Fortran computes this ONCE from the RIGHT-SIDE m0mat - # and reuses it for both left and right vmat (matching Fortran STRIDE). - # When alpha_override is provided (for the left-side call), use that instead. - # Fortran: di = m0(1,1)*m0(2,2) - m0(2,1)*m0(1,2); alpha = sqrt(-di) - # This matches eigenvalues only when tr(m0mat_block) = 0. - alpha = if alpha_override !== nothing - alpha_override - else - # Match Fortran exactly: alpha = sqrt(-det(m0mat_block)) for each resonant mode - [sqrt(-ComplexF64(m0mat[(2*(i-1)+1), (2*(i-1)+1)] * m0mat[(2*i), (2*i)] - - m0mat[(2*i), (2*(i-1)+1)] * m0mat[(2*(i-1)+1), (2*i)])) - for i in eachindex(r1)] - end - - # This is the parameter α but for all modes - α = 0 for non-resonant modes - power[ipert_res] .= -alpha - power[ipert_res .+ intr.numpert_total] .= alpha - - # Zeroth-order non-resonant solutions + # Zeroth-order non-resonant solutions. Nonresonant powers stay 0. for ipert in 1:intr.numpert_total vmat[ipert, ipert, 1, 1] = 1 vmat[ipert, ipert+intr.numpert_total, 2, 1] = 1 end - # Zeroth-order resonant solutions: v_big_ξ' = -(m0(1,1) ± sig·α)/m0(1,2). - # Matches Fortran STRIDE sing_vmat (sig·α sign convention separates left vs right side). + # This side's block: power is sig·λ (big first); the ξ' row is the eigenvector for λ. + alpha = Vector{ComplexF64}(undef, length(r1)) for i in eachindex(r1) m0mat_block = m0mat[(2*(i-1)+1):(2*i), (2*(i-1)+1):(2*i)] + alpha_i, p_big, p_small = resonant_block_exponents(m0mat_block, sig) + alpha[i] = alpha_i + # The assembled block's trace is ideally zero; a large one flags a non-converged block. + block_tr = m0mat_block[1, 1] + m0mat_block[2, 2] + abs(block_tr) / 2 > 1e-3 * abs(alpha_i) && + @warn @sprintf("Resonant M₀ block at ψ=%.6f, m=%d (sig=%+.0f): |tr|/2 = %.3e is not small against |α| = %.3e", + singp.psifac, singp.m[i], sig, abs(block_tr) / 2, abs(alpha_i)) r1_i = r1[i] + power[r1_i] = p_big + power[r1_i+intr.numpert_total] = p_small r2_i = r1_i + intr.numpert_total - alpha_i = alpha[i] vmat[r1_i, r1_i, 1, 1] = 1 vmat[r1_i, r2_i, 1, 1] = 1 - vmat[r1_i, r1_i, 2, 1] = -(m0mat_block[1, 1] + sig * alpha_i) / m0mat_block[1, 2] - vmat[r1_i, r2_i, 2, 1] = -(m0mat_block[1, 1] - sig * alpha_i) / m0mat_block[1, 2] + vmat[r1_i, r1_i, 2, 1] = -(m0mat_block[1, 1] - sig * p_big) / m0mat_block[1, 2] + vmat[r1_i, r2_i, 2, 1] = -(m0mat_block[1, 1] - sig * p_small) / m0mat_block[1, 2] end # Higher order solutions — sig propagates through the recursion (Fortran STRIDE sing_solve). @@ -111,8 +130,10 @@ 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] - msg *= @sprintf(" di= %+.12e, alpha= %+.12e %+.12ei\n", real(di), real(alpha[1]), imag(alpha[1])) + block_det = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2] + block_tr = m0mat[1, 1] + m0mat[2, 2] + msg *= @sprintf(" det= %+.12e, tr= %+.12e %+.12ei\n", real(block_det), real(block_tr), imag(block_tr)) + msg *= @sprintf(" alpha= %+.12e %+.12ei\n", 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])) msg *= @sprintf(" vmat(ip,ip+N,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0+N, 2, 1]), imag(vmat[ipert0, ipert0+N, 2, 1])) @@ -457,9 +478,9 @@ See equation 47 in the Glasser 2016 DCON paper. Identical to the Fortran - `vmat::Array{ComplexF64,4}`: V matrix power series (modified in-place) - `mmat::Array{ComplexF64,4}`: M matrix power series - `m0mat::Matrix{ComplexF64}`: Zeroth order M matrix - - `alpha::Vector{ComplexF64}`: Eigenvalues of M₀ for resonant modes + - `alpha::Vector{ComplexF64}`: Mercier exponent of each resonant block - `r1, r2, n1, n2::Vector{Int}`: Resonant and nonresonant indices - - `power::Vector{ComplexF64}`: α values for all modes (0 for nonresonant) + - `power::Vector{ComplexF64}`: Frobenius power of each solution column (0 for nonresonant) - `k::Int`: The current order in the power series expansion """ @with_pool pool function solve_higher_order_vmat!( @@ -569,7 +590,6 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64) # SingAsymptotics (left vs right vmat built with sig=-1 or sig=+1). # Matches Fortran STRIDE sing_get_ua: sqrtfac=SQRT(dpsi), always positive. sqrtfac = sqrt(dpsi) - pfac_base = dpsi # used for dpsi^alpha below # Compute power series via Horner's method (eq. 45 in Glasser 2016) ua = copy(sing_asymp.vmat[:, :, :, 2*sing_asymp.sing_order+1]) @@ -577,11 +597,10 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64) ua .= ua .* sqrtfac .+ sing_asymp.vmat[:, :, :, iorder+1] end - # Restore powers (unshear v→u) — matches Fortran STRIDE sing_get_ua + # Restore powers (unshear v→u). Resonant columns use their Frobenius powers. for i in eachindex(r1) - 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^α + _apply_column_powers!(@view(ua[:, r2[2*i-1], :]), @view(ua[:, r2[2*i], :]), + sing_asymp.power[r2[2*i-1]], sing_asymp.power[r2[2*i]], dpsi) ua[r1[i], :, 1] ./= sqrtfac # resonant row ξ: /√dpsi ua[r1[i], :, 2] .*= sqrtfac # resonant row ξ': *√dpsi end @@ -600,8 +619,8 @@ outer-region Galerkin solver (`gal_extension`, `sing_matvec`). Direction is carried by the `SingAsymptotics` (left vs right vmat built with sig=∓1), exactly as in `sing_get_ua`, so `dpsi` is always positive and the arithmetic is real (no analytic continuation). The term-by-term derivative is built via a `power` array tracking each term's total exponent (in -half-powers of `dpsi`: the 2*sing_order offset, the ∓1 shearing on resonant rows, and the ∓2α on -resonant columns), then the same shearing/`pfac` restoration and the chain-rule factor `1/(2·dpsi)` +half-powers of `dpsi`: the 2*sing_order offset, the ∓1 shearing on resonant rows, and twice the +column's Frobenius power), then the same shearing restoration and the chain-rule factor `1/(2·dpsi)` are applied. The physical d/dψ sign for the left side (ψ = ψ_res − dpsi) is applied by the caller. ### Arguments @@ -622,8 +641,8 @@ function sing_get_dua(sing_asymp::SingAsymptotics, dpsi::Float64) power[r1, :, 1] .-= 1 power[r1, :, 2] .+= 1 for i in eachindex(r1) - power[:, r2[2*i-1], :] .-= 2 * sing_asymp.alpha[i] - power[:, r2[2*i], :] .+= 2 * sing_asymp.alpha[i] + power[:, r2[2*i-1], :] .+= 2 * sing_asymp.power[r2[2*i-1]] + power[:, r2[2*i], :] .+= 2 * sing_asymp.power[r2[2*i]] end # Power series derivative by Horner's method (Fortran sing_get_dua) @@ -635,9 +654,8 @@ function sing_get_dua(sing_asymp::SingAsymptotics, dpsi::Float64) # Restore shearing and resonant powers, then the chain-rule factor 1/(2·dpsi) (Fortran sing_get_dua) for i in eachindex(r1) - pfac = dpsi^sing_asymp.alpha[i] - dua[:, r2[2*i-1], :] ./= pfac - dua[:, r2[2*i], :] .*= pfac + _apply_column_powers!(@view(dua[:, r2[2*i-1], :]), @view(dua[:, r2[2*i], :]), + sing_asymp.power[r2[2*i-1]], sing_asymp.power[r2[2*i]], dpsi) dua[r1[i], :, 1] ./= sqrtfac dua[r1[i], :, 2] .*= sqrtfac end @@ -685,12 +703,9 @@ function sing_get_ua_res!(out::AbstractArray{ComplexF64,3}, sing_asymp::SingAsym out[ii, 2, k] = as end - # Unshear v→u: big column /dpsi^α, small column *dpsi^α; resonant row /√dpsi (qty1), *√dpsi (qty2). - pfac = dpsi^sing_asymp.alpha[1] - @inbounds for k in 1:2, ii in 1:N - out[ii, 1, k] /= pfac - out[ii, 2, k] *= pfac - end + # Unshear v→u by the column Frobenius powers; resonant row /√dpsi (qty1), *√dpsi (qty2). + _apply_column_powers!(@view(out[:, 1, :]), @view(out[:, 2, :]), + sing_asymp.power[cbig], sing_asymp.power[csml], dpsi) @inbounds out[ρ, 1, 1] /= sqrtfac @inbounds out[ρ, 2, 1] /= sqrtfac @inbounds out[ρ, 1, 2] *= sqrtfac @@ -713,15 +728,15 @@ function sing_get_dua_res!(out::AbstractArray{ComplexF64,3}, sing_asymp::SingAsy ρ = sing_asymp.r1[1] cbig = sing_asymp.r2[1] csml = sing_asymp.r2[2] - α = sing_asymp.alpha[1] + p_big = sing_asymp.power[cbig] + p_small = sing_asymp.power[csml] nterm = 2 * order - # Per-term exponent base (×1, in half-powers of dpsi): 2·order + row offset (∓1 on the resonant row, - # by qty) + column offset (−2α big, +2α small). Tracked as scalars per column instead of a full array. + # Exponent base in half-powers of dpsi: 2·order + resonant-row offset (∓1 by qty) + 2·power. @inbounds for k in 1:2, ii in 1:N roff = ii == ρ ? (k == 1 ? ComplexF64(-1) : ComplexF64(1)) : ComplexF64(0) - pb = nterm + roff - 2 * α - ps = nterm + roff + 2 * α + pb = nterm + roff + 2 * p_big + ps = nterm + roff + 2 * p_small ab = vmat[ii, cbig, k, nterm+1] * pb as = vmat[ii, csml, k, nterm+1] * ps for t in (nterm-1):-1:0 @@ -734,12 +749,8 @@ function sing_get_dua_res!(out::AbstractArray{ComplexF64,3}, sing_asymp::SingAsy out[ii, 2, k] = as end - pfac = dpsi^α + _apply_column_powers!(@view(out[:, 1, :]), @view(out[:, 2, :]), p_big, p_small, dpsi) twodpsi = 2 * dpsi - @inbounds for k in 1:2, ii in 1:N - out[ii, 1, k] /= pfac - out[ii, 2, k] *= pfac - end @inbounds out[ρ, 1, 1] /= sqrtfac @inbounds out[ρ, 2, 1] /= sqrtfac @inbounds out[ρ, 1, 2] *= sqrtfac @@ -782,11 +793,8 @@ function sing_get_ua_res_cut!(out::AbstractArray{ComplexF64,3}, sing_asymp::Sing out[ii, 2, k] = vmat[ii, csml, k, 1] end - pfac = dpsi^sing_asymp.alpha[1] - @inbounds for k in 1:2, ii in 1:N - out[ii, 1, k] /= pfac - out[ii, 2, k] *= pfac - end + _apply_column_powers!(@view(out[:, 1, :]), @view(out[:, 2, :]), + sing_asymp.power[cbig], sing_asymp.power[csml], dpsi) @inbounds out[ρ, 1, 1] /= sqrtfac @inbounds out[ρ, 2, 1] /= sqrtfac @inbounds out[ρ, 1, 2] *= sqrtfac diff --git a/src/ForceFreeStates/Surfaces/Types.jl b/src/ForceFreeStates/Surfaces/Types.jl index 5f995d2a0..eba795a2d 100644 --- a/src/ForceFreeStates/Surfaces/Types.jl +++ b/src/ForceFreeStates/Surfaces/Types.jl @@ -40,12 +40,12 @@ This data is computed on-demand during singular surface crossings in `cross_idea ## Fields - - `alpha::Vector{ComplexF64}` - Resonant matrix eigenvalues + - `alpha::Vector{ComplexF64}` - Mercier exponent α = √(tr²/4 − det) of each resonant block. - `r1::Vector{Int}` - Resonant indices along first index - `r2::Vector{Int}` - Resonant indices along second index - `n1::Vector{Int}` - Nonresonant indices along first index - `n2::Vector{Int}` - Nonresonant indices along second index - - `power::Vector{ComplexF64}` - Power series coefficients + - `power::Vector{ComplexF64}` - Frobenius power of each column in |ψ − ψ_s| (0 off resonance). Resonant powers are sig·λ of that side's M₀ block, big solution (more negative real part) first. - `vmat::Array{ComplexF64,4}` - Power series of V matrix for asymptotic analysis - `mmat::Array{ComplexF64,4}` - Power series of M matrix for asymptotic analysis - `m0mat::Matrix{ComplexF64}` - Zeroth order M matrix projected onto resonant subspace diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 6322d2c79..45e157f66 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -1303,7 +1303,9 @@ function efc_couplings( ψk = kp.xs ω_E = [kp.omegaE_spline(ψ) for ψ in ψk] ω_star_T = [-2π * kp.Ti_deriv(ψ) / (chrg * chi1) for ψ in ψk] - ω_star_n = [(n = kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk] + ω_star_n = [ + (n=kp.ni_spline(ψ); n > 0 ? -2π * kp.Ti_spline(ψ) * kp.ni_deriv(ψ) / (chrg * chi1 * n) : 0.0) for ψ in ψk + ] ω_φ = ω_E .+ ω_star_n .+ ω_star_T # Momentum weight n_i·dV/dψ with ⟨R²⟩ taken as R₀², inside the plasma only. w = [ψ <= ffs.psilim ? max(kp.ni_spline(ψ), 0.0) * equil.profiles.dVdpsi_spline(clamp(ψ, equil.profiles.xs[1], equil.profiles.xs[end])) : 0.0 for ψ in ψk] diff --git a/test/runtests.jl b/test/runtests.jl index e18762e18..1d4dafb9e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,6 +33,7 @@ else include("./runtests_parallel_integration.jl") include("./runtests_result_struct.jl") include("./runtests_solve_api.jl") + include("./runtests_galerkin_cholesky.jl") include("./runtests_dominant_coupling.jl") include("./runtests_error_fields.jl") include("./runtests_tolerance_toml.jl") diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 98fca1a63..a53cfb676 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -640,3 +640,100 @@ end @test !odet.du_store_populated end end + +@testset "resonant block Frobenius exponents" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + + @testset "traceless block reduces to ∓√(−det) on both sides" begin + m0 = ComplexF64[1 2; 3 -1] + α0 = sqrt(-(m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2])) + for sig in (1.0, -1.0) + α, p_big, p_small = FFS.resonant_block_exponents(m0, sig) + @test α ≈ α0 + @test p_big ≈ -α0 + @test p_small ≈ α0 + end + end + + @testset "nonzero trace uses the block eigenvalues" begin + m0 = ComplexF64[2 1; 0.5 4] + tr = m0[1, 1] + m0[2, 2] + block_det = m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2] + α0 = sqrt(tr * tr / 4 - block_det) + α, p_big, p_small = FFS.resonant_block_exponents(m0, 1.0) + @test α === α0 + @test p_big === tr / 2 - α0 + @test p_small === tr / 2 + α0 + @test abs(det(m0 - p_big * I)) < 1e-12 + @test abs(det(m0 - p_small * I)) < 1e-12 + + α, p_big, p_small = FFS.resonant_block_exponents(m0, -1.0) + p_minus = -(tr / 2 - α0) + p_plus = -(tr / 2 + α0) + @test p_big === (real(p_minus) < real(p_plus) ? p_minus : p_plus) + @test p_small === (real(p_minus) < real(p_plus) ? p_plus : p_minus) + @test real(p_big) < real(p_small) + @test abs(det(m0 - (-p_big) * I)) < 1e-12 + @test abs(det(m0 - (-p_small) * I)) < 1e-12 + end + + @testset "equal real parts keep the sig ≥ 0 ordering" begin + # tr = 2, det = 2 → α = √(−1). Both powers have the same real part. + m0 = ComplexF64[1 1; -1 1] + tr = m0[1, 1] + m0[2, 2] + block_det = m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2] + α0 = sqrt(tr * tr / 4 - block_det) + @test real(α0) == 0 + α, p_big, p_small = FFS.resonant_block_exponents(m0, 1.0) + @test p_big === tr / 2 - α0 + @test p_small === tr / 2 + α0 + α, p_big, p_small = FFS.resonant_block_exponents(m0, -1.0) + @test p_big === -(tr / 2 + α0) + @test p_small === -(tr / 2 - α0) + end + + @testset "Hermitian-structured block: shared imaginary shift, real eigenvectors, same powers on both sides" begin + # Shape the assembly produces: [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real and κ → −κ on the left side. + κ, φ, γ = 0.3 + 0.02im, 0.7, 2.5 + m_right = ComplexF64[-κ+0.5 φ; γ conj(κ)-0.5] + m_left = ComplexF64[κ-0.5 φ; γ -conj(κ)+0.5] + α_right, big_right, small_right = FFS.resonant_block_exponents(m_right, 1.0) + α_left, big_left, small_left = FFS.resonant_block_exponents(m_left, -1.0) + @test α_right ≈ sqrt((real(κ) - 0.5)^2 + φ * γ) + @test α_left ≈ α_right + @test imag(big_right) ≈ -imag(κ) + @test imag(small_right) ≈ -imag(κ) + @test big_left ≈ big_right + @test small_left ≈ small_right + # ξ′ entries of both leading eigenvectors are real, so the big/small cross-Wronskian 2α/φ is real. + @test abs(imag(-(m_right[1, 1] - big_right) / φ)) < 1e-14 + @test abs(imag(-(m_right[1, 1] - small_right) / φ)) < 1e-14 + @test -(m_right[1, 1] - small_right) / φ + (m_right[1, 1] - big_right) / φ ≈ 2 * α_right / φ + end + + @testset "Mercier-unstable block: round-off sign of Im(tr²/4 − det) does not flip α or the big root" begin + # tr²/4 − det = −0.36 ± 1e-17i; without canonicalisation α = ±0.6i, so the big root depends on that sign. + for δ in (1e-17, -1e-17) + m0 = ComplexF64[0.1im 1.0; -0.36+δ*im 0.1im] + α, p_big, p_small = FFS.resonant_block_exponents(m0, 1.0) + @test real(α) == 0 && α ≈ 0.6im + @test p_big ≈ -0.5im + @test p_small ≈ 0.7im + end + end + + @testset "column scaling multiplies each column by dpsi^power" begin + dpsi = 1.3e-4 + p_big = ComplexF64(-1.2 + 0.1im) + p_small = ComplexF64(1.8) + big = ones(ComplexF64, 1, 1) + small = ones(ComplexF64, 1, 1) + FFS._apply_column_powers!(big, small, p_big, p_small, dpsi) + big_ref = ones(ComplexF64, 1, 1) + small_ref = ones(ComplexF64, 1, 1) + big_ref .*= dpsi^p_big + small_ref .*= dpsi^p_small + @test big == big_ref + @test small == small_ref + end +end diff --git a/test/runtests_galerkin_cholesky.jl b/test/runtests_galerkin_cholesky.jl new file mode 100644 index 000000000..01bc5deee --- /dev/null +++ b/test/runtests_galerkin_cholesky.jl @@ -0,0 +1,38 @@ +using Test +using LinearAlgebra + +# gal_banded_solve! on a full-band (gbtrf layout) system built like the Galerkin matrix: Hermitian, with +# Dirichlet identity rows (row zeroed, column kept) and nonzero prescribed values, as in the rpec coil columns. +@testset "Galerkin banded solve: Cholesky first, LU fallback" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + n, kl, fixed = 24, 3, [1, 2, 24] + H = zeros(ComplexF64, n, n) + for j in 1:n, i in max(1, j - kl):min(n, j + kl) + H[i, j] = i == j ? 4 + sin(j) : 0.3 * cis(0.3i + 0.7j) + end + H = (H + H') / 2 # diagonally dominant → positive definite + function system(A) + A = copy(A) + A[fixed, :] .= 0 + A[fixed, fixed] = I(length(fixed)) + ab = zeros(ComplexF64, 3kl + 1, n) + for j in 1:n, i in max(1, j - kl):min(n, j + kl) + ab[2kl+1+i-j, j] = A[i, j] + end + return A, ab + end + rhs = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n, _ in 1:2] + rhs[fixed, 2] .= [1.0, -0.5im, 2.0] + + A, ab = system(H) + sol = copy(rhs) + @test FFS.gal_banded_solve!(ab, sol, kl, kl, "auto") == "cholesky" + @test sol ≈ A \ rhs rtol = 1e-12 + + H[5, 5] = -1 + A, ab = system(H) + sol = copy(rhs) + @test FFS.gal_banded_solve!(copy(ab), sol, kl, kl, "auto") == "LU" + @test sol ≈ A \ rhs rtol = 1e-12 + @test_throws ErrorException FFS.gal_banded_solve!(ab, copy(rhs), kl, kl, "cholesky") +end diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index b427c3795..aa2f54f46 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -555,6 +555,15 @@ using TOML @test isapprox(real(dpm[1, 1]), +9.300000e+00; rtol=3e-2) # q=2, converged (ldp2048 9.244) @test isapprox(real(dpm[2, 2]), -5.344199e+00; rtol=1e-1) # q=3 @test isapprox(real(dpm[3, 3]), -1.590034e+01; rtol=1e-1) # q=4 + + # extended_precision_bvp: the same BVP solved in Complex{Double64} + ctrl_ext = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; + (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., extended_precision_bvp=true) + GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_delta_prime_matrix!( + intr, fm_propagators, fm_chunks; + wv=vac.wv, psio=equil.psio, + S_at_surface_left=fm_S_left, ctrl=ctrl_ext, equil=equil, mats=mats) + @test isapprox(real.(diag(intr.delta_prime_matrix)), real.(diag(dpm)); rtol=1e-8) end end