]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert ARKode interfaces to conform to our callback policy.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 15 May 2023 22:07:51 +0000 (16:07 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 18 May 2023 16:35:55 +0000 (10:35 -0600)
include/deal.II/sundials/arkode.h
include/deal.II/sundials/sunlinsol_wrapper.h
source/sundials/arkode.cc

index f69622f06ed4f46169000fbd2f217febe2e71df8..5ad3d25363e588ca29eb94d4fffcf547ca77f653 100644 (file)
@@ -51,6 +51,7 @@
 #  include <sundials/sundials_math.h>
 #  include <sundials/sundials_types.h>
 
+#  include <exception>
 #  include <memory>
 
 
@@ -579,15 +580,15 @@ namespace SUNDIALS
      * provided. According to which one is provided, explicit, implicit, or
      * mixed RK methods are used.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
     std::function<
-      int(const double t, const VectorType &y, VectorType &explicit_f)>
+      void(const double t, const VectorType &y, VectorType &explicit_f)>
       explicit_function;
 
     /**
@@ -599,14 +600,14 @@ namespace SUNDIALS
      * provided. According to which one is provided, explicit, implicit, or
      * mixed RK methods are used.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(const double t, const VectorType &y, VectorType &res)>
+    std::function<void(const double t, const VectorType &y, VectorType &res)>
       implicit_function;
 
     /**
@@ -619,14 +620,14 @@ namespace SUNDIALS
      * A call to this function should store in `Mv` the result of $M$
      * applied to `v`.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(double t, const VectorType &v, VectorType &Mv)>
+    std::function<void(const double t, const VectorType &v, VectorType &Mv)>
       mass_times_vector;
 
     /**
@@ -656,14 +657,14 @@ namespace SUNDIALS
      *
      * @param t The current evaluation time
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(const double t)> mass_times_setup;
+    std::function<void(const double t)> mass_times_setup;
 
     /**
      * A function object that users may supply and that is intended to compute
@@ -684,18 +685,18 @@ namespace SUNDIALS
      * @param[in] fy  The current value of the implicit right-hand side at y,
      *   $f_I (t_n, y)$.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(const VectorType &v,
-                      VectorType &      Jv,
-                      double            t,
-                      const VectorType &y,
-                      const VectorType &fy)>
+    std::function<void(const VectorType &v,
+                       VectorType &      Jv,
+                       const double      t,
+                       const VectorType &y,
+                       const VectorType &fy)>
       jacobian_times_vector;
 
     /**
@@ -725,14 +726,15 @@ namespace SUNDIALS
      * @param fy  The implicit right-hand side function evaluated at the
      *   current time $t$ and state $y$, i.e., $f_I(y,t)$
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(realtype t, const VectorType &y, const VectorType &fy)>
+    std::function<
+      void(const double t, const VectorType &y, const VectorType &fy)>
       jacobian_times_setup;
 
     /**
@@ -753,6 +755,13 @@ namespace SUNDIALS
      * SUNDIALS packaged SPGMR solver with default settings is used.
      *
      * For more details on the function type refer to LinearSolveFunction.
+     *
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
     LinearSolveFunction<VectorType> solve_linearized_system;
 
@@ -770,6 +779,13 @@ namespace SUNDIALS
      * and applied in mass_times_vector().
      *
      * For more details on the function type refer to LinearSolveFunction.
+     *
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
     LinearSolveFunction<VectorType> solve_mass;
 
@@ -797,21 +813,21 @@ namespace SUNDIALS
      *   (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If
      *   used with a custom solve_mass() function this will be set to zero.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(double            t,
-                      const VectorType &y,
-                      const VectorType &fy,
-                      const VectorType &r,
-                      VectorType &      z,
-                      double            gamma,
-                      double            tol,
-                      int               lr)>
+    std::function<void(const double      t,
+                       const VectorType &y,
+                       const VectorType &fy,
+                       const VectorType &r,
+                       VectorType &      z,
+                       const double      gamma,
+                       const double      tol,
+                       const int         lr)>
       jacobian_preconditioner_solve;
 
     /**
@@ -849,19 +865,19 @@ namespace SUNDIALS
      * @param[in] gamma The value $\gamma$ in $M-\gamma J$. The preconditioner
      *   should approximate the inverse of this matrix.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(double            t,
-                      const VectorType &y,
-                      const VectorType &fy,
-                      int               jok,
-                      int &             jcur,
-                      double            gamma)>
+    std::function<void(const double      t,
+                       const VectorType &y,
+                       const VectorType &fy,
+                       const int         jok,
+                       int &             jcur,
+                       const double      gamma)>
       jacobian_preconditioner_setup;
 
     /**
@@ -884,15 +900,18 @@ namespace SUNDIALS
      *   (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If
      *   used with a custom solve_mass() function this will be set to zero.
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<
-      int(double t, const VectorType &r, VectorType &z, double tol, int lr)>
+    std::function<void(const double      t,
+                       const VectorType &r,
+                       VectorType &      z,
+                       const double      tol,
+                       const int         lr)>
       mass_preconditioner_solve;
 
     /**
@@ -912,14 +931,14 @@ namespace SUNDIALS
      *
      * @param[in] t  The current time
      *
-     * This function should return:
-     * - 0: Success
-     * - >0: Recoverable error, ARKode will reattempt the solution and call this
-     *       function again.
-     * - <0: Unrecoverable error, the computation will be aborted and an
-     *       assertion will be thrown.
+     * @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. In particular, ARKode can deal
+     * with "recoverable" errors in some circumstances, so callbacks
+     * can throw exceptions of type RecoverableUserCallbackError.
      */
-    std::function<int(double t)> mass_preconditioner_setup;
+    std::function<void(const double t)> mass_preconditioner_setup;
 
     /**
      * A function object that users may supply and that is intended to
@@ -983,7 +1002,7 @@ namespace SUNDIALS
      * @endcode
      *
      * @note This function will be called at the end of all other set up right
-     *   before the actual time evloution is started or continued with
+     *   before the actual time evolution is started or continued with
      *   solve_ode(). This function is also called when the solver is restarted,
      *   see solver_should_restart(). Consult the SUNDIALS manual to see which
      *   options are still available at this point.
@@ -1068,6 +1087,14 @@ namespace SUNDIALS
     std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
     std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
 
+  public:
+    /**
+     * 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;
+
 #  ifdef DEAL_II_WITH_PETSC
 #    ifdef PETSC_USE_COMPLEX
     static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
index 509596f629f84b6eadfa45e3473b6f5fc59dd888..fa586df60b0bcb7516e2d1064988183fec06a5c7 100644 (file)
@@ -171,7 +171,7 @@ namespace SUNDIALS
   };
 
   /**
-   * Type of function objects to interface with SUNDIALS linear solvers
+   * Type of function objects to interface with SUNDIALS' linear solvers
    *
    * This function type encapsulates the action of solving $P^{-1}Ax=P^{-1}b$.
    * The LinearOperator @p op encapsulates the matrix vector product $Ax$ and
@@ -188,20 +188,21 @@ namespace SUNDIALS
    * @param[in] b The right-hand side
    * @param[in] tol Tolerance for the iterative solver
    *
-   * This function should return:
-   * - 0: Success
-   * - >0: Recoverable error, ARKode will reattempt the solution and call this
-   *       function again.
-   * - <0: Unrecoverable error, the computation will be aborted and an
-   *       assertion will be thrown.
+   *
+   * @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. In particular, ARKode can deal
+   * with "recoverable" errors in some circumstances, so callbacks
+   * can throw exceptions of type RecoverableUserCallbackError.
    */
   template <typename VectorType>
   using LinearSolveFunction =
-    std::function<int(SundialsOperator<VectorType> &      op,
-                      SundialsPreconditioner<VectorType> &prec,
-                      VectorType &                        x,
-                      const VectorType &                  b,
-                      double                              tol)>;
+    std::function<void(SundialsOperator<VectorType> &      op,
+                       SundialsPreconditioner<VectorType> &prec,
+                       VectorType &                        x,
+                       const VectorType &                  b,
+                       double                              tol)>;
 
   namespace internal
   {
index 8554e7d3ffe25b596c9249f11b16d23a4665f7a2..3186c288923c42a8a2d57ad7f9b7f430ef5644a1 100644 (file)
@@ -49,6 +49,78 @@ 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
+
+
   namespace
   {
     template <typename VectorType>
@@ -65,7 +137,11 @@ namespace SUNDIALS
       auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
       auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
 
-      return solver.explicit_function(tt, *src_yy, *dst_yp);
+      return call_and_possibly_capture_exception(solver.explicit_function,
+                                                 solver.pending_exception,
+                                                 tt,
+                                                 *src_yy,
+                                                 *dst_yp);
     }
 
 
@@ -84,7 +160,11 @@ namespace SUNDIALS
       auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
       auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
 
-      return solver.implicit_function(tt, *src_yy, *dst_yp);
+      return call_and_possibly_capture_exception(solver.implicit_function,
+                                                 solver.pending_exception,
+                                                 tt,
+                                                 *src_yy,
+                                                 *dst_yp);
     }
 
 
@@ -109,7 +189,13 @@ namespace SUNDIALS
 
       auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
 
-      return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
+      return call_and_possibly_capture_exception(solver.jacobian_times_vector,
+                                                 solver.pending_exception,
+                                                 *src_v,
+                                                 *dst_Jv,
+                                                 t,
+                                                 *src_y,
+                                                 *src_fy);
     }
 
 
@@ -128,7 +214,11 @@ namespace SUNDIALS
       auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
       auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
 
-      return solver.jacobian_times_setup(t, *src_y, *src_fy);
+      return call_and_possibly_capture_exception(solver.jacobian_times_setup,
+                                                 solver.pending_exception,
+                                                 t,
+                                                 *src_y,
+                                                 *src_fy);
     }
 
 
@@ -155,8 +245,17 @@ namespace SUNDIALS
 
       auto *dst_z = internal::unwrap_nvector<VectorType>(z);
 
-      return solver.jacobian_preconditioner_solve(
-        t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr);
+      return call_and_possibly_capture_exception(
+        solver.jacobian_preconditioner_solve,
+        solver.pending_exception,
+        t,
+        *src_y,
+        *src_fy,
+        *src_r,
+        *dst_z,
+        gamma,
+        delta,
+        lr);
     }
 
 
@@ -178,8 +277,15 @@ namespace SUNDIALS
       auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
       auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
 
-      return solver.jacobian_preconditioner_setup(
-        t, *src_y, *src_fy, jok, *jcurPtr, gamma);
+      return call_and_possibly_capture_exception(
+        solver.jacobian_preconditioner_setup,
+        solver.pending_exception,
+        t,
+        *src_y,
+        *src_fy,
+        jok,
+        *jcurPtr,
+        gamma);
     }
 
 
@@ -198,7 +304,8 @@ namespace SUNDIALS
       auto *src_v  = internal::unwrap_nvector_const<VectorType>(v);
       auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
 
-      return solver.mass_times_vector(t, *src_v, *dst_Mv);
+      return call_and_possibly_capture_exception(
+        solver.mass_times_vector, solver.pending_exception, t, *src_v, *dst_Mv);
     }
 
 
@@ -211,7 +318,9 @@ namespace SUNDIALS
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(mtimes_data);
 
-      return solver.mass_times_setup(t);
+      return call_and_possibly_capture_exception(solver.mass_times_setup,
+                                                 solver.pending_exception,
+                                                 t);
     }
 
 
@@ -232,7 +341,14 @@ namespace SUNDIALS
       auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
       auto *dst_z = internal::unwrap_nvector<VectorType>(z);
 
-      return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
+      return call_and_possibly_capture_exception(
+        solver.mass_preconditioner_solve,
+        solver.pending_exception,
+        t,
+        *src_r,
+        *dst_z,
+        delta,
+        lr);
     }
 
 
@@ -245,7 +361,8 @@ namespace SUNDIALS
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
 
-      return solver.mass_preconditioner_setup(t);
+      return call_and_possibly_capture_exception(
+        solver.mass_preconditioner_setup, solver.pending_exception, t);
     }
   } // namespace
 
@@ -368,14 +485,41 @@ namespace SUNDIALS
       {
         time.set_desired_next_step_size(data.output_period);
 
-        // Let ARKode advance time by one period:
+        // Having set up all of the ancillary things, finally call the main
+        // ARKode 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());
         double     actual_next_time;
         const auto status = ARKStepEvolve(arkode_mem,
                                           time.get_next_time(),
                                           solution_nvector,
                                           &actual_next_time,
                                           ARK_NORMAL);
-        (void)status;
+        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;
+              }
+          }
+
         AssertARKode(status);
 
         // Then reflect this time advancement in our own DiscreteTime object:

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.