From: Luca Heltai Date: Sun, 16 May 2021 09:11:52 +0000 (+0200) Subject: Restructured and renamed functions. X-Git-Tag: v9.3.0-rc1~53^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=64d87afe535256ed6ba460afc3dc2e366a243a8d;p=dealii.git Restructured and renamed functions. --- diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 72dc4e87fd..80aaa6984a 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -72,9 +72,14 @@ namespace SUNDIALS * - reinit_vector; * - residual; * - setup_jacobian; - * - solve_jacobian_system; + * - solve_jacobian_system/solve_jacobian_system_up_to_tolerance; * - * Optionally, also the following functions could be rewritten. By default + * 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 * 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. @@ -595,25 +600,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_iter(const int n_iter) - { - this->n_iter = n_iter; - } + set_n_iterations(const int n_iter); + /** + * Return the number of iterations of the last Jacobian solve. + */ int - get_n_iter() const - { - return this->n_iter; - } + get_n_iterations() const; + /** * Integrate differential-algebraic equations. This function returns the * final number of computed steps. @@ -670,7 +675,9 @@ 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(). + * 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). * * The Jacobian $J$ should be a (possibly inexact) computation of * \f[ @@ -681,13 +688,15 @@ 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() can - * assume that this function has been called at least once. + * 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. * * 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(), to - * obtain a solution $x$ to the system $J x = b$. + * 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$. * * This function should return: * - 0: Success @@ -734,20 +743,22 @@ 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. 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 + * 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 * Newton process. * * The jacobian $J$ should be (an approximation of) the system Jacobian @@ -780,7 +791,7 @@ namespace SUNDIALS VectorType & dst, int & n_iter, const double tolerance)> - solve_with_jacobian; + solve_jacobian_system_up_to_tolerance; /** * Process solution. This function is called by IDA at fixed time steps, @@ -876,7 +887,7 @@ namespace SUNDIALS /** * IDA configuration data. */ - AdditionalData data; + const AdditionalData data; /** * IDA memory object. @@ -884,9 +895,9 @@ namespace SUNDIALS void *ida_mem; /** - * Number of iteration required to solve the Jacobian system + * Number of iteration that were required to solve the last Jacobian system */ - int n_iter; + int n_iterations; /** * MPI communicator. SUNDIALS solver runs happily in @@ -895,8 +906,6 @@ namespace SUNDIALS */ MPI_Comm communicator; - - /** * Memory pool of vectors. */ @@ -915,7 +924,6 @@ 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 e289f97933..4c1b8bc2ce 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 bd207108cf..653c98dcd6 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -173,36 +173,30 @@ namespace SUNDIALS SUNMatrix /*ignored*/, N_Vector x, N_Vector b, - 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) + 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; - 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); + 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()); return err; } @@ -212,7 +206,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 ? @@ -239,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) @@ -247,7 +259,7 @@ namespace SUNDIALS double h = data.initial_step_size; unsigned int step_number = 0; - this->n_iter = 1; + this->n_iterations = 1; int status; (void)status; @@ -399,18 +411,12 @@ namespace SUNDIALS return 0; }; - 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")); - } + 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 // SUNDIALS docs If an iterative method computes the preconditioned initial // residual and returns with a successful solve without performing any @@ -426,7 +432,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_iter(); + 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 @@ -530,8 +536,7 @@ namespace SUNDIALS GrowingVectorMemory mem; typename VectorMemory::Pointer v(mem); reinit_vector(*v); - const unsigned int size = v->size(); - return complete_index_set(size); + return v->locally_owned_elements(); }; } diff --git a/tests/sundials/ida_01.cc b/tests/sundials/ida_01.cc index a44c679751..be104c411e 100644 --- a/tests/sundials/ida_01.cc +++ b/tests/sundials/ida_01.cc @@ -99,17 +99,19 @@ 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; }; - time_stepper.solve_with_jacobian = [&](const VectorType &src, - VectorType & dst, - int & n_iter, - const double) -> int { + // 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 { Jinv.vmult(dst, src); n_iter = 1; return 0; diff --git a/tests/sundials/ida_02.cc b/tests/sundials/ida_02.cc index 2da335b998..700f676be4 100644 --- a/tests/sundials/ida_02.cc +++ b/tests/sundials/ida_02.cc @@ -101,10 +101,11 @@ public: return 0; }; - time_stepper.solve_with_jacobian = [&](const VectorType &src, - VectorType & dst, - int & n_iter, - 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());