From d6405cb90cbf5ff76e3415bd4d13eb8d48188448 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 1 Jun 2023 15:36:37 -0600 Subject: [PATCH] Unify multiple copies of the same function. --- include/deal.II/sundials/utilities.h | 106 +++++++++++++++++++++++ source/sundials/arkode.cc | 124 ++++++--------------------- source/sundials/ida.cc | 119 ++++++------------------- source/sundials/kinsol.cc | 104 ++++------------------ 4 files changed, 175 insertions(+), 278 deletions(-) create mode 100644 include/deal.II/sundials/utilities.h diff --git a/include/deal.II/sundials/utilities.h b/include/deal.II/sundials/utilities.h new file mode 100644 index 0000000000..7a509bb838 --- /dev/null +++ b/include/deal.II/sundials/utilities.h @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_sundials_utilities_h +#define dealii_sundials_utilities_h + +#include + +#ifdef DEAL_II_WITH_SUNDIALS +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace SUNDIALS +{ + namespace Utilities + { + /** + * 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 Utilities +} // namespace SUNDIALS + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SUNDIALS + +#endif diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index dd29cf2188..fdbd600b96 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -38,6 +38,7 @@ # include # include +# include # include # include @@ -49,79 +50,6 @@ 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 ARKode::ARKode(const AdditionalData &data) : ARKode(data, MPI_COMM_SELF) @@ -350,11 +278,12 @@ namespace SUNDIALS auto *src_yy = internal::unwrap_nvector_const(yy); auto *dst_yp = internal::unwrap_nvector(yp); - return call_and_possibly_capture_exception(solver.explicit_function, - solver.pending_exception, - tt, - *src_yy, - *dst_yp); + return Utilities::call_and_possibly_capture_exception( + solver.explicit_function, + solver.pending_exception, + tt, + *src_yy, + *dst_yp); }; @@ -367,11 +296,12 @@ namespace SUNDIALS auto *src_yy = internal::unwrap_nvector_const(yy); auto *dst_yp = internal::unwrap_nvector(yp); - return call_and_possibly_capture_exception(solver.implicit_function, - solver.pending_exception, - tt, - *src_yy, - *dst_yp); + return Utilities::call_and_possibly_capture_exception( + solver.implicit_function, + solver.pending_exception, + tt, + *src_yy, + *dst_yp); }; arkode_mem = ARKStepCreate(explicit_function ? explicit_function_callback : @@ -501,7 +431,7 @@ namespace SUNDIALS auto *dst_Jv = internal::unwrap_nvector(Jv); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.jacobian_times_vector, solver.pending_exception, *src_v, @@ -520,7 +450,7 @@ namespace SUNDIALS auto *src_y = internal::unwrap_nvector_const(y); auto *src_fy = internal::unwrap_nvector_const(fy); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.jacobian_times_setup, solver.pending_exception, t, @@ -554,7 +484,7 @@ namespace SUNDIALS auto *dst_z = internal::unwrap_nvector(z); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.jacobian_preconditioner_solve, solver.pending_exception, t, @@ -581,7 +511,7 @@ namespace SUNDIALS auto *src_y = internal::unwrap_nvector_const(y); auto *src_fy = internal::unwrap_nvector_const(fy); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.jacobian_preconditioner_setup, solver.pending_exception, t, @@ -696,9 +626,8 @@ namespace SUNDIALS ARKode &solver = *static_cast *>(mtimes_data); - return call_and_possibly_capture_exception(solver.mass_times_setup, - solver.pending_exception, - t); + return Utilities::call_and_possibly_capture_exception( + solver.mass_times_setup, solver.pending_exception, t); }; auto mass_matrix_times_vector_callback = @@ -710,11 +639,12 @@ namespace SUNDIALS auto *src_v = internal::unwrap_nvector_const(v); auto *dst_Mv = internal::unwrap_nvector(Mv); - return call_and_possibly_capture_exception(solver.mass_times_vector, - solver.pending_exception, - t, - *src_v, - *dst_Mv); + return Utilities::call_and_possibly_capture_exception( + solver.mass_times_vector, + solver.pending_exception, + t, + *src_v, + *dst_Mv); }; status = ARKStepSetMassTimes(arkode_mem, @@ -733,7 +663,7 @@ namespace SUNDIALS ARKode &solver = *static_cast *>(user_data); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.mass_preconditioner_setup, solver.pending_exception, t); }; @@ -750,7 +680,7 @@ namespace SUNDIALS auto *src_r = internal::unwrap_nvector_const(r); auto *dst_z = internal::unwrap_nvector(z); - return call_and_possibly_capture_exception( + return Utilities::call_and_possibly_capture_exception( solver.mass_preconditioner_solve, solver.pending_exception, t, diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 8b8597a477..ec4a3e6f5a 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -36,6 +36,7 @@ # endif # include +# include # include # include @@ -44,78 +45,6 @@ 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 IDA::IDA(const AdditionalData &data) : IDA(data, MPI_COMM_SELF) @@ -306,12 +235,13 @@ namespace SUNDIALS auto *src_yp = internal::unwrap_nvector_const(yp); auto *residual = internal::unwrap_nvector(rr); - return call_and_possibly_capture_exception(solver.residual, - solver.pending_exception, - tt, - *src_yy, - *src_yp, - *residual); + return Utilities::call_and_possibly_capture_exception( + solver.residual, + solver.pending_exception, + tt, + *src_yy, + *src_yp, + *residual); }, current_time, yy, @@ -421,16 +351,18 @@ namespace SUNDIALS auto *src_b = internal::unwrap_nvector_const(b); auto *dst_x = internal::unwrap_nvector(x); if (solver.solve_with_jacobian) - return call_and_possibly_capture_exception(solver.solve_with_jacobian, - solver.pending_exception, - *src_b, - *dst_x, - tol); + return Utilities::call_and_possibly_capture_exception( + solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); else if (solver.solve_jacobian_system) - return call_and_possibly_capture_exception(solver.solve_jacobian_system, - solver.pending_exception, - *src_b, - *dst_x); + return Utilities::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. @@ -510,12 +442,13 @@ namespace SUNDIALS auto *src_yy = internal::unwrap_nvector_const(yy); auto *src_yp = internal::unwrap_nvector_const(yp); - return call_and_possibly_capture_exception(solver.setup_jacobian, - solver.pending_exception, - tt, - *src_yy, - *src_yp, - cj); + return Utilities::call_and_possibly_capture_exception( + solver.setup_jacobian, + solver.pending_exception, + tt, + *src_yy, + *src_yp, + cj); }); AssertIDA(status); status = IDASetMaxOrd(ida_mem, data.maximum_order); diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index a2bd709493..aaaa7c0103 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -37,6 +37,7 @@ # endif # include +# include // Make sure we #include the SUNDIALS config file... # include @@ -52,77 +53,6 @@ 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 (...) - { - eptr = std::current_exception(); - return -1; - } - } - } // namespace - template KINSOL::AdditionalData::AdditionalData( const SolutionStrategy &strategy, @@ -344,11 +274,11 @@ namespace SUNDIALS Assert(solver.iteration_function, ExcInternalError()); - const int err = - call_and_possibly_capture_exception(solver.iteration_function, - solver.pending_exception, - *src_yy, - *dst_FF); + const int err = Utilities::call_and_possibly_capture_exception( + solver.iteration_function, + solver.pending_exception, + *src_yy, + *dst_FF); return err; }, @@ -366,7 +296,7 @@ namespace SUNDIALS Assert(solver.residual, ExcInternalError()); - const int err = call_and_possibly_capture_exception( + const int err = Utilities::call_and_possibly_capture_exception( solver.residual, solver.pending_exception, *src_yy, *dst_FF); return err; @@ -452,12 +382,12 @@ namespace SUNDIALS auto src_b = internal::unwrap_nvector_const(b); auto dst_x = internal::unwrap_nvector(x); - const int err = - call_and_possibly_capture_exception(solver.solve_with_jacobian, - solver.pending_exception, - *src_b, - *dst_x, - tol); + const int err = Utilities::call_and_possibly_capture_exception( + solver.solve_with_jacobian, + solver.pending_exception, + *src_b, + *dst_x, + tol); return err; } @@ -481,7 +411,7 @@ namespace SUNDIALS // src_ycur and src_fcur, and so we simply pass dummy vector in. // These vectors will have zero lengths because we don't reinit // them above. - const int err = call_and_possibly_capture_exception( + const int err = Utilities::call_and_possibly_capture_exception( solver.solve_jacobian_system, solver.pending_exception, *src_ycur, @@ -552,10 +482,8 @@ namespace SUNDIALS auto fcur = internal::unwrap_nvector(f); // Call the user-provided setup function with these arguments: - return call_and_possibly_capture_exception(solver.setup_jacobian, - solver.pending_exception, - *ycur, - *fcur); + return Utilities::call_and_possibly_capture_exception( + solver.setup_jacobian, solver.pending_exception, *ycur, *fcur); }); AssertKINSOL(status); } -- 2.39.5