Repository navigation
Keep the matrix-free operator on a retained SNES - #4571
Lewin-Shen wants to merge 2 commits into
Conversation
reuse_preconditioner() keeps the SNES alive between solves. PetscNonlinearSolver::solve() registers the SNES matrices on every solve (the single-matrix overload passes the preconditioning matrix as both Amat and Pmat), and on a retained SNES the SNESSetUp() that follows returns at its early exit, so the matrix-free operator that -snes_mf_operator or -snes_mf asked for is never re-created. Replacing the SNES matrices drops the SNES's reference to it, and the options step that follows re-sets the KSP's operators and drops the last one, so the operator is destroyed: from the second solve on, a matrix-free Newton-Krylov solve is silently a Newton iteration on the registered Amat, which is the preconditioning matrix unless a separate operator matrix was passed. Under -snes_mf the Jacobian callback PETSc owns is lost as well, because init() attaches a new DM every solve and SNESSetDM() copies the old DMSNES only onto a DM that has none. When the retained SNES already holds a MATMFFD operator, keep it: pass a null Amat and refresh only the preconditioning matrix and the callback, or, for the ordinary -snes_mf state (one matrix-free matrix as both Amat and Pmat), pass null for both matrices and re-register MatMFFDComputeJacobian (after SNESSetFromOptions(), which may install a -snes_fd or -snes_fd_color callback) and PCNONE. The PCNONE re-application is compiled only for PETSc releases before 3.26: from 3.26 on (PETSc commit f4330cacc5, "Error with -snes_mf and incompatible pc") the matrix-free setup no longer forces PCNONE, SNESSetFromOptions() sets it only when no PC type was chosen, and the type the first solve ended with stays on the retained SNES. A separate operator matrix passed to the two-matrix solve() is not installed on that path either: under -snes_mf_operator the first solve's SNESSetUp() replaced it with the matrix-free operator, which later solves now keep. An explicit change of the matrix-free options read by SNESSetFromOptions() on a solver that retained a matrix-free operator is rejected with an error naming force_new_preconditioner(), unless the SNES has a left nonlinear preconditioner, whose setup rewrites the flags that check compares; callers who change the mode programmatically or from a configuration hook must call force_new_preconditioner() first. The change is a no-op on the first solve of every solve type, whenever the solver is not matrix-free, and whenever the solver is not retained. Adds tests/systems/nonlinear_implicit_system_reuse_test.C, a unit test with seven two-solve cases: the operator type and matrix identity under each matrix-free option, a residual-only -snes_mf solve that asserts convergence, the solution and PCNONE (with -pc_type ilu before PETSc 3.26, where the first solve's setup overrides it), a SolverConfiguration that installs a new shell preconditioner on every solve and a new KSP on the second, a rejected matrix-free mode change, -snes_mf beside -snes_fd, and unchanged options under a left nonlinear preconditioner. Fixes libMesh#4570. Co-Authored-By: Claude <noreply@anthropic.com>
Regenerate tests/Makefile.in with the autotools the tree ships (./bootstrap --build-autotools: autoconf 2.72, automake 1.18.1, libtool 2.5.4) after the previous commit added tests/systems/nonlinear_implicit_system_reuse_test.C to tests/Makefile.am. Only the lines for that source are taken: a full bootstrap of devel today also adds the meshes/quad4_tri3_hourglass.xda.gz line that tests/Makefile.am has carried since libMesh#4544 without a re-bootstrap, which is unrelated to this change and left out here. Co-Authored-By: Claude <noreply@anthropic.com>
|
I'm not seeing any problems here, but I'm just not confident enough in my PETSc-ese to be sure there's not something wrong that I'm missing. If @lindsayad or @jwpeterson find time to review this I'll defer to either of their judgements; if not then I'll probably merge next week. |
|
Can you have the LLM shorten the PR body? If you're using Claude Code and have access to Opus 5.5 I highly recommend it. It's the first LLM I've seen that writes prose resembling what a human would write. @roystgnr is more patient than I am; my eyes glazed over by the third row of your table. Intuitively your bug report makes sense to me |
lindsayad
left a comment
There was a problem hiding this comment.
I appreciate the contribution; it's a genuine bug fix so a very important very first contribution! Some of my comments may seem harsh but it's only because of aggregate frustration with reviewing LLM-generated codecomments. Also I'm not guilt-free of submitting bad LLM-written comments to reviewers, so I'm no better than you are!
| LibmeshPetscCall(SNESSetFunction (_snes, r->vec(), libmesh_petsc_snes_residual, this)); | ||
|
|
||
| // A SNES kept across solves (reuse_preconditioner()) is already set up, so the | ||
| // SNESSetUp() below returns at its early exit and never redoes the matrix-free setup |
There was a problem hiding this comment.
| // SNESSetUp() below returns at its early exit and never redoes the matrix-free setup | |
| // SNESSetUp() below returns early and never redoes the matrix-free setup |
|
|
||
| // A SNES kept across solves (reuse_preconditioner()) is already set up, so the | ||
| // SNESSetUp() below returns at its early exit and never redoes the matrix-free setup | ||
| // -snes_mf_operator or -snes_mf asked for: its operator and its callback are preserved here. |
There was a problem hiding this comment.
| // -snes_mf_operator or -snes_mf asked for: its operator and its callback are preserved here. | |
| // -snes_mf_operator or -snes_mf asked for. Consequently we must preserve the operator and its callback here. |
| bool retained_mf_operator = false; // -snes_mf_operator: the operator is PETSc's, Pmat ours | ||
| bool retained_mf = false; // -snes_mf: operator, Pmat and callback are PETSc's | ||
|
|
||
| #if !PETSC_VERSION_LESS_THAN(3,8,0) // SNESGetUseMatrixFree() exists from PETSc 3.8 on |
There was a problem hiding this comment.
I don't think we need to clutter the code with versions this old
| if (amat && (mf_operator || mf)) | ||
| LibmeshPetscCall(PetscObjectTypeCompare((PetscObject)amat, MATMFFD, &amat_is_mf)); | ||
|
|
||
| // One matrix-free matrix as both Amat and Pmat is the signature of PETSc's ordinary |
There was a problem hiding this comment.
| // One matrix-free matrix as both Amat and Pmat is the signature of PETSc's ordinary | |
| // One matrix-free matrix as both Amat and Pmat is the signature of PETSc's |
| retained_mf = true; | ||
| else if (amat_is_mf) | ||
| retained_mf_operator = true; | ||
| // reinit() calls force_new_preconditioner(), so no retained operator outlives a mesh change. |
There was a problem hiding this comment.
| // reinit() calls force_new_preconditioner(), so no retained operator outlives a mesh change. | |
| // reinit() calls force_new_preconditioner(), so a mesh change always discards operators |
| if (saved.name == name) | ||
| recorded = true; |
| // An identity apply is enough: the system is diagonal. | ||
| PetscErrorCode shell_pc_apply(PC, Vec x, Vec y) { return VecCopy(x, y); } |
There was a problem hiding this comment.
A diagonal system can have strongly varying elements such that the identity is quite a bad preconditioner relative to Jacobi, so I don't understand this comment. This is an instance where an unnecessary comment detracts from what should be a simple unit test that should have little to do with numerics
| void tearDown() {} | ||
|
|
||
| // Nine EDGE2 elements with one unknown each: the problem is diagonal and needs nine | ||
| // unknowns, not geometry, so the cases run on a libMesh built for one dimension. |
| static bool passes_everywhere(bool passes_here) | ||
| { | ||
| unsigned int failed = passes_here ? 0 : 1; | ||
| TestCommWorld->max(failed); | ||
| return failed == 0; |
There was a problem hiding this comment.
| static bool passes_everywhere(bool passes_here) | |
| { | |
| unsigned int failed = passes_here ? 0 : 1; | |
| TestCommWorld->max(failed); | |
| return failed == 0; | |
| static bool passes_everywhere(bool passes) | |
| { | |
| TestCommWorld->min(passes); | |
| return passes; |
Yo dawg I heard you like double negatives
| CPPUNIT_ASSERT_MESSAGE("the solve did not converge", | ||
| passes_everywhere(reason > 0)); |
There was a problem hiding this comment.
so maybe having a non-trivial (non-linear) function is important for the unit testing? That should be motivated up by the introduction of the cubic residual functor
I apologize for these unreadable LLM-written stuff without sufficient human review. I will review and refactor all comments and related code in the next few days |
|
Again, don't feel bad. Thanks so much for the contribution! |
Fixes #4570.
reuse_preconditioner()keeps the SNES alive between solves.PetscNonlinearSolver<T>::solve()registers the SNES matrices on every solve (the single-matrix overload passes the assembled preconditioning matrix as both Amat and Pmat), and on a retained SNES theSNESSetUp()that follows returns at its early exit, so the matrix-free operator that-snes_mf_operatoror-snes_mfasked for is never re-created: from the second solve on a matrix-free Newton-Krylov solve is silently a Newton iteration on the registered Amat, and under-snes_mfthe Jacobian callback PETSc owns is lost as well (the new DM attached on every solve carries its own DMSNES, soSNESSetDMdoes not copy the old one), which crashes inattach_dof_mapwhen only a residual is registered.When the retained SNES already holds a
MATMFFDoperator, keep it:solve()SNESSetJacobian(snes, jac->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this), runs unchanged when a Jacobian callback is registered, and residual-only solves register nothing, as before-snes_mf_operator(Amatmffd, Pmat assembled)SNESSetJacobian(snes, NULL, pre->mat(), libmesh_petsc_snes_jacobian, this): a null Amat leaves the operator untouched; Pmat and the callback are refreshed (the callback because the new DM has a new DMSNES). An operator matrix passed to the two-matrixsolve()is not installed on this path: under-snes_mf_operatorthe first solve'sSNESSetUp()replaced it with the matrix-free operator, and later solves keep that operator-snes_mf(onemffdas both Amat and Pmat)SNESSetJacobian(snes, NULL, NULL, MatMFFDComputeJacobian, NULL)afterSNESSetFromOptions, so that a-snes_fdor-snes_fd_colorcallback the options step installs is superseded again, as the first solve's setup superseded it; and afterSNESSetUp()thePCNONEthat PETSc forces for-snes_mfis re-applied (withPCSHELL, andPCH2OPUSfrom PETSc 3.16, exempt as in PETSc), becauseSNESSetFromOptionsre-reads-pc_typeon every solve and the early-exitSNESSetUpno longer overrides it; the KSP and PC are fetched afteruser_presolveand theSolverConfigurationhook, which may have replaced them. The re-application is compiled only for PETSc releases before 3.26: from 3.26 on (PETSc commit f4330cacc5, "Error with -snes_mf and incompatible pc") the matrix-free setup no longer forcesPCNONE,SNESSetFromOptionssets it only when no PC type was chosen, and the type the first solve ended with stays on the retained SNES, so there is nothing to put backSNESSetFromOptions)libmesh_error_msg_ifwith a message namingforce_new_preconditioner(): a mode change cannot be honored on a retained, set-up SNES; a programmatic change (SNESSetUseMatrixFree, or from a hook) is not detected and needsforce_new_preconditioner()first, as the commit message states. The check is skipped under a left nonlinear preconditioner, where PETSc's setup rewrites the flags and an unchanged-snes_mf_operatorwould read as a changeVersion guards:
SNESGetUseMatrixFreeexists from PETSc 3.8 (below that the change compiles out and the old behavior stays);PCH2OPUSfrom 3.16 (checked on the release tags); thePCNONEre-application compiles only below 3.26, as the table says (PETSC_VERSION_LESS_THANis PETSc'sPETSC_VERSION_LT, false on every development build: a development snapshot older than f4330cacc5 therefore compiles the re-application out while its setup still forcesPCNONE, so development-build support assumes f4330cacc5 or later; 3.26 is the anticipated release boundary, and neithermainnor 3.26 was built here). Everything else exists at 3.6.Three things worth declaring: a retained operator keeps the
-mat_mffd_typeit was created with, so a differencing type changed between solves does not reach it; a mesh change goes throughNonlinearImplicitSystem::reinit(), which callsforce_new_preconditioner()and destroys the SNES, so no retained operator outlives the layout it was created for (no size guard is needed, and a rank-local one would be wrong); and the detection keys on theMATMFFDtype of PETSc's default matrix-free operator, so the deprecated-snes_mf_version 2operator, aMATSHELL, is left exactly as it is today.Test:
tests/systems/nonlinear_implicit_system_reuse_test.C(suiteNonlinearImplicitSystemReuseTest, next to #4537'snonlinear_implicit_system_test.C), seven two-solve cases on a diagonal cubic system (residual and Jacobian as objects: the matrix-free residual callback serves the residual objects, not the bare function pointer): the operator type and matrix identity after each solve under-snes_mf_operatorand-snes_mf; a residual-only-snes_mfcase asserting convergence, the root, andPCNONEafter each solve, with an explicit-pc_type ilubefore PETSc 3.26 (where the first solve's setup overrides it and only the re-application keeps the second solve off ILU) and an explicitly cleared-pc_typefrom 3.26 on; aSolverConfigurationthat installs a new shell PC on every solve and a new KSP on the second; the options-change rejection (underLIBMESH_ENABLE_EXCEPTIONS);-snes_mfbeside-snes_fdwith a Jacobian registered (the options step re-installs PETSc's finite-difference callback on every solve, and the second solve must end up matrix-free as the first did); and unchanged options under a left nonlinear preconditioner, on a residual-only system that converges at its initial iterate (compiled out on a debug PETSc, whoseDMCreateMatrix()check libMesh's DM does not satisfy, a pre-existing condition). Options the cases set or clear are restored on exit; every assertion is reduced over the communicator. With PETSc 3.25.4 (a release build) the seven cases pass on this branch atdevelcf92df06d9: each alone at 1 and 2 ranks, the suite at 1, 2, 4 and 10 ranks (ten ranks on a nine-element mesh), together with the sibling suite, inside the whole unit-test driver at 1 and 2 ranks, and in a debug (dbg) libMesh build; the first five cases also pass at90766057ef. On the unpatched library (at90766057ef, and atdevelwith only this test added) the first five cases fail their second-solve assertion or segfault. Removing any single repair (tried at90766057ef) is caught by at least one case, except two whose removal no case detects: theAmat == Pmatterm that distinguishes ordinary-snes_mffrom PETSc's left-nonlinear-preconditioner branch, and the refresh of Pmat in the operator-keeping arm (every case passes the same preconditioning matrix to both solves).tests/Makefile.amgains the file, and a second commit regeneratestests/Makefile.inwith the autotools the tree ships (./bootstrap --build-autotools: autoconf 2.72, automake 1.18.1, libtool 2.5.4), restricted to the lines for the new source. A full bootstrap ofdeveltoday also adds themeshes/quad4_tri3_hourglass.xda.gzline that #4544 put intests/Makefile.amwithout a re-bootstrap; it is unrelated to this change and left out, but I can add it here if you prefer. The unit-test driver builds from the regenerated file as committed.Application-level evidence (MOOSE on libMesh
90766057ef, 32 runs at 1 and 2 ranks: PJFNK with ilu, gamg, hypre and jacobi, JFNK, NEWTON): the patch turns every matrix-free reuse run from a diverged trajectory or a crash into one that matches its reuse-off reference (to 1e-10 or bit for bit), while the three existing NEWTON reuse regression cases keep reproducing their gold files byte for byte.Disclosure: this change was developed with LLM assistance (Anthropic's Claude, through Claude Code), which is why the commits carry the
Co-Authored-By: Claudetrailer. The code, the unit test, the build-system edits and this description were drafted with it and reviewed by me before posting; the mechanism was traced in the libMesh and PETSc sources, and every result quoted above comes from builds and runs on my machine (PETSc 3.25.4, 1 and 2 ranks). If the repository uses theAI Assistedlabel, please apply it.