]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert PETSc SNES interfaces to callback convention.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 25 May 2023 21:20:32 +0000 (15:20 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 25 May 2023 21:20:32 +0000 (15:20 -0600)
include/deal.II/lac/petsc_snes.h
include/deal.II/lac/petsc_snes.templates.h

index 6167b0d9cdc782959b6ae48158a4ac86cdc17680..3cf07e360b5a93018ce9ced5e0b782e36e894cf9 100644 (file)
@@ -29,6 +29,7 @@
 
 #  include <petscsnes.h>
 
+#  include <exception>
 #  if defined(DEAL_II_HAVE_CXX20)
 #    include <concepts>
 #  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<int(const VectorType &x, VectorType &res)> residual;
+    std::function<void(const VectorType &x, VectorType &res)> residual;
 
     /**
      * Callback for the computation of the Jacobian
      * $\dfrac{\partial F}{\partial x}$.
-     */
-    std::function<int(const VectorType &x, AMatrixType &A, PMatrixType &P)>
+     *
+     * @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<void(const VectorType &x, AMatrixType &A, PMatrixType &P)>
       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<int(const VectorType & x,
-                      const unsigned int step_number,
-                      const real_type    f_norm)>
+     *
+     * @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<void(const VectorType & x,
+                       const unsigned int step_number,
+                       const real_type    f_norm)>
       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<int(const VectorType &x)> setup_jacobian;
+    std::function<void(const VectorType &x)> 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<int(const VectorType &src, VectorType &dst)>
+     *
+     * @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<void(const VectorType &src, VectorType &dst)>
       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<int(const VectorType &x, real_type &energy_value)> energy;
+    std::function<void(const VectorType &x, real_type &energy_value)> 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
index adac937aecc8b62770b6adb943e18754965cad02..db5d8d7965c9da9607552e70047955e2ed62d47c 100644 (file)
@@ -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 <typename F, typename... Args>
+    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>(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 <typename VectorType, typename PMatrixType, typename AMatrixType>
   DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
                           std::constructible_from<VectorType, Vec>))
   NonlinearSolver<VectorType, PMatrixType, AMatrixType>::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<VectorType, PMatrixType, AMatrixType>::~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<const Vec &>(indst));
           const VectorType src(static_cast<const Vec &>(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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.