Custom NVector for GPU

cmake/SetupBOUTThirdParty.cmake

@@ -285,12 +285,8 @@ if(BOUT_USE_ADIOS2)
     find_package(ADIOS2 REQUIRED)
   endif()
 
-
-  foreach(_adios2_candidate
-    adios2::cxx_mpi
-    adios2::cxx20_mpi
-    adios2::cxx17_mpi
-    adios2::cxx11_mpi
+  foreach(_adios2_candidate adios2::cxx_mpi adios2::cxx20_mpi adios2::cxx17_mpi
+                            adios2::cxx11_mpi
   )
     if(TARGET ${_adios2_candidate})
       set(_adios2_target ${_adios2_candidate})
@@ -299,8 +295,7 @@ if(BOUT_USE_ADIOS2)
   endforeach()
 
   if(NOT DEFINED _adios2_target)
-    message(FATAL_ERROR
-      "Could not find a usable ADIOS2 CXX CMake target. ")
+    message(FATAL_ERROR "Could not find a usable ADIOS2 CXX CMake target. ")
   endif()
   message(STATUS "Using ADIOS2 CMake target: ${_adios2_target}")
 
@@ -358,6 +353,7 @@ cmake_dependent_option(
   BOUT_USE_SUNDIALS "Enable support for SUNDIALS time solvers" OFF
   "NOT BOUT_DOWNLOAD_SUNDIALS" ON
 )
+set(BOUT_HAS_SUNDIALS_MANYVECTOR OFF)
 if(BOUT_USE_SUNDIALS)
   enable_language(C)
   if(BOUT_DOWNLOAD_SUNDIALS)
@@ -411,6 +407,14 @@ if(BOUT_USE_SUNDIALS)
     set(SUNDIALS_ROOT "${SUNDIALS_DIR}")
   endif()
   target_link_libraries(bout++ PUBLIC SUNDIALS::nvecparallel)
+  if(TARGET SUNDIALS::nvecmanyvector)
+    target_link_libraries(bout++ PUBLIC SUNDIALS::nvecmanyvector)
+    set(BOUT_HAS_SUNDIALS_MANYVECTOR ON)
+  else()
+    message(
+      STATUS "SUNDIALS ManyVector support not found; custom N_Vector disabled"
+    )
+  endif()
   target_link_libraries(bout++ PUBLIC SUNDIALS::cvode)
   target_link_libraries(bout++ PUBLIC SUNDIALS::ida)
   target_link_libraries(bout++ PUBLIC SUNDIALS::arkode)

cmake_build_defines.hxx.in

@@ -20,6 +20,7 @@
 #cmakedefine01 BOUT_HAS_SCOREP
 #cmakedefine01 BOUT_HAS_SLEPC
 #cmakedefine01 BOUT_HAS_SUNDIALS
+#cmakedefine01 BOUT_HAS_SUNDIALS_MANYVECTOR
 #cmakedefine01 BOUT_HAS_UUID_SYSTEM_GENERATOR
 #cmakedefine01 BOUT_USE_COLOR
 #cmakedefine01 BOUT_USE_OPENMP

include/bout/build_config.hxx

@@ -30,6 +30,7 @@ constexpr auto has_uuid_system_generator =
     static_cast<bool>(BOUT_HAS_UUID_SYSTEM_GENERATOR);
 constexpr auto has_slepc = static_cast<bool>(BOUT_HAS_SLEPC);
 constexpr auto has_sundials = static_cast<bool>(BOUT_HAS_SUNDIALS);
+constexpr auto has_sundials_manyvector = static_cast<bool>(BOUT_HAS_SUNDIALS_MANYVECTOR);
 constexpr auto use_backtrace = static_cast<bool>(BOUT_USE_BACKTRACE);
 constexpr auto use_color = static_cast<bool>(BOUT_USE_COLOR);
 constexpr auto use_openmp = static_cast<bool>(BOUT_USE_OPENMP);

include/bout/field2d.hxx

@@ -274,6 +274,7 @@ public:
   void applyTDerivBoundary() override;
   void setBoundaryTo(const Field2D& f2d); ///< Copy the boundary region
 
+  void swapData(Field2D& other);
   friend void swap(Field2D& first, Field2D& second) noexcept;
 
   int size() const override { return nx * ny; }

include/bout/field3d.hxx

@@ -536,6 +536,7 @@ public:
   void applyParallelBoundary(const std::string& region, const std::string& condition,
                              Field3D* f);
 
+  void swapData(Field3D& other);
   friend void swap(Field3D& first, Field3D& second) noexcept;
 
   int size() const override { return nx * ny * nz; };

include/bout/solver.hxx

@@ -50,6 +50,7 @@
 // C function pointer types
 
 class Solver;
+class SundialsNVectorInterface;
 
 /// RHS function pointer
 using rhsfunc = int (*)(BoutReal);
@@ -358,6 +359,8 @@ public:
   int getIterationOffset() const { return iteration_offset; }
 
 protected:
+  friend class SundialsNVectorInterface;
+
   /// Number of command-line arguments
   static int* pargc;
   /// Command-line arguments

include/bout/sundials_backports.hxx

@@ -32,10 +32,11 @@
 
 // NOLINTBEGIN(cppcoreguidelines-macro-usage)
 #define SUNDIALS_VERSION_AT_LEAST(major, minor, patch) \
-  ((major) < SUNDIALS_VERSION_MAJOR                      \
-   || ((major) == SUNDIALS_VERSION_MAJOR                 \
-       && ((minor) < SUNDIALS_VERSION_MINOR              \
-           || ((minor) == SUNDIALS_VERSION_MINOR && (patch) <= SUNDIALS_VERSION_PATCH))))
+  ((major) < SUNDIALS_VERSION_MAJOR                    \
+   || ((major) == SUNDIALS_VERSION_MAJOR               \
+       && ((minor) < SUNDIALS_VERSION_MINOR            \
+           || ((minor) == SUNDIALS_VERSION_MINOR       \
+               && (patch) <= SUNDIALS_VERSION_PATCH))))
 #define SUNDIALS_VERSION_LESS_THAN(major, minor, patch) \
   (!SUNDIALS_VERSION_AT_LEAST(major, minor, patch))
 // NOLINTEND(cppcoreguidelines-macro-usage)
@@ -77,13 +78,13 @@ inline sundials::Context createSUNContext([[maybe_unused]] MPI_Comm& comm) {
 #endif
 }
 
-template <typename Func, typename... Args>
-inline decltype(auto) callWithSUNContext(Func f, [[maybe_unused]] sundials::Context& ctx,
+template <typename Func, typename Ctx, typename... Args>
+inline decltype(auto) callWithSUNContext(Func f, [[maybe_unused]] Ctx&& ctx,
                                          Args&&... args) {
 #if SUNDIALS_VERSION_LESS_THAN(6, 0, 0)
   return f(std::forward<Args>(args)...);
 #else
-  return f(std::forward<Args>(args)..., ctx);
+  return f(std::forward<Args>(args)..., std::forward<Ctx>(ctx));
 #endif
 }
 

manual/sphinx/user_docs/time_integration.rst

@@ -33,7 +33,7 @@ needed to make the solver available.
 
 .. _tab-solvers:
 .. table:: Available time integration solvers
-	   
+
    +---------------+-----------------------------------------+------------------------+
    | Name          | Description                             | Compile options        |
    +===============+=========================================+========================+
@@ -68,7 +68,7 @@ given in table :numref:`tab-solveropts`.
 
 .. _tab-solveropts:
 .. table:: Time integration solver options
-	   
+
    +--------------------------+--------------------------------------------+-------------------------------------+
    | Option                   | Description                                | Solvers used                        |
    +==========================+============================================+=====================================+
@@ -104,13 +104,40 @@ given in table :numref:`tab-solveropts`.
    +--------------------------+--------------------------------------------+-------------------------------------+
    | diagnose                 | Collect and print additional diagnostics   | cvode, imexbdf2, beuler             |
    +--------------------------+--------------------------------------------+-------------------------------------+
+   | nvector                  | ``N_Vector`` backend for SUNDIALS solvers: | cvode, ida, arkode                  |
+   |                          | ``sundials`` or ``manyvector``             |                                     |
+   +--------------------------+--------------------------------------------+-------------------------------------+
 
 |
 
 The most commonly changed options are the absolute and relative solver
 tolerances, ``atol`` and ``rtol`` which should be varied to check
 convergence.
 
+SUNDIALS ``N_Vector`` backends
+------------------------------
+
+The SUNDIALS-based solvers ``cvode``, ``ida``, and ``arkode`` can select
+the ``N_Vector`` backend at runtime using ``solver:nvector``:
+
+.. code-block:: cfg
+
+    [solver]
+    type = cvode
+    nvector = sundials
+
+Valid values are:
+
+- ``sundials`` uses the standard SUNDIALS parallel ``N_Vector``. This is the
+  default.
+- ``manyvector`` uses the BOUT++ field-backed custom ``N_Vector`` built on top
+  of SUNDIALS ManyVector support.
+
+The ``manyvector`` option is only available when BOUT++ was built with SUNDIALS
+ManyVector support. If ``solver:nvector=manyvector`` is selected in a build
+that does not provide this support, solver initialisation will throw an
+exception.
+
 CVODE
 -----
 

src/field/field2d.cxx

@@ -327,6 +327,8 @@ void Field2D::setBoundaryTo(const Field2D& f2d) {
   }
 }
 
+void Field2D::swapData(Field2D& other) { std::swap(data, other.data); }
+
 ////////////// NON-MEMBER OVERLOADED OPERATORS //////////////
 
 // Unary minus

src/field/field3d.cxx

@@ -670,6 +670,8 @@ void Field3D::applyParallelBoundary(const std::string& region,
   }
 }
 
+void Field3D::swapData(Field3D& other) { std::swap(data, other.data); }
+
 /***************************************************************
  *               NON-MEMBER OVERLOADED OPERATORS
  ***************************************************************/