]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support SUNDIALS 6.0 and later. 13438/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 23 Feb 2022 03:59:07 +0000 (20:59 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 28 Feb 2022 21:24:12 +0000 (14:24 -0700)
include/deal.II/sundials/arkode.h
include/deal.II/sundials/ida.h
include/deal.II/sundials/kinsol.h
include/deal.II/sundials/n_vector.h
include/deal.II/sundials/n_vector.templates.h
include/deal.II/sundials/sunlinsol_wrapper.h
source/sundials/arkode.cc
source/sundials/ida.cc
source/sundials/kinsol.cc
source/sundials/n_vector.inst.in
source/sundials/sunlinsol_wrapper.cc

index ad2ca7d863d543042c7a5f1e4c95b33adf3b1a5e..eb31f2e22fcc5a55aa8a587600f419fefb4c59c2 100644 (file)
@@ -1267,12 +1267,19 @@ namespace SUNDIALS
      */
     void *arkode_mem;
 
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    /**
+     * A context object associated with the ARKode solver.
+     */
+    SUNContext arkode_ctx;
+#  endif
+
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI
      * support, MPI_Comm is aliased as int.
      */
-    MPI_Comm communicator;
+    MPI_Comm mpi_communicator;
 
     /**
      * The final time in the last call to solve_ode().
index 16f6f5070cc5fee0286d749d9120559250a8dc03..bc50c9b7e4ecd5b7c9fbe1445359fe650345c52f 100644 (file)
@@ -882,12 +882,19 @@ namespace SUNDIALS
      */
     void *ida_mem;
 
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    /**
+     * A context object associated with the IDA solver.
+     */
+    SUNContext ida_ctx;
+#  endif
+
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI
      * support, MPI_Comm is aliased as int.
      */
-    MPI_Comm communicator;
+    MPI_Comm mpi_communicator;
 
     /**
      * Memory pool of vectors.
index 89c5a146c51b6522077e9f3124545c649a23dae7..df2f6b16c744e75542c04b9a2aa26a24d3513115 100644 (file)
@@ -684,11 +684,23 @@ namespace SUNDIALS
      */
     AdditionalData data;
 
+    /**
+     * The MPI communicator to be used by this solver, if any.
+     */
+    MPI_Comm mpi_communicator;
+
     /**
      * KINSOL memory object.
      */
     void *kinsol_mem;
 
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    /**
+     * A context object associated with the KINSOL solver.
+     */
+    SUNContext kinsol_ctx;
+#  endif
+
     /**
      * Memory pool of vectors.
      */
index df2c27f1298f6e4795de940f10be2fdd494c39c8..6f959444e693e46886423a0acefb79c9ce59207f 100644 (file)
@@ -56,13 +56,22 @@ namespace SUNDIALS
      * @tparam VectorType Type of the viewed vector. This parameter can be
      *   deduced automatically and will respect a potential const-qualifier.
      * @param vector The vector to view.
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+     * @param nvector_context A SUNDIALS context object to be used for the
+     *   operations we want to do on this vector.
+#  endif
      * @return A NVectorView of the @p vector.
      *
      * @related NVectorView
      */
     template <typename VectorType>
     NVectorView<VectorType>
-    make_nvector_view(VectorType &vector);
+    make_nvector_view(VectorType &vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                      ,
+                      SUNContext nvector_context
+#  endif
+    );
 
     /**
      * Retrieve the underlying vector attached to N_Vector @p v. This call will
@@ -129,7 +138,12 @@ namespace SUNDIALS
       /**
        * Constructor. Create view of @p vector.
        */
-      NVectorView(VectorType &vector);
+      NVectorView(VectorType &vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                  ,
+                  SUNContext nvector_context
+#  endif
+      );
 
       /**
        * Move assignment.
index 002eb5eb66d7e9c74e98dbe628027f4965e694b6..ee5220d023a3cc21e2ad4f9d9682ceae8c91657e 100644 (file)
@@ -125,7 +125,12 @@ namespace SUNDIALS
      */
     template <typename VectorType>
     N_Vector
-    create_nvector(NVectorContent<VectorType> *content);
+    create_nvector(NVectorContent<VectorType> *content
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                   ,
+                   SUNContext nvector_context
+#  endif
+    );
 
     /**
      * Helper to create an empty vector with all operations but no content.
@@ -133,7 +138,11 @@ namespace SUNDIALS
      */
     template <typename VectorType>
     N_Vector
-    create_empty_nvector();
+    create_empty_nvector(
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+      SUNContext nvector_context
+#  endif
+    );
 
     /**
      * Collection of all operations specified by SUNDIALS N_Vector
@@ -356,9 +365,19 @@ namespace SUNDIALS
 
     template <typename VectorType>
     NVectorView<VectorType>
-    make_nvector_view(VectorType &vector)
+    make_nvector_view(VectorType &vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                      ,
+                      SUNContext nvector_context
+#  endif
+    )
     {
-      return NVectorView<VectorType>(vector);
+      return NVectorView<VectorType>(vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                     ,
+                                     nvector_context
+#  endif
+      );
     }
 
 
@@ -390,11 +409,21 @@ namespace SUNDIALS
 
 
     template <typename VectorType>
-    NVectorView<VectorType>::NVectorView(VectorType &vector)
+    NVectorView<VectorType>::NVectorView(VectorType &vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                         ,
+                                         SUNContext nvector_context
+#  endif
+                                         )
       : vector_ptr(
           create_nvector(
             new NVectorContent<typename std::remove_const<VectorType>::type>(
-              &vector)),
+              &vector)
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+              ,
+            nvector_context
+#  endif
+            ),
           [](N_Vector v) { N_VDestroy(v); })
     {}
 
@@ -421,10 +450,19 @@ namespace SUNDIALS
 
     template <typename VectorType>
     N_Vector
-    create_nvector(NVectorContent<VectorType> *content)
+    create_nvector(NVectorContent<VectorType> *content
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                   ,
+                   SUNContext nvector_context
+#  endif
+    )
     {
       // Create an N_Vector with operators attached and empty content
+#  if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
       N_Vector v = create_empty_nvector<VectorType>();
+#  else
+      N_Vector v = create_empty_nvector<VectorType>(nvector_context);
+#  endif
       Assert(v != nullptr, ExcInternalError());
 
       v->content = content;
@@ -447,7 +485,11 @@ namespace SUNDIALS
       clone_empty(N_Vector w)
       {
         Assert(w != nullptr, ExcInternalError());
+#  if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
         N_Vector v = N_VNewEmpty();
+#  else
+        N_Vector v = N_VNewEmpty(w->sunctx);
+#  endif
         Assert(v != nullptr, ExcInternalError());
 
         int status = N_VCopyOps(w, v);
@@ -940,9 +982,17 @@ namespace SUNDIALS
 
     template <typename VectorType>
     N_Vector
-    create_empty_nvector()
+    create_empty_nvector(
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+      SUNContext nvector_context
+#  endif
+    )
     {
+#  if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
       N_Vector v = N_VNewEmpty();
+#  else
+      N_Vector v = N_VNewEmpty(nvector_context);
+#  endif
       Assert(v != nullptr, ExcInternalError());
 
       /* constructors, destructors, and utility operations */
index 4b94e1dc9b21f0f23eb482e93dcfa15d49ecfd27..7ae539e82c295871a3a10e2658abf58e1c5405d4 100644 (file)
@@ -156,7 +156,7 @@ namespace SUNDIALS
 
   namespace internal
   {
-    /*!
+    /**
      * Attach wrapper functions to SUNDIALS' linear solver interface. We pretend
      * that the user-supplied linear solver is matrix-free, even though it can
      * be matrix-based. This way SUNDIALS does not need to understand our matrix
@@ -166,7 +166,12 @@ namespace SUNDIALS
     class LinearSolverWrapper
     {
     public:
-      explicit LinearSolverWrapper(LinearSolveFunction<VectorType> lsolve);
+      explicit LinearSolverWrapper(LinearSolveFunction<VectorType> lsolve
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                   ,
+                                   SUNContext &linsol_ctx
+#    endif
+      );
 
       ~LinearSolverWrapper();
 
index c68e57abd7d3cb92fcf1ee860edd315cd4607129..83cd451abf8152976bc55efcf85349a86d2f15cd 100644 (file)
@@ -368,9 +368,9 @@ namespace SUNDIALS
                              const MPI_Comm &      mpi_comm)
     : data(data)
     , arkode_mem(nullptr)
-    , communicator(is_serial_vector<VectorType>::value ?
-                     MPI_COMM_SELF :
-                     Utilities::MPI::duplicate_communicator(mpi_comm))
+    , mpi_communicator(is_serial_vector<VectorType>::value ?
+                         MPI_COMM_SELF :
+                         Utilities::MPI::duplicate_communicator(mpi_comm))
     , last_end_time(data.initial_time)
   {
     set_functions_to_trigger_an_assert();
@@ -382,15 +382,23 @@ namespace SUNDIALS
   ARKode<VectorType>::~ARKode()
   {
     if (arkode_mem)
+      {
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-      ARKodeFree(&arkode_mem);
+        ARKodeFree(&arkode_mem);
 #  else
-      ARKStepFree(&arkode_mem);
+        ARKStepFree(&arkode_mem);
+#  endif
+
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        const int status = SUNContext_Free(&arkode_ctx);
+        (void)status;
+        AssertARKode(status);
 #  endif
+      }
 
 #  ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
-      Utilities::MPI::free_communicator(communicator);
+      Utilities::MPI::free_communicator(mpi_communicator);
 #  endif
   }
 
@@ -443,7 +451,12 @@ namespace SUNDIALS
                       time.get_step_number());
       }
 
-    auto solution_nvector = internal::make_nvector_view(solution);
+    auto solution_nvector = internal::make_nvector_view(solution
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                        ,
+                                                        arkode_ctx
+#  endif
+    );
 
     while (!time.is_at_end())
       {
@@ -510,7 +523,12 @@ namespace SUNDIALS
 
     // just a view on the memory in solution, all write operations on yy by
     // ARKODE will automatically be mirrored to solution
-    auto initial_condition_nvector = internal::make_nvector_view(solution);
+    auto initial_condition_nvector = internal::make_nvector_view(solution
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                                 ,
+                                                                 arkode_ctx
+#    endif
+    );
 
     status = ARKodeInit(
       arkode_mem,
@@ -522,8 +540,12 @@ namespace SUNDIALS
 
     if (get_local_tolerances)
       {
-        const auto abs_tols =
-          internal::make_nvector_view(get_local_tolerances());
+        const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                            ,
+                                                          arkode_ctx
+#    endif
+        );
         status =
           ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
         AssertARKode(status);
@@ -605,30 +627,58 @@ namespace SUNDIALS
   {
     last_end_time = current_time;
     if (arkode_mem)
-      ARKStepFree(&arkode_mem);
+      {
+        ARKStepFree(&arkode_mem);
+
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        const int status = SUNContext_Free(&arkode_ctx);
+        (void)status;
+        AssertARKode(status);
+#    endif
+      }
 
     int status;
     (void)status;
 
     // just a view on the memory in solution, all write operations on yy by
     // ARKODE will automatically be mirrored to solution
-    auto initial_condition_nvector = internal::make_nvector_view(solution);
+    auto initial_condition_nvector = internal::make_nvector_view(solution
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                                 ,
+                                                                 arkode_ctx
+#    endif
+    );
 
     Assert(explicit_function || implicit_function,
            ExcFunctionNotProvided("explicit_function || implicit_function"));
 
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
     arkode_mem = ARKStepCreate(
       explicit_function ? &explicit_function_callback<VectorType> : nullptr,
       implicit_function ? &implicit_function_callback<VectorType> : nullptr,
       current_time,
       initial_condition_nvector);
+#    else
+    status = SUNContext_Create(&mpi_communicator, &arkode_ctx);
+    AssertARKode(status);
 
+    arkode_mem = ARKStepCreate(
+      explicit_function ? &explicit_function_callback<VectorType> : nullptr,
+      implicit_function ? &implicit_function_callback<VectorType> : nullptr,
+      current_time,
+      initial_condition_nvector,
+      arkode_ctx);
+#    endif
     Assert(arkode_mem != nullptr, ExcInternalError());
 
     if (get_local_tolerances)
       {
-        const auto abs_tols =
-          internal::make_nvector_view(get_local_tolerances());
+        const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                            ,
+                                                          arkode_ctx
+#    endif
+        );
         status =
           ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
         AssertARKode(status);
@@ -683,18 +733,36 @@ namespace SUNDIALS
           {
             linear_solver =
               std::make_unique<internal::LinearSolverWrapper<VectorType>>(
-                solve_linearized_system);
+                solve_linearized_system
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                ,
+                arkode_ctx
+#    endif
+              );
             sun_linear_solver = *linear_solver;
           }
         else
           {
             // use default solver from SUNDIALS
             // TODO give user options
-            auto y_template = internal::make_nvector_view(solution);
+            auto y_template = internal::make_nvector_view(solution
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                          ,
+                                                          arkode_ctx
+#    endif
+            );
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
             sun_linear_solver =
               SUNLinSol_SPGMR(y_template,
                               PREC_NONE,
                               0 /*krylov subvectors, 0 uses default*/);
+#    else
+            sun_linear_solver =
+              SUNLinSol_SPGMR(y_template,
+                              PREC_NONE,
+                              0 /*krylov subvectors, 0 uses default*/,
+                              arkode_ctx);
+#    endif
           }
         status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr);
         AssertARKode(status);
@@ -724,11 +792,23 @@ namespace SUNDIALS
       }
     else
       {
-        auto y_template = internal::make_nvector_view(solution);
+        auto y_template = internal::make_nvector_view(solution
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                      ,
+                                                      arkode_ctx
+#    endif
+        );
 
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
         SUNNonlinearSolver fixed_point_solver =
           SUNNonlinSol_FixedPoint(y_template,
                                   data.anderson_acceleration_subspace);
+#    else
+        SUNNonlinearSolver fixed_point_solver =
+          SUNNonlinSol_FixedPoint(y_template,
+                                  data.anderson_acceleration_subspace,
+                                  arkode_ctx);
+#    endif
 
         status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
         AssertARKode(status);
@@ -755,16 +835,34 @@ namespace SUNDIALS
           {
             mass_solver =
               std::make_unique<internal::LinearSolverWrapper<VectorType>>(
-                solve_mass);
+                solve_mass
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                ,
+                arkode_ctx
+#    endif
+              );
             sun_mass_linear_solver = *mass_solver;
           }
         else
           {
-            auto y_template = internal::make_nvector_view(solution);
+            auto y_template = internal::make_nvector_view(solution
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                          ,
+                                                          arkode_ctx
+#    endif
+            );
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
             sun_mass_linear_solver =
               SUNLinSol_SPGMR(y_template,
                               PREC_NONE,
                               0 /*krylov subvectors, 0 uses default*/);
+#    else
+            sun_mass_linear_solver =
+              SUNLinSol_SPGMR(y_template,
+                              PREC_NONE,
+                              0 /*krylov subvectors, 0 uses default*/,
+                              arkode_ctx);
+#    endif
           }
         booleantype mass_time_dependent =
           data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
index fe212f4e3e9d93efd4b9e46bb6989b68d3ac3d44..8c886bca8d34ba503c5ebd00ddc6fa1ed0782276 100644 (file)
@@ -197,21 +197,32 @@ namespace SUNDIALS
   IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm)
     : data(data)
     , ida_mem(nullptr)
-    , communicator(is_serial_vector<VectorType>::value ?
-                     MPI_COMM_SELF :
-                     Utilities::MPI::duplicate_communicator(mpi_comm))
+    , mpi_communicator(is_serial_vector<VectorType>::value ?
+                         MPI_COMM_SELF :
+                         Utilities::MPI::duplicate_communicator(mpi_comm))
   {
     set_functions_to_trigger_an_assert();
   }
 
+
+
   template <typename VectorType>
   IDA<VectorType>::~IDA()
   {
     if (ida_mem)
-      IDAFree(&ida_mem);
+      {
+        IDAFree(&ida_mem);
+
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        const int status = SUNContext_Free(&ida_ctx);
+        (void)status;
+        AssertIDA(status);
+#  endif
+      }
+
 #  ifdef DEAL_II_WITH_MPI
     if (is_serial_vector<VectorType>::value == false)
-      Utilities::MPI::free_communicator(communicator);
+      Utilities::MPI::free_communicator(mpi_communicator);
 #  endif
   }
 
@@ -242,8 +253,18 @@ namespace SUNDIALS
       {
         next_time += data.output_period;
 
-        auto yy = internal::make_nvector_view(solution);
-        auto yp = internal::make_nvector_view(solution_dot);
+        auto yy = internal::make_nvector_view(solution
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                              ida_ctx
+#  endif
+        );
+        auto yp = internal::make_nvector_view(solution_dot
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                              ida_ctx
+#  endif
+        );
 
         status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
         AssertIDA(status);
@@ -274,23 +295,53 @@ namespace SUNDIALS
     bool first_step = (current_time == data.initial_time);
 
     if (ida_mem)
-      IDAFree(&ida_mem);
+      {
+        IDAFree(&ida_mem);
 
-    ida_mem = IDACreate();
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        const int status = SUNContext_Free(&ida_ctx);
+        (void)status;
+        AssertIDA(status);
+#  endif
+      }
 
     int status;
     (void)status;
 
-    auto yy = internal::make_nvector_view(solution);
-    auto yp = internal::make_nvector_view(solution_dot);
+
+#  if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    ida_mem = IDACreate();
+#  else
+    status = SUNContext_Create(&mpi_communicator, &ida_ctx);
+    AssertIDA(status);
+
+    ida_mem            = IDACreate(ida_ctx);
+#  endif
+
+    auto yy = internal::make_nvector_view(solution
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                          ,
+                                          ida_ctx
+#  endif
+    );
+    auto yp = internal::make_nvector_view(solution_dot
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                          ,
+                                          ida_ctx
+#  endif
+    );
 
     status =
       IDAInit(ida_mem, residual_callback<VectorType>, current_time, yy, yp);
     AssertIDA(status);
     if (get_local_tolerances)
       {
-        const auto abs_tols =
-          internal::make_nvector_view(get_local_tolerances());
+        const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                            ,
+                                                          ida_ctx
+#  endif
+        );
         status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
         AssertIDA(status);
       }
@@ -319,8 +370,13 @@ namespace SUNDIALS
           diff_comp_vector[*i] = 1.0;
         diff_comp_vector.compress(VectorOperation::insert);
 
-        const auto diff_id = internal::make_nvector_view(diff_comp_vector);
-        status             = IDASetId(ida_mem, diff_id);
+        const auto diff_id = internal::make_nvector_view(diff_comp_vector
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                         ,
+                                                         ida_ctx
+#  endif
+        );
+        status = IDASetId(ida_mem, diff_id);
         AssertIDA(status);
       }
 
@@ -352,7 +408,12 @@ namespace SUNDIALS
     // called do not actually receive the IDAMEM object, just the LS
     // object, so we have to store a pointer to the current
     // object in the LS object
-    LS          = SUNLinSolNewEmpty();
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    LS = SUNLinSolNewEmpty();
+#    else
+    LS = SUNLinSolNewEmpty(ida_ctx);
+#    endif
+
     LS->content = this;
 
     LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
@@ -398,7 +459,11 @@ namespace SUNDIALS
     // if we don't set it, it won't call the functions that set up
     // the matrix object (i.e., the argument to the 'IDASetJacFn'
     // function below).
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
     J          = SUNMatNewEmpty();
+#    else
+    J  = SUNMatNewEmpty(ida_ctx);
+#    endif
     J->content = this;
 
     J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
index d26f1bfe12a8f44ce290b117e712e7ab3cde6f32..36cbc51c4b846dfe832f6e31d88edf0098ac1b45 100644 (file)
@@ -303,8 +303,12 @@ namespace SUNDIALS
 
 
   template <typename VectorType>
-  KINSOL<VectorType>::KINSOL(const AdditionalData &data, const MPI_Comm &)
+  KINSOL<VectorType>::KINSOL(const AdditionalData &data,
+                             const MPI_Comm &      mpi_comm)
     : data(data)
+    , mpi_communicator(is_serial_vector<VectorType>::value ?
+                         MPI_COMM_SELF :
+                         Utilities::MPI::duplicate_communicator(mpi_comm))
     , kinsol_mem(nullptr)
   {
     set_functions_to_trigger_an_assert();
@@ -315,8 +319,23 @@ namespace SUNDIALS
   template <typename VectorType>
   KINSOL<VectorType>::~KINSOL()
   {
+    int status = 0;
+    (void)status;
+
     if (kinsol_mem)
-      KINFree(&kinsol_mem);
+      {
+        KINFree(&kinsol_mem);
+
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        status = SUNContext_Free(&kinsol_ctx);
+        AssertKINSOL(status);
+#  endif
+      }
+
+#  ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      Utilities::MPI::free_communicator(mpi_communicator);
+#  endif
   }
 
 
@@ -330,21 +349,41 @@ namespace SUNDIALS
     VectorType u_scale_temp, f_scale_temp;
 
     if (get_solution_scaling)
-      u_scale = internal::make_nvector_view(get_solution_scaling());
+      u_scale = internal::make_nvector_view(get_solution_scaling()
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                            kinsol_ctx
+#  endif
+      );
     else
       {
         reinit_vector(u_scale_temp);
         u_scale_temp = 1.0;
-        u_scale      = internal::make_nvector_view(u_scale_temp);
+        u_scale      = internal::make_nvector_view(u_scale_temp
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                              kinsol_ctx
+#  endif
+        );
       }
 
     if (get_function_scaling)
-      f_scale = internal::make_nvector_view(get_function_scaling());
+      f_scale = internal::make_nvector_view(get_function_scaling()
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                            kinsol_ctx
+#  endif
+      );
     else
       {
         reinit_vector(f_scale_temp);
         f_scale_temp = 1.0;
-        f_scale      = internal::make_nvector_view(f_scale_temp);
+        f_scale      = internal::make_nvector_view(f_scale_temp
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                              ,
+                                              kinsol_ctx
+#  endif
+        );
       }
 
     // Make sure we have what we need
@@ -361,15 +400,34 @@ namespace SUNDIALS
                  "solve_jacobian_system || solve_with_jacobian"));
       }
 
-    auto solution = internal::make_nvector_view(initial_guess_and_solution);
+    auto solution = internal::make_nvector_view(initial_guess_and_solution
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                                ,
+                                                kinsol_ctx
+#  endif
+    );
+
+    int status = 0;
+    (void)status;
 
     if (kinsol_mem)
-      KINFree(&kinsol_mem);
+      {
+        KINFree(&kinsol_mem);
 
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+        status = SUNContext_Free(&kinsol_ctx);
+        AssertKINSOL(status);
+#  endif
+      }
+
+#  if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
     kinsol_mem = KINCreate();
+#  else
+    status = SUNContext_Create(&mpi_communicator, &kinsol_ctx);
+    AssertKINSOL(status);
 
-    int status = 0;
-    (void)status;
+    kinsol_mem = KINCreate(kinsol_ctx);
+#  endif
 
     status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
     AssertKINSOL(status);
@@ -443,7 +501,11 @@ namespace SUNDIALS
         // 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
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
         LS          = SUNLinSolNewEmpty();
+#    else
+        LS = SUNLinSolNewEmpty(kinsol_ctx);
+#    endif
         LS->content = this;
 
         LS->ops->gettype =
@@ -473,7 +535,11 @@ namespace SUNDIALS
         // 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).
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
         J          = SUNMatNewEmpty();
+#    else
+        J  = SUNMatNewEmpty(kinsol_ctx);
+#    endif
         J->content = this;
 
         J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
index 899195b7950c0747205f71c7c97489c9423dda75..53757bffaacb44636d8ea41b991f44c79d1ab3fd 100644 (file)
 for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
   {
     template SUNDIALS::internal::NVectorView<V<S>>
-    SUNDIALS::internal::make_nvector_view<>(V<S> &);
+    SUNDIALS::internal::make_nvector_view<>(V<S> &
+#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#endif
+    );
     template SUNDIALS::internal::NVectorView<const V<S>>
-                         SUNDIALS::internal::make_nvector_view<>(const V<S> &);
+    SUNDIALS::internal::make_nvector_view<>(const V<S> &
+#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#endif
+    );
     template V<S> *      SUNDIALS::internal::unwrap_nvector<V<S>>(N_Vector);
     template const V<S> *SUNDIALS::internal::unwrap_nvector_const<V<S>>(
       N_Vector);
@@ -32,11 +42,21 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
 for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
   {
     template SUNDIALS::internal::NVectorView<LinearAlgebra::distributed::V<S>>
-    SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V<S> &);
+    SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V<S> &
+#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#endif
+    );
     template SUNDIALS::internal::NVectorView<
       const LinearAlgebra::distributed::V<S>>
     SUNDIALS::internal::make_nvector_view<>(
-      const LinearAlgebra::distributed::V<S> &);
+      const LinearAlgebra::distributed::V<S> &
+#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+      ,
+      SUNContext
+#endif
+    );
     template LinearAlgebra::distributed::V<S>
       *SUNDIALS::internal::unwrap_nvector<LinearAlgebra::distributed::V<S>>(
         N_Vector);
@@ -54,9 +74,19 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
 for (V : DEAL_II_VEC_TEMPLATES)
   {
     template SUNDIALS::internal::NVectorView<TrilinosWrappers::MPI::V>
-    SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V &);
+    SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V &
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#  endif
+    );
     template SUNDIALS::internal::NVectorView<const TrilinosWrappers::MPI::V>
-    SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V &);
+    SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V &
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#  endif
+    );
     template TrilinosWrappers::MPI::V
       *SUNDIALS::internal::unwrap_nvector<TrilinosWrappers::MPI::V>(N_Vector);
     template const TrilinosWrappers::MPI::V
@@ -74,9 +104,19 @@ for (V : DEAL_II_VEC_TEMPLATES)
 for (V : DEAL_II_VEC_TEMPLATES)
   {
     template SUNDIALS::internal::NVectorView<PETScWrappers::MPI::V>
-    SUNDIALS::internal::make_nvector_view<>(PETScWrappers::MPI::V &);
+    SUNDIALS::internal::make_nvector_view<>(PETScWrappers::MPI::V &
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#  endif
+    );
     template SUNDIALS::internal::NVectorView<const PETScWrappers::MPI::V>
-    SUNDIALS::internal::make_nvector_view<>(const PETScWrappers::MPI::V &);
+    SUNDIALS::internal::make_nvector_view<>(const PETScWrappers::MPI::V &
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                            ,
+                                            SUNContext
+#  endif
+    );
     template PETScWrappers::MPI::V
       *SUNDIALS::internal::unwrap_nvector<PETScWrappers::MPI::V>(N_Vector);
     template const PETScWrappers::MPI::V *
index 78bc3db2fdaac8e51a7a1f9cd60dfc28537bf213..506be94b7e4546da8bf06a5132173f95f806d1ac 100644 (file)
@@ -176,10 +176,20 @@ namespace SUNDIALS
 
   template <typename VectorType>
   internal::LinearSolverWrapper<VectorType>::LinearSolverWrapper(
-    LinearSolveFunction<VectorType> lsolve)
+    LinearSolveFunction<VectorType> lsolve
+#    if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    ,
+    SUNContext &linsol_ctx
+#    endif
+    )
     : content(std::make_unique<LinearSolverContent<VectorType>>())
   {
-    sun_linear_solver                  = SUNLinSolNewEmpty();
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+    sun_linear_solver = SUNLinSolNewEmpty();
+#    else
+    sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
+#    endif
+
     sun_linear_solver->ops->gettype    = arkode_linsol_get_type;
     sun_linear_solver->ops->solve      = arkode_linsol_solve<VectorType>;
     sun_linear_solver->ops->setup      = arkode_linsol_setup<VectorType>;
@@ -230,11 +240,21 @@ namespace SUNDIALS
   SundialsOperator<VectorType>::vmult(VectorType &      dst,
                                       const VectorType &src) const
   {
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
     auto sun_dst = internal::make_nvector_view(dst);
     auto sun_src = internal::make_nvector_view(src);
     int  status  = a_times_fn(A_data, sun_src, sun_dst);
     (void)status;
     AssertSundialsSolver(status);
+#    else
+    // We don't currently know how to implement this: the
+    // internal::make_nvector_view() function called above requires a
+    // SUNContext object, but we don't actually have access to such an
+    // object here...
+    (void)dst;
+    (void)src;
+    Assert(false, ExcNotImplemented());
+#    endif
   }
 
 
@@ -263,6 +283,7 @@ namespace SUNDIALS
         return;
       }
 
+#    if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
     auto sun_dst = internal::make_nvector_view(dst);
     auto sun_src = internal::make_nvector_view(src);
     // for custom preconditioners no distinction between left and right
@@ -271,6 +292,15 @@ namespace SUNDIALS
       p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
     (void)status;
     AssertSundialsSolver(status);
+#    else
+    // We don't currently know how to implement this: the
+    // internal::make_nvector_view() function called above requires a
+    // SUNContext object, but we don't actually have access to such an
+    // object here...
+    (void)dst;
+    (void)src;
+    Assert(false, ExcNotImplemented());
+#    endif
   }
 
   template struct SundialsOperator<Vector<double>>;

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.