Revert SUNDIALS copy to fix KINSOL, part 1

diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h
index 5b2c50cd194c25df9252a90d5d69e284e28acc3b..9bcdfda83c1b379e6eaceea0b86d5d65b4e95555 100644 (file)
--- a/include/deal.II/sundials/kinsol.h
+++ b/include/deal.II/sundials/kinsol.h
@@ -691,6 +691,21 @@ 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;
+
 #  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
     /**
      * A context object associated with the KINSOL solver.

diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc
index 66b40dfc846dda4c1f111a0e11224a4b9397dac4..921755853f1b98f5426c5012705cab67eea3b6bf 100644 (file)
--- a/source/sundials/kinsol.cc
+++ b/source/sundials/kinsol.cc
@@ -34,6 +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...
@@ -134,9 +135,15 @@ namespace SUNDIALS
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
 
-      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *dst_FF = internal::unwrap_nvector<VectorType>(FF);
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer dst_FF(mem);
+      solver.reinit_vector(*dst_FF);
+
+      internal::copy(*src_yy, yy);
 
       int err = 0;
       if (solver.residual)
@@ -144,6 +151,8 @@ namespace SUNDIALS
       else
         Assert(false, ExcInternalError());
 
+      internal::copy(FF, *dst_FF);
+
       return err;
     }
 
@@ -155,9 +164,15 @@ namespace SUNDIALS
     {
       KINSOL<VectorType> &solver =
         *static_cast<KINSOL<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
 
-      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *dst_FF = internal::unwrap_nvector<VectorType>(FF);
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer dst_FF(mem);
+      solver.reinit_vector(*dst_FF);
+
+      internal::copy(*src_yy, yy);
 
       int err = 0;
       if (solver.iteration_function)
@@ -165,6 +180,8 @@ namespace SUNDIALS
       else
         Assert(false, ExcInternalError());
 
+      internal::copy(FF, *dst_FF);
+
       return err;
     }
 
@@ -185,8 +202,16 @@ namespace SUNDIALS
       const KINSOL<VectorType> &solver =
         *static_cast<const KINSOL<VectorType> *>(user_data);
 
-      auto *ycur = internal::unwrap_nvector_const<VectorType>(u);
-      auto *fcur = internal::unwrap_nvector_const<VectorType>(f);
+      // 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);
+
+      internal::copy(*ycur, u);
+      internal::copy(*fcur, f);
 
       // Call the user-provided setup function with these arguments:
       solver.setup_jacobian(*ycur, *fcur);
@@ -214,11 +239,21 @@ namespace SUNDIALS
       // signals to call. Let's first check the more modern one:
       if (solver.solve_with_jacobian)
         {
-          auto *src_b = internal::unwrap_nvector<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);
+
+          internal::copy(*src_b, b);
 
           const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
 
+          internal::copy(x, *dst_x);
+
           return err;
         }
       else
@@ -233,9 +268,13 @@ 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);
 
-          auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
-          auto *dst_x = internal::unwrap_nvector<VectorType>(x);
+          internal::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
@@ -245,6 +284,8 @@ namespace SUNDIALS
           const int err =
             solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
 
+          internal::copy(x, *dst_x);
+
           return err;
         }
     }
@@ -260,6 +301,9 @@ namespace SUNDIALS
                          MPI_COMM_SELF :
                          Utilities::MPI::duplicate_communicator(mpi_comm))
     , kinsol_mem(nullptr)
+    , solution(nullptr)
+    , u_scale(nullptr)
+    , f_scale(nullptr)
   {
     set_functions_to_trigger_an_assert();
 
@@ -299,48 +343,6 @@ namespace SUNDIALS
   unsigned int
   KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
   {
-    internal::NVectorView<VectorType> u_scale, f_scale;
-
-    VectorType u_scale_temp, f_scale_temp;
-
-    if (get_solution_scaling)
-      u_scale = internal::make_nvector_view(get_solution_scaling()
-#  if DEAL_II_SUNDIALS_VERSION_GTE(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
-#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
-                                              ,
-                                              kinsol_ctx
-#  endif
-        );
-      }
-
-    if (get_function_scaling)
-      f_scale = internal::make_nvector_view(get_function_scaling()
-#  if DEAL_II_SUNDIALS_VERSION_GTE(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
-#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
-                                              ,
-                                              kinsol_ctx
-#  endif
-        );
-      }
-
     // Make sure we have what we need
     if (data.strategy == AdditionalData::fixed_point)
       {
@@ -355,13 +357,7 @@ namespace SUNDIALS
                  "solve_jacobian_system || solve_with_jacobian"));
       }
 
-    auto solution = internal::make_nvector_view(initial_guess_and_solution
-#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
-                                                ,
-                                                kinsol_ctx
-#  endif
-    );
-
+    // Create a new solver object:
     int status = 0;
     (void)status;
 
@@ -387,6 +383,79 @@ namespace SUNDIALS
     status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
     AssertKINSOL(status);
 
+    const auto 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 auto local_system_size = is.n_elements();
+
+        solution = N_VNew_Parallel(mpi_communicator,
+                                   local_system_size,
+                                   system_size
+#    if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                   ,
+                                   kinsol_ctx
+#    endif
+        );
+
+        u_scale = N_VNew_Parallel(mpi_communicator,
+                                  local_system_size,
+                                  system_size
+#    if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                  ,
+                                  kinsol_ctx
+#    endif
+        );
+        N_VConst_Parallel(1.e0, u_scale);
+
+        f_scale = N_VNew_Parallel(mpi_communicator,
+                                  local_system_size,
+                                  system_size
+#    if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                  ,
+                                  kinsol_ctx
+#    endif
+        );
+        N_VConst_Parallel(1.e0, f_scale);
+      }
+    else
+#  endif
+      {
+        Assert(is_serial_vector<VectorType>::value, ExcInternalError());
+        solution = N_VNew_Serial(system_size
+#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                 ,
+                                 kinsol_ctx
+#  endif
+        );
+        u_scale = N_VNew_Serial(system_size
+#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                ,
+                                kinsol_ctx
+#  endif
+        );
+        N_VConst_Serial(1.e0, u_scale);
+        f_scale = N_VNew_Serial(system_size
+#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                ,
+                                kinsol_ctx
+#  endif
+        );
+        N_VConst_Serial(1.e0, f_scale);
+      }
+
+    if (get_solution_scaling)
+      internal::copy(u_scale, get_solution_scaling());
+
+    if (get_function_scaling)
+      internal::copy(f_scale, get_function_scaling());
+
+    internal::copy(solution, initial_guess_and_solution);
+
     // This must be called before KINSetMAA
     status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
     AssertKINSOL(status);
@@ -512,6 +581,24 @@ namespace SUNDIALS
     status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
     AssertKINSOL(status);
 
+    internal::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);