From fd22ffa278dad25662bfb59450060217e3d5626a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 2 May 2023 17:49:57 -0600 Subject: [PATCH] Convert KINSOL callbacks to use exceptions. --- include/deal.II/sundials/kinsol.h | 25 ++++-- source/sundials/kinsol.cc | 137 ++++++++++++++++++++++++++++-- 2 files changed, 145 insertions(+), 17 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index ee17684d4e..df4177b520 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -38,6 +38,7 @@ # include # include +# include # include @@ -430,7 +431,7 @@ namespace SUNDIALS * - <0: Unrecoverable error the computation will be aborted and an * assertion will be thrown. */ - std::function residual; + std::function residual; /** * A function object that users should supply and that is intended to @@ -445,7 +446,7 @@ namespace SUNDIALS * - <0: Unrecoverable error; the computation will be aborted and an * assertion will be thrown. */ - std::function + std::function iteration_function; /** @@ -493,7 +494,8 @@ namespace SUNDIALS * - <0: Unrecoverable error the computation will be aborted and an * assertion will be thrown. */ - std::function + std::function setup_jacobian; /** @@ -555,10 +557,10 @@ namespace SUNDIALS * that the Jacobian should be re-computed in every iteration. */ DEAL_II_DEPRECATED - std::function + std::function solve_jacobian_system; /** @@ -601,7 +603,7 @@ namespace SUNDIALS * 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; /** @@ -670,6 +672,13 @@ namespace SUNDIALS << ". 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 diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index fc1d9adbd5..13a25fc3bb 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -52,6 +52,77 @@ 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 + template KINSOL::AdditionalData::AdditionalData( const SolutionStrategy &strategy, @@ -148,7 +219,10 @@ namespace SUNDIALS 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()); @@ -177,7 +251,10 @@ namespace SUNDIALS 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()); @@ -215,9 +292,10 @@ namespace SUNDIALS 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); } @@ -251,7 +329,12 @@ namespace SUNDIALS 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); @@ -283,7 +366,12 @@ namespace SUNDIALS // 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); @@ -304,7 +392,8 @@ namespace SUNDIALS template KINSOL::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) @@ -342,6 +431,8 @@ namespace SUNDIALS (void)status; AssertKINSOL(status); # endif + + Assert(pending_exception == nullptr, ExcInternalError()); } @@ -584,8 +675,36 @@ namespace SUNDIALS 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); -- 2.39.5