Revert "Moved to new interface in sundials."

This reverts commit b108fe75d06f08db71b33b031f2ea57daa7f4251.

diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h
index 5d006d969f4be1bee7ba442a2dc72106b0786ad0..43d52e16af8dacfeaa8bb643237285035252d7e9 100644 (file)
--- a/include/deal.II/sundials/ida.h
+++ b/include/deal.II/sundials/ida.h
@@ -883,11 +883,32 @@ namespace SUNDIALS
      */
     void *ida_mem;
 
+    /**
+     * IDA solution vector.
+     */
+    N_Vector yy;
+
+    /**
+     * IDA solution derivative vector.
+     */
+    N_Vector yp;
+
+    /**
+     * IDA absolute tolerances vector.
+     */
+    N_Vector abs_tolls;
+
+    /**
+     * IDA differential components vector.
+     */
+    N_Vector diff_id;
+
     /**
      * Number of iteration required to solve the Jacobian system
      */
     int n_iter;
 
+
     /**
      * MPI communicator. SUNDIALS solver runs happily in
      * parallel. Note that if the library is compiled without MPI

diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h
index 027fd017f1c52f02fe0d3b2a12b0917ac572c7fb..0d2acde788ff03bf31686d7127046483dc27835b 100644 (file)
--- a/include/deal.II/sundials/n_vector.templates.h
+++ b/include/deal.II/sundials/n_vector.templates.h
@@ -220,10 +220,6 @@ namespace SUNDIALS
       realtype
       weighted_rms_norm(N_Vector x, N_Vector w);
 
-      template <typename VectorType>
-      realtype
-      weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask);
-
       template <typename VectorType>
       realtype
       max_norm(N_Vector x);
@@ -669,24 +665,6 @@ SUNDIALS::internal::NVectorOperations::weighted_rms_norm(N_Vector x, N_Vector w)
 
 
 
-template <typename VectorType>
-realtype
-SUNDIALS::internal::NVectorOperations::weighted_rms_norm_mask(N_Vector x,
-                                                              N_Vector w,
-                                                              N_Vector mask)
-{
-  // TODO copy can be avoided by a custom kernel
-  VectorType tmp         = *unwrap_nvector_const<VectorType>(x);
-  auto *     w_dealii    = unwrap_nvector_const<VectorType>(w);
-  auto *     mask_dealii = unwrap_nvector_const<VectorType>(mask);
-  const auto n           = tmp.size();
-  tmp.scale(*w_dealii);
-  tmp.scale(*mask_dealii);
-  return tmp.l2_norm() / std::sqrt(n);
-}
-
-
-
 template <typename VectorType>
 realtype
 SUNDIALS::internal::NVectorOperations::max_norm(N_Vector x)
@@ -962,11 +940,11 @@ SUNDIALS::internal::create_empty_nvector()
   v->ops->nvdotprod   = NVectorOperations::dot_product<VectorType>;
   v->ops->nvmaxnorm   = NVectorOperations::max_norm<VectorType>;
   v->ops->nvwrmsnorm  = NVectorOperations::weighted_rms_norm<VectorType>;
+  //  v->ops->nvwrmsnormmask = undef;
   v->ops->nvmin     = NVectorOperations::min_element<VectorType>;
   v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm<VectorType>;
   v->ops->nvl1norm  = NVectorOperations::l1_norm<VectorType>;
-  v->ops->nvwrmsnormmask =
-    NVectorOperations::weighted_rms_norm_mask<VectorType>;
+  ;
   //  v->ops->nvcompare      = undef;
   //  v->ops->nvinvtest      = undef;
   //  v->ops->nvconstrmask   = undef;

diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc
index bd207108cfa8e635d9844a6ea9885e3fc7378bcb..6788c62fca4b3d9cbec7f1f9e8481625d2895e4d 100644 (file)
--- a/source/sundials/ida.cc
+++ b/source/sundials/ida.cc
@@ -14,17 +14,12 @@
 //-----------------------------------------------------------
 #include <deal.II/base/config.h>
 
-#include <deal.II/lac/vector_operation.h>
-
 #include <deal.II/sundials/ida.h>
 
 #ifdef DEAL_II_WITH_SUNDIALS
 #  include <deal.II/base/utilities.h>
 
 #  include <deal.II/lac/block_vector.h>
-
-#  include <deal.II/sundials/n_vector.h>
-#  include <deal.II/sundials/sunlinsol_wrapper.h>
 #  ifdef DEAL_II_WITH_TRILINOS
 #    include <deal.II/lac/trilinos_parallel_block_vector.h>
 #    include <deal.II/lac/trilinos_vector.h>
@@ -34,6 +29,7 @@
 #    include <deal.II/lac/petsc_vector.h>
 #  endif
 
+#  include <deal.II/sundials/copy.h>
 #  include <deal.II/sundials/n_vector.h>
 
 #  include <sundials/sundials_config.h>
@@ -67,13 +63,24 @@ namespace SUNDIALS
                    void *   user_data)
     {
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
 
-      auto *src_yy   = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *src_yp   = internal::unwrap_nvector_const<VectorType>(yp);
-      auto *residual = internal::unwrap_nvector<VectorType>(rr);
+      typename VectorMemory<VectorType>::Pointer src_yp(mem);
+      solver.reinit_vector(*src_yp);
+
+      typename VectorMemory<VectorType>::Pointer residual(mem);
+      solver.reinit_vector(*residual);
+
+      copy(*src_yy, yy);
+      copy(*src_yp, yp);
 
       int err = solver.residual(tt, *src_yy, *src_yp, *residual);
 
+      copy(rr, *residual);
+
       return err;
     }
 
@@ -96,9 +103,16 @@ namespace SUNDIALS
       (void)resp;
       IDA<VectorType> &solver =
         *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
+      GrowingVectorMemory<VectorType> mem;
 
-      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer src_yp(mem);
+      solver.reinit_vector(*src_yp);
+
+      copy(*src_yy, yy);
+      copy(*src_yp, yp);
 
       int err = solver.setup_jacobian(IDA_mem->ida_tn,
                                       *src_yy,
@@ -127,13 +141,16 @@ namespace SUNDIALS
         *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
       GrowingVectorMemory<VectorType> mem;
 
+      typename VectorMemory<VectorType>::Pointer src(mem);
+      solver.reinit_vector(*src);
+
       typename VectorMemory<VectorType>::Pointer dst(mem);
       solver.reinit_vector(*dst);
 
-      auto *src = internal::unwrap_nvector<VectorType>(b);
+      copy(*src, b);
 
       int err = solver.solve_jacobian_system(*src, *dst);
-      *src    = *dst;
+      copy(b, *dst);
 
       return err;
     }
@@ -156,12 +173,21 @@ namespace SUNDIALS
     {
       Assert(user_data != nullptr, ExcInternalError());
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer src_yp(mem);
+      solver.reinit_vector(*src_yp);
+
+      copy(*src_yy, yy);
+      copy(*src_yp, yp);
 
-      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
 
       int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj);
 
+
       return err;
     }
 
@@ -178,11 +204,21 @@ namespace SUNDIALS
       const IDA<VectorType> &solver =
         *static_cast<const IDA<VectorType> *>(LS->content);
 
-      auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
-      auto *dst_x = internal::unwrap_nvector<VectorType>(x);
+      // Allocate temporary (deal.II-type) vectors into which to copy the
+      // N_vectors
+      GrowingVectorMemory<VectorType>            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);
 
       const int err = solver.solve_jacobian_system(*src_b, *dst_x);
 
+      copy(x, *dst_x);
+
       return err;
     }
 
@@ -197,12 +233,20 @@ namespace SUNDIALS
     {
       IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
 
-      auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
-      auto *dst_x = internal::unwrap_nvector<VectorType>(x);
+      // Allocate temporary (deal.II-type) vectors into which to copy the
+      // N_vectors
+      GrowingVectorMemory<VectorType>            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);
       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);
+      copy(x, *dst_x);
 
       return err;
     }
@@ -215,6 +259,10 @@ namespace SUNDIALS
   IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm mpi_comm)
     : data(data)
     , ida_mem(nullptr)
+    , yy(nullptr)
+    , yp(nullptr)
+    , abs_tolls(nullptr)
+    , diff_id(nullptr)
     , communicator(is_serial_vector<VectorType>::value ?
                      MPI_COMM_SELF :
                      Utilities::MPI::duplicate_communicator(mpi_comm))
@@ -243,6 +291,8 @@ namespace SUNDIALS
   unsigned int
   IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
   {
+    unsigned int system_size = solution.size();
+
     double       t           = data.initial_time;
     double       h           = data.initial_step_size;
     unsigned int step_number = 0;
@@ -251,11 +301,36 @@ namespace SUNDIALS
     int status;
     (void)status;
 
-    reset(data.initial_time, data.initial_step_size, solution, solution_dot);
-
     // The solution is stored in
     // solution. Here we take only a
     // view of it.
+#  ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        const IndexSet    is                = solution.locally_owned_elements();
+        const std::size_t local_system_size = is.n_elements();
+
+        yy = N_VNew_Parallel(communicator, local_system_size, system_size);
+
+        yp = N_VNew_Parallel(communicator, local_system_size, system_size);
+
+        diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
+
+        abs_tolls =
+          N_VNew_Parallel(communicator, local_system_size, system_size);
+      }
+    else
+#  endif
+      {
+        Assert(is_serial_vector<VectorType>::value,
+               ExcInternalError(
+                 "Trying to use a serial code with a parallel vector."));
+        yy        = N_VNew_Serial(system_size);
+        yp        = N_VNew_Serial(system_size);
+        diff_id   = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+    reset(data.initial_time, data.initial_step_size, solution, solution_dot);
 
     double next_time = data.initial_time;
 
@@ -265,15 +340,15 @@ namespace SUNDIALS
       {
         next_time += data.output_period;
 
-        auto yy = internal::make_nvector_view(solution);
-        auto yp = internal::make_nvector_view(solution_dot);
-
         status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
         AssertIDA(status);
 
         status = IDAGetLastStep(ida_mem, &h);
         AssertIDA(status);
 
+        copy(solution, yy);
+        copy(solution_dot, yp);
+
         while (solver_should_restart(t, solution, solution_dot))
           reset(t, h, solution, solution_dot);
 
@@ -282,6 +357,24 @@ namespace SUNDIALS
         output_step(t, solution, solution_dot, step_number);
       }
 
+      // Free the vectors which are no longer used.
+#  ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        N_VDestroy_Parallel(yy);
+        N_VDestroy_Parallel(yp);
+        N_VDestroy_Parallel(abs_tolls);
+        N_VDestroy_Parallel(diff_id);
+      }
+    else
+#  endif
+      {
+        N_VDestroy_Serial(yy);
+        N_VDestroy_Serial(yp);
+        N_VDestroy_Serial(abs_tolls);
+        N_VDestroy_Serial(diff_id);
+      }
+
     return step_number;
   }
 
@@ -294,26 +387,71 @@ namespace SUNDIALS
                          VectorType & solution,
                          VectorType & solution_dot)
   {
-    bool first_step = (current_time == data.initial_time);
+    unsigned int system_size;
+    bool         first_step = (current_time == data.initial_time);
 
     if (ida_mem)
       IDAFree(&ida_mem);
 
     ida_mem = IDACreate();
 
+    // Free the vectors which are no longer used.
+    if (yy)
+      {
+#  ifdef DEAL_II_WITH_MPI
+        if (is_serial_vector<VectorType>::value == false)
+          {
+            N_VDestroy_Parallel(yy);
+            N_VDestroy_Parallel(yp);
+            N_VDestroy_Parallel(abs_tolls);
+            N_VDestroy_Parallel(diff_id);
+          }
+        else
+#  endif
+          {
+            N_VDestroy_Serial(yy);
+            N_VDestroy_Serial(yp);
+            N_VDestroy_Serial(abs_tolls);
+            N_VDestroy_Serial(diff_id);
+          }
+      }
+
     int status;
     (void)status;
+    system_size = solution.size();
+#  ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        const IndexSet    is                = solution.locally_owned_elements();
+        const std::size_t local_system_size = is.n_elements();
+
+        yy = N_VNew_Parallel(communicator, local_system_size, system_size);
 
-    auto yy = internal::make_nvector_view(solution);
-    auto yp = internal::make_nvector_view(solution_dot);
+        yp = N_VNew_Parallel(communicator, local_system_size, system_size);
+
+        diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
+
+        abs_tolls =
+          N_VNew_Parallel(communicator, local_system_size, system_size);
+      }
+    else
+#  endif
+      {
+        yy        = N_VNew_Serial(system_size);
+        yp        = N_VNew_Serial(system_size);
+        diff_id   = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+
+    copy(yy, solution);
+    copy(yp, solution_dot);
 
     status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
     AssertIDA(status);
     if (get_local_tolerances)
       {
-        const auto abs_tols =
-          internal::make_nvector_view(get_local_tolerances());
-        status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
+        copy(abs_tolls, get_local_tolerances());
+        status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls);
         AssertIDA(status);
       }
     else
@@ -339,10 +477,9 @@ namespace SUNDIALS
         auto dc          = differential_components();
         for (auto i = dc.begin(); i != dc.end(); ++i)
           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);
+        copy(diff_id, diff_comp_vector);
+        status = IDASetId(ida_mem, diff_id);
         AssertIDA(status);
       }
 
@@ -488,6 +625,9 @@ namespace SUNDIALS
 
         status = IDAGetConsistentIC(ida_mem, yy, yp);
         AssertIDA(status);
+
+        copy(solution, yy);
+        copy(solution_dot, yp);
       }
     else if (type == AdditionalData::use_y_diff)
       {
@@ -497,6 +637,9 @@ namespace SUNDIALS
 
         status = IDAGetConsistentIC(ida_mem, yy, yp);
         AssertIDA(status);
+
+        copy(solution, yy);
+        copy(solution_dot, yp);
       }
   }