SET(${var} TRUE)
#
- # We require at least sundials 3.0.0
+ # We require at least sundials 5.4.0
#
- SET(_version_required 3.0.0)
+ SET(_version_required 5.4.0)
IF(SUNDIALS_VERSION VERSION_LESS ${_version_required})
MESSAGE(STATUS "Could not find a sufficient Sundials installation: "
"deal.II requires at least version ${_version_required}, "
# include <deal.II/lac/vector_memory.h>
# include <arkode/arkode.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-# include <arkode/arkode_impl.h>
-# endif
# include <nvector/nvector_serial.h>
# ifdef DEAL_II_WITH_MPI
# include <nvector/nvector_parallel.h>
std::function<int(const double t, const VectorType &y, VectorType &res)>
implicit_function;
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- /**
- * A function object that users may supply and that is intended to
- * prepare the linear solver for subsequent calls to
- * solve_jacobian_system().
- *
- * Make sure that after a call to this function, we know how to compute
- * solutions of systems $A x = b$, where $A$ is some approximation to the
- * Newton matrix, $M - \gamma \partial f_I/\partial y$. This function is
- * optional. If the user does not provide it, then solve_jacobian_system()
- * is assumed to also perform the setup internally.
- *
- * The setup_jacobian() function may call a user-supplied function to
- * compute needed data related to the Jacobian matrix. Alternatively, it may
- * choose to retrieve and use stored values of this data. In either case,
- * setup_jacobian() may also preprocess that data as needed for
- * solve_jacobian_system(), which may involve calling a generic function
- * (such as for LU factorization).
- *
- * This data may be intended either for direct use (in a direct linear
- * solver) or for use in a preconditioner (in a preconditioned iterative
- * linear solver). The setup_jacobian() function is not called at every
- * stage solve (or even every time step), but only as frequently as the
- * solver determines that it is appropriate to perform the setup task. In
- * this way, Jacobian-related data generated by setup_jacobian() is
- * expected to be used over a number of time steps.
- *
- * If the user uses a matrix based computation of the Jacobian, then this
- * is the right place where an assembly routine shoulde be called to
- * assemble both a matrix and a preconditioner for the Jacobian system.
- * Subsequent calls (possibly more than one) to solve_jacobian_system() can
- * assume that this function has been called at least once.
- *
- * Notice that no assumption is made by this interface on what the user
- * should do in this function. ARKode only assumes that after a call to
- * setup_jacobian() it is possible to call solve_jacobian_system(), to
- * obtain a solution $x$ to the system $J x = b$. If this function is not
- * provided, then it is never called.
- *
- * Arguments to the function are
- *
- * @param[in] t the current time
- * @param[in] gamma the current factor to use in the jacobian computation
- * @param[in] ypred is the predicted $y$ vector for the current ARKode
- * internal step
- * @param[in] fpred is the value of the implicit right-hand side at ypred,
- * $f_I (t_n, ypred)$.
- *
- * @param[in] convfail Input flag used to indicate any problem that
- * occurred during the solution of the nonlinear equation on the current
- * time step for which the linear solver is being used. This flag can be
- * used to help decide whether the Jacobian data kept by a linear solver
- * needs to be updated or not. Its possible values are:
- *
- * - ARK_NO_FAILURES: this value is passed if either this is the first
- * call for this step, or the local error test failed on the previous
- * attempt at this step (but the Newton iteration converged).
- *
- * - ARK_FAIL_BAD_J: this value is passed if (a) the previous Newton
- * corrector iteration did not converge and the linear solver's setup
- * function indicated that its Jacobian-related data is not current, or
- * (b) during the previous Newton corrector iteration, the linear solver's
- * solve function failed in a recoverable manner and the linear solver's
- * setup function indicated that its Jacobian-related data is not
- * current.
- *
- * - ARK_FAIL_OTHER: this value is passed if during the current internal
- * step try, the previous Newton iteration failed to converge even
- * though the linear solver was using current Jacobian-related data.
- *
- * @param[out] j_is_current: a boolean to be filled in by setup_jacobian().
- * The value should be set to `true` if the Jacobian data is current after
- * the call, and should be set to `false` if its Jacobian data is not
- * current. If setup_jacobian() calls for re-evaluation of Jacobian data
- * (based on convfail and ARKode state data), then it should set
- * `j_is_current` to `true` unconditionally, otherwise an infinite loop can
- * result.
- *
- * This function should return:
- * - 0: Success
- * - >0: Recoverable error (ARKodeReinit will be called if this happens, and
- * then last function will be attempted again
- * - <0: Unrecoverable error the computation will be aborted and an
- * assertion will be thrown.
- */
- std::function<int(const int convfail,
- const double t,
- const double gamma,
- const VectorType &ypred,
- const VectorType &fpred,
- bool & j_is_current)>
- setup_jacobian;
-
- /**
- * A function object that users may supply and that is intended to solve
- * the Jacobian linear system. This function will be called by ARKode
- * (possibly several times) after setup_jacobian() has been called at least
- * once. ARKode tries to do its best to call setup_jacobian() the minimum
- * amount of times. If convergence can be achieved without updating the
- * Jacobian, then ARKode does not call setup_jacobian() again. If, on the
- * contrary, internal ARKode convergence tests fail, then ARKode calls
- * again setup_jacobian() with updated vectors and coefficients so that
- * successive calls to solve_jacobian_systems() lead to better convergence
- * in the Newton process.
- *
- * If you do not specify a solve_jacobian_system() function, then a fixed
- * point iteration is used instead of a Newton method. Notice that this may
- * not converge, or may converge very slowly.
- *
- * The jacobian $J$ should be (an approximation of) the system Jacobian
- * \f[
- * J = M - \gamma \frac{\partial f_I}{\partial y}
- * \f]
- * evaluated at `t`, `ycur`. `fcur` is $f_I(t,ycur)$.
- *
- * A call to this function should store in `dst` the result of $J^{-1}$
- * applied to `src`, i.e., `J*dst = src`. It is the users responsibility to
- * set up proper solvers and preconditioners inside this function.
- *
- *
- * Arguments to the function are
- *
- * @param[in] t the current time
- * @param[in] gamma the current factor to use in the jacobian computation
- * @param[in] ycur is the current $y$ vector for the current ARKode
- * internal step
- * @param[in] fcur is the current value of the implicit right-hand side at
- * ycur, $f_I (t_n, ypred)$.
- *
- *
- * This function should return:
- * - 0: Success
- * - >0: Recoverable error (ARKodeReinit will be called if this happens, and
- * then last function will be attempted again
- * - <0: Unrecoverable error the computation will be aborted and an
- * assertion will be thrown.
- */
- std::function<int(const double t,
- const double gamma,
- const VectorType &ycur,
- const VectorType &fcur,
- const VectorType &rhs,
- VectorType & dst)>
- solve_jacobian_system;
-
-
- /**
- * A function object that users may supply and that is intended to set up
- * the mass matrix. This function is called by ARKode any time a mass
- * matrix update is required. The user should compute the mass matrix (or
- * update all the variables that allow the application of the mass matrix).
- * This function is called by ARKode once, before any call to
- * solve_mass_system().
- *
- * ARKode supports the case where the mass matrix may depend on time, but
- * not the case where the mass matrix depends on the solution itself.
- *
- * If the user does not provide a solve_mass_matrix() function, then the
- * identity is used. If the setup_mass() function is not provided, then
- * solve_mass_system() should do all the work by itself.
- *
- * If the user uses a matrix based computation of the mass matrix, then
- * this is the right place where an assembly routine shoulde be called to
- * assemble both a matrix and a preconditioner. Subsequent calls (possibly
- * more than one) to solve_mass_system() can assume that this function
- * has been called at least once.
- *
- * Notice that no assumption is made by this interface on what the user
- * should do in this function. ARKode only assumes that after a call to
- * setup_mass() it is possible to call solve_mass_system(), to
- * obtain a solution $x$ to the system $M x = b$.
- *
- * This function should return:
- * - 0: Success
- * - >0: Recoverable error (ARKodeReinit will be called if this happens, and
- * then last function will be attempted again
- * - <0: Unrecoverable error the computation will be aborted and an
- * assertion will be thrown.
- */
- std::function<int(const double t)> setup_mass;
-
- /**
- * A function object that users may supply and that is intended to solve
- * the mass matrix linear system. This function will be called by ARKode
- * (possibly several times) after setup_mass() has been called at least
- * once. ARKode tries to do its best to call setup_mass() the minimum
- * amount of times.
- *
- * A call to this function should store in `dst` the result of $M^{-1}$
- * applied to `src`, i.e., `M*dst = src`. It is the users responsibility to
- * set up proper solvers and preconditioners inside this function.
- *
- * This function should return:
- * - 0: Success
- * - >0: Recoverable error (ARKodeReinit will be called if this happens, and
- * then last function will be attempted again
- * - <0: Unrecoverable error the computation will be aborted and an
- * assertion will be thrown.
- */
- std::function<int(const VectorType &rhs, VectorType &dst)>
- solve_mass_system;
-# else
-
/**
* A function object that users may supply and that is intended to compute
* the product of the mass matrix with a given vector `v`. This function
* assertion will be thrown.
*/
std::function<int(double t)> mass_preconditioner_setup;
-# endif
/**
* A function object that users may supply and that is intended to
dealii::DiscreteTime &time,
const bool do_reset);
-# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
-
/**
* Set up the (non)linear solver and preconditioners in the ARKODE memory
* object based on the user-specified functions.
void
setup_mass_solver(const VectorType &solution);
-# endif
-
/**
* This function is executed at construction time to set the
* std::function above to trigger an assert if they are not
*/
double last_end_time;
-# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
-# endif
# ifdef DEAL_II_WITH_PETSC
# ifdef PETSC_USE_COMPLEX
# include <ida/ida.h>
# endif
-# include <sundials/sundials_config.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
-# include <ida/ida_spbcgs.h>
-# include <ida/ida_spgmr.h>
-# include <ida/ida_sptfqmr.h>
-# endif
# include <deal.II/sundials/sunlinsol_wrapper.h>
# include <boost/signals2.hpp>
# include <nvector/nvector_serial.h>
+# include <sundials/sundials_config.h>
# include <sundials/sundials_math.h>
# include <sundials/sundials_types.h>
# include <boost/signals2.hpp>
# include <kinsol/kinsol.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
-# include <kinsol/kinsol_impl.h>
-# endif
# include <nvector/nvector_serial.h>
# include <sundials/sundials_math.h>
# include <sundials/sundials_types.h>
# include <deal.II/lac/trilinos_vector.h>
# include <deal.II/lac/vector_memory.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-# include <deal.II/sundials/sundials_backport.h>
-# endif
-
DEAL_II_NAMESPACE_OPEN
namespace SUNDIALS
v->ops->nvcloneempty = &NVectorOperations::clone_empty;
v->ops->nvdestroy = &NVectorOperations::destroy<VectorType>;
// v->ops->nvspace = undef;
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
v->ops->nvgetcommunicator =
&NVectorOperations::get_communicator_as_void_ptr<VectorType>;
v->ops->nvgetlength = &NVectorOperations::get_global_length<VectorType>;
-# endif
/* standard vector operations */
v->ops->nvlinearsum = &NVectorOperations::linear_sum<VectorType>;
+++ /dev/null
-//-----------------------------------------------------------
-//
-// Copyright (C) 2020 - 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.
-//
-//-----------------------------------------------------------
-
-/*
- * The functions in this file are based on an implementation distributed within
- * the SUNDIALS package, see the license here:
- * https://computing.llnl.gov/projects/sundials/license.
- * -----------------------------------------------------------------
- * Programmer(s): Daniel Reynolds @ SMU
- * David J. Gardner, Carol S. Woodward, and
- * Slaven Peles @ LLNL
- * -----------------------------------------------------------------*/
-
-#ifndef dealii_sundials_backport_h
-#define dealii_sundials_backport_h
-
-#include <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_SUNDIALS
-
-# include <sundials/sundials_nvector.h>
-
-DEAL_II_NAMESPACE_OPEN
-namespace SUNDIALS
-{
- namespace internal
- {
- N_Vector
- N_VNewEmpty()
- {
- N_Vector v = new _generic_N_Vector;
- N_Vector_Ops ops = new _generic_N_Vector_Ops;
-
- /* initialize operations to nullptr */
-
- /* constructors, destructors, and utility operations */
- ops->nvgetvectorid = nullptr;
- ops->nvclone = nullptr;
- ops->nvcloneempty = nullptr;
- ops->nvdestroy = nullptr;
- ops->nvspace = nullptr;
- ops->nvgetarraypointer = nullptr;
- ops->nvsetarraypointer = nullptr;
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
- ops->nvgetcommunicator = nullptr;
- ops->nvgetlength = nullptr;
-# endif
-
- /* standard vector operations */
- ops->nvlinearsum = nullptr;
- ops->nvconst = nullptr;
- ops->nvprod = nullptr;
- ops->nvdiv = nullptr;
- ops->nvscale = nullptr;
- ops->nvabs = nullptr;
- ops->nvinv = nullptr;
- ops->nvaddconst = nullptr;
- ops->nvdotprod = nullptr;
- ops->nvmaxnorm = nullptr;
- ops->nvwrmsnormmask = nullptr;
- ops->nvwrmsnorm = nullptr;
- ops->nvmin = nullptr;
- ops->nvwl2norm = nullptr;
- ops->nvl1norm = nullptr;
- ops->nvcompare = nullptr;
- ops->nvinvtest = nullptr;
- ops->nvconstrmask = nullptr;
- ops->nvminquotient = nullptr;
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
- /* fused vector operations (optional) */
- ops->nvlinearcombination = nullptr;
- ops->nvscaleaddmulti = nullptr;
- ops->nvdotprodmulti = nullptr;
-
- /* vector array operations (optional) */
- ops->nvlinearsumvectorarray = nullptr;
- ops->nvscalevectorarray = nullptr;
- ops->nvconstvectorarray = nullptr;
- ops->nvwrmsnormvectorarray = nullptr;
- ops->nvwrmsnormmaskvectorarray = nullptr;
- ops->nvscaleaddmultivectorarray = nullptr;
- ops->nvlinearcombinationvectorarray = nullptr;
-
- /* local reduction operations (optional) */
- ops->nvdotprodlocal = nullptr;
- ops->nvmaxnormlocal = nullptr;
- ops->nvminlocal = nullptr;
- ops->nvl1normlocal = nullptr;
- ops->nvinvtestlocal = nullptr;
- ops->nvconstrmasklocal = nullptr;
- ops->nvminquotientlocal = nullptr;
- ops->nvwsqrsumlocal = nullptr;
- ops->nvwsqrsummasklocal = nullptr;
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0)
- /* XBraid interface operations */
- ops->nvbufsize = nullptr;
- ops->nvbufpack = nullptr;
- ops->nvbufunpack = nullptr;
-# endif
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 3, 0)
- /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */
- ops->nvprint = nullptr;
- ops->nvprintfile = nullptr;
-# endif
-# endif
-
- /* attach ops and initialize content to nullptr */
- v->ops = ops;
- v->content = nullptr;
-
- return v;
- }
-
-
-
- void
- N_VFreeEmpty(N_Vector v)
- {
- if (v == nullptr)
- return;
-
- /* free non-nullptr ops structure */
- if (v->ops)
- delete v->ops;
- v->ops = nullptr;
-
- /* free overall N_Vector object and return */
- delete v;
- }
-
-
-
- int
- N_VCopyOps(N_Vector w, N_Vector v)
- {
- /* Check that ops structures exist */
- if (w == nullptr || v == nullptr)
- return (-1);
- if (w->ops == nullptr || v->ops == nullptr)
- return (-1);
-
- /* Copy ops from w to v */
-
- /* constructors, destructors, and utility operations */
- v->ops->nvgetvectorid = w->ops->nvgetvectorid;
- v->ops->nvclone = w->ops->nvclone;
- v->ops->nvcloneempty = w->ops->nvcloneempty;
- v->ops->nvdestroy = w->ops->nvdestroy;
- v->ops->nvspace = w->ops->nvspace;
- v->ops->nvgetarraypointer = w->ops->nvgetarraypointer;
- v->ops->nvsetarraypointer = w->ops->nvsetarraypointer;
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
- v->ops->nvgetcommunicator = w->ops->nvgetcommunicator;
- v->ops->nvgetlength = w->ops->nvgetlength;
-# endif
-
- /* standard vector operations */
- v->ops->nvlinearsum = w->ops->nvlinearsum;
- v->ops->nvconst = w->ops->nvconst;
- v->ops->nvprod = w->ops->nvprod;
- v->ops->nvdiv = w->ops->nvdiv;
- v->ops->nvscale = w->ops->nvscale;
- v->ops->nvabs = w->ops->nvabs;
- v->ops->nvinv = w->ops->nvinv;
- v->ops->nvaddconst = w->ops->nvaddconst;
- v->ops->nvdotprod = w->ops->nvdotprod;
- v->ops->nvmaxnorm = w->ops->nvmaxnorm;
- v->ops->nvwrmsnormmask = w->ops->nvwrmsnormmask;
- v->ops->nvwrmsnorm = w->ops->nvwrmsnorm;
- v->ops->nvmin = w->ops->nvmin;
- v->ops->nvwl2norm = w->ops->nvwl2norm;
- v->ops->nvl1norm = w->ops->nvl1norm;
- v->ops->nvcompare = w->ops->nvcompare;
- v->ops->nvinvtest = w->ops->nvinvtest;
- v->ops->nvconstrmask = w->ops->nvconstrmask;
- v->ops->nvminquotient = w->ops->nvminquotient;
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
- /* fused vector operations */
- v->ops->nvlinearcombination = w->ops->nvlinearcombination;
- v->ops->nvscaleaddmulti = w->ops->nvscaleaddmulti;
- v->ops->nvdotprodmulti = w->ops->nvdotprodmulti;
-
- /* vector array operations */
- v->ops->nvlinearsumvectorarray = w->ops->nvlinearsumvectorarray;
- v->ops->nvscalevectorarray = w->ops->nvscalevectorarray;
- v->ops->nvconstvectorarray = w->ops->nvconstvectorarray;
- v->ops->nvwrmsnormvectorarray = w->ops->nvwrmsnormvectorarray;
- v->ops->nvwrmsnormmaskvectorarray = w->ops->nvwrmsnormmaskvectorarray;
- v->ops->nvscaleaddmultivectorarray = w->ops->nvscaleaddmultivectorarray;
- v->ops->nvlinearcombinationvectorarray =
- w->ops->nvlinearcombinationvectorarray;
-
- /* local reduction operations */
- v->ops->nvdotprodlocal = w->ops->nvdotprodlocal;
- v->ops->nvmaxnormlocal = w->ops->nvmaxnormlocal;
- v->ops->nvminlocal = w->ops->nvminlocal;
- v->ops->nvl1normlocal = w->ops->nvl1normlocal;
- v->ops->nvinvtestlocal = w->ops->nvinvtestlocal;
- v->ops->nvconstrmasklocal = w->ops->nvconstrmasklocal;
- v->ops->nvminquotientlocal = w->ops->nvminquotientlocal;
- v->ops->nvwsqrsumlocal = w->ops->nvwsqrsumlocal;
- v->ops->nvwsqrsummasklocal = w->ops->nvwsqrsummasklocal;
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0)
- /* XBraid interface operations */
- v->ops->nvbufsize = w->ops->nvbufsize;
- v->ops->nvbufpack = w->ops->nvbufpack;
- v->ops->nvbufunpack = w->ops->nvbufunpack;
-# endif
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(5, 3, 0)
- /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */
- v->ops->nvprint = w->ops->nvprint;
- v->ops->nvprintfile = w->ops->nvprintfile;
-# endif
-# endif
-
- return (0);
- }
- } // namespace internal
-} // namespace SUNDIALS
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_WITH_SUNDIALS
-#endif // dealii_sundials_sunlinsol_newempty_h
+++ /dev/null
-//-----------------------------------------------------------
-//
-// Copyright (C) 2020 - 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.
-//
-//-----------------------------------------------------------
-
-/*
- * The functions in this file are based on an implementation distributed within
- * the SUNDIALS package, see the license here:
- * https://computing.llnl.gov/projects/sundials/license.
- * -----------------------------------------------------------------
- * Programmer(s): Daniel Reynolds @ SMU
- * David J. Gardner, Carol S. Woodward, and
- * Slaven Peles @ LLNL
- * -----------------------------------------------------------------*/
-
-#ifndef dealii_sundials_sunlinsol_newempty_h
-#define dealii_sundials_sunlinsol_newempty_h
-
-#include <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_SUNDIALS
-
-# include <sundials/sundials_linearsolver.h>
-
-DEAL_II_NAMESPACE_OPEN
-namespace SUNDIALS
-{
- namespace internal
- {
- /**
- * Create a new SUNLinearSolver structure without any content and
- * operations set to `nullptr`.
- */
- inline SUNLinearSolver
- SUNLinSolNewEmpty()
- {
- /* create linear solver object */
- SUNLinearSolver LS = new _generic_SUNLinearSolver;
-
- /* create linear solver ops structure */
- SUNLinearSolver_Ops ops = new _generic_SUNLinearSolver_Ops;
-
- /* initialize operations to nullptr */
- ops->gettype = nullptr;
- ops->setatimes = nullptr;
- ops->setpreconditioner = nullptr;
- ops->setscalingvectors = nullptr;
- ops->initialize = nullptr;
- ops->setup = nullptr;
- ops->solve = nullptr;
- ops->numiters = nullptr;
- ops->resnorm = nullptr;
- ops->resid = nullptr;
- ops->lastflag = nullptr;
- ops->space = nullptr;
- ops->free = nullptr;
-
- /* attach ops and initialize content to nullptr */
- LS->ops = ops;
- LS->content = nullptr;
-
- return (LS);
- }
-
- /**
- * Free the memory associated with @p solver which was previously allocated
- * with a call to SUNLinSolNewEmpty().
- *
- * @note A call to this function does not deallocate the `content` field.
- *
- * @param solver The solver memory to free
- */
- inline void
- SUNLinSolFreeEmpty(SUNLinearSolver solver)
- {
- if (solver == nullptr)
- return;
-
- /* free non-nullptr ops structure */
- if (solver->ops)
- delete solver->ops;
- solver->ops = nullptr;
-
- /* free overall linear solver object */
- delete solver;
- return;
- }
-
- } // namespace internal
-} // namespace SUNDIALS
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_WITH_SUNDIALS
-#endif // dealii_sundials_sunlinsol_newempty_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>
+# include <sundials/sundials_linearsolver.h>
-# include <functional>
-# include <memory>
+# include <functional>
+# include <memory>
DEAL_II_NAMESPACE_OPEN
namespace SUNDIALS
{
-# ifndef DOXYGEN
+# ifndef DOXYGEN
// forward declarations
namespace internal
{
template <typename VectorType>
struct LinearSolverContent;
}
-# endif
+# endif
/**
* A linear operator that wraps SUNDIALS functionality.
{
public:
explicit LinearSolverWrapper(LinearSolveFunction<VectorType> lsolve
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
SUNContext &linsol_ctx
-# endif
+# endif
);
~LinearSolverWrapper();
DEAL_II_NAMESPACE_CLOSE
-# endif
#endif
#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>
-# include <sundials/sundials_config.h>
-# else
-# include <arkode/arkode_arkstep.h>
-# include <sunlinsol/sunlinsol_spgmr.h>
-# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-# include <deal.II/sundials/sunlinsol_newempty.h>
-# endif
-# endif
+# include <arkode/arkode_arkstep.h>
+# include <sunlinsol/sunlinsol_spgmr.h>
+# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
# include <iostream>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-// Make sure we know how to call sundials own ARKode() function
-const auto &SundialsARKode = ARKode;
-# endif
-
DEAL_II_NAMESPACE_OPEN
namespace SUNDIALS
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- template <typename VectorType>
- int
- setup_jacobian_callback(ARKodeMem arkode_mem,
- int convfail,
- N_Vector ypred,
- N_Vector fpred,
- booleantype *jcurPtr,
- N_Vector,
- N_Vector,
- N_Vector)
- {
- Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
- ARKode<VectorType> &solver =
- *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
- auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
- auto *src_fpred = internal::unwrap_nvector_const<VectorType>(fpred);
- // avoid reinterpret_cast
- bool jcurPtr_tmp = false;
- int err = solver.setup_jacobian(convfail,
- arkode_mem->ark_tn,
- arkode_mem->ark_gamma,
- *src_ypred,
- *src_fpred,
- jcurPtr_tmp);
- *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE;
-
- return err;
- }
-
-
-
- template <typename VectorType>
- int
- solve_with_jacobian_callback(ARKodeMem arkode_mem,
- N_Vector b,
- N_Vector ycur,
- N_Vector fcur)
- {
- Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
- ARKode<VectorType> &solver =
- *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
- auto *dst = internal::unwrap_nvector<VectorType>(b);
- auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
- auto *src_fcur = internal::unwrap_nvector_const<VectorType>(fcur);
-
- // make a temporary copy to work on in the user call
- VectorType src = *dst;
-
- int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
- arkode_mem->ark_gamma,
- *src_ycur,
- *src_fcur,
- src,
- *dst);
-
-
- return err;
- }
-
-
-
- template <typename VectorType>
- int
- setup_mass_matrix_callback(ARKodeMem arkode_mem,
- N_Vector,
- N_Vector,
- N_Vector)
- {
- Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
- ARKode<VectorType> &solver =
- *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
- int err = solver.setup_mass(arkode_mem->ark_tn);
- return err;
- }
-
-
-
- template <typename VectorType>
- int
- solve_with_mass_matrix_callback(ARKodeMem arkode_mem, N_Vector b)
- {
- Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
- ARKode<VectorType> &solver =
- *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
- auto *dst = internal::unwrap_nvector<VectorType>(b);
-
- // make a temporary copy to work on in the user call
- VectorType src = *dst;
-
- int err = solver.solve_mass_system(src, *dst);
-
- return err;
- }
-# else
-
template <typename VectorType>
int
jacobian_times_vector_callback(N_Vector v,
return solver.mass_preconditioner_setup(t);
}
-# endif
-
} // namespace
{
if (arkode_mem)
{
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- ARKodeFree(&arkode_mem);
-# else
ARKStepFree(&arkode_mem);
-# endif
# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
const int status = SUNContext_Free(&arkode_ctx);
time.set_desired_next_step_size(data.output_period);
// Let ARKode advance time by one period:
- double actual_next_time;
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- const auto status = SundialsARKode(arkode_mem,
- time.get_next_time(),
- solution_nvector,
- &actual_next_time,
- ARK_NORMAL);
-# else
+ double actual_next_time;
const auto status = ARKStepEvolve(arkode_mem,
time.get_next_time(),
solution_nvector,
&actual_next_time,
ARK_NORMAL);
-# endif
(void)status;
AssertARKode(status);
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
template <typename VectorType>
void
ARKode<VectorType>::reset(const double current_time,
const VectorType &solution)
{
last_end_time = current_time;
- if (arkode_mem)
- ARKodeFree(&arkode_mem);
-
- arkode_mem = ARKodeCreate();
-
int status;
(void)status;
- Assert(explicit_function || implicit_function,
- ExcFunctionNotProvided("explicit_function || implicit_function"));
-
- // just a view on the memory in solution, all write operations on yy by
- // ARKODE will automatically be mirrored to solution
- auto initial_condition_nvector = internal::make_nvector_view(solution
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- arkode_ctx
-# endif
- );
-
- status = ARKodeInit(
- arkode_mem,
- explicit_function ? &explicit_function_callback<VectorType> : nullptr,
- implicit_function ? &implicit_function_callback<VectorType> : nullptr,
- current_time,
- initial_condition_nvector);
- AssertARKode(status);
-
- if (get_local_tolerances)
- {
- const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- arkode_ctx
-# endif
- );
- status =
- ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
- AssertARKode(status);
- }
- else
- {
- status = ARKodeSStolerances(arkode_mem,
- data.relative_tolerance,
- data.absolute_tolerance);
- AssertARKode(status);
- }
-
- status = ARKodeSetInitStep(arkode_mem, current_time_step);
- AssertARKode(status);
-
- status = ARKodeSetUserData(arkode_mem, this);
- AssertARKode(status);
-
- status = ARKodeSetStopTime(arkode_mem, data.final_time);
- AssertARKode(status);
-
- status =
- ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
- AssertARKode(status);
-
- // Initialize solver
- auto ARKode_mem = static_cast<ARKodeMem>(arkode_mem);
-
- if (solve_jacobian_system)
- {
- status = ARKodeSetNewton(arkode_mem);
- AssertARKode(status);
- if (data.implicit_function_is_linear)
- {
- status = ARKodeSetLinear(
- arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
- AssertARKode(status);
- }
-
-
- ARKode_mem->ark_lsolve = solve_with_jacobian_callback<VectorType>;
- if (setup_jacobian)
- {
- ARKode_mem->ark_lsetup = setup_jacobian_callback<VectorType>;
- }
- }
- else
- {
- status =
- ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations);
- AssertARKode(status);
- }
-
-
- if (solve_mass_system)
- {
- ARKode_mem->ark_msolve = solve_with_mass_matrix_callback<VectorType>;
-
- if (setup_mass)
- {
- ARKode_mem->ark_msetup = setup_mass_matrix_callback<VectorType>;
- }
- }
-
- status = ARKodeSetOrder(arkode_mem, data.maximum_order);
- AssertARKode(status);
-
- if (custom_setup)
- custom_setup(arkode_mem);
- }
-
-# else
-
- template <typename VectorType>
- void
- ARKode<VectorType>::reset(const double current_time,
- const double current_time_step,
- const VectorType &solution)
- {
- last_end_time = current_time;
- int status;
- (void)status;
-
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
if (arkode_ctx)
{
status = SUNContext_Free(&arkode_ctx);
}
status = SUNContext_Create(&mpi_communicator, &arkode_ctx);
AssertARKode(status);
-# endif
+# endif
if (arkode_mem)
{
// just a view on the memory in solution, all write operations on yy by
// ARKODE will automatically be mirrored to solution
auto initial_condition_nvector = internal::make_nvector_view(solution
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
Assert(explicit_function || implicit_function,
implicit_function ? &implicit_function_callback<VectorType> : nullptr,
current_time,
initial_condition_nvector
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
Assert(arkode_mem != nullptr, ExcInternalError());
if (get_local_tolerances)
{
const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
status =
ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
linear_solver =
std::make_unique<internal::LinearSolverWrapper<VectorType>>(
solve_linearized_system
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
sun_linear_solver = *linear_solver;
}
// use default solver from SUNDIALS
// TODO give user options
auto y_template = internal::make_nvector_view(solution
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
sun_linear_solver =
SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/);
-# else
+# else
sun_linear_solver =
SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/,
arkode_ctx);
-# endif
+# endif
}
status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr);
AssertARKode(status);
else
{
auto y_template = internal::make_nvector_view(solution
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
SUNNonlinearSolver fixed_point_solver =
SUNNonlinSol_FixedPoint(y_template,
data.anderson_acceleration_subspace);
-# else
+# else
SUNNonlinearSolver fixed_point_solver =
SUNNonlinSol_FixedPoint(y_template,
data.anderson_acceleration_subspace,
arkode_ctx);
-# endif
+# endif
status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
AssertARKode(status);
mass_solver =
std::make_unique<internal::LinearSolverWrapper<VectorType>>(
solve_mass
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
sun_mass_linear_solver = *mass_solver;
}
else
{
auto y_template = internal::make_nvector_view(solution
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
arkode_ctx
-# endif
+# endif
);
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
sun_mass_linear_solver =
SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/);
-# else
+# else
sun_mass_linear_solver =
SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/,
arkode_ctx);
-# endif
+# endif
}
booleantype mass_time_dependent =
data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
}
}
}
-# endif
# include <deal.II/sundials/n_vector.h>
-# include <sundials/sundials_config.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-# ifdef DEAL_II_SUNDIALS_WITH_IDAS
-# include <idas/idas_impl.h>
-# else
-# include <ida/ida_impl.h>
-# endif
-# endif
-# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-# include <deal.II/sundials/sunlinsol_newempty.h>
-# endif
# include <iomanip>
# include <iostream>
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- template <typename VectorType>
- int
- t_dae_lsetup(IDAMem IDA_mem,
- N_Vector yy,
- N_Vector yp,
- N_Vector resp,
- N_Vector tmp1,
- N_Vector tmp2,
- N_Vector tmp3)
- {
- (void)tmp1;
- (void)tmp2;
- (void)tmp3;
- (void)resp;
- IDA<VectorType> &solver =
- *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
-
- auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
- auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
-
- int err = solver.setup_jacobian(IDA_mem->ida_tn,
- *src_yy,
- *src_yp,
- IDA_mem->ida_cj);
-
- return err;
- }
-
-
-
- template <typename VectorType>
- int
- solve_with_jacobian_callback(IDAMem IDA_mem,
- N_Vector b,
- N_Vector weight,
- N_Vector yy,
- N_Vector yp,
- N_Vector resp)
- {
- (void)weight;
- (void)yy;
- (void)yp;
- (void)resp;
- IDA<VectorType> &solver =
- *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer dst(mem);
- solver.reinit_vector(*dst);
-
- auto *src = internal::unwrap_nvector<VectorType>(b);
-
- int err = solver.solve_jacobian_system(*src, *dst);
- *src = *dst;
-
- return err;
- }
-
-
-
-# else
template <typename VectorType>
int
setup_jacobian_callback(realtype tt,
return err;
}
-# endif
} // namespace
status = SUNContext_Create(&mpi_communicator, &ida_ctx);
AssertIDA(status);
- ida_mem = IDACreate(ida_ctx);
+ ida_mem = IDACreate(ida_ctx);
# endif
auto yy = internal::make_nvector_view(solution
AssertIDA(status);
// Initialize solver
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- auto IDA_mem = static_cast<IDAMem>(ida_mem);
-
- IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
-
- if (solve_jacobian_system)
- IDA_mem->ida_lsolve = solve_with_jacobian_callback<VectorType>;
- else
- AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
-# else
SUNMatrix J = nullptr;
SUNLinearSolver LS = nullptr;
// called do not actually receive the IDAMEM object, just the LS
// object, so we have to store a pointer to the current
// object in the LS object
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
LS = SUNLinSolNewEmpty();
-# else
- LS = SUNLinSolNewEmpty(ida_ctx);
-# endif
+# else
+ LS = SUNLinSolNewEmpty(ida_ctx);
+# endif
LS->content = this;
// if we don't set it, it won't call the functions that set up
// the matrix object (i.e., the argument to the 'IDASetJacFn'
// function below).
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
- J = SUNMatNewEmpty();
-# else
- J = SUNMatNewEmpty(ida_ctx);
-# endif
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+ J = SUNMatNewEmpty();
+# else
+ J = SUNMatNewEmpty(ida_ctx);
+# endif
J->content = this;
J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
// calling IDASetLinearSolver
status = IDASetJacFn(ida_mem, &setup_jacobian_callback<VectorType>);
AssertIDA(status);
-# endif
status = IDASetMaxOrd(ida_mem, data.maximum_order);
AssertIDA(status);
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
- template <typename VectorType>
- int
- setup_jacobian_callback(KINMem kinsol_mem)
- {
- KINSOL<VectorType> &solver =
- *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-
- auto *src_ycur =
- internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
- auto *src_fcur =
- internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
-
- int err = solver.setup_jacobian(*src_ycur, *src_fcur);
- return err;
- }
-
-
-
- template <typename VectorType>
- int
- solve_with_jacobian_callback(KINMem kinsol_mem,
- N_Vector x,
- N_Vector b,
- realtype *sJpnorm,
- realtype *sFdotJp)
- {
- KINSOL<VectorType> &solver =
- *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-
- auto *src_ycur =
- internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
- auto *src_fcur =
- internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
- auto *src = internal::unwrap_nvector_const<VectorType>(b);
- auto *dst = internal::unwrap_nvector<VectorType>(x);
-
- int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst);
-
- *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale);
- N_VProd(b, kinsol_mem->kin_fscale, b);
- N_VProd(b, kinsol_mem->kin_fscale, b);
- *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b);
-
- return err;
- }
-
-# else // SUNDIALS 5.0 or later
-
template <typename VectorType>
int
setup_jacobian_callback(N_Vector u,
return err;
}
}
-# endif
} // namespace
solve_with_jacobian) // user assigned a function object to the solver
// slot
{
-/* interface up to and including 4.0 */
-# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
- auto KIN_mem = static_cast<KINMem>(kinsol_mem);
- // Old version only works with solve_jacobian_system
- Assert(solve_jacobian_system,
- ExcFunctionNotProvided("solve_jacobian_system"))
- KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
- if (setup_jacobian) // user assigned a function object to the Jacobian
- // set-up slot
- KIN_mem->kin_lsetup = setup_jacobian_callback<VectorType>;
-
-/* interface up to and including 4.1 */
-# elif DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-
- // deal.II does not currently have support for KINSOL in
- // SUNDIALS 4.1. One could write this and update this section,
- // but it does not seem worthwhile spending the time to
- // interface with an old version of SUNDIAL given that the
- // code below supports modern SUNDIAL versions just fine.
- Assert(false, ExcNotImplemented());
-
-# else /* interface starting with SUNDIALS 5.0 */
// Set the operations we care for in the sun_linear_solver object
// and attach it to the KINSOL object. The functions that will get
// called do not actually receive the KINSOL object, just the LS
// object, so we have to store a pointer to the current
// object in the LS object
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
- LS = SUNLinSolNewEmpty();
-# else
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+ LS = SUNLinSolNewEmpty();
+# else
LS = SUNLinSolNewEmpty(kinsol_ctx);
-# endif
+# endif
LS->content = this;
LS->ops->gettype =
// if we don't set it, it won't call the functions that set up
// the matrix object (i.e., the argument to the 'KINSetJacFn'
// function below).
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
- J = SUNMatNewEmpty();
-# else
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+ J = SUNMatNewEmpty();
+# else
J = SUNMatNewEmpty(kinsol_ctx);
-# endif
+# endif
J->content = this;
J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
};
status = KINSetJacFn(kinsol_mem, &setup_jacobian_callback<VectorType>);
AssertKINSOL(status);
-# endif
}
// call to KINSol
#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/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>
+# 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
-# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-# include <deal.II/sundials/sunlinsol_newempty.h>
-# endif
+# include <deal.II/sundials/n_vector.h>
DEAL_II_NAMESPACE_OPEN
<< " functions returned a negative error code: " << arg1
<< ". Please consult SUNDIALS manual.");
-# define AssertSundialsSolver(code) \
- Assert(code >= 0, ExcSundialsSolverError(code))
+# define AssertSundialsSolver(code) \
+ Assert(code >= 0, ExcSundialsSolverError(code))
namespace internal
{
template <typename VectorType>
internal::LinearSolverWrapper<VectorType>::LinearSolverWrapper(
LinearSolveFunction<VectorType> lsolve
-# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
,
SUNContext &linsol_ctx
-# endif
+# endif
)
: content(std::make_unique<LinearSolverContent<VectorType>>())
{
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
sun_linear_solver = SUNLinSolNewEmpty();
-# else
+# else
sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
-# endif
+# endif
sun_linear_solver->ops->gettype = arkode_linsol_get_type;
sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
SundialsOperator<VectorType>::vmult(VectorType & dst,
const VectorType &src) const
{
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
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);
-# else
+# else
// We don't currently know how to implement this: the
// internal::make_nvector_view() function called above requires a
// SUNContext object, but we don't actually have access to such an
(void)dst;
(void)src;
Assert(false, ExcNotImplemented());
-# endif
+# endif
}
return;
}
-# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
+# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
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
p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
(void)status;
AssertSundialsSolver(status);
-# else
+# else
// We don't currently know how to implement this: the
// internal::make_nvector_view() function called above requires a
// SUNContext object, but we don't actually have access to such an
(void)dst;
(void)src;
Assert(false, ExcNotImplemented());
-# endif
+# endif
}
template struct SundialsOperator<Vector<double>>;
template class internal::LinearSolverWrapper<
LinearAlgebra::distributed::BlockVector<double>>;
-# ifdef DEAL_II_WITH_MPI
-# ifdef DEAL_II_WITH_TRILINOS
+# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_TRILINOS
template struct SundialsOperator<TrilinosWrappers::MPI::Vector>;
template struct SundialsOperator<TrilinosWrappers::MPI::BlockVector>;
template class internal::LinearSolverWrapper<TrilinosWrappers::MPI::Vector>;
template class internal::LinearSolverWrapper<
TrilinosWrappers::MPI::BlockVector>;
-# endif // DEAL_II_WITH_TRILINOS
+# endif // DEAL_II_WITH_TRILINOS
-# ifdef DEAL_II_WITH_PETSC
-# ifndef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifndef PETSC_USE_COMPLEX
template struct SundialsOperator<PETScWrappers::MPI::Vector>;
template struct SundialsOperator<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 // PETSC_USE_COMPLEX
+# endif // DEAL_II_WITH_PETSC
-# endif // DEAL_II_WITH_MPI
+# endif // DEAL_II_WITH_MPI
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
-# endif
#endif
projects / dealii.git / commitdiff