]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make KINSOL available for SUNDIALS versions >4.x.
authorWolfgang Bangerth <bangerth@colostate.edu>
Sun, 28 Feb 2021 16:53:56 +0000 (17:53 +0100)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sun, 21 Mar 2021 20:38:38 +0000 (21:38 +0100)
include/deal.II/sundials/kinsol.h
source/sundials/kinsol.cc

index 82ba7f519c08aac9bbe5afbf5f292ad40c912c09..71dd9ad5f45672163a3c26c8b9c21a0f805c84de 100644 (file)
 #include <deal.II/base/config.h>
 
 #ifdef DEAL_II_WITH_SUNDIALS
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
 
 
-#    include <deal.II/base/conditional_ostream.h>
-#    include <deal.II/base/exceptions.h>
-#    include <deal.II/base/logstream.h>
-#    include <deal.II/base/mpi.h>
-#    include <deal.II/base/parameter_handler.h>
+#  include <deal.II/base/conditional_ostream.h>
+#  include <deal.II/base/exceptions.h>
+#  include <deal.II/base/logstream.h>
+#  include <deal.II/base/mpi.h>
+#  include <deal.II/base/parameter_handler.h>
 
-#    include <deal.II/lac/vector.h>
-#    include <deal.II/lac/vector_memory.h>
+#  include <deal.II/lac/vector.h>
+#  include <deal.II/lac/vector_memory.h>
 
-#    include <boost/signals2.hpp>
+#  include <boost/signals2.hpp>
 
-#    include <kinsol/kinsol.h>
+#  include <kinsol/kinsol.h>
+#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
 #    include <kinsol/kinsol_impl.h>
-#    include <nvector/nvector_serial.h>
-#    include <sundials/sundials_math.h>
-#    include <sundials/sundials_types.h>
+#  endif
+#  include <nvector/nvector_serial.h>
+#  include <sundials/sundials_math.h>
+#  include <sundials/sundials_types.h>
 
-#    include <memory>
+#  include <memory>
 
 
 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<int(const VectorType &src, VectorType &dst)>
@@ -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<int(const VectorType &ycur,
                       const VectorType &fcur,
@@ -635,7 +652,7 @@ namespace SUNDIALS
 
 
 DEAL_II_NAMESPACE_CLOSE
-#  endif
+
 #endif
 
 #endif
index 9e7d553242f15abcb0e76f743bfcbb0d375c54c7..6738c0cc01a0a62d64714310ae7f83732b586635 100644 (file)
 
 #ifdef DEAL_II_WITH_SUNDIALS
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+#  include <deal.II/base/utilities.h>
 
-#    include <deal.II/base/utilities.h>
-
-#    include <deal.II/lac/block_vector.h>
-#    ifdef DEAL_II_WITH_TRILINOS
-#      include <deal.II/lac/trilinos_parallel_block_vector.h>
-#      include <deal.II/lac/trilinos_vector.h>
-#    endif
-#    ifdef DEAL_II_WITH_PETSC
-#      include <deal.II/lac/petsc_block_vector.h>
-#      include <deal.II/lac/petsc_vector.h>
-#    endif
+#  include <deal.II/lac/block_vector.h>
+#  ifdef DEAL_II_WITH_TRILINOS
+#    include <deal.II/lac/trilinos_parallel_block_vector.h>
+#    include <deal.II/lac/trilinos_vector.h>
+#  endif
+#  ifdef DEAL_II_WITH_PETSC
+#    include <deal.II/lac/petsc_block_vector.h>
+#    include <deal.II/lac/petsc_vector.h>
+#  endif
 
-#    include <deal.II/sundials/copy.h>
+#  include <deal.II/sundials/copy.h>
 
 // Make sure we #include the SUNDIALS config file...
-#    include <sundials/sundials_config.h>
+#  include <sundials/sundials_config.h>
 // ...before the rest of the SUNDIALS files:
-#    include <kinsol/kinsol_direct.h>
-#    include <sunlinsol/sunlinsol_dense.h>
-#    include <sunmatrix/sunmatrix_dense.h>
+#  include <kinsol/kinsol_direct.h>
+#  include <sunlinsol/sunlinsol_dense.h>
+#  include <sunmatrix/sunmatrix_dense.h>
 
-#    include <iomanip>
-#    include <iostream>
+#  include <iomanip>
+#  include <iostream>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -162,6 +160,7 @@ namespace SUNDIALS
 
 
 
+#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
     template <typename VectorType>
     int
     t_kinsol_setup_jacobian(KINMem kinsol_mem)
@@ -224,8 +223,88 @@ namespace SUNDIALS
 
       return err;
     }
+
+#  else // SUNDIALS 5.0 or later
+
+    template <typename VectorType>
+    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<VectorType> &solver =
+        *static_cast<const KINSOL<VectorType> *>(user_data);
+
+      // Allocate temporary (deal.II-type) vectors into which to copy the
+      // N_vectors
+      GrowingVectorMemory<VectorType>            mem;
+      typename VectorMemory<VectorType>::Pointer ycur(mem);
+      typename VectorMemory<VectorType>::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 <typename VectorType>
+    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<VectorType> &solver =
+        *static_cast<const KINSOL<VectorType> *>(LS->content);
+
+      // Allocate temporary (deal.II-type) vectors into which to copy the
+      // N_vectors
+      GrowingVectorMemory<VectorType>            mem;
+      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
+      typename VectorMemory<VectorType>::Pointer src_fcur(mem);
+      typename VectorMemory<VectorType>::Pointer src_b(mem);
+      typename VectorMemory<VectorType>::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 <typename VectorType>
   KINSOL<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::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<KINMem>(kinsol_mem);
         KIN_mem->kin_lsolve = t_kinsol_solve_jacobian<VectorType>;
-        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<VectorType>;
+#  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<VectorType>;
+
+        // 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<VectorType>);
+            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<VectorType>::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<Vector<double>>;
   template class KINSOL<BlockVector<double>>;
 
-#    ifdef DEAL_II_WITH_MPI
+#  ifdef DEAL_II_WITH_MPI
 
-#      ifdef DEAL_II_WITH_TRILINOS
+#    ifdef DEAL_II_WITH_TRILINOS
   template class KINSOL<TrilinosWrappers::MPI::Vector>;
   template class KINSOL<TrilinosWrappers::MPI::BlockVector>;
-#      endif
+#    endif
 
-#      ifdef DEAL_II_WITH_PETSC
-#        ifndef PETSC_USE_COMPLEX
+#    ifdef DEAL_II_WITH_PETSC
+#      ifndef PETSC_USE_COMPLEX
   template class KINSOL<PETScWrappers::MPI::Vector>;
   template class KINSOL<PETScWrappers::MPI::BlockVector>;
-#        endif
 #      endif
-
 #    endif
 
+#  endif
+
 } // namespace SUNDIALS
 
 DEAL_II_NAMESPACE_CLOSE
-#  endif
+
 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.