From 610e5714ef1d3db938d1dbf2e78a890b1eaf4ec0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 15 May 2023 16:09:21 -0600 Subject: [PATCH] Adjust tests to our new interface. --- include/deal.II/sundials/sunlinsol_wrapper.h | 9 +- source/sundials/arkode.cc | 9 +- source/sundials/sunlinsol_wrapper.cc | 106 +++++++++++++++++-- tests/sundials/arkode_01.cc | 15 +-- tests/sundials/arkode_02.cc | 23 ++-- tests/sundials/arkode_03.cc | 30 ++---- tests/sundials/arkode_04.cc | 25 ++--- tests/sundials/arkode_05.cc | 29 ++--- tests/sundials/arkode_06.cc | 48 ++++----- tests/sundials/arkode_07.cc | 52 ++++----- tests/sundials/arkode_08.cc | 62 ++++------- tests/sundials/arkode_repeated_solve.cc | 18 ++-- 12 files changed, 218 insertions(+), 208 deletions(-) diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h index fa586df60b..0dc0d84282 100644 --- a/include/deal.II/sundials/sunlinsol_wrapper.h +++ b/include/deal.II/sundials/sunlinsol_wrapper.h @@ -23,6 +23,7 @@ # include +# include # include # include @@ -216,10 +217,12 @@ namespace SUNDIALS class LinearSolverWrapper { public: - explicit LinearSolverWrapper(LinearSolveFunction lsolve + explicit LinearSolverWrapper( + const LinearSolveFunction &lsolve, + std::exception_ptr & pending_exception # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - SUNContext linsol_ctx + , + SUNContext linsol_ctx # endif ); diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index 3186c28892..824899c4e7 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -383,6 +383,7 @@ namespace SUNDIALS # endif , mpi_communicator(mpi_comm) , last_end_time(data.initial_time) + , pending_exception(nullptr) { set_functions_to_trigger_an_assert(); @@ -413,6 +414,8 @@ namespace SUNDIALS (void)status; AssertARKode(status); # endif + + Assert(pending_exception == nullptr, ExcInternalError()); } @@ -662,7 +665,8 @@ namespace SUNDIALS , arkode_ctx # endif - ); + , + pending_exception); sun_linear_solver = *linear_solver; } else @@ -764,7 +768,8 @@ namespace SUNDIALS , arkode_ctx # endif - ); + , + pending_exception); sun_mass_linear_solver = *mass_solver; } else diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc index 579ff6fd2f..a3e07c9acb 100644 --- a/source/sundials/sunlinsol_wrapper.cc +++ b/source/sundials/sunlinsol_wrapper.cc @@ -52,6 +52,78 @@ namespace SUNDIALS # define AssertSundialsSolver(code) \ Assert(code >= 0, ExcSundialsSolverError(code)) + 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 internal { /** @@ -60,7 +132,7 @@ namespace SUNDIALS template struct LinearSolverContent { - LinearSolverContent() + LinearSolverContent(std::exception_ptr &pending_exception) : a_times_fn(nullptr) , preconditioner_setup(nullptr) , preconditioner_solve(nullptr) @@ -69,6 +141,7 @@ namespace SUNDIALS # endif , P_data(nullptr) , A_data(nullptr) + , pending_exception(pending_exception) {} ATimesFn a_times_fn; @@ -83,6 +156,12 @@ namespace SUNDIALS void *P_data; void *A_data; + + /** + * A reference to a location where we can store exceptions, should they + * be thrown by a linear solver object. + */ + std::exception_ptr &pending_exception; }; } // namespace internal @@ -119,7 +198,7 @@ namespace SUNDIALS N_Vector b, realtype tol) { - auto content = access_content(LS); + LinearSolverContent *content = access_content(LS); auto *src_b = unwrap_nvector_const(b); auto *dst_x = unwrap_nvector(x); @@ -140,7 +219,13 @@ namespace SUNDIALS # endif tol); - return content->lsolve(op, preconditioner, *dst_x, *src_b, tol); + return call_and_possibly_capture_exception(content->lsolve, + content->pending_exception, + op, + preconditioner, + *dst_x, + *src_b, + tol); } @@ -149,7 +234,8 @@ namespace SUNDIALS int arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/) { - auto content = access_content(LS); + LinearSolverContent *content = access_content(LS); + if (content->preconditioner_setup) return content->preconditioner_setup(content->P_data); return 0; @@ -171,7 +257,8 @@ namespace SUNDIALS int arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes) { - auto content = access_content(LS); + LinearSolverContent *content = access_content(LS); + content->A_data = A_data; content->a_times_fn = ATimes; return 0; @@ -186,7 +273,8 @@ namespace SUNDIALS PSetupFn p_setup, PSolveFn p_solve) { - auto content = access_content(LS); + LinearSolverContent *content = access_content(LS); + content->P_data = P_data; content->preconditioner_setup = p_setup; content->preconditioner_solve = p_solve; @@ -198,13 +286,15 @@ namespace SUNDIALS template internal::LinearSolverWrapper::LinearSolverWrapper( - LinearSolveFunction lsolve + const LinearSolveFunction &lsolve, + std::exception_ptr & pending_exception # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , SUNContext linsol_ctx # endif ) - : content(std::make_unique>()) + : content( + std::make_unique>(pending_exception)) { # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx); diff --git a/tests/sundials/arkode_01.cc b/tests/sundials/arkode_01.cc index e559e97ea1..b2e38d370c 100644 --- a/tests/sundials/arkode_01.cc +++ b/tests/sundials/arkode_01.cc @@ -75,23 +75,18 @@ main() double kappa = 1.0; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = y[1]; ydot[1] = -kappa * kappa * y[0]; - return 0; }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl; + }; Vector y(2); y[0] = 0; y[1] = kappa; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_02.cc b/tests/sundials/arkode_02.cc index 3059f38e07..f01a6d831c 100644 --- a/tests/sundials/arkode_02.cc +++ b/tests/sundials/arkode_02.cc @@ -68,27 +68,22 @@ main() double kappa = 1.0; - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = y[1]; ydot[1] = -kappa * kappa * y[0]; - return 0; }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - // limit the output to every 10th step and increase the precision to make - // the test more robust - if (step_number % 10 == 0) - deallog << t << ' ' << std::setprecision(7) << sol[0] << ' ' << sol[1] - << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + // limit the output to every 10th step and increase the precision to make + // the test more robust + if (step_number % 10 == 0) + deallog << t << ' ' << std::setprecision(7) << sol[0] << ' ' << sol[1] + << std::endl; + }; Vector y(2); y[0] = 0; y[1] = kappa; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_03.cc b/tests/sundials/arkode_03.cc index ba63b093b0..35c525b7a0 100644 --- a/tests/sundials/arkode_03.cc +++ b/tests/sundials/arkode_03.cc @@ -74,45 +74,38 @@ main() // Parameters double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5; - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 0; ydot[1] = 0; ydot[2] = -y[2] / eps; - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; ydot[2] = b / eps - y[2] * y[0]; - return 0; }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - // limit the output to every 10th step and increase the precision to make - // the test more robust - if (step_number % 10 == 0) - deallog << t << ' ' << std::setprecision(10) << sol[0] << ' ' << sol[1] - << ' ' << sol[2] << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + // limit the output to every 10th step and increase the precision to make + // the test more robust + if (step_number % 10 == 0) + deallog << t << ' ' << std::setprecision(10) << sol[0] << ' ' << sol[1] + << ' ' << sol[2] << std::endl; + }; // This test, for reasons I don't fully understand, generates some output // which varies between environments much more than the other ARKODE // tests. Work around it by setting a fairly stringent maximum time step. - ode.custom_setup = [&](void *arkode_mem) -> int { + ode.custom_setup = [&](void *arkode_mem) { int ierr = ARKStepSetMinStep(arkode_mem, 1e-8); AssertThrow(ierr == 0, ExcInternalError()); ierr = ARKStepSetMaxStep(arkode_mem, 1e-4); AssertThrow(ierr == 0, ExcInternalError()); ierr = ARKStepSetMaxNumSteps(arkode_mem, 5000); AssertThrow(ierr == 0, ExcInternalError()); - return 0; }; VectorType y(3); @@ -120,5 +113,4 @@ main() y[1] = v0; y[2] = w0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_04.cc b/tests/sundials/arkode_04.cc index 1f5a284e1d..3c74f7ca13 100644 --- a/tests/sundials/arkode_04.cc +++ b/tests/sundials/arkode_04.cc @@ -75,45 +75,36 @@ main() FullMatrix J(3, 3); J(2, 2) = -1.0 / eps; - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 0; ydot[1] = 0; ydot[2] = -y[2] / eps; - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; ydot[2] = b / eps - y[2] * y[0]; - return 0; }; ode.jacobian_times_vector = [&](const VectorType &src, VectorType & dst, double t, const VectorType & /*y*/, - const VectorType & /*fy*/) -> int { + const VectorType & /*fy*/) { J.vmult(dst, src); - - return 0; }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1] - << ' ' << sol[2] << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1] + << ' ' << sol[2] << std::endl; + }; Vector y(3); y[0] = u0; y[1] = v0; y[2] = w0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_05.cc b/tests/sundials/arkode_05.cc index 40d320f6d8..76b8909551 100644 --- a/tests/sundials/arkode_05.cc +++ b/tests/sundials/arkode_05.cc @@ -74,21 +74,17 @@ main() // Explicit jacobian. FullMatrix J(3, 3); - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 0; ydot[1] = 0; ydot[2] = -y[2] / eps; - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; ydot[2] = b / eps - y[2] * y[0]; - return 0; }; @@ -97,14 +93,13 @@ main() const double gamma, const VectorType &, const VectorType &, - bool &j_is_current) -> int { + bool &j_is_current) { J = 0; J(0, 0) = 1; J(1, 1) = 1; J(2, 2) = 1 + gamma / eps; J.gauss_jordan(); j_is_current = true; - return 0; }; ode.solve_jacobian_system = [&](const double t, @@ -112,23 +107,17 @@ main() const VectorType &, const VectorType &, const VectorType &src, - VectorType & dst) -> int { - J.vmult(dst, src); - return 0; - }; + VectorType & dst) { J.vmult(dst, src); }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] - << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] + << std::endl; + }; Vector y(3); y[0] = u0; y[1] = v0; y[2] = w0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_06.cc b/tests/sundials/arkode_06.cc index 35908ec72c..05221b7559 100644 --- a/tests/sundials/arkode_06.cc +++ b/tests/sundials/arkode_06.cc @@ -81,59 +81,48 @@ main() // Explicit jacobian. FullMatrix J(3, 3); - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 0; ydot[1] = 0; ydot[2] = -y[2] / eps; - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; ydot[2] = b / eps - y[2] * y[0]; - return 0; }; ode.jacobian_times_setup = - [&](realtype t, const VectorType &y, const VectorType &fy) -> int { - J = 0; - J(2, 2) = -1.0 / eps; - return 0; - }; + [&](realtype t, const VectorType &y, const VectorType &fy) { + J = 0; + J(2, 2) = -1.0 / eps; + }; ode.jacobian_times_vector = [&](const VectorType &v, VectorType & Jv, double t, const VectorType &y, - const VectorType &fy) -> int { - J.vmult(Jv, v); - return 0; - }; + const VectorType &fy) { J.vmult(Jv, v); }; ode.solve_linearized_system = [&](SUNDIALS::SundialsOperator &op, SUNDIALS::SundialsPreconditioner &, VectorType & x, const VectorType &b, - double tol) -> int { - ReductionControl control; - SolverCG solver_cg(control); - solver_cg.solve(op, x, b, PreconditionIdentity()); - return 0; - }; - - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] - << std::endl; - return 0; - }; + double tol) { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, PreconditionIdentity()); + }; + + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] + << std::endl; + }; // after 5.2.0 a special interpolation mode should be used for stiff problems #if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0) @@ -147,5 +136,4 @@ main() y[1] = v0; y[2] = w0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_07.cc b/tests/sundials/arkode_07.cc index b9cd3ad97a..5aa4fd662b 100644 --- a/tests/sundials/arkode_07.cc +++ b/tests/sundials/arkode_07.cc @@ -81,60 +81,50 @@ main() // Explicit jacobian. FullMatrix J(3, 3); - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 0; ydot[1] = 0; ydot[2] = -y[2] / eps; - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; ydot[2] = b / eps - y[2] * y[0]; - return 0; }; ode.jacobian_times_setup = - [&](realtype t, const VectorType &y, const VectorType &fy) -> int { - J = 0; - J(2, 2) = -1.0 / eps; - return 0; - }; + [&](realtype t, const VectorType &y, const VectorType &fy) { + J = 0; + J(2, 2) = -1.0 / eps; + }; ode.jacobian_times_vector = [&](const VectorType &v, VectorType & Jv, double t, const VectorType &y, - const VectorType &fy) -> int { - J.vmult(Jv, v); - return 0; - }; + const VectorType &fy) { J.vmult(Jv, v); }; ode.solve_linearized_system = [&](SUNDIALS::SundialsOperator & op, SUNDIALS::SundialsPreconditioner &prec, VectorType & x, const VectorType & b, - double tol) -> int { - ReductionControl control; - SolverCG solver_cg(control); - solver_cg.solve(op, x, b, prec); - return 0; - }; + double tol) { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, prec); + }; ode.jacobian_preconditioner_setup = [&](double t, const VectorType &y, const VectorType &fy, int jok, int & jcur, - double gamma) -> int { + double gamma) { deallog << "jacobian_preconditioner_setup called\n"; - return 0; }; ode.jacobian_preconditioner_solve = [&](double t, @@ -144,20 +134,17 @@ main() VectorType & z, double gamma, double delta, - int lr) -> int { + int lr) { deallog << "jacobian_preconditioner_solve called\n"; z = r; - return 0; }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] - << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2] + << std::endl; + }; // after 5.2.0 a special interpolation mode should be used for stiff problems #if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0) @@ -171,5 +158,4 @@ main() y[1] = v0; y[2] = w0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_08.cc b/tests/sundials/arkode_08.cc index 5fb8af7db5..c70c943eba 100644 --- a/tests/sundials/arkode_08.cc +++ b/tests/sundials/arkode_08.cc @@ -79,40 +79,32 @@ main() M(0, 0) = M(1, 1) = 2.0 / 3; M(1, 0) = M(0, 1) = 1.0 / 3; - ode.implicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) { K.vmult(ydot, y); - return 0; }; - ode.explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = 1; ydot[1] = 2; - return 0; }; ode.jacobian_times_vector = [&](const VectorType &v, VectorType & Jv, double t, const VectorType &y, - const VectorType &fy) -> int { - K.vmult(Jv, v); - return 0; - }; + const VectorType &fy) { K.vmult(Jv, v); }; const auto solve_function = [&](SUNDIALS::SundialsOperator & op, SUNDIALS::SundialsPreconditioner &prec, VectorType & x, const VectorType & b, - double tol) -> int { - ReductionControl control; - SolverCG solver_cg(control); - solver_cg.solve(op, x, b, prec); - return 0; - }; + double tol) { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, prec); + }; ode.solve_linearized_system = solve_function; @@ -120,41 +112,31 @@ main() FullMatrix M_inv(2, 2); - ode.mass_preconditioner_solve = [&](double t, - const VectorType &r, - VectorType & z, - double gamma, - int lr) -> int { - LogStream::Prefix prefix("mass_preconditioner_solve"); - deallog << "applied" << std::endl; - M_inv.vmult(z, r); - return 0; - }; + ode.mass_preconditioner_solve = + [&](double t, const VectorType &r, VectorType &z, double gamma, int lr) { + LogStream::Prefix prefix("mass_preconditioner_solve"); + deallog << "applied" << std::endl; + M_inv.vmult(z, r); + }; - ode.mass_preconditioner_setup = [&](double t) -> int { + ode.mass_preconditioner_setup = [&](double t) { LogStream::Prefix prefix("mass_preconditioner_setup"); deallog << "applied" << std::endl; M_inv.invert(M); - return 0; }; - ode.mass_times_vector = - [&](const double t, const VectorType &v, VectorType &Mv) -> int { - M.vmult(Mv, v); - return 0; - }; + ode.mass_times_vector = [&](const double t, + const VectorType &v, + VectorType & Mv) { M.vmult(Mv, v); }; - ode.output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl; - return 0; - }; + ode.output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl; + }; Vector y(2); y[0] = 1; y[1] = 0; ode.solve_ode(y); - return 0; } diff --git a/tests/sundials/arkode_repeated_solve.cc b/tests/sundials/arkode_repeated_solve.cc index a2649ba339..a096b886ff 100644 --- a/tests/sundials/arkode_repeated_solve.cc +++ b/tests/sundials/arkode_repeated_solve.cc @@ -41,20 +41,16 @@ create_solver() auto ode = std::make_unique>(data); // will yield analytic solution y[0] = sin(kappa*t); y[1] = kappa*cos(kappa*t) - ode->explicit_function = - [&](double, const VectorType &y, VectorType &ydot) -> int { + ode->explicit_function = [&](double, const VectorType &y, VectorType &ydot) { ydot[0] = y[1]; ydot[1] = -kappa * kappa * y[0]; - return 0; }; - ode->output_step = [&](const double t, - const VectorType & sol, - const unsigned int step_number) -> int { - deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1] - << std::endl; - return 0; - }; + ode->output_step = + [&](const double t, const VectorType &sol, const unsigned int step_number) { + deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1] + << std::endl; + }; return ode; } @@ -122,6 +118,4 @@ main() ode->reset(0.0, 0.01, y0); ode->solve_ode_incrementally(y, 2.0); } - - return 0; } -- 2.39.5