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 5544b0e60..9e0831265 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -192,7 +192,7 @@ gpec.toml. extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). # --- 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/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..6da63d2f1 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -60,7 +60,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 @@ -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) @@ -361,3 +352,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/test/runtests.jl b/test/runtests.jl index e1256e56d..f3f4dde76 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_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