Repository navigation
ForceFreeStates.Galerkin - FEATURE - Add the banded Cholesky solve for gal_solver = "cholesky" - #513
StuartBenjamin wants to merge 4 commits into
Conversation
…r gal_solver = "cholesky"
gal_solver = "cholesky" called LinearAlgebra.LAPACK.pbtrf!/pbtrs!, which do not exist, so the
option died with an UndefVarError. Call LAPACK zpbtrf('L')/zpbtrs('L') directly through
libblastrampoline, unscaled, on the lower band the assembly already builds for Cholesky. Unlike
the reference implementation, a non-positive-definite matrix (zpbtrf info > 0) raises an error
naming gal_solver = "LU" instead of continuing with a partial factor.
Add examples/LAR_gal_cholesky_example (LAR analytic equilibrium, positive-definite Galerkin
matrix), regression case lar_gal_cholesky, test/runtests_galerkin_cholesky.jl, and a solver note
in docs/src/galerkin.md.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
matt-pharr
left a comment
matt-pharr
left a comment
There was a problem hiding this comment.
Couple cleanup items and one main question, being whether we should do away with the flag altogether and just make every galerkin solve try to use cholesky first and then use LU as a fallback. See the relevant comment.
|
|
||
| # --- banded solve (gal.f) --- | ||
| ws.sol .= ws.rhs | ||
| if ws.solver == "LU" |
There was a problem hiding this comment.
So we cannot always use Cholesky. However, we could remove this flag by making the code default to Cholesky since it is always the superior solve method, and then if we are not positive definite fall back to LU. This would increase the runtime by cases that are not LU, but the galerkin matrix solve is not the biggest piece in the rdcon runtime (I believe it is the computation of the weak form integrals which we have slowed down here by switching to Gauss–Kronrod quadrature from LSODE). How do you feel about this @StuartBenjamin ? Perhaps we should just make it, in pseudocode:
try:
sol = choelsky(galmat,b)
except MatrixIsNotPositiveDefinite:
sol = LU(galmat, b)
| x = ComplexF64[cis(0.1j) * (1 + j / n) for j in 1:n] | ||
| B = A * hcat(x, -3x) | ||
| ab = lowerband(A) | ||
| FFS.gal_zpbtrf!(ab, kl) |
There was a problem hiding this comment.
is the point of this test just to make sure that gal_zpbtrf! works and does a cholesky solve?
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…to enable Cholesky solver for ncoil > 0
This branch has not been deployed
2 participants
Release note
diiid_n1,gal_resistive_diiidanddiiid_n1_riccatiare bit-identical. The new caselar_gal_choleskysets its first values. (harness @ b5b430f)gal_solver = "cholesky"now works for the Galerkin Δ′ solver. Before, it crashed withUndefVarError: pbtrf!. It factors the Galerkin matrix with LAPACK's banded Hermitian Cholesky (zpbtrf/zpbtrs, lower band, unscaled), the same method as the original RDCON implementation. When the matrix is not positive definite, it stops with an error that points togal_solver = "LU". The new exampleexamples/LAR_gal_cholesky_exampleshows a case where Cholesky applies.Regression report
Notes for reviewers
Reading the harness report.
lar_gal_choleskyshowsFAILEDfororigin/develop, which makes the harness exit non-zero. That is expected: the case's example deck is new in this PR, sodevelophas nothing to run. On this branch it records the values below.gal_resistive_peis N/A on both refs. This is pre-existing ondevelop: with the Galerkin integrator, PerturbedEquilibrium skips the plasma-response and singular-coupling stages ("free_boundarywas not produced by the galerkin integrator").What changed (
src/ForceFreeStates/Galerkin/GalerkinSolve.jl, Cholesky branch only):LinearAlgebra.LAPACKhas nopbtrf!/pbtrs!, so the option could never run. Two smallccallwrappers,gal_zpbtrf!andgal_zpbtrs!, call LAPACKzpbtrf('L')/zpbtrs('L')throughlibblastrampoline. That is the same library the LU path'sgbtrf!uses, so no new dependency is added.ldab = kl + 1). There is no scaling, matching the reference implementation.zpbtrf'sinfocode, so on a matrix that is not positive definite it continues with a partial factor and returns a wrong Δ′ silently. Hereinfo > 0raises an error naminggal_solver = "LU".GalerkinStructs.jl, a solver note indocs/src/galerkin.md, the new example, the new harness case and a new unit test.bugfix/galerkin-scaled-lu-prchanges only the LU branch, so the two PRs edit disjoint lines and merge cleanly in either order (checked withgit merge-tree).Where Cholesky applies. It needs a Hermitian positive-definite Galerkin matrix.
gal_nx256 and 1024gal_nx256, 1024 and 2048gal_nxdoes not help.Evidence:
examples/LAR_gal_cholesky_example(n = 1), Re Δ′. All runs on this PR @b5b430f57, same deck.gal_nx = 128gal_nx = 256(example as committed)gal_nx = 512gal_nx = 1024gal_nx = 256LAR_ideal_match_testODE settings (eulerlagrange_tolerance = 1e-12, …)gal_nx: 2/1 doesn't change, and 3/1 moves 0.04% from 128 to 1024.LAR_ideal_match_testsettings.develop, the same deck fails withUndefVarError: pbtrf! not defined in LinearAlgebra.LAPACK.Tests. The new
test/runtests_galerkin_cholesky.jlchecks four things:choleskyfactor;At
b5b430f57(1 CPU job):Galerkin banded Choleskysolve APIForceFreeStatesResultRunner: Control + run_slayer + HDF5 outputgpec.h5 schema namingReproducing. All commands run from a GPEC checkout of this branch.
Unit test:
The example:
julia --project=. -e 'using GeneralizedPerturbedEquilibrium; GeneralizedPerturbedEquilibrium.main(["examples/LAR_gal_cholesky_example"])'Δ′ is written to
SingularSurfaces/Delta_prime_matrixinexamples/LAR_gal_cholesky_example/gpec.h5.Regression report:
julia --project=regression-harness regression-harness/regress.jl \ --cases diiid_n1,gal_resistive_diiid,gal_resistive_pe,diiid_n1_riccati,lar_gal_cholesky \ --refs origin/develop,feature/galerkin-cholesky --forceFull table above (LU, Riccati and
gal_nxvariants, the DIII-D refusal, and thedevelopfailure). The attachedcholesky_pr_package.zip
contains
reproduce/:demo.jlandrun_demo.sh;The expected logs of the runs cited here are in
results_b5b430f/.New regression case.
lar_gal_choleskyrunsexamples/LAR_gal_cholesky_exampleand tracks:D_Iand α.On
developthis case has no example deck and Cholesky fails, so the harness reports nodevelopvalues for it. The report above records its first values.