From 939cc98f48a0f85655753d83ab3e3e8d1ef3b38e Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Mon, 10 May 2021 10:46:53 +0200
Subject: [PATCH] Enable SUNDIALS::KINSOL for
 LinearAlgebra::distributed::Vector

---
 include/deal.II/sundials/copy.h | 14 +++++++
 source/sundials/copy.cc         | 71 +++++++++++++++++++++++++++++++++
 source/sundials/kinsol.cc       |  5 +++
 3 files changed, 90 insertions(+)

diff --git a/include/deal.II/sundials/copy.h b/include/deal.II/sundials/copy.h
index 4a14ae9a5b..afc2ae3b52 100644
--- a/include/deal.II/sundials/copy.h
+++ b/include/deal.II/sundials/copy.h
@@ -24,6 +24,8 @@
 #    include <nvector/nvector_parallel.h>
 #  endif
 #  include <deal.II/lac/block_vector.h>
+#  include <deal.II/lac/la_parallel_block_vector.h>
+#  include <deal.II/lac/la_parallel_vector.h>
 #  include <deal.II/lac/vector.h>
 
 #  include <nvector/nvector_serial.h>
@@ -82,6 +84,18 @@ namespace SUNDIALS
     copy(Vector<double> &dst, const N_Vector &src);
     void
     copy(N_Vector &dst, const Vector<double> &src);
+
+    void
+    copy(LinearAlgebra::distributed::BlockVector<double> &dst,
+         const N_Vector &                                 src);
+    void
+    copy(N_Vector &                                             dst,
+         const LinearAlgebra::distributed::BlockVector<double> &src);
+
+    void
+    copy(LinearAlgebra::distributed::Vector<double> &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const LinearAlgebra::distributed::Vector<double> &src);
   } // namespace internal
 } // namespace SUNDIALS
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/sundials/copy.cc b/source/sundials/copy.cc
index 1c18eef17f..b11b279cc0 100644
--- a/source/sundials/copy.cc
+++ b/source/sundials/copy.cc
@@ -56,6 +56,77 @@ namespace SUNDIALS
       return static_cast<std::size_t>(length);
     }
 
+
+    void
+    copy(LinearAlgebra::distributed::Vector<double> &dst, const N_Vector &src)
+    {
+      const IndexSet    is = dst.locally_owned_elements();
+      const std::size_t N  = is.n_elements();
+      AssertDimension(N, N_Vector_length(src));
+      dst.zero_out_ghost_values();
+      for (std::size_t i = 0; i < N; ++i)
+        {
+#  ifdef DEAL_II_WITH_MPI
+          dst.local_element(i) = NV_Ith_P(src, i);
+#  else
+          dst.local_element(i)        = NV_Ith_S(src, i);
+#  endif
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &dst, const LinearAlgebra::distributed::Vector<double> &src)
+    {
+      const IndexSet    is = src.locally_owned_elements();
+      const std::size_t N  = is.n_elements();
+      AssertDimension(N, N_Vector_length(dst));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+#  ifdef DEAL_II_WITH_MPI
+          NV_Ith_P(dst, i) = src.local_element(i);
+#  else
+          NV_Ith_S(dst, i)            = src.local_element(i);
+#  endif
+        }
+    }
+
+    void
+    copy(LinearAlgebra::distributed::BlockVector<double> &dst,
+         const N_Vector &                                 src)
+    {
+      const IndexSet    is = dst.locally_owned_elements();
+      const std::size_t N  = is.n_elements();
+      AssertDimension(N, N_Vector_length(src));
+      dst.zero_out_ghost_values();
+      for (std::size_t i = 0; i < N; ++i)
+        {
+#  ifdef DEAL_II_WITH_MPI
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+#  else
+          dst[is.nth_index_in_set(i)] = NV_Ith_S(src, i);
+#  endif
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &                                             dst,
+         const LinearAlgebra::distributed::BlockVector<double> &src)
+    {
+      IndexSet          is = src.locally_owned_elements();
+      const std::size_t N  = is.n_elements();
+      AssertDimension(N, N_Vector_length(dst));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+#  ifdef DEAL_II_WITH_MPI
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+#  else
+          NV_Ith_S(dst, i)            = src[is.nth_index_in_set(i)];
+#  endif
+        }
+    }
+
 #  ifdef DEAL_II_WITH_MPI
 
 #    ifdef DEAL_II_WITH_TRILINOS
diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc
index 83327e55f2..ce3cb44077 100644
--- a/source/sundials/kinsol.cc
+++ b/source/sundials/kinsol.cc
@@ -23,6 +23,8 @@
 #  include <deal.II/base/utilities.h>
 
 #  include <deal.II/lac/block_vector.h>
+#  include <deal.II/lac/la_parallel_block_vector.h>
+#  include <deal.II/lac/la_parallel_vector.h>
 #  ifdef DEAL_II_WITH_TRILINOS
 #    include <deal.II/lac/trilinos_parallel_block_vector.h>
 #    include <deal.II/lac/trilinos_vector.h>
@@ -612,6 +614,9 @@ namespace SUNDIALS
   template class KINSOL<Vector<double>>;
   template class KINSOL<BlockVector<double>>;
 
+  template class KINSOL<LinearAlgebra::distributed::Vector<double>>;
+  template class KINSOL<LinearAlgebra::distributed::BlockVector<double>>;
+
 #  ifdef DEAL_II_WITH_MPI
 
 #    ifdef DEAL_II_WITH_TRILINOS
-- 
2.39.5