]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert "Remove sudials/copy.h"
authorDavid Wells <drwells@email.unc.edu>
Fri, 3 Jun 2022 15:03:56 +0000 (11:03 -0400)
committerDavid Wells <drwells@email.unc.edu>
Fri, 3 Jun 2022 17:09:07 +0000 (13:09 -0400)
This reverts commit 51b6f10f36a4dd61ed4daf20c137cc48c19de5f0.

While we are here - patch the copy test to work with SUNDIALS 6.

include/deal.II/sundials/copy.h [new file with mode: 0644]
source/sundials/CMakeLists.txt
source/sundials/copy.cc [new file with mode: 0644]
tests/sundials/copy_01.cc [new file with mode: 0644]
tests/sundials/copy_01.with_mpi=true.with_petsc=true.with_petsc_with_complex=false.mpirun=4.output [new file with mode: 0644]

diff --git a/include/deal.II/sundials/copy.h b/include/deal.II/sundials/copy.h
new file mode 100644 (file)
index 0000000..afc2ae3
--- /dev/null
@@ -0,0 +1,104 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2017 - 2018 by the deal.II authors
+//
+//    This file is part of the deal.II library.
+//
+//    The deal.II library is free software; you can use it, redistribute
+//    it, and/or modify it under the terms of the GNU Lesser General
+//    Public License as published by the Free Software Foundation; either
+//    version 2.1 of the License, or (at your option) any later version.
+//    The full text of the license can be found in the file LICENSE.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#ifndef dealii_sundials_copy_h
+#define dealii_sundials_copy_h
+
+#include <deal.II/base/config.h>
+#ifdef DEAL_II_WITH_SUNDIALS
+
+#  include <sundials/sundials_nvector.h>
+#  ifdef DEAL_II_WITH_MPI
+#    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>
+
+#  ifdef DEAL_II_WITH_TRILINOS
+#    include <deal.II/lac/trilinos_parallel_block_vector.h>
+#    include <deal.II/lac/trilinos_vector.h>
+#  endif
+
+#  ifdef DEAL_II_WITH_PETSC
+#    include <deal.II/lac/petsc_block_vector.h>
+#    include <deal.II/lac/petsc_vector.h>
+#  endif
+
+DEAL_II_NAMESPACE_OPEN
+namespace SUNDIALS
+{
+  namespace internal
+  {
+    // The following internal functions are used by SUNDIALS wrappers to copy
+    // to and from deal.II vector types.
+#  ifdef DEAL_II_WITH_MPI
+
+#    ifdef DEAL_II_WITH_TRILINOS
+    void
+    copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const TrilinosWrappers::MPI::Vector &src);
+    void
+    copy(TrilinosWrappers::MPI::BlockVector &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const TrilinosWrappers::MPI::BlockVector &src);
+#    endif // DEAL_II_WITH_TRILINOS
+
+#    ifdef DEAL_II_WITH_PETSC
+#      ifndef PETSC_USE_COMPLEX
+    void
+    copy(PETScWrappers::MPI::Vector &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const PETScWrappers::MPI::Vector &src);
+    void
+    copy(PETScWrappers::MPI::BlockVector &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const PETScWrappers::MPI::BlockVector &src);
+#      endif // PETSC_USE_COMPLEX
+#    endif   // DEAL_II_WITH_PETSC
+
+#  endif
+
+    void
+    copy(BlockVector<double> &dst, const N_Vector &src);
+    void
+    copy(N_Vector &dst, const BlockVector<double> &src);
+
+    void
+    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
+
+#endif // DEAL_II_WITH_SUNDIALS
+#endif // dealii_sundials_copy_h
index 2c48fbc2a0f3fd8cbeb9a647adb8214f270aa81d..35458673a3c20ba2095c4710a866e14aac43b467 100644 (file)
@@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 SET(_src
   arkode.cc
   ida.cc
+  copy.cc
   kinsol.cc
   n_vector.cc
   sunlinsol_wrapper.cc
diff --git a/source/sundials/copy.cc b/source/sundials/copy.cc
new file mode 100644 (file)
index 0000000..b11b279
--- /dev/null
@@ -0,0 +1,292 @@
+//-----------------------------------------------------------
+//
+// Copyright (C) 2017 - 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/sundials/copy.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+
+DEAL_II_NAMESPACE_OPEN
+namespace SUNDIALS
+{
+  namespace internal
+  {
+    /**
+     * SUNDIALS provides different macros for getting the local length of a
+     * vector for serial and parallel vectors (as well as various parallel
+     * vectors that are not yet supported by deal.II). This function provides
+     * a generic interface to both and does a (checked) conversion from long
+     * int (the type SUNDIALS uses for lengths) to std::size_t.
+     */
+    inline std::size_t
+    N_Vector_length(const N_Vector &vec)
+    {
+      const N_Vector_ID id     = N_VGetVectorID(vec);
+      long int          length = -1;
+      switch (id)
+        {
+          case SUNDIALS_NVEC_SERIAL:
+            {
+              length = NV_LENGTH_S(vec);
+              break;
+            }
+#  ifdef DEAL_II_WITH_MPI
+          case SUNDIALS_NVEC_PARALLEL:
+            {
+              length = NV_LOCLENGTH_P(vec);
+              break;
+            }
+#  endif
+          default:
+            Assert(false, ExcNotImplemented());
+        }
+
+      Assert(length >= 0, ExcInternalError());
+      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
+
+
+    void
+    copy(TrilinosWrappers::MPI::Vector &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));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &dst, const TrilinosWrappers::MPI::Vector &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)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+    void
+    copy(TrilinosWrappers::MPI::BlockVector &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));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &dst, const TrilinosWrappers::MPI::BlockVector &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)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+#    endif // DEAL_II_WITH_TRILINOS
+
+#    ifdef DEAL_II_WITH_PETSC
+#      ifndef PETSC_USE_COMPLEX
+
+    void
+    copy(PETScWrappers::MPI::Vector &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));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &dst, const PETScWrappers::MPI::Vector &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)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+    void
+    copy(PETScWrappers::MPI::BlockVector &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));
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[is.nth_index_in_set(i)] = NV_Ith_P(src, i);
+        }
+      dst.compress(VectorOperation::insert);
+    }
+
+    void
+    copy(N_Vector &dst, const PETScWrappers::MPI::BlockVector &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)
+        {
+          NV_Ith_P(dst, i) = src[is.nth_index_in_set(i)];
+        }
+    }
+
+#      endif // PETSC_USE_COMPLEX
+#    endif   // DEAL_II_WITH_PETSC
+
+#  endif // mpi
+
+    void
+    copy(BlockVector<double> &dst, const N_Vector &src)
+    {
+      const std::size_t N = dst.size();
+      AssertDimension(N_Vector_length(src), N);
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[i] = NV_Ith_S(src, i);
+        }
+    }
+
+    void
+    copy(N_Vector &dst, const BlockVector<double> &src)
+    {
+      const std::size_t N = src.size();
+      AssertDimension(N_Vector_length(dst), N);
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          NV_Ith_S(dst, i) = src[i];
+        }
+    }
+
+    void
+    copy(Vector<double> &dst, const N_Vector &src)
+    {
+      const std::size_t N = dst.size();
+      AssertDimension(N_Vector_length(src), N);
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          dst[i] = NV_Ith_S(src, i);
+        }
+    }
+
+    void
+    copy(N_Vector &dst, const Vector<double> &src)
+    {
+      const std::size_t N = src.size();
+      AssertDimension(N_Vector_length(dst), N);
+      for (std::size_t i = 0; i < N; ++i)
+        {
+          NV_Ith_S(dst, i) = src[i];
+        }
+    }
+  } // namespace internal
+} // namespace SUNDIALS
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/sundials/copy_01.cc b/tests/sundials/copy_01.cc
new file mode 100644 (file)
index 0000000..eff6ccf
--- /dev/null
@@ -0,0 +1,97 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Verify that we can correctly copy a PETSc vector to a parallel (non-PETSc)
+// SUNDIALS vector.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/petsc_vector.h>
+
+#include <deal.II/sundials/copy.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log_all;
+
+  MPI_Comm communicator = MPI_COMM_WORLD;
+
+  const unsigned int current_process =
+    Utilities::MPI::this_mpi_process(communicator);
+  const unsigned int n_processes =
+    Utilities::MPI::n_mpi_processes(communicator);
+
+  // give the zeroth process 1 dof, the first 2, etc
+  const types::global_dof_index n_local_dofs = 1 + current_process;
+  const types::global_dof_index begin_dof =
+    (current_process * (current_process + 1)) / 2;
+  const types::global_dof_index end_dof = begin_dof + n_local_dofs;
+
+  const types::global_dof_index n_global_dofs =
+    ((n_processes - 1) * (n_processes)) / 2 + n_processes;
+  IndexSet local_dofs(n_global_dofs);
+  local_dofs.add_range(begin_dof, end_dof);
+  local_dofs.compress();
+  AssertDimension(local_dofs.n_elements(), n_local_dofs);
+
+
+  PETScWrappers::MPI::Vector vec(local_dofs, communicator);
+  for (const types::global_dof_index local_dof : local_dofs)
+    {
+      vec[local_dof] = 2 * local_dof - current_process;
+    }
+  vec.compress(VectorOperation::insert);
+  const std::size_t n_dofs = vec.size();
+
+#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+  SUNContext context;
+  int        status = SUNContext_Create(&communicator, &context);
+  AssertThrow(status == 0, ExcInternalError());
+#endif
+
+  N_Vector sundials_vector = N_VNew_Parallel(communicator,
+                                             n_local_dofs,
+                                             n_dofs
+#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+                                             ,
+                                             context
+#endif
+  );
+  SUNDIALS::internal::copy(sundials_vector, vec);
+
+  PETScWrappers::MPI::Vector vec2(local_dofs, communicator);
+  SUNDIALS::internal::copy(vec2, sundials_vector);
+
+  AssertThrow(vec2 == vec,
+              ExcMessage(
+                "The two PETSc vectors should be equal since they are "
+                "copies (by means of an intermediate SUNDIALs vector)"));
+
+  deallog << "n_local_dofs: " << n_local_dofs << std::endl;
+  deallog << "OK" << std::endl;
+
+  N_VDestroy_Parallel(sundials_vector);
+#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+  status = SUNContext_Free(&context);
+  AssertThrow(status == 0, ExcInternalError());
+#endif
+}
diff --git a/tests/sundials/copy_01.with_mpi=true.with_petsc=true.with_petsc_with_complex=false.mpirun=4.output b/tests/sundials/copy_01.with_mpi=true.with_petsc=true.with_petsc_with_complex=false.mpirun=4.output
new file mode 100644 (file)
index 0000000..ea6a78c
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL:0::n_local_dofs: 1
+DEAL:0::OK
+
+DEAL:1::n_local_dofs: 2
+DEAL:1::OK
+
+
+DEAL:2::n_local_dofs: 3
+DEAL:2::OK
+
+
+DEAL:3::n_local_dofs: 4
+DEAL:3::OK
+

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.