Port SNES wrapper lessons to the Python bindings - #4445
Merged
Merged
Conversation
Ported two things learnt while writing the C++ SNESSolver wrapper: NonlinearProblem.set_update() gives a convenience wrapper for SNESSetUpdate matching SNESSolver::set_update, and the assemble_residual/assemble_jacobian docstrings now warn that the b/J/P_mat passed to the callback are not always the cached problem vectors/matrices (e.g. line search work vectors). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
jorgensd
approved these changes
Aug 24, 2026
pull Bot
pushed a commit
to gnikit/dolfinx
that referenced
this pull request
Aug 29, 2026
* Add nls::petsc::NonlinearProblem, a C++ wrapper for PETSc SNES Adapts C++ callables to the SNES residual and Jacobian callback interface, and owns the SNES object along with references to the vector and matrices assembled into. Configuration of the solve is left to the user via the options database or the wrapped SNES object. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add SNES update hook, propagate callback exceptions, drop const solve Expose SNESSetUpdate via set_update. As SNESSetUpdate takes no context argument, the problem is recovered in the callback from a PetscContainer composed on the SNES object. Exceptions thrown by a callback are stored and re-thrown by solve rather than being reduced to a PETSc error code. PETSc does not restore its state as the aborted solve unwinds, so this is documented as terminal for the problem. KrylovSolver::solve and NonlinearProblem::solve mutate the solver state, so neither is const. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use CHECK_ERROR and reinterpret_cast in the PETSc wrappers Replace the hand-written PETSc error checks in NonlinearProblem with the CHECK_ERROR macro used by la/petsc.cpp and fem/petsc.h, and replace C-style casts of PETSc handles with reinterpret_cast, here and in la/petsc.cpp. Complete the Doxygen comments on the NonlinearProblem accessors. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use nls::petsc::NonlinearProblem in the C++ hyperelasticity demo Replaces the hand-written Newton loop and its Krylov solver with a SNES driven by NonlinearProblem, configured through the PETSc options database under the problem's options prefix. The residual is now assembled into the vector that the solver passes to the callback, rather than into a vector owned by the demo class: PETSc line searches evaluate the residual in a work vector of their own, which is not the vector registered with set_F. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Clarify that callbacks assemble into the objects they are passed The vector and matrices registered with set_F and set_J define the layout of the objects the solver passes to the callbacks; they are not necessarily those objects. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add assemble_residual and assemble_jacobian for SNES callbacks C++ counterparts of the functions of the same name in dolfinx.fem.petsc, for use as the bodies of the nls::petsc::NonlinearProblem callbacks. They copy the current iterate into the solution function, then assemble the residual (with lifting) or the Jacobian and preconditioner into the objects the solver passes in. Use them in the hyperelasticity demo, which no longer spells out the zero/assemble/ghost-update/set_bc sequence by hand. It now converges in three Newton iterations rather than four, as the hand-written residual did not lift the Dirichlet conditions. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Drop the HyperElasticProblem class from the hyperelasticity demo With assembly in fem::petsc::assemble_residual/assemble_jacobian and the SNES owned by nls::petsc::NonlinearProblem, the class held nothing but the objects the callbacks need. Set the solver up in main instead. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Test PETSc reference counting in NonlinearProblem Check that the problem and the SNES each hold one reference to the residual vector and to each matrix, that re-setting with the same object does not accumulate references, that a Jacobian used as its own preconditioner is referenced twice and released twice, and that everything is released when the problem is destroyed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Explain the expected reference counts in the tests Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Document the lifetime contract of the SNES wrapping constructor The callbacks hold a pointer to the problem, which is not reference counted, so a SNES that outlives the problem carries dead callbacks. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix undefined behaviour when reading an unset PETSc options prefix PETSc reports an unset prefix as a null pointer, which was passed straight to the std::string constructor. Affects get_options_prefix on la::petsc::Vector, la::petsc::Matrix, la::petsc::KrylovSolver and nls::petsc::NonlinearProblem, and crashes rather than returning "". Point set_F and set_J at the fem::petsc functions that assemble the residual and Jacobian of a form. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Extend the NonlinearProblem tests Cover the update hook across a move, which reaches the problem through the container composed on the SNES rather than a context pointer; exceptions from the Jacobian callback and the update hook, each of which has its own trampoline; that the vector passed to solve is used as the initial guess, by converging to the negative root from a negative guess; and that a problem that has converged can be solved again. The move test re-registers the hook on the moved-to problem with a distinct target. A moved-from std::function is left in a valid but unspecified state and remains callable on libc++, so a hook sharing state with the moved-from one cannot show which problem was reached. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Only insert the unit diagonal for square forms Matches dolfinx.fem.petsc.assemble_matrix, which guards the flush and the diagonal insertion on the test and trial spaces being the same. It makes no difference for a Jacobian, but a preconditioner form need not be square. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Rename nls::petsc::NonlinearProblem to SNESSolver The class wraps a SNES and holds callbacks; it does not hold the forms, boundary conditions or solution function that make up a problem, so the Python class of the same name sets the wrong expectation. SNESSolver says what it is, and leaves NonlinearProblem free for a form-based class mirroring dolfinx.fem.petsc.NonlinearProblem. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Test SNESSolver with nest matrices and vectors The solver holds Vec and Mat, so MATNEST and VECNEST pass through it untouched. Solve two decoupled blocks with different roots, so that a mix-up between them would show up in the solution, preconditioned with fieldsplit as a nest matrix cannot be factored directly. Under MPI the blocks are MPIAIJ, and the suite already runs on one, two and three ranks. Also restore the test case description, which a rename had turned into "Solve nonlinear solver with SNES". Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Recover the solver from the residual context in the update hook SNESSetUpdate passes no context, which a PetscContainer composed on the SNES worked around. SNESGetFunction returns the context registered by set_F, which is the same pointer, so the container and its lifetime, naming and re-composition on move all go away. Guard the recovery on the residual callback being this class's own: a caller holding the SNES can register their own function, and its context must not be cast to a solver. Test that solving through the SNES object directly, rather than through solve(), works and still runs the update hook. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Name the demo vectors and matrices for what they hold The locals named x and b shadowed the callback parameters of the same name, which are the point the solver evaluates at and the object it assembles into. The registered objects set the layout of the ones the solver passes in, so they are named A_layout and b_layout, and the vector sharing the degrees-of-freedom of u is u_vec. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Hand edit docstrings and comments. * Correct the assemble_residual and assemble_jacobian documentation The assemble_jacobian example predated moving the preconditioner form to a trailing defaulted argument, and would not compile. The assemble_residual example registered a vector named b, shadowing the callback parameter of the same name. Say that the vector and matrices are the ones the solver passed to the callback rather than those registered with it, that x can be a line search trial point, and that u must be the function the forms hold as a coefficient. Qualify the unit diagonal, which is only set for forms whose test and trial spaces are the same. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Pass the evaluation point first in the assembly functions assemble_residual and assemble_jacobian took the object to assemble into before the point to evaluate at, the reverse of the callback signatures they are written for. Both arguments are Vec or Mat, so transposing them compiles and silently assembles into the solution vector. Match the callbacks, and dolfinx.fem.petsc, by taking x first. Also name the callback arguments in the documentation, in the order the solver passes them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Distinguish the registered objects from the callback arguments set_F and set_J named their vector and matrices b, Jmat and Pmat, the same as the parameters of the callbacks they take, so the documentation could not refer to either without ambiguity. Name the registered objects b_layout, J_layout and P_layout, as the hyperelasticity demo does, and name the callback arguments in the documentation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Explain the two re-throws in the exception handling Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use PetscInt for iteration counts and callback indices These come straight from PETSc, so they take its integer type rather than being narrowed to int on the way out: the value returned by KrylovSolver::solve and SNESSolver::solve, and the step index passed to the SNESSolver update hook. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Say which direction the PETSc index rule applies in The rule covered values passed into a PETSc call, leaving the type of a count returned from one, or handed to a callback, undecided. Name the nls::petsc wrapper alongside the la ones. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix the single precision and documentation build failures The tolerance on the computed root was a fixed 1e-6, which the float32 and complex64 builds miss by 2e-6. Take it from the solver instead: for f(x) = x^2 - c the root error is bounded by the residual norm SNES stopped at divided by |f'(root)|, so the check follows whatever precision and tolerances are in force. It is clamped below by round-off and above so that a solve stopping far from the root fails rather than widening its own tolerance. Doxygen matches a definition to its declaration including the parameter names inside a std::function type, so the set_F, set_J and set_update definitions now spell them as the header does. Leaving them out failed the documentation build, which treats the warning as an error. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Put options adjustment in a more natural place. * Check the return code of every raw PETSc call in snessolver.cpp Wrap every PETSc call, including inside the callback helper functions, with CHECK(... == 0) for strict error checking, matching the convention applied to la::petsc's own test suite. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017uXpk1iyjECr5Eei7Sq1dH * SonarCloud * clang-format * Use CHECK_ERROR consistently in la::petsc Convert remaining explicit if (ierr != 0) petsc::error(...) sites to the CHECK_ERROR macro, and add missing error checking to Matrix::set_options_prefix, get_options_prefix, and set_from_options. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Merge la::petsc::Operator into Matrix; remove apply and norm Operator had exactly one subclass (Matrix), and nothing consumed it polymorphically (KrylovSolver::set_operator(s) takes a raw Mat, not Operator&), so the base class was pure indirection. Fold its members (mat(), size(), create_vector(dim), Mat lifetime/ref-counting) directly into Matrix. Matrix::apply(AssemblyType) and Matrix::norm(Norm) are unused throughout the codebase (demos and Python assemble/norm operate on the raw PETSc Mat directly) and are dropped along with the now-unused Norm forward declaration. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Replace CHECK_ERROR macros with a check() helper using std::source_location Add la::petsc::check(ierr, petsc_function), a thin wrapper around error() that throws on a non-zero PETSc error code. Both take a std::source_location defaulted to the call site, so error() no longer needs a caller-supplied __FILE__ and reports the exact line, not just the file. Being a plain namespaced function rather than a macro, check() is usable from header-inline code, so it also replaces the private CHECK_ERROR macros duplicated across la/petsc.cpp, la/slepc.cpp, and fem/petsc.h, and the leftover explicit if (ierr != 0) sites in la/petsc.h and fem/petsc.cpp that those macros couldn't reach. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Replace CHECK_ERROR macros with a check() helper using std::source_location Add la::petsc::check(ierr, petsc_function), a thin wrapper around error() that throws on a non-zero PETSc error code. Both take a std::source_location defaulted to the call site, so error() no longer needs a caller-supplied __FILE__ and reports the exact line, not just the file. Being a plain namespaced function rather than a macro, check() is usable from header-inline code, so it also replaces the private CHECK_ERROR macros duplicated across la/petsc.cpp, la/slepc.cpp, fem/petsc.h, and nls/SNESSolver.cpp, and the leftover explicit if (ierr != 0) sites in la/petsc.h and fem/petsc.cpp that those macros couldn't reach. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Fix CI: convert last la::petsc::error(ierr, __FILE__, ...) call site python/dolfinx/wrappers/petsc.cpp called the old 3-argument la::petsc::error(ierr, __FILE__, name) directly, missed when error()'s signature was changed to take a std::source_location. Every CI leg that builds the Python bindings failed on this one site. Switch it to la::petsc::check(ierr, name), matching every other call site. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Fix docstring conventions in la/petsc.h and fem/petsc.h Add @brief to every non-trivial doc comment (has @PARAM, @note, or @return) that was missing it, per AGENTS.md. Also fix Matrix's deleted copy-assignment operator, which used a Doxygen /// comment instead of the plain // required for deleted members (a pre-existing bug on both the old Operator and Matrix classes, carried over by the Operator merge), and a stray space in create_nullspace's @PARAM [in] tag. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Move PETSc error handling from la::petsc to common::petsc Consolidates error()/check() into a single common::petsc namespace so fem::petsc and la::petsc/slepc share one implementation instead of la owning it on fem's behalf. * Use inline common::petsc::check(PetscCall(...), name) style Replaces the two-line "PetscErrorCode ierr = PetscCall(...); check(ierr, name);" pattern with a single inline call across la::petsc, la::slepc, fem::petsc, and the petsc4py wrapper, dropping the intermediate ierr variable wherever it wasn't needed afterward. * Assert that common::petsc::error() is only called with a nonzero code error() is only reachable from check() after ierr != 0, so a caller passing 0 is a bug, not a condition to silently no-op on. * Move la::petsc::options into flat common::petsc functions Relocates the PETSc options-database helpers to common::petsc, since they are already used outside la (nls, demos) and only depend on common::petsc::check. Flattens options::set/clear into set_option/set_option<T>/clear_option/clear_options rather than keeping a third-level options namespace, matching the rest of the library's dolfinx::<module>::<submodule> convention. * Apply gersemi formatting to common/CMakeLists.txt * Update copyright headers for PETSc error-handling work Bump end year to 2026 across touched files, and add Jack S. Hale to files with substantial contributions in this branch. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Check remaining unchecked PETSc calls in demos, tests, and wrapper Bring the demos onto common::petsc::check for all Mat/Vec calls, wrap the previously bare slepc.cpp test assertions in CHECK(... == 0) to match cpp/test/petsc.cpp, and check the unchecked MatSetOption calls in the Python discrete-operator wrappers. Destructors and catch-block cleanup that destroy PETSc objects are now checked too rather than silently discarding the error code, in anticipation of --with-strict-petscerrorcode. A thrown error in these paths calls std::terminate instead of propagating, which each site now notes explicitly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Include call site and function name in PETSc error messages Merge the two spdlog::error calls into one and fold the same file, line, and enclosing function name into the thrown std::runtime_error message, not just the log line. Callers that only see e.what() (e.g. via the Python bindings) previously lost the call site entirely. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Check the remaining unchecked PETSc calls in the hyperelasticity demo VecGhostUpdateBegin/End and SNESGetLinearSolveIterations were left unchecked after the merge, unlike the other demos. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Make SNESSolver::solve() return void, not the iteration count The Python NonlinearProblem (fem/petsc.py), the SNES-based analogue of this wrapper, returns only the solution from solve() and leaves iteration counts and convergence status to be queried from the SNES object directly. Follow the same convention in C++: the iteration count is equally available via SNESGetIterationNumber(solver.snes(), ...), so returning it from solve() was redundant. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Revert "Make SNESSolver::solve() return void, not the iteration count" This reverts commit e9dd730. * Simplify SNESSolver's exception handling Three call sites duplicated the same try/catch-and-store-exception logic; factor it into one invoke() helper. store_exception() used to re-throw and re-catch its own exception_ptr just to read the message for logging -- unnecessary when each catch clause already has the exception by type. Also recover the solver in the update hook via the SNES application context (SNESSetApplicationContext/SNESGetApplicationContext) rather than by inspecting SNESGetFunction and checking it points at residual(). This is more direct and does not depend on set_F having been called. Restore the reinterpret_cast<PetscObject> casts that the more-petsc-work merge had turned into C-style casts. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * List captures explicitly in the SNES callback lambdas [&] captured everything in scope; list only what each lambda actually uses. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Throw on a null SNES in the wrapping constructor Match the KrylovSolver sibling constructor: an O(1) API-boundary check should throw rather than assert, so a null SNES fails cleanly in Release builds instead of proceeding into undefined behaviour. * Make invoke() a free function It only touches _exception, passed explicitly, so it does not need to be a SNESSolver member. This drops it from the class's private interface in the header. * Add SNES update-step hook and document callback vector/matrix aliasing (FEniCS#4445) --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
pull Bot
pushed a commit
to gnikit/dolfinx
that referenced
this pull request
Aug 30, 2026
…niCS#4448) * Add nls::petsc::NonlinearProblem, a C++ wrapper for PETSc SNES Adapts C++ callables to the SNES residual and Jacobian callback interface, and owns the SNES object along with references to the vector and matrices assembled into. Configuration of the solve is left to the user via the options database or the wrapped SNES object. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add SNES update hook, propagate callback exceptions, drop const solve Expose SNESSetUpdate via set_update. As SNESSetUpdate takes no context argument, the problem is recovered in the callback from a PetscContainer composed on the SNES object. Exceptions thrown by a callback are stored and re-thrown by solve rather than being reduced to a PETSc error code. PETSc does not restore its state as the aborted solve unwinds, so this is documented as terminal for the problem. KrylovSolver::solve and NonlinearProblem::solve mutate the solver state, so neither is const. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use CHECK_ERROR and reinterpret_cast in the PETSc wrappers Replace the hand-written PETSc error checks in NonlinearProblem with the CHECK_ERROR macro used by la/petsc.cpp and fem/petsc.h, and replace C-style casts of PETSc handles with reinterpret_cast, here and in la/petsc.cpp. Complete the Doxygen comments on the NonlinearProblem accessors. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use nls::petsc::NonlinearProblem in the C++ hyperelasticity demo Replaces the hand-written Newton loop and its Krylov solver with a SNES driven by NonlinearProblem, configured through the PETSc options database under the problem's options prefix. The residual is now assembled into the vector that the solver passes to the callback, rather than into a vector owned by the demo class: PETSc line searches evaluate the residual in a work vector of their own, which is not the vector registered with set_F. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Clarify that callbacks assemble into the objects they are passed The vector and matrices registered with set_F and set_J define the layout of the objects the solver passes to the callbacks; they are not necessarily those objects. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add assemble_residual and assemble_jacobian for SNES callbacks C++ counterparts of the functions of the same name in dolfinx.fem.petsc, for use as the bodies of the nls::petsc::NonlinearProblem callbacks. They copy the current iterate into the solution function, then assemble the residual (with lifting) or the Jacobian and preconditioner into the objects the solver passes in. Use them in the hyperelasticity demo, which no longer spells out the zero/assemble/ghost-update/set_bc sequence by hand. It now converges in three Newton iterations rather than four, as the hand-written residual did not lift the Dirichlet conditions. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Drop the HyperElasticProblem class from the hyperelasticity demo With assembly in fem::petsc::assemble_residual/assemble_jacobian and the SNES owned by nls::petsc::NonlinearProblem, the class held nothing but the objects the callbacks need. Set the solver up in main instead. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Test PETSc reference counting in NonlinearProblem Check that the problem and the SNES each hold one reference to the residual vector and to each matrix, that re-setting with the same object does not accumulate references, that a Jacobian used as its own preconditioner is referenced twice and released twice, and that everything is released when the problem is destroyed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Explain the expected reference counts in the tests Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Document the lifetime contract of the SNES wrapping constructor The callbacks hold a pointer to the problem, which is not reference counted, so a SNES that outlives the problem carries dead callbacks. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix undefined behaviour when reading an unset PETSc options prefix PETSc reports an unset prefix as a null pointer, which was passed straight to the std::string constructor. Affects get_options_prefix on la::petsc::Vector, la::petsc::Matrix, la::petsc::KrylovSolver and nls::petsc::NonlinearProblem, and crashes rather than returning "". Point set_F and set_J at the fem::petsc functions that assemble the residual and Jacobian of a form. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Extend the NonlinearProblem tests Cover the update hook across a move, which reaches the problem through the container composed on the SNES rather than a context pointer; exceptions from the Jacobian callback and the update hook, each of which has its own trampoline; that the vector passed to solve is used as the initial guess, by converging to the negative root from a negative guess; and that a problem that has converged can be solved again. The move test re-registers the hook on the moved-to problem with a distinct target. A moved-from std::function is left in a valid but unspecified state and remains callable on libc++, so a hook sharing state with the moved-from one cannot show which problem was reached. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Only insert the unit diagonal for square forms Matches dolfinx.fem.petsc.assemble_matrix, which guards the flush and the diagonal insertion on the test and trial spaces being the same. It makes no difference for a Jacobian, but a preconditioner form need not be square. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Rename nls::petsc::NonlinearProblem to SNESSolver The class wraps a SNES and holds callbacks; it does not hold the forms, boundary conditions or solution function that make up a problem, so the Python class of the same name sets the wrong expectation. SNESSolver says what it is, and leaves NonlinearProblem free for a form-based class mirroring dolfinx.fem.petsc.NonlinearProblem. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Test SNESSolver with nest matrices and vectors The solver holds Vec and Mat, so MATNEST and VECNEST pass through it untouched. Solve two decoupled blocks with different roots, so that a mix-up between them would show up in the solution, preconditioned with fieldsplit as a nest matrix cannot be factored directly. Under MPI the blocks are MPIAIJ, and the suite already runs on one, two and three ranks. Also restore the test case description, which a rename had turned into "Solve nonlinear solver with SNES". Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Recover the solver from the residual context in the update hook SNESSetUpdate passes no context, which a PetscContainer composed on the SNES worked around. SNESGetFunction returns the context registered by set_F, which is the same pointer, so the container and its lifetime, naming and re-composition on move all go away. Guard the recovery on the residual callback being this class's own: a caller holding the SNES can register their own function, and its context must not be cast to a solver. Test that solving through the SNES object directly, rather than through solve(), works and still runs the update hook. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Name the demo vectors and matrices for what they hold The locals named x and b shadowed the callback parameters of the same name, which are the point the solver evaluates at and the object it assembles into. The registered objects set the layout of the ones the solver passes in, so they are named A_layout and b_layout, and the vector sharing the degrees-of-freedom of u is u_vec. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Hand edit docstrings and comments. * Correct the assemble_residual and assemble_jacobian documentation The assemble_jacobian example predated moving the preconditioner form to a trailing defaulted argument, and would not compile. The assemble_residual example registered a vector named b, shadowing the callback parameter of the same name. Say that the vector and matrices are the ones the solver passed to the callback rather than those registered with it, that x can be a line search trial point, and that u must be the function the forms hold as a coefficient. Qualify the unit diagonal, which is only set for forms whose test and trial spaces are the same. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Pass the evaluation point first in the assembly functions assemble_residual and assemble_jacobian took the object to assemble into before the point to evaluate at, the reverse of the callback signatures they are written for. Both arguments are Vec or Mat, so transposing them compiles and silently assembles into the solution vector. Match the callbacks, and dolfinx.fem.petsc, by taking x first. Also name the callback arguments in the documentation, in the order the solver passes them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Distinguish the registered objects from the callback arguments set_F and set_J named their vector and matrices b, Jmat and Pmat, the same as the parameters of the callbacks they take, so the documentation could not refer to either without ambiguity. Name the registered objects b_layout, J_layout and P_layout, as the hyperelasticity demo does, and name the callback arguments in the documentation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Explain the two re-throws in the exception handling Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Use PetscInt for iteration counts and callback indices These come straight from PETSc, so they take its integer type rather than being narrowed to int on the way out: the value returned by KrylovSolver::solve and SNESSolver::solve, and the step index passed to the SNESSolver update hook. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Say which direction the PETSc index rule applies in The rule covered values passed into a PETSc call, leaving the type of a count returned from one, or handed to a callback, undecided. Name the nls::petsc wrapper alongside the la ones. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Fix the single precision and documentation build failures The tolerance on the computed root was a fixed 1e-6, which the float32 and complex64 builds miss by 2e-6. Take it from the solver instead: for f(x) = x^2 - c the root error is bounded by the residual norm SNES stopped at divided by |f'(root)|, so the check follows whatever precision and tolerances are in force. It is clamped below by round-off and above so that a solve stopping far from the root fails rather than widening its own tolerance. Doxygen matches a definition to its declaration including the parameter names inside a std::function type, so the set_F, set_J and set_update definitions now spell them as the header does. Leaving them out failed the documentation build, which treats the warning as an error. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Put options adjustment in a more natural place. * Check the return code of every raw PETSc call in snessolver.cpp Wrap every PETSc call, including inside the callback helper functions, with CHECK(... == 0) for strict error checking, matching the convention applied to la::petsc's own test suite. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017uXpk1iyjECr5Eei7Sq1dH * SonarCloud * clang-format * Use CHECK_ERROR consistently in la::petsc Convert remaining explicit if (ierr != 0) petsc::error(...) sites to the CHECK_ERROR macro, and add missing error checking to Matrix::set_options_prefix, get_options_prefix, and set_from_options. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Merge la::petsc::Operator into Matrix; remove apply and norm Operator had exactly one subclass (Matrix), and nothing consumed it polymorphically (KrylovSolver::set_operator(s) takes a raw Mat, not Operator&), so the base class was pure indirection. Fold its members (mat(), size(), create_vector(dim), Mat lifetime/ref-counting) directly into Matrix. Matrix::apply(AssemblyType) and Matrix::norm(Norm) are unused throughout the codebase (demos and Python assemble/norm operate on the raw PETSc Mat directly) and are dropped along with the now-unused Norm forward declaration. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Replace CHECK_ERROR macros with a check() helper using std::source_location Add la::petsc::check(ierr, petsc_function), a thin wrapper around error() that throws on a non-zero PETSc error code. Both take a std::source_location defaulted to the call site, so error() no longer needs a caller-supplied __FILE__ and reports the exact line, not just the file. Being a plain namespaced function rather than a macro, check() is usable from header-inline code, so it also replaces the private CHECK_ERROR macros duplicated across la/petsc.cpp, la/slepc.cpp, and fem/petsc.h, and the leftover explicit if (ierr != 0) sites in la/petsc.h and fem/petsc.cpp that those macros couldn't reach. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Replace CHECK_ERROR macros with a check() helper using std::source_location Add la::petsc::check(ierr, petsc_function), a thin wrapper around error() that throws on a non-zero PETSc error code. Both take a std::source_location defaulted to the call site, so error() no longer needs a caller-supplied __FILE__ and reports the exact line, not just the file. Being a plain namespaced function rather than a macro, check() is usable from header-inline code, so it also replaces the private CHECK_ERROR macros duplicated across la/petsc.cpp, la/slepc.cpp, fem/petsc.h, and nls/SNESSolver.cpp, and the leftover explicit if (ierr != 0) sites in la/petsc.h and fem/petsc.cpp that those macros couldn't reach. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Fix CI: convert last la::petsc::error(ierr, __FILE__, ...) call site python/dolfinx/wrappers/petsc.cpp called the old 3-argument la::petsc::error(ierr, __FILE__, name) directly, missed when error()'s signature was changed to take a std::source_location. Every CI leg that builds the Python bindings failed on this one site. Switch it to la::petsc::check(ierr, name), matching every other call site. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Fix docstring conventions in la/petsc.h and fem/petsc.h Add @brief to every non-trivial doc comment (has @PARAM, @note, or @return) that was missing it, per AGENTS.md. Also fix Matrix's deleted copy-assignment operator, which used a Doxygen /// comment instead of the plain // required for deleted members (a pre-existing bug on both the old Operator and Matrix classes, carried over by the Operator merge), and a stray space in create_nullspace's @PARAM [in] tag. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Move PETSc error handling from la::petsc to common::petsc Consolidates error()/check() into a single common::petsc namespace so fem::petsc and la::petsc/slepc share one implementation instead of la owning it on fem's behalf. * Use inline common::petsc::check(PetscCall(...), name) style Replaces the two-line "PetscErrorCode ierr = PetscCall(...); check(ierr, name);" pattern with a single inline call across la::petsc, la::slepc, fem::petsc, and the petsc4py wrapper, dropping the intermediate ierr variable wherever it wasn't needed afterward. * Assert that common::petsc::error() is only called with a nonzero code error() is only reachable from check() after ierr != 0, so a caller passing 0 is a bug, not a condition to silently no-op on. * Move la::petsc::options into flat common::petsc functions Relocates the PETSc options-database helpers to common::petsc, since they are already used outside la (nls, demos) and only depend on common::petsc::check. Flattens options::set/clear into set_option/set_option<T>/clear_option/clear_options rather than keeping a third-level options namespace, matching the rest of the library's dolfinx::<module>::<submodule> convention. * Apply gersemi formatting to common/CMakeLists.txt * Update copyright headers for PETSc error-handling work Bump end year to 2026 across touched files, and add Jack S. Hale to files with substantial contributions in this branch. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Check remaining unchecked PETSc calls in demos, tests, and wrapper Bring the demos onto common::petsc::check for all Mat/Vec calls, wrap the previously bare slepc.cpp test assertions in CHECK(... == 0) to match cpp/test/petsc.cpp, and check the unchecked MatSetOption calls in the Python discrete-operator wrappers. Destructors and catch-block cleanup that destroy PETSc objects are now checked too rather than silently discarding the error code, in anticipation of --with-strict-petscerrorcode. A thrown error in these paths calls std::terminate instead of propagating, which each site now notes explicitly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Include call site and function name in PETSc error messages Merge the two spdlog::error calls into one and fold the same file, line, and enclosing function name into the thrown std::runtime_error message, not just the log line. Callers that only see e.what() (e.g. via the Python bindings) previously lost the call site entirely. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Check the remaining unchecked PETSc calls in the hyperelasticity demo VecGhostUpdateBegin/End and SNESGetLinearSolveIterations were left unchecked after the merge, unlike the other demos. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Make SNESSolver::solve() return void, not the iteration count The Python NonlinearProblem (fem/petsc.py), the SNES-based analogue of this wrapper, returns only the solution from solve() and leaves iteration counts and convergence status to be queried from the SNES object directly. Follow the same convention in C++: the iteration count is equally available via SNESGetIterationNumber(solver.snes(), ...), so returning it from solve() was redundant. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Revert "Make SNESSolver::solve() return void, not the iteration count" This reverts commit e9dd730. * Simplify SNESSolver's exception handling Three call sites duplicated the same try/catch-and-store-exception logic; factor it into one invoke() helper. store_exception() used to re-throw and re-catch its own exception_ptr just to read the message for logging -- unnecessary when each catch clause already has the exception by type. Also recover the solver in the update hook via the SNES application context (SNESSetApplicationContext/SNESGetApplicationContext) rather than by inspecting SNESGetFunction and checking it points at residual(). This is more direct and does not depend on set_F having been called. Restore the reinterpret_cast<PetscObject> casts that the more-petsc-work merge had turned into C-style casts. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * List captures explicitly in the SNES callback lambdas [&] captured everything in scope; list only what each lambda actually uses. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Throw on a null SNES in the wrapping constructor Match the KrylovSolver sibling constructor: an O(1) API-boundary check should throw rather than assert, so a null SNES fails cleanly in Release builds instead of proceeding into undefined behaviour. * Make invoke() a free function It only touches _exception, passed explicitly, so it does not need to be a SNESSolver member. This drops it from the class's private interface in the header. * Add SNES update-step hook and document callback vector/matrix aliasing (FEniCS#4445) * Make SNES/KSP/EPS solve() return [[nodiscard]] convergence reason Force callers to check the outcome of la::petsc::KrylovSolver::solve, nls::petsc::SNESSolver::solve, and la::SLEPcEigenSolver::solve rather than silently ignoring non-convergence. Iteration counts are no longer returned; recover them via the raw PETSc API (SNESGetIterationNumber/KSPGetIterationNumber) when needed, as shown in the updated demos. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * Give [[nodiscard]] on solve() a reason message State the positive/negative convergence-reason rule directly in the nodiscard message so it surfaces in compiler diagnostics. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com> Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
NonlinearProblem.set_update(), a thin wrapper aroundSNES.setUpdatemirroringnls::petsc::SNESSolver::set_updateadded in the C++ SNES wrapper on Add nls::petsc::SNESSolver, a C++ wrapper for PETSc SNES #4433.assemble_residual/assemble_jacobianthat theb/J/P_matpassed to the callback are not always the cachedproblem.b/problem.A/problem.P_mat(e.g. line-search work vectors) — the same aliasing footgun already called out inSNESSolver.h.🤖 Generated with Claude Code