From: Wolfgang Bangerth Date: Thu, 8 Apr 2021 19:30:13 +0000 (-0600) Subject: Update the KINSOL callbacks to what KINSOL current provides. X-Git-Tag: v9.3.0-rc1~215^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bf0a6fba004aa2f07054f276e1bca66e0f2a2f52;p=dealii.git Update the KINSOL callbacks to what KINSOL current provides. --- diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index f37bc2723c..f0d29b970d 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -510,24 +510,29 @@ namespace SUNDIALS setup_jacobian; /** + * @deprecated Versions of SUNDIALS after 4.0 no longer provide all + * of the information necessary for this callback (see below). Use the + * `solve_with_jacobian` callback described below. + * * A function object that users may supply and that is intended to solve - * the Jacobian linear system. This function will be called by KINSOL - * (possibly several times) after setup_jacobian() has been called at least - * once. KINSOL tries to do its best to call setup_jacobian() the minimum - * number of times. If convergence can be achieved without updating the - * Jacobian, then KINSOL does not call setup_jacobian() again. If, on the - * contrary, internal KINSOL convergence tests fail, then KINSOL calls + * a linear system with the Jacobian matrix. This function will be called by + * KINSOL (possibly several times) after setup_jacobian() has been called at + * least once. KINSOL tries to do its best to call setup_jacobian() the + * minimum number of times. If convergence can be achieved without updating + * the Jacobian, then KINSOL does not call setup_jacobian() again. If, on + * the contrary, internal KINSOL convergence tests fail, then KINSOL calls * setup_jacobian() again with updated vectors and coefficients so that * successive calls to solve_jacobian_systems() lead to better convergence * in the Newton process. * - * If you do not specify a solve_jacobian_system() function, then only a - * fixed point iteration strategy can be used. Notice that this may not - * converge, or may converge very slowly. + * If you do not specify a `solve_jacobian_system` or `solve_with_jacobian` + * function, then only a fixed point iteration strategy can be used. Notice + * that this may not converge, or may converge very slowly. * * A call to this function should store in `dst` the result of $J^{-1}$ * applied to `rhs`, i.e., `J*dst = rhs`. It is the user's responsibility - * to set up proper solvers and preconditioners inside this function. + * to set up proper solvers and preconditioners inside this function + * (or in the `setup_jacobian` callback above). * * * Arguments to the function are: @@ -547,7 +552,7 @@ namespace SUNDIALS * - <0: Unrecoverable error the computation will be aborted and an * assertion will be thrown. * - * @warning Starting with SUNDIALS 4.0, SUNDIALS no longer provides the + * @warning Starting with SUNDIALS 4.1, SUNDIALS no longer provides the * `ycur` and `fcur` variables -- only `rhs` is provided and `dst` * needs to be returned. The first two arguments will therefore be * empty vectors in that case. In practice, that means that one @@ -561,12 +566,56 @@ namespace SUNDIALS * AdditionalData::maximum_newton_step variable to one, indicating * that the Jacobian should be re-computed in every iteration. */ + DEAL_II_DEPRECATED_EARLY std::function solve_jacobian_system; + /** + * A function object that users may supply and that is intended to solve + * a linear system with the Jacobian matrix. This function will be called by + * KINSOL (possibly several times) after setup_jacobian() has been called at + * least once. KINSOL tries to do its best to call setup_jacobian() the + * minimum number of times. If convergence can be achieved without updating + * the Jacobian, then KINSOL does not call setup_jacobian() again. If, on + * the contrary, internal KINSOL convergence tests fail, then KINSOL calls + * setup_jacobian() again with updated vectors and coefficients so that + * successive calls to solve_jacobian_systems() lead to better convergence + * in the Newton process. + * + * If you do not specify a `solve_with_jacobian` function, then only a + * fixed point iteration strategy can be used. Notice that this may not + * converge, or may converge very slowly. + * + * A call to this function should store in `dst` the result of $J^{-1}$ + * applied to `rhs`, i.e., `J*dst = rhs`. It is the user's responsibility + * to set up proper solvers and preconditioners inside this function + * (or in the `setup_jacobian` callback above). The function attached + * to this callback is also provided with a tolerance to the linear solver, + * indicating that it is not necessary to solve the linear system with + * the Jacobian matrix exactly, but only to a tolerance that KINSOL will + * adapt over time. + * + * Arguments to the function are: + * + * @param[in] rhs The system right hand side to solve for. + * @param[out] dst The solution of $J^{-1} * src$. + * @param[in] tolerance The tolerance with which to solve the linear system + * of equations. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error (KINSOL will try to change its internal + * parameters and attempt a new solution step) + * - <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; + /** * A function object that users may supply and that is intended to return a * vector whose components are the weights used by KINSOL to compute the diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index 0c27dc2369..d8009ce468 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -266,7 +266,7 @@ namespace SUNDIALS SUNMatrix /*ignored*/, N_Vector x, N_Vector b, - realtype /*tol*/) + realtype tol) { // Receive the object that describes the linear solver and // unpack the pointer to the KINSOL object from which we can then @@ -274,30 +274,59 @@ namespace SUNDIALS const KINSOL &solver = *static_cast *>(LS->content); - // Allocate temporary (deal.II-type) vectors into which to copy the - // N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_ycur(mem); - typename VectorMemory::Pointer src_fcur(mem); - typename VectorMemory::Pointer src_b(mem); - typename VectorMemory::Pointer dst_x(mem); + // This is where we have to make a decision about which of the two + // signals to call. Let's first check the more modern one: + if (solver.solve_with_jacobian) + { + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); - solver.reinit_vector(*src_b); - solver.reinit_vector(*dst_x); + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); - copy(*src_b, b); + copy(*src_b, b); - // Call the user-provided setup function with these arguments. Note - // that Sundials 4.x and later no longer provide values for - // 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 = - solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x); + const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol); - copy(x, *dst_x); + copy(x, *dst_x); - return err; + return err; + } + else + { + // User has not provided the modern callback, so the fact that we are + // here means that they must have given us something for the old + // signal. Check this. + Assert(solver.solve_jacobian_system, ExcInternalError()); + + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_ycur(mem); + typename VectorMemory::Pointer src_fcur(mem); + typename VectorMemory::Pointer src_b(mem); + typename VectorMemory::Pointer dst_x(mem); + + solver.reinit_vector(*src_b); + solver.reinit_vector(*dst_x); + + copy(*src_b, b); + + // Call the user-provided setup function with these arguments. Note + // that Sundials 4.x and later no longer provide values for + // 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 = + solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x); + + copy(x, *dst_x); + + return err; + } } # endif @@ -428,8 +457,9 @@ namespace SUNDIALS SUNMatrix J = nullptr; SUNLinearSolver LS = nullptr; - if (solve_jacobian_system) // user assigned a function object to the solver - // slot + if (solve_jacobian_system || + solve_with_jacobian) // user assigned a function object to the solver + // slot { # if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) auto KIN_mem = static_cast(kinsol_mem);