Migrate timestepping to direct SUNDIALS wrappers

CMakeLists.txt

@@ -332,7 +332,7 @@ option(BUILD_SUNDIALS
 set(DiFfRG_BOOST_MIN_VERSION 1.81)
 set(DiFfRG_TBB_MIN_VERSION 2021)
 set(DiFfRG_HDF5_MIN_VERSION 1.12.0)
-set(DiFfRG_SUNDIALS_MIN_VERSION 5.4.0) # deal.II's own SUNDIALS floor
+set(DiFfRG_SUNDIALS_MIN_VERSION 7.7.0) # direct SUNDIALS API target
 set(DiFfRG_SUNDIALS_KNOWN_GOOD_MAX 7.7) # newest version validated with bundled
                                         # deal.II
 
@@ -516,10 +516,8 @@ else()
 endif()
 
 # ----------------------------------------------------------------------------
-# SUNDIALS — needed only by deal.II (DiFfRG links it transitively via
-# dealii::dealii). Replicate deal.II's own discovery: the
-# arkode/kinsol/nvecserial libraries plus version >= 5.4.0; otherwise build the
-# bundled copy.
+# SUNDIALS — DiFfRG links directly to IDA/KINSOL. Require the direct solver
+# libraries plus version >= 7.7.0; otherwise build the bundled copy.
 # ----------------------------------------------------------------------------
 set(_build_sundials OFF)
 set(_sundials_prefix "${DEPEND_INSTALL_PREFIX}")
@@ -536,8 +534,8 @@ endif()
 
 if(NOT _build_sundials)
   find_library(
-    _SUN_ARKODE
-    NAMES sundials_arkode
+    _SUN_IDA
+    NAMES sundials_ida
     HINTS ${_sun_hint}
     PATH_SUFFIXES lib lib64)
   find_library(
@@ -561,7 +559,7 @@ if(NOT _build_sundials)
   set(_sun_ok FALSE)
   set(_sun_ver "")
   set(_sun_reject_reason "")
-  if(_SUN_ARKODE
+  if(_SUN_IDA
      AND _SUN_KINSOL
      AND _SUN_NVEC
      AND _SUN_INC)
@@ -590,12 +588,9 @@ if(NOT _build_sundials)
       string(REGEX REPLACE "^.*PATCH[ \t]+([0-9]+).*" "\\1" _sun_pt
                            "${_sun_pt}")
       set(_sun_ver "${_sun_mj}.${_sun_mn}.${_sun_pt}")
-      # deal.II requires deal.II and SUNDIALS to agree on MPI support; otherwise
-      # it silently sets DEAL_II_WITH_SUNDIALS=OFF, which leaves the
-      # dealii::SUNDIALS namespace undefined and breaks our KINSOL/IDA wrappers
-      # at compile time. Probe the system SUNDIALS MPI setting
-      # (sundials_config.h) and require it to match the MPI option this build
-      # passes to deal.II.
+      # Probe the system SUNDIALS MPI setting (sundials_config.h) and require
+      # it to match this build's MPI option. The direct wrapper implemented in
+      # DiFfRG currently supports the serial N_Vector path only.
       set(_sun_mpi OFF)
       if(EXISTS "${_SUN_INC}/sundials/sundials_config.h")
         file(STRINGS "${_SUN_INC}/sundials/sundials_config.h" _sun_mpi_probe
@@ -629,8 +624,8 @@ if(NOT _build_sundials)
       message(
         WARNING
           "Using system SUNDIALS ${_sun_ver}, newer than the last version "
-          "(${DiFfRG_SUNDIALS_KNOWN_GOOD_MAX}) validated with the bundled deal.II. "
-          "If the deal.II build fails, configure with -DBUILD_SUNDIALS=ON.")
+          "(${DiFfRG_SUNDIALS_KNOWN_GOOD_MAX}) validated with DiFfRG's direct SUNDIALS interface. "
+          "If the DiFfRG build fails, configure with -DBUILD_SUNDIALS=ON.")
     endif()
     set(_sundials_status
         "  ${BoldGreen}[SUNDIALS] Using system SUNDIALS ${_sun_ver} (${_SUN_INC}); skipping bundled build.${ColourReset}"
@@ -642,19 +637,19 @@ if(NOT _build_sundials)
         FATAL_ERROR
           "SUNDIALS_DIR=${SUNDIALS_DIR} provides SUNDIALS ${_sun_ver} built with "
           "SUNDIALS_MPI_ENABLED=${_sun_mpi}, but DiFfRG is configured with MPI=${MPI}. "
-          "deal.II requires deal.II and SUNDIALS to agree on MPI support.\n"
+          "The direct SUNDIALS wrapper currently requires matching MPI support.\n"
           "  Use -DBUILD_SUNDIALS=ON to build a matching bundled SUNDIALS, or point "
           "SUNDIALS_DIR at a SUNDIALS whose MPI setting matches MPI=${MPI}.")
     else()
       set(_sundials_status
-          "  ${BoldYellow}[SUNDIALS] System SUNDIALS ${_sun_ver} (${_SUN_INC}) MPI setting (SUNDIALS_MPI_ENABLED=${_sun_mpi}) does not match MPI=${MPI}; building the bundled SUNDIALS instead. deal.II requires both to agree on MPI.${ColourReset}"
+          "  ${BoldYellow}[SUNDIALS] System SUNDIALS ${_sun_ver} (${_SUN_INC}) MPI setting (SUNDIALS_MPI_ENABLED=${_sun_mpi}) does not match MPI=${MPI}; building the bundled SUNDIALS instead.${ColourReset}"
       )
       set(_build_sundials ON)
     endif()
   elseif(_sun_explicit)
     if(NOT _sun_reject_reason)
       set(_sun_reject_reason
-          "the required libraries (arkode, kinsol, nvecserial) or version header were not found"
+          "the required libraries (ida, kinsol, nvecserial) or version header were not found"
       )
     endif()
     message(
@@ -676,21 +671,6 @@ if(_build_sundials)
   set(_sundials_dep sundials_dep)
 endif()
 
-# deal.II's FindSUNDIALS, when handed an explicit SUNDIALS_DIR, searches only
-# <DIR>/lib and <DIR>/lib64 -- it misses Debian/Ubuntu multiarch paths such as
-# /usr/lib/x86_64-linux-gnu, so it finds the headers but not the libraries and
-# then *silently* disables SUNDIALS (it is auto/optional). DiFfRG, however,
-# requires it (SUNDIALS::KINSOL/IDA), so that failure only surfaces later when
-# kinsol.cc fails to compile. Therefore: only hand deal.II the hint for a
-# bundled build or an explicit -DSUNDIALS_DIR (both have a plain lib/); for an
-# auto-detected system install, omit the hint and let deal.II's standard,
-# multiarch-aware search locate the libraries.
-if(_build_sundials OR _sun_explicit)
-  set(_dealii_sundials_arg "-DSUNDIALS_DIR:STRING=${_sundials_prefix}")
-else()
-  set(_dealii_sundials_arg "")
-endif()
-
 message("${BoldCyan}[1/8] Configuring HDF5...${ColourReset}")
 if(_build_hdf5)
   ExternalProject_Add(
@@ -1003,10 +983,7 @@ ExternalProject_Add(
     -DCMAKE_CXX_EXTENSIONS:BOOL=OFF
     -DKOKKOS_DIR:STRING=${DEPEND_INSTALL_PREFIX}
     -DTBB_DIR:STRING=${_tbb_dir}
-    # Required by DiFfRG -- force ON so an unusable SUNDIALS fails here at
-    # configure rather than silently disabling and breaking kinsol.cc later.
-    -DDEAL_II_WITH_SUNDIALS:BOOL=ON
-    ${_dealii_sundials_arg}
+    -DDEAL_II_WITH_SUNDIALS:BOOL=OFF
     -DBOOST_ROOT:PATH=${_boost_root}
     -DBOOST_DIR:STRING=${_boost_root}
     -DBOOST_INCLUDEDIR:PATH=${_boost_incdir}

DiFfRG/cmake/setup_build_system.cmake

@@ -204,6 +204,20 @@ endif()
 diffrg_find_package(Kokkos HINTS ${BUNDLED_DIR})
 message(STATUS "Found Kokkos in ${Kokkos_DIR}")
 
+# Find SUNDIALS directly. DiFfRG uses the C API for IDA/KINSOL and no longer
+# relies on deal.II's optional SUNDIALS wrappers.
+diffrg_find_package(
+  SUNDIALS
+  VERSION
+  7.7.0
+  HINTS
+  ${BUNDLED_DIR}
+  COMPONENTS
+  ida
+  kinsol
+  nvecserial)
+message(STATUS "Found SUNDIALS in ${SUNDIALS_DIR}")
+
 # Find Boost. find_package also honors BOOST_ROOT/Boost_DIR and standard system
 # paths, so a system Boost (selected via BOOST_DIR/BUILD_BOOST in the top-level
 # build) is picked up here when BUNDLED_DIR does not contain one. Use Boost's own
@@ -389,6 +403,7 @@ endfunction()
 _diffrg_summary_row("deal.II"  "${deal.II_VERSION}"  "${deal.II_DIR}")
 _diffrg_summary_row("TBB"      "${TBB_VERSION}"      "${TBB_DIR}")
 _diffrg_summary_row("Kokkos"   "${Kokkos_VERSION}"   "${Kokkos_DIR}")
+_diffrg_summary_row("SUNDIALS" "${SUNDIALS_VERSION}" "${SUNDIALS_DIR}")
 _diffrg_summary_row("Boost"    "${Boost_VERSION}"    "${_boost_path}")
 _diffrg_summary_row("Eigen3"   "${Eigen3_VERSION}"   "${Eigen3_DIR}")
 _diffrg_summary_row("GSL"      "${GSL_VERSION}"      "${GSL_INCLUDE_DIRS}")
@@ -459,6 +474,7 @@ function(setup_target TARGET)
   target_link_libraries(${TARGET} PUBLIC ${Boost_LIBRARIES})
   target_link_libraries(${TARGET} PUBLIC TBB::tbb)
   target_link_libraries(${TARGET} PUBLIC Kokkos::kokkos)
+  target_link_libraries(${TARGET} PUBLIC SUNDIALS::ida SUNDIALS::kinsol SUNDIALS::nvecserial)
   # target_link_libraries(${TARGET} PUBLIC petsc)
 
   if(${DiFfRG_MPI})

DiFfRG/include/DiFfRG/timestepping/solver/kinsol.hh

@@ -2,13 +2,21 @@
 
 // standard library
 #include <functional>
+#include <memory>
+#include <vector>
 
 // external libraries
-#include <deal.II/sundials/kinsol.h>
+#include <deal.II/base/config.h>
 
 namespace DiFfRG
 {
   using namespace dealii;
+
+  namespace sundials
+  {
+    template <typename VectorType> class KINSOL;
+  }
+
   /**
    * @brief A newton solver, using local error estimates for each vector component.
    *
@@ -21,6 +29,7 @@ namespace DiFfRG
 
     KINSOL(double abstol_ = 1e-10, double reltol_ = 1e-6, double assemble_threshold_ = 0., unsigned int max_steps_ = 21,
            unsigned int n_stepsize_iterations_ = 21);
+    ~KINSOL();
 
     /**
      * @brief Should it be necessary, one can force a recalculation of the jacobian hereby.
@@ -105,6 +114,6 @@ namespace DiFfRG
      */
     bool check();
 
-    std::shared_ptr<SUNDIALS::KINSOL<VectorType>> kinsol;
+    std::unique_ptr<sundials::KINSOL<VectorType>> kinsol;
   };
-} // namespace DiFfRG
+} // namespace DiFfRG

DiFfRG/include/DiFfRG/timestepping/sundials/common.hh

@@ -0,0 +1,94 @@
+#pragma once
+
+#include <sundials/sundials_context.h>
+#include <sundials/sundials_types.h>
+
+#include <exception>
+#include <stdexcept>
+#include <string>
+
+namespace DiFfRG::sundials
+{
+  class SolveFailure : public std::runtime_error
+  {
+  public:
+    SolveFailure(const int flag_, const char *function)
+        : std::runtime_error(std::string(function) + " failed with SUNDIALS flag " + std::to_string(flag_)),
+          flag_value(flag_)
+    {}
+
+    int flag() const noexcept { return flag_value; }
+
+  private:
+    int flag_value;
+  };
+
+  inline void check_flag(const int flag, const char *function)
+  {
+    if (flag < 0)
+      throw std::runtime_error(std::string(function) + " failed with SUNDIALS flag " + std::to_string(flag));
+  }
+
+  inline void check_solve_flag(const int flag, const char *function)
+  {
+    if (flag < 0) throw SolveFailure(flag, function);
+  }
+
+  class Context
+  {
+  public:
+    Context()
+    {
+      check_flag(SUNContext_Create(SUN_COMM_NULL, &context), "SUNContext_Create");
+    }
+
+    Context(const Context &) = delete;
+    Context &operator=(const Context &) = delete;
+
+    Context(Context &&other) noexcept : context(other.context) { other.context = nullptr; }
+
+    Context &operator=(Context &&other) noexcept
+    {
+      if (this != &other) {
+        reset();
+        context = other.context;
+        other.context = nullptr;
+      }
+      return *this;
+    }
+
+    ~Context() { reset(); }
+
+    SUNContext get() const { return context; }
+
+  private:
+    void reset() noexcept
+    {
+      if (context != nullptr) {
+        SUNContext_Free(&context);
+        context = nullptr;
+      }
+    }
+
+    SUNContext context = nullptr;
+  };
+
+  inline int callback_failed(std::exception_ptr &pending_exception)
+  {
+    try {
+      throw;
+    } catch (...) {
+      pending_exception = std::current_exception();
+      return -1;
+    }
+  }
+
+  inline void rethrow_pending(std::exception_ptr &pending_exception)
+  {
+    if (pending_exception) {
+      auto exception = pending_exception;
+      pending_exception = nullptr;
+      std::rethrow_exception(exception);
+    }
+  }
+} // namespace DiFfRG::sundials