From 94dfd734e0091c30b68ab1c4b8c885fb0c01d2c0 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 5 Nov 2021 11:24:56 -0400 Subject: [PATCH] Revert "Restructured and renamed functions." This reverts commit 64d87afe535256ed6ba460afc3dc2e366a243a8d. --- include/deal.II/sundials/ida.h | 74 +++++++--------- include/deal.II/sundials/n_vector.templates.h | 6 +- source/sundials/ida.cc | 87 +++++++++---------- tests/sundials/ida_01.cc | 12 ++- tests/sundials/ida_02.cc | 9 +- 5 files changed, 86 insertions(+), 102 deletions(-) diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 5c549b7956..5d006d969f 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -72,14 +72,9 @@ namespace SUNDIALS * - reinit_vector; * - residual; * - setup_jacobian; - * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance; + * - solve_jacobian_system; * - * The function `solve_jacobian_system` should be implemented for SUNDIALS - * < 4.0.0. For later versions, you should use - * `solve_jacobian_system_up_to_tolerance` to leverage better non-linear - * algorithms. - * - * Optionally, also the following functions could be provided. By default + * Optionally, also the following functions could be rewritten. By default * they do nothing, or are not required. If you call the constructor in a way * that requires a not-implemented function, an Assertion will be * thrown. @@ -600,25 +595,25 @@ namespace SUNDIALS * @param mpi_comm MPI communicator */ IDA(const AdditionalData &data = AdditionalData(), - const MPI_Comm & mpi_comm = MPI_COMM_WORLD); + const MPI_Comm mpi_comm = MPI_COMM_WORLD); /** * Destructor. */ ~IDA(); - /** - * Save the number of iterations of the last Jacobian solve. - */ + void - set_n_iterations(const int n_iter); + set_n_iter(const int n_iter) + { + this->n_iter = n_iter; + } - /** - * Return the number of iterations of the last Jacobian solve. - */ int - get_n_iterations() const; - + get_n_iter() const + { + return this->n_iter; + } /** * Integrate differential-algebraic equations. This function returns the * final number of computed steps. @@ -675,9 +670,7 @@ namespace SUNDIALS * Compute Jacobian. This function is called by IDA any time a Jacobian * 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_jacobian_system_up_to_tolerance() (for - * SUNDIALS >= 4.0.0). + * is called by IDA once, before any call to solve_jacobian_system(). * * The Jacobian $J$ should be a (possibly inexact) computation of * \f[ @@ -688,15 +681,13 @@ namespace SUNDIALS * If the user uses a matrix based computation of the Jacobian, than 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 - * solve_jacobian_system_up_to_tolerance() can assume that this function has - * been called at least once. + * Subsequent calls (possibly more than one) to solve_jacobian_system() 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_jacobian_system_up_to_tolerance() to obtain a solution $x$ to the - * system $J x = b$. + * setup_jacobian() it is possible to call solve_jacobian_system(), to + * obtain a solution $x$ to the system $J x = b$. * * This function should return: * - 0: Success @@ -743,22 +734,20 @@ 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; + /** - * Solve the Jacobian linear system up to a specified tolerance. This - * function will be called by IDA (possibly several times) after - * setup_jacobian() has been called at least once. IDA tries to do its best - * to call setup_jacobian() the minimum amount of times. If convergence can - * be achieved without updating the Jacobian, then IDA does not call - * 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_jacobian_system_up_to_tolerance() lead to better convergence in the + * Solve the Jacobian linear system. This function will be called by IDA + * (possibly several times) after setup_jacobian() has been called at least + * once. IDA tries to do its best to call setup_jacobian() the minimum + * amount of times. If convergence can be achieved without updating the + * Jacobian, then IDA does not call 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_jacobian_systems() lead to better convergence in the * Newton process. * * The jacobian $J$ should be (an approximation of) the system Jacobian @@ -791,7 +780,7 @@ namespace SUNDIALS VectorType & dst, int & n_iter, const double tolerance)> - solve_jacobian_system_up_to_tolerance; + solve_with_jacobian; /** * Process solution. This function is called by IDA at fixed time steps, @@ -887,7 +876,7 @@ namespace SUNDIALS /** * IDA configuration data. */ - const AdditionalData data; + AdditionalData data; /** * IDA memory object. @@ -895,9 +884,9 @@ namespace SUNDIALS void *ida_mem; /** - * Number of iteration that were required to solve the last Jacobian system + * Number of iteration required to solve the Jacobian system */ - int n_iterations; + int n_iter; /** * MPI communicator. SUNDIALS solver runs happily in @@ -906,6 +895,8 @@ namespace SUNDIALS */ MPI_Comm communicator; + + /** * Memory pool of vectors. */ @@ -924,6 +915,7 @@ namespace SUNDIALS # endif // PETSC_USE_COMPLEX # endif // DEAL_II_WITH_PETSC }; + } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index ee6294cd6d..027fd017f1 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -962,9 +962,9 @@ SUNDIALS::internal::create_empty_nvector() v->ops->nvdotprod = NVectorOperations::dot_product; v->ops->nvmaxnorm = NVectorOperations::max_norm; v->ops->nvwrmsnorm = NVectorOperations::weighted_rms_norm; - v->ops->nvmin = NVectorOperations::min_element; - v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm; - v->ops->nvl1norm = NVectorOperations::l1_norm; + v->ops->nvmin = NVectorOperations::min_element; + v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm; + v->ops->nvl1norm = NVectorOperations::l1_norm; v->ops->nvwrmsnormmask = NVectorOperations::weighted_rms_norm_mask; // v->ops->nvcompare = undef; diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 653c98dcd6..bd207108cf 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -173,30 +173,36 @@ namespace SUNDIALS SUNMatrix /*ignored*/, N_Vector x, N_Vector b, - realtype tol) + realtype /*tol*/) + { + const IDA &solver = + *static_cast *>(LS->content); + + auto *src_b = internal::unwrap_nvector_const(b); + auto *dst_x = internal::unwrap_nvector(x); + + const int err = solver.solve_jacobian_system(*src_b, *dst_x); + + return err; + } + + + template + int + t_dae_solve_with_jacobian(SUNLinearSolver LS, + SUNMatrix /*ignored*/, + N_Vector x, + N_Vector b, + realtype tol) { IDA &solver = *static_cast *>(LS->content); auto *src_b = internal::unwrap_nvector_const(b); auto *dst_x = internal::unwrap_nvector(x); - int err = 0; - 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); - } - else - // We have already checked this outside, so we should never get here. - Assert(false, ExcInternalError()); + int n_iter; + const int err = solver.solve_with_jacobian(*src_b, *dst_x, n_iter, tol); + solver.set_n_iter(n_iter > 0 ? n_iter : 1); return err; } @@ -206,7 +212,7 @@ namespace SUNDIALS template - IDA::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm) + IDA::IDA(const AdditionalData &data, const MPI_Comm mpi_comm) : data(data) , ida_mem(nullptr) , communicator(is_serial_vector::value ? @@ -233,24 +239,6 @@ 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) @@ -259,7 +247,7 @@ namespace SUNDIALS double h = data.initial_step_size; unsigned int step_number = 0; - this->n_iterations = 1; + this->n_iter = 1; int status; (void)status; @@ -411,12 +399,18 @@ namespace SUNDIALS return 0; }; - 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; - + if (solve_with_jacobian) + { + LS->ops->solve = t_dae_solve_with_jacobian; + } + else if (solve_jacobian_system) + { + LS->ops->solve = t_dae_solve_jacobian_system; + } + else + { + AssertThrow(false, ExcFunctionNotProvided("solve_with_jacobian")); + } // When we set an iterative solver IDA requires that resid is provided. From // SUNDIALS docs If an iterative method computes the preconditioned initial // residual and returns with a successful solve without performing any @@ -432,7 +426,7 @@ namespace SUNDIALS // has provided we set 1. This is clearly suboptimal. LS->ops->numiters = [](SUNLinearSolver LS) -> int { IDA &solver = *static_cast *>(LS->content); - return solver.get_n_iterations(); + return solver.get_n_iter(); }; // 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 @@ -536,7 +530,8 @@ namespace SUNDIALS GrowingVectorMemory mem; typename VectorMemory::Pointer v(mem); reinit_vector(*v); - return v->locally_owned_elements(); + const unsigned int size = v->size(); + return complete_index_set(size); }; } diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index be104c411e..a44c679751 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -99,19 +99,17 @@ public: return 0; }; - // Used only in ver < 4.0.0 + time_stepper.solve_jacobian_system = [&](const VectorType &src, VectorType & dst) -> int { Jinv.vmult(dst, src); return 0; }; - // Used in ver >= 4.0.0 - time_stepper.solve_jacobian_system_up_to_tolerance = - [&](const VectorType &src, - VectorType & dst, - int & n_iter, - const double) -> int { + time_stepper.solve_with_jacobian = [&](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 700f676be4..2da335b998 100644 --- a/tests/sundials/ida_02.cc +++ b/tests/sundials/ida_02.cc @@ -101,11 +101,10 @@ public: return 0; }; - time_stepper.solve_jacobian_system_up_to_tolerance = - [&](const VectorType &src, - VectorType & dst, - int & n_iter, - const double tolerance) -> int { + time_stepper.solve_with_jacobian = [&](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()); -- 2.39.5