# include <sundials/sundials_math.h>
# include <sundials/sundials_types.h>
+# include <exception>
# include <memory>
* 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;
/**
* 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;
/**
* 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;
/**
*
* @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
* @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;
/**
* @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;
/**
* 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;
* 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;
* (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;
/**
* @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;
/**
* (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;
/**
*
* @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
* @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.
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,
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>
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);
}
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);
}
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);
}
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);
}
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);
}
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);
}
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);
}
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);
}
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);
}
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
{
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:
projects / dealii.git / commitdiff