From 7ffab7a3664a77db1f8cb9f95623bc7609e378dc Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 25 May 2023 15:20:32 -0600 Subject: [PATCH] Convert PETSc SNES interfaces to callback convention. --- include/deal.II/lac/petsc_snes.h | 78 ++++++++++++--- include/deal.II/lac/petsc_snes.templates.h | 106 ++++++++++++++++++--- 2 files changed, 160 insertions(+), 24 deletions(-) diff --git a/include/deal.II/lac/petsc_snes.h b/include/deal.II/lac/petsc_snes.h index 6167b0d9cd..3cf07e360b 100644 --- a/include/deal.II/lac/petsc_snes.h +++ b/include/deal.II/lac/petsc_snes.h @@ -29,6 +29,7 @@ # include +# include # if defined(DEAL_II_HAVE_CXX20) # include # endif @@ -350,14 +351,30 @@ namespace PETScWrappers /** * Callback for the computation of the nonlinear residual $F(x)$. + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. */ - std::function residual; + std::function residual; /** * Callback for the computation of the Jacobian * $\dfrac{\partial F}{\partial x}$. - */ - std::function + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. + */ + std::function jacobian; /** @@ -366,10 +383,18 @@ namespace PETScWrappers * This function is called by NonlinearSolver at the beginning * of each time step. Input arguments are the current step number * and the current value for ||F(x)||. - */ - std::function + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. + */ + std::function monitor; /** @@ -379,16 +404,32 @@ namespace PETScWrappers * operator $\dfrac{\partial F}{\partial x}$. * * Solvers must be provided via NonlinearSolver::solve_with_jacobian. + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. */ - std::function setup_jacobian; + std::function setup_jacobian; /** * Callback for the solution of the tangent system set up with * NonlinearSolver::setup_jacobian. * * This is used as a preconditioner inside the Krylov process. - */ - std::function + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. + */ + std::function solve_with_jacobian; /** @@ -402,8 +443,16 @@ namespace PETScWrappers * the reduction in a trust region step. * * The value of the energy function must be returned in @p energy_value. + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. PETSc's SNES package can not deal + * with "recoverable" errors, so if a callback + * throws an exception of type RecoverableUserCallbackError, then this + * exception is treated like any other exception. */ - std::function energy; + std::function energy; private: /** @@ -418,6 +467,13 @@ namespace PETScWrappers * This flag is used to support versions of PETSc older than 3.13. */ bool need_dummy_assemble; + + /** + * A pointer to any exception that may have been thrown in user-defined + * call-backs and that we have to deal after the KINSOL function we call + * has returned. + */ + mutable std::exception_ptr pending_exception; }; } // namespace PETScWrappers diff --git a/include/deal.II/lac/petsc_snes.templates.h b/include/deal.II/lac/petsc_snes.templates.h index adac937aec..db5d8d7965 100644 --- a/include/deal.II/lac/petsc_snes.templates.h +++ b/include/deal.II/lac/petsc_snes.templates.h @@ -50,12 +50,54 @@ DEAL_II_NAMESPACE_OPEN namespace PETScWrappers { + namespace + { + /** + * A function that calls the function object given by its first argument + * with the set of arguments following at the end. If the call returns + * regularly, the current function returns zero to indicate success. If + * the call fails with an exception, then the current function returns with + * an error code of -1. In that case, the exception thrown by `f` is + * captured and `eptr` is set to the exception. In case of success, + * `eptr` is set to `nullptr`. + */ + template + int + call_and_possibly_capture_exception(const F & f, + std::exception_ptr &eptr, + Args &&...args) + { + // See whether there is already something in the exception pointer + // variable. There is no reason why this should be so, and + // we should probably bail out: + AssertThrow(eptr == nullptr, ExcInternalError()); + + // Call the function and if that succeeds, return zero: + try + { + f(std::forward(args)...); + eptr = nullptr; + return 0; + } + // In case of an exception, capture the exception and + // return -1: + catch (...) + { + eptr = std::current_exception(); + return -1; + } + } + } // namespace + + + template DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type || std::constructible_from)) NonlinearSolver::NonlinearSolver( const NonlinearSolverData &data, const MPI_Comm mpi_comm) + : pending_exception(nullptr) { AssertPETSc(SNESCreate(mpi_comm, &snes)); AssertPETSc(SNESSetApplicationContext(snes, this)); @@ -101,6 +143,8 @@ namespace PETScWrappers NonlinearSolver::~NonlinearSolver() { AssertPETSc(SNESDestroy(&snes)); + + Assert(pending_exception == nullptr, ExcInternalError()); } @@ -201,9 +245,10 @@ namespace PETScWrappers VectorType xdealii(x); VectorType fdealii(f); - AssertUser(user->residual(xdealii, fdealii), "residual"); + const int err = call_and_possibly_capture_exception( + user->residual, user->pending_exception, xdealii, fdealii); petsc_increment_state_counter(f); - PetscFunctionReturn(0); + PetscFunctionReturn(err); }; const auto snes_jacobian = @@ -214,7 +259,9 @@ namespace PETScWrappers VectorType xdealii(x); AMatrixType Adealii(A); PMatrixType Pdealii(P); - AssertUser(user->jacobian(xdealii, Adealii, Pdealii), "jacobian"); + + const int err = call_and_possibly_capture_exception( + user->jacobian, user->pending_exception, xdealii, Adealii, Pdealii); petsc_increment_state_counter(P); // Handle the Jacobian-free case @@ -229,7 +276,7 @@ namespace PETScWrappers } else petsc_increment_state_counter(A); - PetscFunctionReturn(0); + PetscFunctionReturn(err); }; const auto snes_jacobian_with_setup = @@ -243,7 +290,10 @@ namespace PETScWrappers user->A = &Adealii; user->P = &Pdealii; - AssertUser(user->setup_jacobian(xdealii), "setup_jacobian"); + const int err = + call_and_possibly_capture_exception(user->setup_jacobian, + user->pending_exception, + xdealii); petsc_increment_state_counter(P); // Handle older versions of PETSc for which we cannot pass a MATSHELL @@ -267,7 +317,8 @@ namespace PETScWrappers } else petsc_increment_state_counter(A); - PetscFunctionReturn(0); + + PetscFunctionReturn(err); }; const auto snes_monitor = @@ -278,8 +329,9 @@ namespace PETScWrappers Vec x; AssertPETSc(SNESGetSolution(snes, &x)); VectorType xdealii(x); - AssertUser(user->monitor(xdealii, it, f), "monitor"); - PetscFunctionReturn(0); + const int err = call_and_possibly_capture_exception( + user->monitor, user->pending_exception, xdealii, it, f); + PetscFunctionReturn(err); }; const auto snes_objective = @@ -289,9 +341,10 @@ namespace PETScWrappers real_type v; VectorType xdealii(x); - AssertUser(user->energy(xdealii, v), "energy"); + const int err = call_and_possibly_capture_exception( + user->energy, user->pending_exception, xdealii, v); *f = v; - PetscFunctionReturn(0); + PetscFunctionReturn(err); }; AssertThrow(residual, @@ -365,7 +418,10 @@ namespace PETScWrappers precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int { VectorType dst(static_cast(indst)); const VectorType src(static_cast(insrc)); - return solve_with_jacobian(src, dst); + return call_and_possibly_capture_exception(solve_with_jacobian, + pending_exception, + src, + dst); }; // Default Krylov solver (preconditioner only) @@ -385,8 +441,32 @@ namespace PETScWrappers // Having set everything up, now do the actual work // and let PETSc solve the system. // Older versions of PETSc requires the solution vector specified even - // if we specified SNESSetSolution upfront - AssertPETSc(SNESSolve(snes, nullptr, x.petsc_vector())); + // if we specified SNESSetSolution upfront. + // + // If there is a pending exception, then one of the user callbacks + // threw an exception we didn't know how to deal with + // at the time. It is possible that PETSc managed to + // recover anyway and in that case would have returned + // a zero error code -- if so, just eat the exception and + // continue on; otherwise, just rethrow the exception + // and get outta here. + const int status = SNESSolve(snes, nullptr, x.petsc_vector()); + if (pending_exception) + { + try + { + std::rethrow_exception(pending_exception); + } + catch (...) + { + pending_exception = nullptr; + if (status == 0) + /* just eat the exception */; + else + throw; + } + } + AssertPETSc(status); // Get the number of steps taken. PetscInt nt; -- 2.39.5