# include <petscts.h>
+# include <exception>
+
# if defined(DEAL_II_HAVE_CXX20)
# include <concepts>
# endif
/**
* Callback for the computation of the implicit residual $F(t, y, \dot y)$.
- */
- std::function<int(const real_type t,
- const VectorType &y,
- const VectorType &y_dot,
- VectorType & res)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ VectorType & res)>
implicit_function;
/**
* All implicit solvers implementations are recast to use the above
* linearization. The $\alpha$ parameter is time-step and solver-type
* specific.
- */
- std::function<int(const real_type t,
- const VectorType &y,
- const VectorType &y_dot,
- const real_type alpha,
- AMatrixType & A,
- PMatrixType & P)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t,
+ const VectorType &y,
+ const VectorType &y_dot,
+ const real_type alpha,
+ AMatrixType & A,
+ PMatrixType & P)>
implicit_jacobian;
/**
* Callback for the computation of the explicit residual $G(t, y)$.
- */
- std::function<int(const real_type t, const VectorType &y, VectorType &res)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t, const VectorType &y, VectorType &res)>
explicit_function;
/**
* Callback for the computation of the explicit Jacobian $\dfrac{\partial
* G}{\partial y}$.
- */
- std::function<int(const real_type t,
- const VectorType &y,
- AMatrixType & A,
- PMatrixType & P)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t,
+ const VectorType &y,
+ AMatrixType & A,
+ PMatrixType & P)>
explicit_jacobian;
/**
*
* This function is called by TimeStepper at the beginning
* of each time step.
- */
- std::function<int(const real_type t,
- const VectorType & y,
- const unsigned int step_number)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t,
+ const VectorType & y,
+ const unsigned int step_number)>
monitor;
/**
* specific.
*
* Solvers must be provided via TimeStepper::solve_with_jacobian.
- */
- std::function<int(const real_type t,
- const VectorType &y,
- const VectorType &ydot,
- const real_type alpha)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const real_type t,
+ const VectorType &y,
+ const VectorType &ydot,
+ const real_type alpha)>
setup_jacobian;
/**
* TimeStepper::setup_jacobian.
*
* This is used as a preconditioner inside the Krylov process.
- */
- std::function<int(const VectorType &src, VectorType &dst)>
+ *
+ * @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. PETSc's TS package can not deal
+ * with "recoverable" errors, so if a callback
+ * throws an exception of type RecoverableUserCallbackError, then this
+ * exception is treated like any other exception.
+ */
+ std::function<void(const VectorType &src, VectorType &dst)>
solve_with_jacobian;
/**
* This flag is used to support versions of PETSc older than 3.13.
*/
bool need_dummy_assemble;
+
+ /**
+ * 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;
};
} // namespace PETScWrappers
} \
while (0)
-// Shorthand notation for User error codes.
-# define AssertUser(code, name) \
- do \
- { \
- int ierr = (code); \
- AssertThrow(ierr == 0, \
- StandardExceptions::ExcFunctionNonzeroReturn(name, ierr)); \
- } \
- while (0)
-
namespace PETScWrappers
{
+ 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, then the current function returns with
+ * an error code of -1. In that case, 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. There is no reason why this should be so, and
+ // we should probably bail out:
+ AssertThrow(eptr == nullptr, ExcInternalError());
+
+ // Call the function and if that succeeds, return zero:
+ try
+ {
+ f(std::forward<Args>(args)...);
+ eptr = nullptr;
+ return 0;
+ }
+ // In case of an exception, capture the exception and
+ // return -1:
+ catch (...)
+ {
+ eptr = std::current_exception();
+ return -1;
+ }
+ }
+ } // namespace
+
+
+
template <typename VectorType, typename PMatrixType, typename AMatrixType>
DEAL_II_CXX20_REQUIRES((concepts::is_dealii_petsc_vector_type<VectorType> ||
std::constructible_from<VectorType, Vec>))
TimeStepper<VectorType, PMatrixType, AMatrixType>::TimeStepper(
const TimeStepperData &data,
const MPI_Comm mpi_comm)
+ : pending_exception(nullptr)
{
AssertPETSc(TSCreate(mpi_comm, &ts));
AssertPETSc(TSSetApplicationContext(ts, this));
TimeStepper<VectorType, PMatrixType, AMatrixType>::~TimeStepper()
{
AssertPETSc(TSDestroy(&ts));
+
+ Assert(pending_exception == nullptr, ExcInternalError());
}
VectorType xdealii(x);
VectorType xdotdealii(xdot);
VectorType fdealii(f);
- AssertUser(user->implicit_function(t, xdealii, xdotdealii, fdealii),
- "implicit_function");
+ const int err =
+ call_and_possibly_capture_exception(user->implicit_function,
+ user->pending_exception,
+ t,
+ xdealii,
+ xdotdealii,
+ fdealii);
petsc_increment_state_counter(f);
- PetscFunctionReturn(0);
+ PetscFunctionReturn(err);
};
const auto ts_ijacobian =
AMatrixType Adealii(A);
PMatrixType Pdealii(P);
- AssertUser(
- user->implicit_jacobian(t, xdealii, xdotdealii, s, Adealii, Pdealii),
- "implicit_jacobian");
+ const int err =
+ call_and_possibly_capture_exception(user->implicit_jacobian,
+ user->pending_exception,
+ t,
+ xdealii,
+ xdotdealii,
+ s,
+ Adealii,
+ Pdealii);
+
petsc_increment_state_counter(P);
// Handle the Jacobian-free case
}
else
petsc_increment_state_counter(A);
- PetscFunctionReturn(0);
+
+ PetscFunctionReturn(err);
};
const auto ts_ijacobian_with_setup =
user->A = &Adealii;
user->P = &Pdealii;
- AssertUser(user->setup_jacobian(t, xdealii, xdotdealii, s),
- "setup_jacobian");
+ const int err =
+ call_and_possibly_capture_exception(user->setup_jacobian,
+ user->pending_exception,
+ t,
+ xdealii,
+ xdotdealii,
+ s);
+
petsc_increment_state_counter(P);
// Handle older versions of PETSc for which we cannot pass a MATSHELL
}
else
petsc_increment_state_counter(A);
- PetscFunctionReturn(0);
+
+ PetscFunctionReturn(err);
};
const auto ts_rhsfunction =
VectorType xdealii(x);
VectorType fdealii(f);
- AssertUser(user->explicit_function(t, xdealii, fdealii),
- "explicit_function");
+ const int err = call_and_possibly_capture_exception(
+ user->explicit_function, user->pending_exception, t, xdealii, fdealii);
petsc_increment_state_counter(f);
- PetscFunctionReturn(0);
+ PetscFunctionReturn(err);
};
const auto ts_rhsjacobian =
AMatrixType Adealii(A);
PMatrixType Pdealii(P);
- AssertUser(user->explicit_jacobian(t, xdealii, Adealii, Pdealii),
- "explicit_jacobian");
+ const int err =
+ call_and_possibly_capture_exception(user->explicit_jacobian,
+ user->pending_exception,
+ t,
+ xdealii,
+ Adealii,
+ Pdealii);
+
petsc_increment_state_counter(P);
// Handle older versions of PETSc for which we cannot pass a MATSHELL
}
else
petsc_increment_state_counter(A);
- PetscFunctionReturn(0);
+
+ PetscFunctionReturn(err);
};
const auto ts_monitor =
auto user = static_cast<TimeStepper *>(ctx);
VectorType xdealii(x);
- AssertUser(user->monitor(t, xdealii, it), "monitor");
- PetscFunctionReturn(0);
+ const int err = call_and_possibly_capture_exception(
+ user->monitor, user->pending_exception, t, xdealii, it);
+ PetscFunctionReturn(err);
};
AssertThrow(explicit_function || implicit_function,
precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
VectorType dst(static_cast<const Vec &>(indst));
const VectorType src(static_cast<const Vec &>(insrc));
- return solve_with_jacobian(src, dst);
+ return call_and_possibly_capture_exception(solve_with_jacobian,
+ pending_exception,
+ src,
+ dst);
};
// Default Krylov solver (preconditioner only)
}
// Having set everything up, now do the actual work
- // and let PETSc do the time stepping.
- AssertPETSc(TSSolve(ts, nullptr));
+ // and let PETSc do the time stepping. If there is
+ // a pending exception, then one of the user callbacks
+ // threw an exception we didn't know how to deal with
+ // at the time. It is possible that PETSc managed to
+ // recover anyway and in that case would have returned
+ // a zero error code -- if so, just eat the exception and
+ // continue on; otherwise, just rethrow the exception
+ // and get outta here.
+ const int status = TSSolve(ts, nullptr);
+ if (pending_exception)
+ {
+ try
+ {
+ std::rethrow_exception(pending_exception);
+ }
+ catch (...)
+ {
+ pending_exception = nullptr;
+ if (status == 0)
+ /* just eat the exception */;
+ else
+ throw;
+ }
+ }
+ AssertPETSc(status);
// Get the number of steps taken.
auto nt = ts_get_step_number(ts);
projects / dealii.git / commitdiff