From 232e7b2d3cc2c6560b0c28ea7fcb214aba43462c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 28 Feb 2021 17:53:56 +0100 Subject: [PATCH] Make KINSOL available for SUNDIALS versions >4.x. --- include/deal.II/sundials/kinsol.h | 61 +++++--- source/sundials/kinsol.cc | 224 +++++++++++++++++++++++++----- 2 files changed, 228 insertions(+), 57 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 82ba7f519c..71dd9ad5f4 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -19,33 +19,34 @@ #include #ifdef DEAL_II_WITH_SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# include -# include -# include -# include +# include +# include +# include +# include +# include -# include -# include +# include +# include -# include +# include -# include +# include +# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) # include -# include -# include -# include +# endif +# include +# include +# include -# include +# include DEAL_II_NAMESPACE_OPEN // Shorthand notation for KINSOL error codes. -# define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code)) +# define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code)) namespace SUNDIALS { @@ -357,7 +358,9 @@ namespace SUNDIALS * The maximum number of nonlinear iterations that can be * performed between calls to the setup_jacobian() function. * - * If set to zero, default values provided by KINSOL will be used. + * If set to zero, default values provided by KINSOL will be used, + * and in practice this often means that KINSOL will re-use a + * Jacobian matrix computed in one iteration for later iterations. */ unsigned int maximum_setup_calls; @@ -452,7 +455,7 @@ namespace SUNDIALS * - 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 + * - <0: Unrecoverable error; the computation will be aborted and an * assertion will be thrown. */ std::function @@ -529,13 +532,13 @@ namespace SUNDIALS * * Arguments to the function are: * - * @param[in] ycur is the current $y$ vector for the current KINSOL + * @param[in] ycur The current $y$ vector for the current KINSOL * internal step. In the documentation above, this $y$ vector is generally * denoted by $u$. - * @param[in] fcur is the current value of the implicit right-hand side at + * @param[in] fcur The current value of the implicit right-hand side at * `ycur`, $f_I (t_n, ypred)$. - * @param[in] rhs the system right hand side to solve for - * @param[out] dst the solution of $A^{-1} * src$ + * @param[in] rhs The system right hand side to solve for + * @param[out] dst The solution of $J^{-1} * src$ * * This function should return: * - 0: Success @@ -543,6 +546,20 @@ namespace SUNDIALS * parameters and attempt a new solution step) * - <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 + * `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 + * can no longer compute a Jacobian matrix for the current iterate + * within this function. Rather, this has to happen inside the + * `setup_jacobian` function above that receives this information. + * If it is important that the Jacobian corresponds to the *current* + * iterate (rather than a re-used Jacobian matrix that had been + * computed in a previous iteration and that therefore corresponds + * to a *previous* iterate), then you will also have to set the + * AdditionalData::maximum_newton_step variable to one, indicating + * that the Jacobian should be re-computed in every iteration. */ std::function -# include - -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif -# include +# include // Make sure we #include the SUNDIALS config file... -# include +# include // ...before the rest of the SUNDIALS files: -# include -# include -# include +# include +# include +# include -# include -# include +# include +# include DEAL_II_NAMESPACE_OPEN @@ -162,6 +160,7 @@ namespace SUNDIALS +# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) template int t_kinsol_setup_jacobian(KINMem kinsol_mem) @@ -224,8 +223,88 @@ namespace SUNDIALS return err; } + +# else // SUNDIALS 5.0 or later + + template + int + t_kinsol_setup_jacobian(N_Vector u, + N_Vector f, + SUNMatrix /* ignored */, + void *user_data, + N_Vector /* tmp1 */, + N_Vector /* tmp2 */) + { + // Receive the object that describes the linear solver and + // unpack the pointer to the KINSOL object from which we can then + // get the 'setup' function. + const KINSOL &solver = + *static_cast *>(user_data); + + // Allocate temporary (deal.II-type) vectors into which to copy the + // N_vectors + GrowingVectorMemory mem; + typename VectorMemory::Pointer ycur(mem); + typename VectorMemory::Pointer fcur(mem); + solver.reinit_vector(*ycur); + solver.reinit_vector(*fcur); + + copy(*ycur, u); + copy(*fcur, f); + + // Call the user-provided setup function with these arguments: + solver.setup_jacobian(*ycur, *fcur); + + return 0; + } + + + + template + int + t_kinsol_solve_jacobian(SUNLinearSolver LS, + SUNMatrix /*ignored*/, + N_Vector x, + N_Vector b, + realtype /*tol*/) + { + // Receive the object that describes the linear solver and + // unpack the pointer to the KINSOL object from which we can then + // get the 'reinit' and 'solve' functions. + 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); + + 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 } // namespace + + template KINSOL::KINSOL(const AdditionalData &data, const MPI_Comm & mpi_comm) @@ -248,14 +327,15 @@ namespace SUNDIALS { if (kinsol_mem) KINFree(&kinsol_mem); -# ifdef DEAL_II_WITH_MPI + +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const int ierr = MPI_Comm_free(&communicator); (void)ierr; AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr)); } -# endif +# endif } @@ -269,7 +349,7 @@ namespace SUNDIALS // The solution is stored in // solution. Here we take only a // view of it. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = initial_guess_and_solution.locally_owned_elements(); @@ -285,7 +365,7 @@ namespace SUNDIALS N_VConst_Parallel(1.e0, f_scale); } else -# endif +# endif { Assert(is_serial_vector::value, ExcInternalError( @@ -347,12 +427,86 @@ namespace SUNDIALS SUNMatrix J = nullptr; SUNLinearSolver LS = nullptr; - if (solve_jacobian_system) + if (solve_jacobian_system) // 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); KIN_mem->kin_lsolve = t_kinsol_solve_jacobian; - if (setup_jacobian) + if (setup_jacobian) // user assigned a function object to the Jacobian + // set-up slot KIN_mem->kin_lsetup = t_kinsol_setup_jacobian; +# else + // Set the operations we care for in the sun_linear_solver object + // and attach it to the KINSOL object. The functions that will get + // called do not actually receive the KINSOL object, just the LS + // object, so we have to store a pointer to the current + // object in the LS object + LS = SUNLinSolNewEmpty(); + LS->content = this; + + LS->ops->gettype = + [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type { + return SUNLINEARSOLVER_MATRIX_ITERATIVE; + }; + + LS->ops->free = [](SUNLinearSolver LS) -> int { + if (LS->content) + { + LS->content = nullptr; + } + if (LS->ops) + { + free(LS->ops); + LS->ops = nullptr; + } + free(LS); + LS = nullptr; + return 0; + }; + + LS->ops->solve = t_kinsol_solve_jacobian; + + // Even though we don't use it, KINSOL 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 + // the matrix object (i.e., the argument to the 'KINSetJacFn' + // function below). + J = SUNMatNewEmpty(); + J->content = this; + + J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { + return SUNMATRIX_CUSTOM; + }; + + J->ops->destroy = [](SUNMatrix A) { + if (A->content) + { + A->content = nullptr; + } + if (A->ops) + { + free(A->ops); + A->ops = nullptr; + } + free(A); + A = nullptr; + }; + + // Now set the linear system and Jacobian objects in the solver: + status = KINSetLinearSolver(kinsol_mem, LS, J); + AssertKINSOL(status); + + // Finally, if we were given a set-up function, tell KINSOL about + // it as well. The manual says that this must happen *after* + // calling KINSetLinearSolver + if (setup_jacobian) + { + status = + KINSetJacFn(kinsol_mem, &t_kinsol_setup_jacobian); + AssertKINSOL(status); + } +# endif } else { @@ -377,7 +531,7 @@ namespace SUNDIALS copy(initial_guess_and_solution, solution); // Free the vectors which are no longer used. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(solution); @@ -385,7 +539,7 @@ namespace SUNDIALS N_VDestroy_Parallel(f_scale); } else -# endif +# endif { N_VDestroy_Serial(solution); N_VDestroy_Serial(u_scale); @@ -415,24 +569,24 @@ namespace SUNDIALS template class KINSOL>; template class KINSOL>; -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI -# ifdef DEAL_II_WITH_TRILINOS +# ifdef DEAL_II_WITH_TRILINOS template class KINSOL; template class KINSOL; -# endif +# endif -# ifdef DEAL_II_WITH_PETSC -# ifndef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifndef PETSC_USE_COMPLEX template class KINSOL; template class KINSOL; -# endif # endif - # endif +# endif + } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -# endif + #endif -- 2.39.5