# include <petscsnes.h>
+# include <exception>
# if defined(DEAL_II_HAVE_CXX20)
# include <concepts>
# endif
/**
* Callback for the computation of the nonlinear residual $F(x)$.
+ *
+ * @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 SNES 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<int(const VectorType &x, VectorType &res)> residual;
+ std::function<void(const VectorType &x, VectorType &res)> residual;
/**
* Callback for the computation of the Jacobian
* $\dfrac{\partial F}{\partial x}$.
- */
- std::function<int(const VectorType &x, 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 SNES 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 &x, AMatrixType &A, PMatrixType &P)>
jacobian;
/**
* This function is called by NonlinearSolver at the beginning
* of each time step. Input arguments are the current step number
* and the current value for ||F(x)||.
- */
- std::function<int(const VectorType & x,
- const unsigned int step_number,
- const real_type f_norm)>
+ *
+ * @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 SNES 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 & x,
+ const unsigned int step_number,
+ const real_type f_norm)>
monitor;
/**
* operator $\dfrac{\partial F}{\partial x}$.
*
* Solvers must be provided via NonlinearSolver::solve_with_jacobian.
+ *
+ * @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 SNES 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<int(const VectorType &x)> setup_jacobian;
+ std::function<void(const VectorType &x)> setup_jacobian;
/**
* Callback for the solution of the tangent system set up with
* NonlinearSolver::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 SNES 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;
/**
* the reduction in a trust region step.
*
* The value of the energy function must be returned in @p energy_value.
+ *
+ * @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 SNES 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<int(const VectorType &x, real_type &energy_value)> energy;
+ std::function<void(const VectorType &x, real_type &energy_value)> energy;
private:
/**
* 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
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>))
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::NonlinearSolver(
const NonlinearSolverData &data,
const MPI_Comm mpi_comm)
+ : pending_exception(nullptr)
{
AssertPETSc(SNESCreate(mpi_comm, &snes));
AssertPETSc(SNESSetApplicationContext(snes, this));
NonlinearSolver<VectorType, PMatrixType, AMatrixType>::~NonlinearSolver()
{
AssertPETSc(SNESDestroy(&snes));
+
+ Assert(pending_exception == nullptr, ExcInternalError());
}
VectorType xdealii(x);
VectorType fdealii(f);
- AssertUser(user->residual(xdealii, fdealii), "residual");
+ const int err = call_and_possibly_capture_exception(
+ user->residual, user->pending_exception, xdealii, fdealii);
petsc_increment_state_counter(f);
- PetscFunctionReturn(0);
+ PetscFunctionReturn(err);
};
const auto snes_jacobian =
VectorType xdealii(x);
AMatrixType Adealii(A);
PMatrixType Pdealii(P);
- AssertUser(user->jacobian(xdealii, Adealii, Pdealii), "jacobian");
+
+ const int err = call_and_possibly_capture_exception(
+ user->jacobian, user->pending_exception, xdealii, Adealii, Pdealii);
petsc_increment_state_counter(P);
// Handle the Jacobian-free case
}
else
petsc_increment_state_counter(A);
- PetscFunctionReturn(0);
+ PetscFunctionReturn(err);
};
const auto snes_jacobian_with_setup =
user->A = &Adealii;
user->P = &Pdealii;
- AssertUser(user->setup_jacobian(xdealii), "setup_jacobian");
+ const int err =
+ call_and_possibly_capture_exception(user->setup_jacobian,
+ user->pending_exception,
+ xdealii);
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 snes_monitor =
Vec x;
AssertPETSc(SNESGetSolution(snes, &x));
VectorType xdealii(x);
- AssertUser(user->monitor(xdealii, it, f), "monitor");
- PetscFunctionReturn(0);
+ const int err = call_and_possibly_capture_exception(
+ user->monitor, user->pending_exception, xdealii, it, f);
+ PetscFunctionReturn(err);
};
const auto snes_objective =
real_type v;
VectorType xdealii(x);
- AssertUser(user->energy(xdealii, v), "energy");
+ const int err = call_and_possibly_capture_exception(
+ user->energy, user->pending_exception, xdealii, v);
*f = v;
- PetscFunctionReturn(0);
+ PetscFunctionReturn(err);
};
AssertThrow(residual,
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 solve the system.
// Older versions of PETSc requires the solution vector specified even
- // if we specified SNESSetSolution upfront
- AssertPETSc(SNESSolve(snes, nullptr, x.petsc_vector()));
+ // if we specified SNESSetSolution upfront.
+ //
+ // 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 = SNESSolve(snes, nullptr, x.petsc_vector());
+ 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.
PetscInt nt;
projects / dealii.git / commitdiff