# 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))
*/
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).
*
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
# 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.
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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
copy.cc
kinsol.cc
n_vector.cc
+ sunlinsol_wrapper.cc
)
SET(_inst
# 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>
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,
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
{
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
{
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
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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