]> https://gitweb.dealii.org/ - dealii.git/commitdiff
KINSOL: do not copy vectors internally 15194/head
authorSebastian Proell <sebastian.proell@tum.de>
Thu, 11 May 2023 14:19:48 +0000 (16:19 +0200)
committerSebastian Proell <sebastian.proell@tum.de>
Sun, 14 May 2023 14:04:43 +0000 (16:04 +0200)
include/deal.II/sundials/kinsol.h
source/sundials/kinsol.cc

index 994f2b225ad7ebbf26969045460fd6bb54ce5bcc..56ce7b6994f71b25604c22f1bbbff11fa2b66e3c 100644 (file)
@@ -739,20 +739,6 @@ namespace SUNDIALS
     SUNContext kinsol_ctx;
 #  endif
 
-    /**
-     * KINSOL solution vector.
-     */
-    N_Vector solution;
-
-    /**
-     * KINSOL solution scale.
-     */
-    N_Vector u_scale;
-
-    /**
-     * KINSOL f scale.
-     */
-    N_Vector f_scale;
 
     /**
      * Memory pool of vectors.
index 1a89c03e6eb633ca7f901be24de6e6cd80b488b4..a2bd709493b999e9ee201c052188175cfdaf1084 100644 (file)
@@ -36,7 +36,6 @@
 #    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...
@@ -217,9 +216,6 @@ namespace SUNDIALS
 #  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
     , kinsol_ctx(nullptr)
 #  endif
-    , solution(nullptr)
-    , u_scale(nullptr)
-    , f_scale(nullptr)
     , pending_exception(nullptr)
   {
     set_functions_to_trigger_an_assert();
@@ -300,78 +296,32 @@ 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
+    // helper function to create N_Vectors compatible with different versions
+    // of SUNDIALS
+    const auto make_compatible_nvector_view = [this](auto &v) {
+      return internal::make_nvector_view(v
+#  if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+                                         ,
+                                         kinsol_ctx
 #  endif
+      );
+    };
+
+
+    VectorType ones;
+    // Prepare constant vector for scaling f or u only if we need it
+    if (!get_function_scaling || !get_solution_scaling)
       {
-        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);
+        reinit_vector(ones);
+        ones = 1.0;
       }
 
-    if (get_solution_scaling)
-      internal::copy(u_scale, get_solution_scaling());
-
-    if (get_function_scaling)
-      internal::copy(f_scale, get_function_scaling());
+    auto u_scale = make_compatible_nvector_view(
+      get_solution_scaling ? get_solution_scaling() : ones);
+    auto f_scale = make_compatible_nvector_view(
+      get_function_scaling ? get_function_scaling() : ones);
 
-    internal::copy(solution, initial_guess_and_solution);
+    auto solution = make_compatible_nvector_view(initial_guess_and_solution);
 
     // This must be called before KINSetMAA
     status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
@@ -388,15 +338,9 @@ namespace SUNDIALS
         [](N_Vector yy, N_Vector FF, void *user_data) -> int {
           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);
 
-          internal::copy(*src_yy, yy);
+          auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+          auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
 
           Assert(solver.iteration_function, ExcInternalError());
 
@@ -406,8 +350,6 @@ namespace SUNDIALS
                                                 *src_yy,
                                                 *dst_FF);
 
-          internal::copy(FF, *dst_FF);
-
           return err;
         },
         solution);
@@ -418,23 +360,15 @@ namespace SUNDIALS
         [](N_Vector yy, N_Vector FF, void *user_data) -> int {
           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);
 
-          internal::copy(*src_yy, yy);
+          auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+          auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
 
           Assert(solver.residual, ExcInternalError());
 
           const int err = call_and_possibly_capture_exception(
             solver.residual, solver.pending_exception, *src_yy, *dst_FF);
 
-          internal::copy(FF, *dst_FF);
-
           return err;
         },
         solution);
@@ -465,8 +399,8 @@ namespace SUNDIALS
     SUNLinearSolver LS = nullptr;
 
     if (solve_jacobian_system ||
-        solve_with_jacobian) // user assigned a function object to the solver
-                             // slot
+        solve_with_jacobian) // user assigned a function
+                             // object to the solver slot
       {
         // Set the operations we care for in the sun_linear_solver object
         // and attach it to the KINSOL object. The functions that will get
@@ -515,16 +449,8 @@ 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);
-
-              internal::copy(*src_b, b);
+              auto src_b = internal::unwrap_nvector_const<VectorType>(b);
+              auto dst_x = internal::unwrap_nvector<VectorType>(x);
 
               const int err =
                 call_and_possibly_capture_exception(solver.solve_with_jacobian,
@@ -533,8 +459,6 @@ namespace SUNDIALS
                                                     *dst_x,
                                                     tol);
 
-              internal::copy(x, *dst_x);
-
               return err;
             }
           else
@@ -544,18 +468,13 @@ namespace SUNDIALS
               // the old signal. Check this.
               Assert(solver.solve_jacobian_system, ExcInternalError());
 
-              // Allocate temporary (deal.II-type) vectors into which to copy
-              // the N_vectors
+              // Allocate temporary (deal.II-type) dummy 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);
-
-              internal::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
@@ -570,8 +489,6 @@ namespace SUNDIALS
                 *src_b,
                 *dst_x);
 
-              internal::copy(x, *dst_x);
-
               return err;
             }
         };
@@ -631,16 +548,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);
-
-            internal::copy(*ycur, u);
-            internal::copy(*fcur, f);
+            auto ycur = internal::unwrap_nvector_const<VectorType>(u);
+            auto fcur = internal::unwrap_nvector<VectorType>(f);
 
             // Call the user-provided setup function with these arguments:
             return call_and_possibly_capture_exception(solver.setup_jacobian,
@@ -669,22 +578,6 @@ namespace SUNDIALS
     status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
 
     ScopeExit upon_exit([this, &J, &LS]() mutable {
-    // 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);
-        }
-
       if (J != nullptr)
         SUNMatDestroy(J);
       if (LS != nullptr)
@@ -714,8 +607,6 @@ namespace SUNDIALS
       }
     AssertKINSOL(status);
 
-    internal::copy(initial_guess_and_solution, solution);
-
     long nniters;
     status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
     AssertKINSOL(status);

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.