From: Wolfgang Bangerth Date: Mon, 15 May 2023 22:07:51 +0000 (-0600) Subject: Convert ARKode interfaces to conform to our callback policy. X-Git-Tag: v9.5.0-rc1~214^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f1501d5790a3d5c19824d14209670ca5a7cd92cf;p=dealii.git Convert ARKode interfaces to conform to our callback policy. --- diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index f69622f06e..5ad3d25363 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -51,6 +51,7 @@ # include # include +# include # include @@ -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 + std::function 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 + std::function 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 mass_times_setup; + std::function 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 + std::function 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 + 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 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 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 + std::function 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 + std::function 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 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 mass_preconditioner_setup; + std::function 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> linear_solver; std::unique_ptr> 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::value, diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h index 509596f629..fa586df60b 100644 --- a/include/deal.II/sundials/sunlinsol_wrapper.h +++ b/include/deal.II/sundials/sunlinsol_wrapper.h @@ -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 using LinearSolveFunction = - std::function & op, - SundialsPreconditioner &prec, - VectorType & x, - const VectorType & b, - double tol)>; + std::function & op, + SundialsPreconditioner &prec, + VectorType & x, + const VectorType & b, + double tol)>; namespace internal { diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index 8554e7d3ff..3186c28892 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -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 + 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)...); + 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 @@ -65,7 +137,11 @@ namespace SUNDIALS auto *src_yy = internal::unwrap_nvector_const(yy); auto *dst_yp = internal::unwrap_nvector(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(yy); auto *dst_yp = internal::unwrap_nvector(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(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(y); auto *src_fy = internal::unwrap_nvector_const(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(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(y); auto *src_fy = internal::unwrap_nvector_const(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(v); auto *dst_Mv = internal::unwrap_nvector(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 &solver = *static_cast *>(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(r); auto *dst_z = internal::unwrap_nvector(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 &solver = *static_cast *>(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: