# include <sundials/sundials_math.h>
# include <sundials/sundials_types.h>
+# include <exception>
# include <memory>
* - <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
* - <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;
/**
* - <0: Unrecoverable error the computation will be aborted and an
* assertion will be thrown.
*/
- std::function<int(const VectorType ¤t_u, const VectorType ¤t_f)>
+ std::function<void(const VectorType ¤t_u,
+ const VectorType ¤t_f)>
setup_jacobian;
/**
* 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;
/**
* 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;
/**
<< ". 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
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,
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());
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());
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);
}
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);
// 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);
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)
(void)status;
AssertKINSOL(status);
# endif
+
+ Assert(pending_exception == nullptr, ExcInternalError());
}
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);