]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert KINSOL callbacks to use exceptions.
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 2 May 2023 23:49:57 +0000 (17:49 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 4 May 2023 23:43:06 +0000 (17:43 -0600)
include/deal.II/sundials/kinsol.h
source/sundials/kinsol.cc

index ee17684d4e7a8b2205e15994d735d70dfd099620..df4177b5207662c7f8520de88d4c2dbf22967a0a 100644 (file)
@@ -38,6 +38,7 @@
 #  include <sundials/sundials_math.h>
 #  include <sundials/sundials_types.h>
 
+#  include <exception>
 #  include <memory>
 
 
@@ -430,7 +431,7 @@ namespace SUNDIALS
      * - <0: Unrecoverable error the computation will be aborted and an
      * assertion will be thrown.
      */
-    std::function<int(const VectorType &src, VectorType &dst)> residual;
+    std::function<void(const VectorType &src, VectorType &dst)> residual;
 
     /**
      * A function object that users should supply and that is intended to
@@ -445,7 +446,7 @@ namespace SUNDIALS
      * - <0: Unrecoverable error; the computation will be aborted and an
      * assertion will be thrown.
      */
-    std::function<int(const VectorType &src, VectorType &dst)>
+    std::function<void(const VectorType &src, VectorType &dst)>
       iteration_function;
 
     /**
@@ -493,7 +494,8 @@ namespace SUNDIALS
      * - <0: Unrecoverable error the computation will be aborted and an
      * assertion will be thrown.
      */
-    std::function<int(const VectorType &current_u, const VectorType &current_f)>
+    std::function<void(const VectorType &current_u,
+                       const VectorType &current_f)>
       setup_jacobian;
 
     /**
@@ -555,10 +557,10 @@ namespace SUNDIALS
      *   that the Jacobian should be re-computed in every iteration.
      */
     DEAL_II_DEPRECATED
-    std::function<int(const VectorType &ycur,
-                      const VectorType &fcur,
-                      const VectorType &rhs,
-                      VectorType &      dst)>
+    std::function<void(const VectorType &ycur,
+                       const VectorType &fcur,
+                       const VectorType &rhs,
+                       VectorType &      dst)>
       solve_jacobian_system;
 
     /**
@@ -601,7 +603,7 @@ namespace SUNDIALS
      * assertion will be thrown.
      */
     std::function<
-      int(const VectorType &rhs, VectorType &dst, const double tolerance)>
+      void(const VectorType &rhs, VectorType &dst, const double tolerance)>
       solve_with_jacobian;
 
     /**
@@ -670,6 +672,13 @@ namespace SUNDIALS
                    << ". Please consult SUNDIALS manual.");
 
 
+    /**
+     * 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;
+
   private:
     /**
      * Throw an exception when a function with the given name is not
index fc1d9adbd53071b0849d0540b869f85b27498147..13a25fc3bb09742309576a78740dd8ca18466b4f 100644 (file)
@@ -52,6 +52,77 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace SUNDIALS
 {
+  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 of type RecoverableUserCallbackError, then
+     * the current function returns 1 to indicate that the called function
+     * object thought the error it encountered is recoverable. If the call fails
+     * with any other exception, then the current function returns with an error
+     * code of -1. In each of the last two cases, 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. This can only happen if we had previously had
+      // a recoverable exception, and the underlying library actually
+      // did recover successfully. In that case, we can abandon the
+      // exception previously thrown. If eptr contains anything other,
+      // then we really don't know how that could have happened, and
+      // should probably bail out:
+      if (eptr)
+        {
+          try
+            {
+              std::rethrow_exception(eptr);
+            }
+          catch (const RecoverableUserCallbackError &)
+            {
+              // ok, ignore, but reset the pointer
+              eptr = nullptr;
+            }
+          catch (...)
+            {
+              // uh oh:
+              AssertThrow(false, ExcInternalError());
+            }
+        }
+
+      // Call the function and if that succeeds, return zero:
+      try
+        {
+          f(std::forward<Args>(args)...);
+          eptr = nullptr;
+          return 0;
+        }
+      // If the call failed with a recoverable error, then
+      // ignore the exception for now (but store a pointer to it)
+      // and return a positive return value (+1). If the underlying
+      // implementation manages to recover
+      catch (const RecoverableUserCallbackError &)
+        {
+          eptr = std::current_exception();
+          return 1;
+        }
+      // For any other exception, capture the exception and
+      // return -1:
+      catch (const std::exception &)
+        {
+          eptr = std::current_exception();
+          return -1;
+        }
+    }
+  } // namespace
+
   template <typename VectorType>
   KINSOL<VectorType>::AdditionalData::AdditionalData(
     const SolutionStrategy &strategy,
@@ -148,7 +219,10 @@ namespace SUNDIALS
 
       int err = 0;
       if (solver.residual)
-        err = solver.residual(*src_yy, *dst_FF);
+        err = call_and_possibly_capture_exception(solver.residual,
+                                                  solver.pending_exception,
+                                                  *src_yy,
+                                                  *dst_FF);
       else
         Assert(false, ExcInternalError());
 
@@ -177,7 +251,10 @@ namespace SUNDIALS
 
       int err = 0;
       if (solver.iteration_function)
-        err = solver.iteration_function(*src_yy, *dst_FF);
+        err = call_and_possibly_capture_exception(solver.iteration_function,
+                                                  solver.pending_exception,
+                                                  *src_yy,
+                                                  *dst_FF);
       else
         Assert(false, ExcInternalError());
 
@@ -215,9 +292,10 @@ namespace SUNDIALS
       internal::copy(*fcur, f);
 
       // Call the user-provided setup function with these arguments:
-      solver.setup_jacobian(*ycur, *fcur);
-
-      return 0;
+      return call_and_possibly_capture_exception(solver.setup_jacobian,
+                                                 solver.pending_exception,
+                                                 *ycur,
+                                                 *fcur);
     }
 
 
@@ -251,7 +329,12 @@ namespace SUNDIALS
 
           internal::copy(*src_b, b);
 
-          const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
+          const int err =
+            call_and_possibly_capture_exception(solver.solve_with_jacobian,
+                                                solver.pending_exception,
+                                                *src_b,
+                                                *dst_x,
+                                                tol);
 
           internal::copy(x, *dst_x);
 
@@ -283,7 +366,12 @@ namespace SUNDIALS
           // These vectors will have zero lengths because we don't reinit them
           // above.
           const int err =
-            solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
+            call_and_possibly_capture_exception(solver.solve_jacobian_system,
+                                                solver.pending_exception,
+                                                *src_ycur,
+                                                *src_fcur,
+                                                *src_b,
+                                                *dst_x);
 
           internal::copy(x, *dst_x);
 
@@ -304,7 +392,8 @@ namespace SUNDIALS
   template <typename VectorType>
   KINSOL<VectorType>::KINSOL(const AdditionalData &data,
                              const MPI_Comm &      mpi_comm)
-    : data(data)
+    : pending_exception(nullptr)
+    , data(data)
     , mpi_communicator(mpi_comm)
     , kinsol_mem(nullptr)
 #  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
@@ -342,6 +431,8 @@ namespace SUNDIALS
     (void)status;
     AssertKINSOL(status);
 #  endif
+
+    Assert(pending_exception == nullptr, ExcInternalError());
   }
 
 
@@ -584,8 +675,36 @@ namespace SUNDIALS
         AssertKINSOL(status);
       }
 
-    // call to KINSol
+    // Having set up all of the ancillary things, finally call the main KINSol
+    // function. One we return, check that what happened:
+    // - If we have a pending recoverable exception, ignore it if SUNDIAL's
+    //   return code was zero -- in that case, SUNDIALS managed to indeed
+    //   recover and we no longer need the exception
+    // - If we have any other exception, rethrow it
+    // - If no exception, test that SUNDIALS really did successfully return
+    Assert(pending_exception == nullptr, ExcInternalError());
     status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
+
+    if (pending_exception)
+      {
+        try
+          {
+            std::rethrow_exception(pending_exception);
+          }
+        catch (const RecoverableUserCallbackError &exc)
+          {
+            pending_exception = nullptr;
+            if (status == 0)
+              /* just eat the exception */;
+            else
+              throw;
+          }
+        catch (...)
+          {
+            pending_exception = nullptr;
+            throw;
+          }
+      }
     AssertKINSOL(status);
 
     internal::copy(initial_guess_and_solution, solution);

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.