From: Wolfgang Bangerth Date: Thu, 18 May 2023 05:50:55 +0000 (-0600) Subject: Convert IDA to use our new callback standard. X-Git-Tag: v9.5.0-rc1~188^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=31e9d4fbb2ecc1d759a9a4f0ab36cc3c5ec8361c;p=dealii.git Convert IDA to use our new callback standard. --- diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 99757a21bc..ae65a5f43e 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -198,17 +198,16 @@ namespace SUNDIALS * time_stepper.residual = [&](const double t, * const VectorType &y, * const VectorType &y_dot, - * VectorType &res) ->int + * VectorType &res) * { * res = y_dot; * A.vmult_add(res, y); - * return 0; * }; * * time_stepper.setup_jacobian = [&](const double , * const VectorType &, * const VectorType &, - * const double alpha) ->int + * const double alpha) * { * J = A; * @@ -216,14 +215,12 @@ namespace SUNDIALS * J(1,1) = alpha; * * Jinv.invert(J); - * return 0; * }; * * time_stepper.solve_jacobian_system = [&](const VectorType &src, - * VectorType &dst) ->int + * VectorType &dst) * { * Jinv.vmult(dst,src); - * return 0; * }; * * y[1] = kappa; @@ -651,17 +648,17 @@ namespace SUNDIALS /** * Compute residual. Return $F(t, y, \dot y)$. * - * This function should return: - * - 0: Success - * - >0: Recoverable error (IDAReinit will be called if this happens, and - * then last function will be attempted 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, IDA can deal + * with "recoverable" errors in some circumstances, so callbacks + * can throw exceptions of type RecoverableUserCallbackError. */ - std::function + std::function residual; /** @@ -677,7 +674,7 @@ namespace SUNDIALS * \alpha \dfrac{\partial F}{\partial \dot y}. * \f] * - * If the user uses a matrix based computation of the Jacobian, than this + * If the user uses a matrix based computation of the Jacobian, then this * is the right place where an assembly routine should be called to * assemble both a matrix and a preconditioner for the Jacobian system. * Subsequent calls (possibly more than one) to solve_jacobian_system() or @@ -690,17 +687,17 @@ namespace SUNDIALS * solve_with_jacobian() to obtain a solution $x$ to the * system $J x = b$. * - * This function should return: - * - 0: Success - * - >0: Recoverable error (IDAReinit will be called if this happens, and - * then last function will be attempted 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, IDA can deal + * with "recoverable" errors in some circumstances, so callbacks + * can throw exceptions of type RecoverableUserCallbackError. */ - std::function + std::function setup_jacobian; /** @@ -724,18 +721,18 @@ namespace SUNDIALS * applied to `src`, i.e., `J*dst = src`. It is the users responsibility * to set up proper solvers and preconditioners inside this function. * - * This function should return: - * - 0: Success - * - >0: Recoverable error (IDAReinit will be called if this happens, and - * then last function will be attempted 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, IDA can deal + * with "recoverable" errors in some circumstances, so callbacks + * can throw exceptions of type RecoverableUserCallbackError. * * @deprecated Use solve_with_jacobian() instead which also uses a numerical * tolerance. */ DEAL_II_DEPRECATED - std::function + std::function solve_jacobian_system; /** @@ -772,15 +769,15 @@ namespace SUNDIALS * `setup_jacobian()`, given that that function is called far less often * than the current one.) * - * This function should return: - * - 0: Success - * - >0: Recoverable error (IDAReinit will be called if this happens, and - * then the last function will be attempted 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, IDA can deal + * with "recoverable" errors in some circumstances, so callbacks + * can throw exceptions of type RecoverableUserCallbackError. */ std::function< - int(const VectorType &rhs, VectorType &dst, const double tolerance)> + void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian; /** @@ -818,6 +815,13 @@ namespace SUNDIALS * * The default implementation simply returns `false`, i.e., no restart is * performed during the evolution. + * + * @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, IDA can deal + * with "recoverable" errors in some circumstances, so callbacks + * can throw exceptions of type RecoverableUserCallbackError. */ std::function solver_should_restart; @@ -902,6 +906,14 @@ namespace SUNDIALS */ GrowingVectorMemory mem; + 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::value, diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 39d36ca214..c149169af0 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -44,6 +44,78 @@ 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 + + namespace { template @@ -60,9 +132,12 @@ namespace SUNDIALS auto *src_yp = internal::unwrap_nvector_const(yp); auto *residual = internal::unwrap_nvector(rr); - int err = solver.residual(tt, *src_yy, *src_yp, *residual); - - return err; + return call_and_possibly_capture_exception(solver.residual, + solver.pending_exception, + tt, + *src_yy, + *src_yp, + *residual); } @@ -86,9 +161,12 @@ namespace SUNDIALS auto *src_yy = internal::unwrap_nvector_const(yy); auto *src_yp = internal::unwrap_nvector_const(yp); - int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj); - - return err; + return call_and_possibly_capture_exception(solver.setup_jacobian, + solver.pending_exception, + tt, + *src_yy, + *src_yp, + cj); } @@ -105,16 +183,23 @@ namespace SUNDIALS auto *src_b = internal::unwrap_nvector_const(b); auto *dst_x = internal::unwrap_nvector(x); - int err = 0; if (solver.solve_with_jacobian) - err = solver.solve_with_jacobian(*src_b, *dst_x, tol); + return call_and_possibly_capture_exception(solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); else if (solver.solve_jacobian_system) - err = solver.solve_jacobian_system(*src_b, *dst_x); + return call_and_possibly_capture_exception(solver.solve_jacobian_system, + solver.pending_exception, + *src_b, + *dst_x); else - // We have already checked this outside, so we should never get here. - Assert(false, ExcInternalError()); - - return err; + { + // We have already checked this outside, so we should never get here. + Assert(false, ExcInternalError()); + return -1; + } } } // namespace @@ -135,6 +220,7 @@ namespace SUNDIALS , ida_ctx(nullptr) # endif , mpi_communicator(mpi_comm) + , pending_exception(nullptr) { // SUNDIALS will always duplicate communicators if we provide them. This // can cause problems if SUNDIALS is configured with MPI and we pass along @@ -163,6 +249,8 @@ namespace SUNDIALS (void)status; AssertIDA(status); # endif + + Assert(pending_exception == nullptr, ExcInternalError()); } @@ -205,7 +293,31 @@ namespace SUNDIALS # endif ); + // Execute time steps. If we ended up with a pending exception, + // see if it was recoverable; if it was, and if IDA recovered, + // just continue on. If IDA did not recover, rethrow the exception. + // Do the same if the exception was not recoverable. status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); + if (pending_exception) + { + try + { + std::rethrow_exception(pending_exception); + } + catch (const RecoverableUserCallbackError &exc) + { + pending_exception = nullptr; + if (status == 0) + /* just eat the exception and continue */; + else + throw; + } + catch (...) + { + pending_exception = nullptr; + throw; + } + } AssertIDA(status); status = IDAGetLastStep(ida_mem, &h);