Remove copy from SUNDIALS::KINSOL

diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h
index 8a33c0ea2741478dd8b21b39132708e4e9be05a5..939215b05f3db1f4f8d6498c24da64a4bc37a542 100644 (file)
--- a/include/deal.II/sundials/kinsol.h
+++ b/include/deal.II/sundials/kinsol.h
@@ -711,21 +711,6 @@ namespace SUNDIALS
      */
     void *kinsol_mem;
 
-    /**
-     * KINSOL solution vector.
-     */
-    N_Vector solution;
-
-    /**
-     * KINSOL solution scale.
-     */
-    N_Vector u_scale;
-
-    /**
-     * KINSOL f scale.
-     */
-    N_Vector f_scale;
-
     /**
      * MPI communicator. SUNDIALS solver runs happily in parallel.
      */

diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h
index a6d6562ff69f0cce8d1cd21d60ad7b642fe137fd..49ea702018c8ca62e08a841a789643b964ab760c 100644 (file)
--- a/include/deal.II/sundials/n_vector.templates.h
+++ b/include/deal.II/sundials/n_vector.templates.h
@@ -168,6 +168,14 @@ namespace SUNDIALS
       realtype
       dot_product(N_Vector x, N_Vector y);
 
+      template <typename VectorType>
+      realtype
+      weighted_l2_norm(N_Vector x, N_Vector y);
+
+      template <typename VectorType>
+      realtype
+      l1_norm(N_Vector x);
+
       template <typename VectorType>
       void
       elementwise_product(N_Vector x, N_Vector y, N_Vector z);
@@ -569,6 +577,28 @@ SUNDIALS::internal::NVectorOperations::dot_product(N_Vector x, N_Vector y)
 
 
 
+template <typename VectorType>
+realtype
+SUNDIALS::internal::NVectorOperations::weighted_l2_norm(N_Vector x, N_Vector w)
+{
+  // TODO copy can be avoided by a custom kernel
+  VectorType tmp      = *unwrap_nvector_const<VectorType>(x);
+  auto *     w_dealii = unwrap_nvector_const<VectorType>(w);
+  tmp.scale(*w_dealii);
+  return tmp.l2_norm();
+}
+
+
+
+template <typename VectorType>
+realtype
+SUNDIALS::internal::NVectorOperations::l1_norm(N_Vector x)
+{
+  return unwrap_nvector_const<VectorType>(x)->l1_norm();
+}
+
+
+
 template <typename VectorType>
 void
 SUNDIALS::internal::NVectorOperations::set_constant(realtype c, N_Vector v)
@@ -911,9 +941,10 @@ SUNDIALS::internal::create_empty_nvector()
   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      = undef;
-  //  v->ops->nvl1norm       = 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->nvcompare      = undef;
   //  v->ops->nvinvtest      = undef;
   //  v->ops->nvconstrmask   = undef;

diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc
index ce3cb440778308ea1b1b910b86a0c804300f7069..cb12781b0eba0984fb457eea491e5bc53964129b 100644 (file)
--- a/source/sundials/kinsol.cc
+++ b/source/sundials/kinsol.cc
@@ -34,7 +34,7 @@
 #    include <deal.II/lac/petsc_vector.h>
 #  endif
 
-#  include <deal.II/sundials/copy.h>
+#  include <deal.II/sundials/n_vector.h>
 
 // Make sure we #include the SUNDIALS config file...
 #  include <sundials/sundials_config.h>
@@ -137,15 +137,9 @@ namespace SUNDIALS
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType> mem;
 
-      typename VectorMemory<VectorType>::Pointer src_yy(mem);
-      solver.reinit_vector(*src_yy);
-
-      typename VectorMemory<VectorType>::Pointer dst_FF(mem);
-      solver.reinit_vector(*dst_FF);
-
-      copy(*src_yy, yy);
+      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+      auto *dst_FF = internal::unwrap_nvector<VectorType>(FF);
 
       int err = 0;
       if (solver.residual)
@@ -155,8 +149,6 @@ namespace SUNDIALS
       else
         Assert(false, ExcInternalError());
 
-      copy(FF, *dst_FF);
-
       return err;
     }
 
@@ -169,16 +161,11 @@ namespace SUNDIALS
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-      GrowingVectorMemory<VectorType> mem;
 
-      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
-      solver.reinit_vector(*src_ycur);
-
-      typename VectorMemory<VectorType>::Pointer src_fcur(mem);
-      solver.reinit_vector(*src_fcur);
-
-      copy(*src_ycur, kinsol_mem->kin_uu);
-      copy(*src_fcur, kinsol_mem->kin_fval);
+      auto *src_ycur =
+        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
+      auto *src_fcur =
+        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
 
       int err = solver.setup_jacobian(*src_ycur, *src_fcur);
       return err;
@@ -196,27 +183,15 @@ namespace SUNDIALS
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
-      solver.reinit_vector(*src_ycur);
 
-      typename VectorMemory<VectorType>::Pointer src_fcur(mem);
-      solver.reinit_vector(*src_fcur);
-
-      copy(*src_ycur, kinsol_mem->kin_uu);
-      copy(*src_fcur, kinsol_mem->kin_fval);
-
-      typename VectorMemory<VectorType>::Pointer src(mem);
-      solver.reinit_vector(*src);
-
-      typename VectorMemory<VectorType>::Pointer dst(mem);
-      solver.reinit_vector(*dst);
-
-      copy(*src, b);
+      auto *src_ycur =
+        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
+      auto *src_fcur =
+        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
+      auto *src = internal::unwrap_nvector_const<VectorType>(b);
+      auto *dst = internal::unwrap_nvector<VectorType>(x);
 
       int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst);
-      copy(x, *dst);
 
       *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale);
       N_VProd(b, kinsol_mem->kin_fscale, b);
@@ -243,16 +218,8 @@ namespace SUNDIALS
       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);
+      auto *ycur = internal::unwrap_nvector_const<VectorType>(u);
+      auto *fcur = internal::unwrap_nvector_const<VectorType>(f);
 
       // Call the user-provided setup function with these arguments:
       solver.setup_jacobian(*ycur, *fcur);
@@ -280,21 +247,11 @@ namespace SUNDIALS
       // 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<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);
+          auto *src_b = internal::unwrap_nvector<VectorType>(b);
+          auto *dst_x = internal::unwrap_nvector<VectorType>(x);
 
           const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
 
-          copy(x, *dst_x);
-
           return err;
         }
       else
@@ -309,13 +266,9 @@ namespace SUNDIALS
           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);
+          auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
+          auto *dst_x = internal::unwrap_nvector<VectorType>(x);
 
           // Call the user-provided setup function with these arguments. Note
           // that Sundials 4.x and later no longer provide values for
@@ -325,8 +278,6 @@ namespace SUNDIALS
           const int err =
             solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
 
-          copy(x, *dst_x);
-
           return err;
         }
     }
@@ -341,9 +292,6 @@ namespace SUNDIALS
                              const MPI_Comm &      mpi_comm)
     : data(data)
     , kinsol_mem(nullptr)
-    , solution(nullptr)
-    , u_scale(nullptr)
-    , f_scale(nullptr)
     , communicator(is_serial_vector<VectorType>::value ?
                      MPI_COMM_SELF :
                      Utilities::MPI::duplicate_communicator(mpi_comm))
@@ -377,44 +325,29 @@ namespace SUNDIALS
   {
     unsigned int system_size = initial_guess_and_solution.size();
 
-    // 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 = initial_guess_and_solution.locally_owned_elements();
-        const unsigned int local_system_size = is.n_elements();
-
-        solution =
-          N_VNew_Parallel(communicator, local_system_size, system_size);
+    NVectorView<VectorType> u_scale, f_scale;
 
-        u_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
-        N_VConst_Parallel(1.e0, u_scale);
+    VectorType u_scale_temp, f_scale_temp;
 
-        f_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
-        N_VConst_Parallel(1.e0, f_scale);
-      }
+    if (get_solution_scaling)
+      u_scale = internal::make_nvector_view(get_solution_scaling());
     else
-#  endif
       {
-        Assert(is_serial_vector<VectorType>::value,
-               ExcInternalError(
-                 "Trying to use a serial code with a parallel vector."));
-        solution = N_VNew_Serial(system_size);
-        u_scale  = N_VNew_Serial(system_size);
-        N_VConst_Serial(1.e0, u_scale);
-        f_scale = N_VNew_Serial(system_size);
-        N_VConst_Serial(1.e0, f_scale);
+        reinit_vector(u_scale_temp);
+        u_scale_temp = 1.0;
+        u_scale      = internal::make_nvector_view(u_scale_temp);
       }
 
-    if (get_solution_scaling)
-      copy(u_scale, get_solution_scaling());
-
     if (get_function_scaling)
-      copy(f_scale, get_function_scaling());
+      f_scale = internal::make_nvector_view(get_function_scaling());
+    else
+      {
+        reinit_vector(f_scale_temp);
+        f_scale_temp = 1.0;
+        f_scale      = internal::make_nvector_view(f_scale_temp);
+      }
 
-    copy(solution, initial_guess_and_solution);
+    auto solution = internal::make_nvector_view(initial_guess_and_solution);
 
     if (kinsol_mem)
       KINFree(&kinsol_mem);
@@ -573,24 +506,6 @@ namespace SUNDIALS
     status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
     AssertKINSOL(status);
 
-    copy(initial_guess_and_solution, solution);
-
-    // Free the vectors which are no longer used.
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        N_VDestroy_Parallel(solution);
-        N_VDestroy_Parallel(u_scale);
-        N_VDestroy_Parallel(f_scale);
-      }
-    else
-#  endif
-      {
-        N_VDestroy_Serial(solution);
-        N_VDestroy_Serial(u_scale);
-        N_VDestroy_Serial(f_scale);
-      }
-
     long nniters;
     status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
     AssertKINSOL(status);