]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Extract SUNLinSolver wrapper from ARKode 11724/head
authorSebastian Proell <proell@lnm.mw.tum.de>
Wed, 10 Feb 2021 13:11:22 +0000 (14:11 +0100)
committerSebastian Proell <proell@lnm.mw.tum.de>
Wed, 10 Feb 2021 17:19:35 +0000 (18:19 +0100)
include/deal.II/sundials/arkode.h
include/deal.II/sundials/sunlinsol_wrapper.h [new file with mode: 0644]
source/sundials/CMakeLists.txt
source/sundials/arkode.cc
source/sundials/sunlinsol_wrapper.cc [new file with mode: 0644]

index 0196c33a81ba17875d92b6fc8e0099e87284d212..a6b59183123ed9983a1629392391c0b8c9885cc0 100644 (file)
@@ -42,6 +42,7 @@
 #    include <nvector/nvector_parallel.h>
 #  endif
 #  include <deal.II/sundials/n_vector.h>
+#  include <deal.II/sundials/sunlinsol_wrapper.h>
 
 #  include <boost/signals2.hpp>
 
 
 DEAL_II_NAMESPACE_OPEN
 
-// Forward declarations
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-#    ifndef DOXYGEN
-namespace SUNDIALS
-{
-  // Forward declaration
-  template <typename VectorType>
-  struct SundialsOperator;
-
-  template <typename VectorType>
-  struct SundialsPreconditioner;
-
-  template <typename VectorType>
-  class SundialsLinearSolverWrapper;
-} // namespace SUNDIALS
-#    endif
-#  endif
 
 // Shorthand notation for ARKODE error codes.
 #  define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
@@ -80,41 +64,6 @@ namespace SUNDIALS
  */
 namespace SUNDIALS
 {
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-  /**
-   * Type of function objects to interface with SUNDIALS linear solvers
-   *
-   * This function type encapsulates the action of solving $P^{-1}Ax=P^{-1}b$.
-   * The LinearOperator @p op encapsulates the matrix vector product $Ax$ and
-   * the LinearOperator @p prec encapsulates the application of the
-   * preconditioner $P^{-1}z$.
-   * The user can specify function objects of this type to attach custom linear
-   * solver routines to SUNDIALS. The two LinearOperators @p op and @p prec are
-   * built internally by SUNDIALS based on user settings. The parameters are
-   * interpreted as follows:
-   *
-   * @param[in] op A LinearOperator that applies the matrix vector product
-   * @param[in] prec A LinearOperator that applies the preconditioner
-   * @param[out] x The output solution vector
-   * @param[in] b The right-hand side
-   * @param[in] tol Tolerance for the iterative solver
-   *
-   * This function should return:
-   * - 0: Success
-   * - >0: Recoverable error, ARKode will reattempt the solution and call this
-   *       function again.
-   * - <0: Unrecoverable error, the computation will be aborted and an
-   *       assertion will be thrown.
-   */
-  template <typename VectorType>
-  using LinearSolveFunction =
-    std::function<int(SundialsOperator<VectorType> &      op,
-                      SundialsPreconditioner<VectorType> &prec,
-                      VectorType &                        x,
-                      const VectorType &                  b,
-                      double                              tol)>;
-#  endif
-
   /**
    * Interface to SUNDIALS additive Runge-Kutta methods (ARKode).
    *
@@ -1336,8 +1285,8 @@ namespace SUNDIALS
     GrowingVectorMemory<VectorType> mem;
 
 #  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-    std::unique_ptr<SundialsLinearSolverWrapper<VectorType>> linear_solver;
-    std::unique_ptr<SundialsLinearSolverWrapper<VectorType>> mass_solver;
+    std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
+    std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
 #  endif
 
 #  ifdef DEAL_II_WITH_PETSC
@@ -1354,91 +1303,6 @@ namespace SUNDIALS
 #  endif   // DEAL_II_WITH_PETSC
   };
 
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-
-  /**
-   * A linear operator that wraps SUNDIALS functionality.
-   */
-  template <typename VectorType>
-  struct SundialsOperator
-  {
-    /**
-     * Apply this LinearOperator to @p src and store the result in @p dst.
-     */
-    void
-    vmult(VectorType &dst, const VectorType &src) const;
-
-    /**
-     * Constructor.
-     *
-     * @param A_data Data required by @p a_times_fn
-     * @param a_times_fn A function pointer to the function that computes A*v
-     */
-    SundialsOperator(void *A_data, ATimesFn a_times_fn);
-
-  private:
-    /**
-     * Data necessary to evaluate a_times_fn.
-     */
-    void *A_data;
-
-    /**
-     * Function pointer declared by SUNDIALS to evaluate the matrix vector
-     * product.
-     */
-    ATimesFn a_times_fn;
-  };
-
-
-
-  /**
-   * A linear operator that wraps preconditioner functionality as specified by
-   * SUNDIALS. The vmult() function solves the preconditioner equation $Px=b$,
-   * i.e., it computes $x=P^{-1}b$.
-   */
-  template <typename VectorType>
-  struct SundialsPreconditioner
-  {
-    /**
-     * Apply the wrapped preconditioner, i.e., solve $Px=b$ where $x$ is the
-     * @p dst vector and $b$ the @p src vector.
-     *
-     * @param dst Result vector of the preconditioner application
-     * @param src Target vector of the preconditioner application
-     */
-    void
-    vmult(VectorType &dst, const VectorType &src) const;
-
-    /**
-     * Constructor.
-     *
-     * @param P_data Data required by @p p_solve_fn
-     * @param p_solve_fn A function pointer to the function that computes A*v
-     * @param tol Tolerance, that an iterative solver should use to judge
-     *   convergence
-     */
-    SundialsPreconditioner(void *P_data, PSolveFn p_solve_fn, double tol);
-
-  private:
-    /**
-     * Data necessary to calls p_solve_fn
-     */
-    void *P_data;
-
-    /**
-     * Function pointer to a function that computes the preconditioner
-     * application.
-     */
-    PSolveFn p_solve_fn;
-
-    /**
-     * Potential tolerance to use in the internal solve of the preconditioner
-     * equation.
-     */
-    double tol;
-  };
-
-#  endif
 
   /**
    * Handle ARKode exceptions.
diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h
new file mode 100644 (file)
index 0000000..d90dc94
--- /dev/null
@@ -0,0 +1,186 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2021 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_sunlinsol_wrapper_h
+#define dealii_sundials_sunlinsol_wrapper_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+
+#    include <sundials/sundials_linearsolver.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+#    ifndef DOXYGEN
+  // forward declarations
+  namespace internal
+  {
+    template <typename VectorType>
+    struct LinearSolverContent;
+  }
+#    endif
+
+  /**
+   * A linear operator that wraps SUNDIALS functionality.
+   */
+  template <typename VectorType>
+  struct SundialsOperator
+  {
+    /**
+     * Apply this LinearOperator to @p src and store the result in @p dst.
+     */
+    void
+    vmult(VectorType &dst, const VectorType &src) const;
+
+    /**
+     * Constructor.
+     *
+     * @param A_data Data required by @p a_times_fn
+     * @param a_times_fn A function pointer to the function that computes A*v
+     */
+    SundialsOperator(void *A_data, ATimesFn a_times_fn);
+
+  private:
+    /**
+     * Data necessary to evaluate a_times_fn.
+     */
+    void *A_data;
+
+    /**
+     * Function pointer declared by SUNDIALS to evaluate the matrix vector
+     * product.
+     */
+    ATimesFn a_times_fn;
+  };
+
+
+
+  /**
+   * A linear operator that wraps preconditioner functionality as specified by
+   * SUNDIALS. The vmult() function solves the preconditioner equation $Px=b$,
+   * i.e., it computes $x=P^{-1}b$.
+   */
+  template <typename VectorType>
+  struct SundialsPreconditioner
+  {
+    /**
+     * Apply the wrapped preconditioner, i.e., solve $Px=b$ where $x$ is the
+     * @p dst vector and $b$ the @p src vector.
+     *
+     * @param dst Result vector of the preconditioner application
+     * @param src Target vector of the preconditioner application
+     */
+    void
+    vmult(VectorType &dst, const VectorType &src) const;
+
+    /**
+     * Constructor.
+     *
+     * @param P_data Data required by @p p_solve_fn
+     * @param p_solve_fn A function pointer to the function that computes A*v
+     * @param tol Tolerance, that an iterative solver should use to judge
+     *   convergence
+     */
+    SundialsPreconditioner(void *P_data, PSolveFn p_solve_fn, double tol);
+
+  private:
+    /**
+     * Data necessary to calls p_solve_fn
+     */
+    void *P_data;
+
+    /**
+     * Function pointer to a function that computes the preconditioner
+     * application.
+     */
+    PSolveFn p_solve_fn;
+
+    /**
+     * Potential tolerance to use in the internal solve of the preconditioner
+     * equation.
+     */
+    double tol;
+  };
+
+  /**
+   * Type of function objects to interface with SUNDIALS linear solvers
+   *
+   * This function type encapsulates the action of solving $P^{-1}Ax=P^{-1}b$.
+   * The LinearOperator @p op encapsulates the matrix vector product $Ax$ and
+   * the LinearOperator @p prec encapsulates the application of the
+   * preconditioner $P^{-1}z$.
+   * The user can specify function objects of this type to attach custom linear
+   * solver routines to SUNDIALS. The two LinearOperators @p op and @p prec are
+   * built internally by SUNDIALS based on user settings. The parameters are
+   * interpreted as follows:
+   *
+   * @param[in] op A LinearOperator that applies the matrix vector product
+   * @param[in] prec A LinearOperator that applies the preconditioner
+   * @param[out] x The output solution vector
+   * @param[in] b The right-hand side
+   * @param[in] tol Tolerance for the iterative solver
+   *
+   * This function should return:
+   * - 0: Success
+   * - >0: Recoverable error, ARKode will reattempt the solution and call this
+   *       function again.
+   * - <0: Unrecoverable error, the computation will be aborted and an
+   *       assertion will be thrown.
+   */
+  template <typename VectorType>
+  using LinearSolveFunction =
+    std::function<int(SundialsOperator<VectorType> &      op,
+                      SundialsPreconditioner<VectorType> &prec,
+                      VectorType &                        x,
+                      const VectorType &                  b,
+                      double                              tol)>;
+
+  namespace internal
+  {
+    /*!
+     * Attach wrapper functions to SUNDIALS' linear solver interface. We pretend
+     * that the user-supplied linear solver is matrix-free, even though it can
+     * be matrix-based. This way SUNDIALS does not need to understand our matrix
+     * types.
+     */
+    template <typename VectorType>
+    class LinearSolverWrapper
+    {
+    public:
+      explicit LinearSolverWrapper(LinearSolveFunction<VectorType> lsolve);
+
+      ~LinearSolverWrapper();
+
+      /**
+       * Implicit conversion to SUNLinearSolver.
+       */
+      operator SUNLinearSolver();
+
+    private:
+      SUNLinearSolver                                  sun_linear_solver;
+      std::unique_ptr<LinearSolverContent<VectorType>> content;
+    };
+  } // namespace internal
+} // namespace SUNDIALS
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif
+#endif
+#endif
index d01e607689f630a4d7d6a53cc9d42b0b7178b872..0b1d3ff82fc01ced65d98d4d2fd3f815a1af8c9a 100644 (file)
@@ -21,6 +21,7 @@ SET(_src
   copy.cc
   kinsol.cc
   n_vector.cc
+  sunlinsol_wrapper.cc
   )
 
 SET(_inst
index 4bb0e2f6926fea8a6acd83f812c7894b346da9e1..31d970927c4249a55932ea3fb76b323d6030fb62 100644 (file)
@@ -36,6 +36,7 @@
 #  endif
 
 #  include <deal.II/sundials/n_vector.h>
+#  include <deal.II/sundials/sunlinsol_wrapper.h>
 
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
 #    include <arkode/arkode_impl.h>
@@ -371,166 +372,10 @@ namespace SUNDIALS
 
       return solver.mass_preconditioner_setup(t);
     }
-
-
-
-    /**
-     * storage for internal content of the linear solver wrapper
-     */
-    template <typename VectorType>
-    struct LinearSolverContent
-    {
-      ATimesFn a_times_fn;
-      PSetupFn preconditioner_setup;
-      PSolveFn preconditioner_solve;
-
-      LinearSolveFunction<VectorType> lsolve;
-
-      void *P_data;
-      void *A_data;
-    };
-
-
-
-    /**
-     * Access our LinearSolverContent from the generic content of the
-     * SUNLinearSolver @p ls.
-     */
-    template <typename VectorType>
-    LinearSolverContent<VectorType> *
-    access_content(SUNLinearSolver ls)
-    {
-      Assert(ls->content != nullptr, ExcInternalError());
-      return static_cast<LinearSolverContent<VectorType> *>(ls->content);
-    }
-
-
-
-    SUNLinearSolver_Type arkode_linsol_get_type(SUNLinearSolver)
-    {
-      return SUNLINEARSOLVER_ITERATIVE;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    arkode_linsol_solve(SUNLinearSolver LS,
-                        SUNMatrix,
-                        N_Vector x,
-                        N_Vector b,
-                        realtype tol)
-    {
-      auto content = access_content<VectorType>(LS);
-
-      auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
-      auto *dst_x = internal::unwrap_nvector<VectorType>(x);
-
-      SundialsOperator<VectorType> op(content->A_data, content->a_times_fn);
-
-      SundialsPreconditioner<VectorType> preconditioner(
-        content->P_data, content->preconditioner_solve, tol);
-
-      return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
-    }
-
-
-
-    template <typename VectorType>
-    int
-    arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix)
-    {
-      auto content = access_content<VectorType>(LS);
-      if (content->preconditioner_setup)
-        return content->preconditioner_setup(content->P_data);
-      return 0;
-    }
-
-
-
-    template <typename VectorType>
-    int arkode_linsol_initialize(SUNLinearSolver)
-    {
-      // this method is currently only provided because SUNDIALS 4.0.0 requires
-      // it - no user-set action is implemented so far
-      return 0;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
-    {
-      auto content        = access_content<VectorType>(LS);
-      content->A_data     = A_data;
-      content->a_times_fn = ATimes;
-      return 0;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    arkode_linsol_set_preconditioner(SUNLinearSolver LS,
-                                     void *          P_data,
-                                     PSetupFn        p_setup,
-                                     PSolveFn        p_solve)
-    {
-      auto content                  = access_content<VectorType>(LS);
-      content->P_data               = P_data;
-      content->preconditioner_setup = p_setup;
-      content->preconditioner_solve = p_solve;
-      return 0;
-    }
 #  endif
 
   } // namespace
 
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-
-  /*!
-   * Attach wrapper functions to SUNDIALS' linear solver interface. We pretend
-   * that the user-supplied linear solver is matrix-free, even though it can
-   * be matrix-based. This way SUNDIALS does not need to understand our matrix
-   * types.
-   */
-  template <typename VectorType>
-  class SundialsLinearSolverWrapper
-  {
-  public:
-    SundialsLinearSolverWrapper(LinearSolveFunction<VectorType> lsolve)
-    {
-      sun_linear_solver                  = SUNLinSolNewEmpty();
-      sun_linear_solver->ops->gettype    = arkode_linsol_get_type;
-      sun_linear_solver->ops->solve      = arkode_linsol_solve<VectorType>;
-      sun_linear_solver->ops->setup      = arkode_linsol_setup<VectorType>;
-      sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
-      sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
-      sun_linear_solver->ops->setpreconditioner =
-        arkode_linsol_set_preconditioner<VectorType>;
-
-      content.lsolve             = lsolve;
-      sun_linear_solver->content = &content;
-    }
-
-    ~SundialsLinearSolverWrapper()
-    {
-      SUNLinSolFreeEmpty(sun_linear_solver);
-    }
-
-    SUNLinearSolver
-    get_wrapped_solver()
-    {
-      return sun_linear_solver;
-    }
-
-  private:
-    SUNLinearSolver                 sun_linear_solver;
-    LinearSolverContent<VectorType> content;
-  };
-
-#  endif
 
   template <typename VectorType>
   ARKode<VectorType>::ARKode(const AdditionalData &data,
@@ -804,9 +649,9 @@ namespace SUNDIALS
         if (solve_linearized_system)
           {
             linear_solver =
-              std::make_unique<SundialsLinearSolverWrapper<VectorType>>(
+              std::make_unique<internal::LinearSolverWrapper<VectorType>>(
                 solve_linearized_system);
-            sun_linear_solver = linear_solver->get_wrapped_solver();
+            sun_linear_solver = *linear_solver;
           }
         else
           {
@@ -875,9 +720,9 @@ namespace SUNDIALS
         if (solve_mass)
           {
             mass_solver =
-              std::make_unique<SundialsLinearSolverWrapper<VectorType>>(
+              std::make_unique<internal::LinearSolverWrapper<VectorType>>(
                 solve_mass);
-            sun_mass_linear_solver = mass_solver->get_wrapped_solver();
+            sun_mass_linear_solver = *mass_solver;
           }
         else
           {
@@ -940,117 +785,24 @@ namespace SUNDIALS
     return arkode_mem;
   }
 
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-
-  template <typename VectorType>
-  SundialsOperator<VectorType>::SundialsOperator(void *   A_data,
-                                                 ATimesFn a_times_fn)
-    : A_data(A_data)
-    , a_times_fn(a_times_fn)
-
-  {
-    Assert(a_times_fn != nullptr, ExcInternalError());
-  }
-
-
-
-  template <typename VectorType>
-  void
-  SundialsOperator<VectorType>::vmult(VectorType &      dst,
-                                      const VectorType &src) const
-  {
-    auto sun_dst = internal::make_nvector_view(dst);
-    auto sun_src = internal::make_nvector_view(src);
-    int  status  = a_times_fn(A_data, sun_src, sun_dst);
-    (void)status;
-    AssertARKode(status);
-  }
-
-
-
-  template <typename VectorType>
-  SundialsPreconditioner<VectorType>::SundialsPreconditioner(
-    void *   P_data,
-    PSolveFn p_solve_fn,
-    double   tol)
-    : P_data(P_data)
-    , p_solve_fn(p_solve_fn)
-    , tol(tol)
-  {}
-
-
-
-  template <typename VectorType>
-  void
-  SundialsPreconditioner<VectorType>::vmult(VectorType &      dst,
-                                            const VectorType &src) const
-  {
-    // apply identity preconditioner if nothing else specified
-    if (!p_solve_fn)
-      {
-        dst = src;
-        return;
-      }
-
-    auto sun_dst = internal::make_nvector_view(dst);
-    auto sun_src = internal::make_nvector_view(src);
-    // for custom preconditioners no distinction between left and right
-    // preconditioning is made
-    int status =
-      p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
-    (void)status;
-    AssertARKode(status);
-  }
-#  endif
-
   template class ARKode<Vector<double>>;
   template class ARKode<BlockVector<double>>;
 
   template class ARKode<LinearAlgebra::distributed::Vector<double>>;
   template class ARKode<LinearAlgebra::distributed::BlockVector<double>>;
 
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-  template struct SundialsOperator<Vector<double>>;
-  template struct SundialsOperator<BlockVector<double>>;
-  template struct SundialsOperator<LinearAlgebra::distributed::Vector<double>>;
-  template struct SundialsOperator<
-    LinearAlgebra::distributed::BlockVector<double>>;
-
-  template struct SundialsPreconditioner<Vector<double>>;
-  template struct SundialsPreconditioner<BlockVector<double>>;
-  template struct SundialsPreconditioner<
-    LinearAlgebra::distributed::Vector<double>>;
-  template struct SundialsPreconditioner<
-    LinearAlgebra::distributed::BlockVector<double>>;
-#  endif
-
 #  ifdef DEAL_II_WITH_MPI
 
 #    ifdef DEAL_II_WITH_TRILINOS
   template class ARKode<TrilinosWrappers::MPI::Vector>;
   template class ARKode<TrilinosWrappers::MPI::BlockVector>;
 
-#      if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-  template struct SundialsOperator<TrilinosWrappers::MPI::Vector>;
-  template struct SundialsOperator<TrilinosWrappers::MPI::BlockVector>;
-
-  template struct SundialsPreconditioner<TrilinosWrappers::MPI::Vector>;
-  template struct SundialsPreconditioner<TrilinosWrappers::MPI::BlockVector>;
-#      endif
 #    endif // DEAL_II_WITH_TRILINOS
 
 #    ifdef DEAL_II_WITH_PETSC
 #      ifndef PETSC_USE_COMPLEX
   template class ARKode<PETScWrappers::MPI::Vector>;
   template class ARKode<PETScWrappers::MPI::BlockVector>;
-
-#        if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-  template struct SundialsOperator<PETScWrappers::MPI::Vector>;
-  template struct SundialsOperator<PETScWrappers::MPI::BlockVector>;
-
-  template struct SundialsPreconditioner<PETScWrappers::MPI::Vector>;
-  template struct SundialsPreconditioner<PETScWrappers::MPI::BlockVector>;
-#        endif
 #      endif // PETSC_USE_COMPLEX
 #    endif   // DEAL_II_WITH_PETSC
 
diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc
new file mode 100644 (file)
index 0000000..9c6010d
--- /dev/null
@@ -0,0 +1,323 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2021 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/base/config.h>
+
+#include <deal.II/sundials/sunlinsol_wrapper.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+
+#    include <deal.II/base/exceptions.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>
+#    include <deal.II/lac/vector.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
+
+#    include <deal.II/sundials/n_vector.h>
+
+#    if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
+#      include <deal.II/sundials/sunlinsol_newempty.h>
+#    endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace SUNDIALS
+{
+  DeclException1(ExcSundialsSolverError,
+                 int,
+                 << "One of the SUNDIALS linear solver internal"
+                 << " functions returned a negative error code: " << arg1
+                 << ". Please consult SUNDIALS manual.");
+
+#    define AssertSundialsSolver(code) \
+      Assert(code >= 0, ExcSundialsSolverError(code))
+
+  namespace internal
+  {
+    /**
+     * storage for internal content of the linear solver wrapper
+     */
+    template <typename VectorType>
+    struct LinearSolverContent
+    {
+      ATimesFn a_times_fn;
+      PSetupFn preconditioner_setup;
+      PSolveFn preconditioner_solve;
+
+      LinearSolveFunction<VectorType> lsolve;
+
+      void *P_data;
+      void *A_data;
+    };
+  } // namespace internal
+
+  namespace
+  {
+    using namespace SUNDIALS::internal;
+
+    /**
+     * Access our LinearSolverContent from the generic content of the
+     * SUNLinearSolver @p ls.
+     */
+    template <typename VectorType>
+    LinearSolverContent<VectorType> *
+    access_content(SUNLinearSolver ls)
+    {
+      Assert(ls->content != nullptr, ExcInternalError());
+      return static_cast<LinearSolverContent<VectorType> *>(ls->content);
+    }
+
+
+
+    SUNLinearSolver_Type arkode_linsol_get_type(SUNLinearSolver)
+    {
+      return SUNLINEARSOLVER_ITERATIVE;
+    }
+
+
+
+    template <typename VectorType>
+    int
+    arkode_linsol_solve(SUNLinearSolver LS,
+                        SUNMatrix,
+                        N_Vector x,
+                        N_Vector b,
+                        realtype tol)
+    {
+      auto content = access_content<VectorType>(LS);
+
+      auto *src_b = unwrap_nvector_const<VectorType>(b);
+      auto *dst_x = unwrap_nvector<VectorType>(x);
+
+      SundialsOperator<VectorType> op(content->A_data, content->a_times_fn);
+
+      SundialsPreconditioner<VectorType> preconditioner(
+        content->P_data, content->preconditioner_solve, tol);
+
+      return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
+    }
+
+
+
+    template <typename VectorType>
+    int
+    arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix)
+    {
+      auto content = access_content<VectorType>(LS);
+      if (content->preconditioner_setup)
+        return content->preconditioner_setup(content->P_data);
+      return 0;
+    }
+
+
+
+    template <typename VectorType>
+    int arkode_linsol_initialize(SUNLinearSolver)
+    {
+      // this method is currently only provided because SUNDIALS 4.0.0 requires
+      // it - no user-set action is implemented so far
+      return 0;
+    }
+
+
+
+    template <typename VectorType>
+    int
+    arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
+    {
+      auto content        = access_content<VectorType>(LS);
+      content->A_data     = A_data;
+      content->a_times_fn = ATimes;
+      return 0;
+    }
+
+
+
+    template <typename VectorType>
+    int
+    arkode_linsol_set_preconditioner(SUNLinearSolver LS,
+                                     void *          P_data,
+                                     PSetupFn        p_setup,
+                                     PSolveFn        p_solve)
+    {
+      auto content                  = access_content<VectorType>(LS);
+      content->P_data               = P_data;
+      content->preconditioner_setup = p_setup;
+      content->preconditioner_solve = p_solve;
+      return 0;
+    }
+  } // namespace
+
+  template <typename VectorType>
+  internal::LinearSolverWrapper<VectorType>::LinearSolverWrapper(
+    LinearSolveFunction<VectorType> lsolve)
+    : content(std::make_unique<LinearSolverContent<VectorType>>())
+  {
+    sun_linear_solver                  = SUNLinSolNewEmpty();
+    sun_linear_solver->ops->gettype    = arkode_linsol_get_type;
+    sun_linear_solver->ops->solve      = arkode_linsol_solve<VectorType>;
+    sun_linear_solver->ops->setup      = arkode_linsol_setup<VectorType>;
+    sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
+    sun_linear_solver->ops->setatimes  = arkode_linsol_set_a_times<VectorType>;
+    sun_linear_solver->ops->setpreconditioner =
+      arkode_linsol_set_preconditioner<VectorType>;
+
+    content->lsolve            = lsolve;
+    sun_linear_solver->content = content.get();
+  }
+
+
+
+  template <typename VectorType>
+  internal::LinearSolverWrapper<VectorType>::~LinearSolverWrapper()
+  {
+    SUNLinSolFreeEmpty(sun_linear_solver);
+  }
+
+
+
+  template <typename VectorType>
+  internal::LinearSolverWrapper<VectorType>::operator SUNLinearSolver()
+  {
+    return sun_linear_solver;
+  }
+
+
+
+  template <typename VectorType>
+  SundialsOperator<VectorType>::SundialsOperator(void *   A_data,
+                                                 ATimesFn a_times_fn)
+    : A_data(A_data)
+    , a_times_fn(a_times_fn)
+
+  {
+    Assert(a_times_fn != nullptr, ExcInternalError());
+  }
+
+
+
+  template <typename VectorType>
+  void
+  SundialsOperator<VectorType>::vmult(VectorType &      dst,
+                                      const VectorType &src) const
+  {
+    auto sun_dst = internal::make_nvector_view(dst);
+    auto sun_src = internal::make_nvector_view(src);
+    int  status  = a_times_fn(A_data, sun_src, sun_dst);
+    (void)status;
+    AssertSundialsSolver(status);
+  }
+
+
+
+  template <typename VectorType>
+  SundialsPreconditioner<VectorType>::SundialsPreconditioner(
+    void *   P_data,
+    PSolveFn p_solve_fn,
+    double   tol)
+    : P_data(P_data)
+    , p_solve_fn(p_solve_fn)
+    , tol(tol)
+  {}
+
+
+
+  template <typename VectorType>
+  void
+  SundialsPreconditioner<VectorType>::vmult(VectorType &      dst,
+                                            const VectorType &src) const
+  {
+    // apply identity preconditioner if nothing else specified
+    if (!p_solve_fn)
+      {
+        dst = src;
+        return;
+      }
+
+    auto sun_dst = internal::make_nvector_view(dst);
+    auto sun_src = internal::make_nvector_view(src);
+    // for custom preconditioners no distinction between left and right
+    // preconditioning is made
+    int status =
+      p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
+    (void)status;
+    AssertSundialsSolver(status);
+  }
+
+  template struct SundialsOperator<Vector<double>>;
+  template struct SundialsOperator<BlockVector<double>>;
+  template struct SundialsOperator<LinearAlgebra::distributed::Vector<double>>;
+  template struct SundialsOperator<
+    LinearAlgebra::distributed::BlockVector<double>>;
+
+  template struct SundialsPreconditioner<Vector<double>>;
+  template struct SundialsPreconditioner<BlockVector<double>>;
+  template struct SundialsPreconditioner<
+    LinearAlgebra::distributed::Vector<double>>;
+  template struct SundialsPreconditioner<
+    LinearAlgebra::distributed::BlockVector<double>>;
+
+  template class internal::LinearSolverWrapper<Vector<double>>;
+  template class internal::LinearSolverWrapper<BlockVector<double>>;
+  template class internal::LinearSolverWrapper<
+    LinearAlgebra::distributed::Vector<double>>;
+  template class internal::LinearSolverWrapper<
+    LinearAlgebra::distributed::BlockVector<double>>;
+
+#    ifdef DEAL_II_WITH_MPI
+#      ifdef DEAL_II_WITH_TRILINOS
+  template struct SundialsOperator<TrilinosWrappers::MPI::Vector>;
+  template struct SundialsOperator<TrilinosWrappers::MPI::BlockVector>;
+
+  template struct SundialsPreconditioner<TrilinosWrappers::MPI::Vector>;
+  template struct SundialsPreconditioner<TrilinosWrappers::MPI::BlockVector>;
+
+  template class internal::LinearSolverWrapper<TrilinosWrappers::MPI::Vector>;
+  template class internal::LinearSolverWrapper<
+    TrilinosWrappers::MPI::BlockVector>;
+#      endif // DEAL_II_WITH_TRILINOS
+
+#      ifdef DEAL_II_WITH_PETSC
+#        ifndef PETSC_USE_COMPLEX
+
+  template struct SundialsOperator<PETScWrappers::MPI::Vector>;
+  template struct SundialsOperator<PETScWrappers::MPI::BlockVector>;
+
+  template struct SundialsPreconditioner<PETScWrappers::MPI::Vector>;
+  template struct SundialsPreconditioner<PETScWrappers::MPI::BlockVector>;
+
+  template class internal::LinearSolverWrapper<PETScWrappers::MPI::Vector>;
+  template class internal::LinearSolverWrapper<PETScWrappers::MPI::BlockVector>;
+#        endif // PETSC_USE_COMPLEX
+#      endif   // DEAL_II_WITH_PETSC
+
+#    endif // DEAL_II_WITH_MPI
+} // namespace SUNDIALS
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif
+#endif

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.