From 7031f0fffa025ce3c3e7dd0de44084c91fdd1463 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 5 Nov 2021 11:24:34 -0400 Subject: [PATCH] Revert "Fixed naming in IDA." This reverts commit 0b157b0d0b680a258bc51ffafe8bc8c4df5038ee. --- include/deal.II/sundials/ida.h | 49 +++++++++++++++++++++++++++------- source/sundials/ida.cc | 47 +++++++++++++++++++++++++++----- tests/sundials/ida_01.cc | 8 ++++-- tests/sundials/ida_02.cc | 9 ++++--- 4 files changed, 92 insertions(+), 21 deletions(-) diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index b8c2ef02de..8726e9b157 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -72,11 +72,11 @@ namespace SUNDIALS * - reinit_vector; * - residual; * - setup_jacobian; - * - solve_jacobian_system/solve_with_jacobian; + * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance; * * The function `solve_jacobian_system` should be implemented for SUNDIALS * < 4.0.0. For later versions, you should use - * `solve_with_jacobian` to leverage better non-linear + * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear * algorithms. * * Optionally, also the following functions could be provided. By default @@ -607,6 +607,21 @@ namespace SUNDIALS */ ~IDA(); + /** + * Save the number of iterations of the last Jacobian solve. + */ + void + set_n_iterations(const int n_iter); + + /** + * Return the number of iterations of the last Jacobian solve. This + * piece of information corresponds to the what the + * solve_jacobian_system_up_to_tolerance() function returns through + * its third argument. + */ + int + get_n_iterations() const; + /** * Integrate differential-algebraic equations. This function returns the * final number of computed steps. @@ -664,7 +679,7 @@ namespace SUNDIALS * update is required. The user should compute the Jacobian (or update all * the variables that allow the application of the Jacobian). This function * is called by IDA once, before any call to solve_jacobian_system() (for - * SUNDIALS < 4.0.0) or solve_with_jacobian() (for + * SUNDIALS < 4.0.0) or solve_jacobian_system_up_to_tolerance() (for * SUNDIALS >= 4.0.0). * * The Jacobian $J$ should be a (possibly inexact) computation of @@ -677,13 +692,13 @@ namespace SUNDIALS * 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 - * solve_with_jacobian() can assume that this function has + * solve_jacobian_system_up_to_tolerance() can assume that this function has * been called at least once. * * Notice that no assumption is made by this interface on what the user * should do in this function. IDA only assumes that after a call to * setup_jacobian() it is possible to call solve_jacobian_system() or - * solve_with_jacobian() to obtain a solution $x$ to the + * solve_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the * system $J x = b$. * * This function should return: @@ -731,7 +746,9 @@ namespace SUNDIALS * specifying the tolerance for the resolution. A part from the tolerance * only `rhs` is provided and `dst` needs to be returned. */ +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) DEAL_II_DEPRECATED_EARLY +# endif std::function solve_jacobian_system; @@ -744,7 +761,7 @@ namespace SUNDIALS * setup_jacobian() again. If, on the contrary, internal IDA convergence * tests fail, then IDA calls again setup_jacobian() with updated vectors * and coefficients so that successive calls to - * solve_with_jacobian() lead to better convergence in the + * solve_jacobian_system_up_to_tolerance() lead to better convergence in the * Newton process. * * The Jacobian $J$ should be (an approximation of) the system Jacobian @@ -757,6 +774,11 @@ namespace SUNDIALS * * @param[in] rhs The system right hand side to solve for. * @param[out] dst The solution of $J^{-1} * src$. + * @param[out] n_iter the number of iterations required to solve the + * Jacobian system. This is an output argument through which the + * function can communicate how many iterations it took to solve the + * linear system, and that can then be queried from the outside using the + * get_n_iterations() function. * @param[in] tolerance The tolerance with which to solve the linear system * of equations. * @@ -776,9 +798,11 @@ namespace SUNDIALS * - <0: Unrecoverable error the computation will be aborted and an * assertion will be thrown. */ - std::function< - int(const VectorType &rhs, VectorType &dst, const double tolerance)> - solve_with_jacobian; + std::function + solve_jacobian_system_up_to_tolerance; /** * Process solution. This function is called by IDA at fixed time steps, @@ -881,6 +905,13 @@ namespace SUNDIALS */ void *ida_mem; + /** + * Number of iteration that were required to solve the last + * Jacobian system. This variable is set by the wrapper that calls + * `solve_jacobian_system_up_to_tolerance`. + */ + int n_iterations; + /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index d03a402876..653c98dcd6 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -180,10 +180,20 @@ 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); + if (solver.solve_jacobian_system_up_to_tolerance) + { + int n_iter = 0; + + err = solver.solve_jacobian_system_up_to_tolerance(*src_b, + *dst_x, + n_iter, + tol); + solver.set_n_iterations(n_iter > 0 ? n_iter : 1); + } else if (solver.solve_jacobian_system) - err = solver.solve_jacobian_system(*src_b, *dst_x); + { + err = solver.solve_jacobian_system(*src_b, *dst_x); + } else // We have already checked this outside, so we should never get here. Assert(false, ExcInternalError()); @@ -223,6 +233,24 @@ namespace SUNDIALS + template + void + IDA::set_n_iterations(const int n_iter) + { + n_iterations = n_iter; + } + + + + template + int + IDA::get_n_iterations() const + { + return n_iterations; + } + + + template unsigned int IDA::solve_dae(VectorType &solution, VectorType &solution_dot) @@ -231,6 +259,7 @@ namespace SUNDIALS double h = data.initial_step_size; unsigned int step_number = 0; + this->n_iterations = 1; int status; (void)status; @@ -382,9 +411,10 @@ namespace SUNDIALS return 0; }; - AssertThrow(solve_jacobian_system || solve_with_jacobian, - ExcFunctionNotProvided( - "solve_jacobian_system or solve_with_jacobian")); + AssertThrow( + solve_jacobian_system || solve_jacobian_system_up_to_tolerance, + ExcFunctionNotProvided( + "solve_jacobian_system or solve_jacobian_system_up_to_tolerance")); LS->ops->solve = t_dae_solve_jacobian_system; // When we set an iterative solver IDA requires that resid is provided. From @@ -400,7 +430,10 @@ namespace SUNDIALS // When we set an iterative solver IDA requires that last number of // iteration is provided. Since we can't know what kind of solver the user // has provided we set 1. This is clearly suboptimal. - LS->ops->numiters = [](SUNLinearSolver /*ignored*/) -> int { return 1; }; + LS->ops->numiters = [](SUNLinearSolver LS) -> int { + IDA &solver = *static_cast *>(LS->content); + return solver.get_n_iterations(); + }; // Even though we don't use it, IDA still wants us to set some // kind of matrix object for the nonlinear solver. This is because // if we don't set it, it won't call the functions that set up diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index 7d87619d3c..be104c411e 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -107,9 +107,13 @@ public: }; // Used in ver >= 4.0.0 - time_stepper.solve_with_jacobian = - [&](const VectorType &src, VectorType &dst, const double) -> int { + time_stepper.solve_jacobian_system_up_to_tolerance = + [&](const VectorType &src, + VectorType & dst, + int & n_iter, + const double) -> int { Jinv.vmult(dst, src); + n_iter = 1; return 0; }; diff --git a/tests/sundials/ida_02.cc b/tests/sundials/ida_02.cc index 73e556ff98..700f676be4 100644 --- a/tests/sundials/ida_02.cc +++ b/tests/sundials/ida_02.cc @@ -101,12 +101,15 @@ public: return 0; }; - time_stepper.solve_with_jacobian = [&](const VectorType &src, - VectorType & dst, - const double tolerance) -> int { + time_stepper.solve_jacobian_system_up_to_tolerance = + [&](const VectorType &src, + VectorType & dst, + int & n_iter, + const double tolerance) -> int { SolverControl solver_control(1000, tolerance); SolverGMRES> solver(solver_control); solver.solve(J, dst, src, PreconditionIdentity()); + n_iter = solver_control.last_step() > 0 ? solver_control.last_step() : 1; return 0; }; -- 2.39.5