# include <deal.II/base/discrete_time.h>
# include <deal.II/sundials/n_vector.h>
+# include <deal.II/sundials/sundials_types.h>
# include <deal.II/sundials/sunlinsol_wrapper.h>
# include <boost/signals2.hpp>
# include <sundials/sundials_linearsolver.h>
# include <sundials/sundials_math.h>
-# include <sundials/sundials_types.h>
# include <exception>
# include <memory>
# include <ida/ida.h>
# endif
+# include <deal.II/sundials/sundials_types.h>
# 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 <memory>
# include <deal.II/lac/vector.h>
# include <deal.II/lac/vector_memory.h>
+# include <deal.II/sundials/sundials_types.h>
+
# include <boost/signals2.hpp>
# include <kinsol/kinsol.h>
# include <nvector/nvector_serial.h>
# include <sundials/sundials_math.h>
-# include <sundials/sundials_types.h>
# include <exception>
# include <memory>
#include <deal.II/base/config.h>
#include <deal.II/sundials/n_vector.h>
+#include <deal.II/sundials/sundials_types.h>
#ifdef DEAL_II_WITH_SUNDIALS
*/
namespace NVectorOperations
{
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- using realtype = ::sunrealtype;
-# else
- using realtype = ::realtype;
-# endif
-
N_Vector_ID
get_vector_id(N_Vector v);
template <typename VectorType>
void
- linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z);
+ linear_sum(SUNDIALS::realtype a,
+ N_Vector x,
+ SUNDIALS::realtype b,
+ N_Vector y,
+ N_Vector z);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
dot_product(N_Vector x, N_Vector y);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_l2_norm(N_Vector x, N_Vector y);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
l1_norm(N_Vector x);
template <typename VectorType>
elementwise_abs(N_Vector x, N_Vector z);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_rms_norm(N_Vector x, N_Vector w);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
max_norm(N_Vector x);
template <typename VectorType,
std::enable_if_t<is_serial_vector<VectorType>::value, int> = 0>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x);
template <typename VectorType,
std::enable_if_t<!is_serial_vector<VectorType>::value &&
!IsBlockVector<VectorType>::value,
int> = 0>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x);
template <typename VectorType,
std::enable_if_t<!is_serial_vector<VectorType>::value &&
IsBlockVector<VectorType>::value,
int> = 0>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x);
template <typename VectorType>
void
- scale(realtype c, N_Vector x, N_Vector z);
+ scale(SUNDIALS::realtype c, N_Vector x, N_Vector z);
template <typename VectorType>
void
- set_constant(realtype c, N_Vector v);
+ set_constant(SUNDIALS::realtype c, N_Vector v);
template <typename VectorType>
void
- add_constant(N_Vector x, realtype b, N_Vector z);
+ add_constant(N_Vector x, SUNDIALS::realtype b, N_Vector z);
template <typename VectorType>
const MPI_Comm &
template <typename VectorType>
void
- linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
+ linear_sum(SUNDIALS::realtype a,
+ N_Vector x,
+ SUNDIALS::realtype b,
+ N_Vector y,
+ N_Vector z)
{
auto *x_dealii = unwrap_nvector_const<VectorType>(x);
auto *y_dealii = unwrap_nvector_const<VectorType>(y);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
dot_product(N_Vector x, N_Vector y)
{
return *unwrap_nvector_const<VectorType>(x) *
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_l2_norm(N_Vector x, N_Vector w)
{
// TODO copy can be avoided by a custom kernel
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
l1_norm(N_Vector x)
{
return unwrap_nvector_const<VectorType>(x)->l1_norm();
template <typename VectorType>
void
- set_constant(realtype c, N_Vector v)
+ set_constant(SUNDIALS::realtype c, N_Vector v)
{
auto *v_dealii = unwrap_nvector<VectorType>(v);
*v_dealii = c;
template <typename VectorType>
void
- add_constant(N_Vector x, realtype b, N_Vector z)
+ add_constant(N_Vector x, SUNDIALS::realtype b, N_Vector z)
{
auto *x_dealii = unwrap_nvector_const<VectorType>(x);
auto *z_dealii = unwrap_nvector<VectorType>(z);
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_rms_norm(N_Vector x, N_Vector w)
{
// TODO copy can be avoided by a custom kernel
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask)
{
// TODO copy can be avoided by a custom kernel
template <typename VectorType>
- realtype
+ SUNDIALS::realtype
max_norm(N_Vector x)
{
return unwrap_nvector_const<VectorType>(x)->linfty_norm();
template <typename VectorType,
std::enable_if_t<is_serial_vector<VectorType>::value, int>>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x)
{
auto *vector = unwrap_nvector_const<VectorType>(x);
std::enable_if_t<!is_serial_vector<VectorType>::value &&
!IsBlockVector<VectorType>::value,
int>>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x)
{
auto *vector = unwrap_nvector_const<VectorType>(x);
std::enable_if_t<!is_serial_vector<VectorType>::value &&
IsBlockVector<VectorType>::value,
int>>
- realtype
+ SUNDIALS::realtype
min_element(N_Vector x)
{
auto *vector = unwrap_nvector_const<VectorType>(x);
template <typename VectorType>
void
- scale(realtype c, N_Vector x, N_Vector z)
+ scale(SUNDIALS::realtype c, N_Vector x, N_Vector z)
{
auto *x_dealii = unwrap_nvector_const<VectorType>(x);
auto *z_dealii = unwrap_nvector<VectorType>(z);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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_types_h
+#define dealii_sundials_types_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+# include <sundials/sundials_types.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+/**
+ * Alias for the real type used by SUNDIALS.
+ */
+# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
+ using realtype = ::sunrealtype;
+# else
+ using realtype = ::realtype;
+# endif
+} // namespace SUNDIALS
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SUNDIALS
+#endif // dealii_sundials_types_h
Assert(explicit_function || implicit_function,
ExcFunctionNotProvided("explicit_function || implicit_function"));
- auto explicit_function_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tt,
-# else
- realtype tt,
-# endif
- N_Vector yy,
- N_Vector yp,
- void *user_data) -> int {
+ auto explicit_function_callback = [](SUNDIALS::realtype tt,
+ N_Vector yy,
+ N_Vector yp,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
};
- auto implicit_function_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tt,
-# else
- realtype tt,
-# endif
- N_Vector yy,
- N_Vector yp,
- void *user_data) -> int {
+ auto implicit_function_callback = [](SUNDIALS::realtype tt,
+ N_Vector yy,
+ N_Vector yp,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr);
AssertARKode(status);
- auto jacobian_times_vector_callback = [](N_Vector v,
- N_Vector Jv,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- N_Vector y,
- N_Vector fy,
- void *user_data,
+ auto jacobian_times_vector_callback = [](N_Vector v,
+ N_Vector Jv,
+ SUNDIALS::realtype t,
+ N_Vector y,
+ N_Vector fy,
+ void *user_data,
N_Vector) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*src_fy);
};
- auto jacobian_times_vector_setup_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- N_Vector y,
- N_Vector fy,
- void *user_data) -> int {
+ auto jacobian_times_vector_setup_callback = [](SUNDIALS::realtype t,
+ N_Vector y,
+ N_Vector fy,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
AssertARKode(status);
if (jacobian_preconditioner_solve)
{
- auto solve_with_jacobian_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- N_Vector y,
- N_Vector fy,
- N_Vector r,
- N_Vector z,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype gamma,
- sunrealtype delta,
-# else
- realtype gamma,
- realtype delta,
-# endif
- int lr,
- void *user_data) -> int {
+ auto solve_with_jacobian_callback = [](SUNDIALS::realtype t,
+ N_Vector y,
+ N_Vector fy,
+ N_Vector r,
+ N_Vector z,
+ SUNDIALS::realtype gamma,
+ SUNDIALS::realtype delta,
+ int lr,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
lr);
};
- auto jacobian_solver_setup_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- N_Vector y,
- N_Vector fy,
- booleantype jok,
- booleantype *jcurPtr,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype gamma,
-# else
- realtype gamma,
-# endif
- void *user_data) -> int {
+ auto jacobian_solver_setup_callback = [](SUNDIALS::realtype t,
+ N_Vector y,
+ N_Vector fy,
+ booleantype jok,
+ booleantype *jcurPtr,
+ SUNDIALS::realtype gamma,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
AssertARKode(status);
auto mass_matrix_times_vector_setup_callback =
- [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- void *mtimes_data) -> int {
+ [](SUNDIALS::realtype t, void *mtimes_data) -> int {
Assert(mtimes_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(mtimes_data);
solver.mass_times_setup, solver.pending_exception, t);
};
- auto mass_matrix_times_vector_callback = [](N_Vector v,
- N_Vector Mv,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
+ auto mass_matrix_times_vector_callback = [](N_Vector v,
+ N_Vector Mv,
+ SUNDIALS::realtype t,
void *mtimes_data) -> int {
Assert(mtimes_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
if (mass_preconditioner_solve)
{
- auto mass_matrix_solver_setup_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- void *user_data) -> int {
+ auto mass_matrix_solver_setup_callback =
+ [](SUNDIALS::realtype t, void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
solver.mass_preconditioner_setup, solver.pending_exception, t);
};
- auto solve_with_mass_matrix_callback = [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-# else
- realtype t,
-# endif
- N_Vector r,
- N_Vector z,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype delta,
-# else
- realtype delta,
-# endif
- int lr,
- void *user_data) -> int {
+ auto solve_with_mass_matrix_callback = [](SUNDIALS::realtype t,
+ N_Vector r,
+ N_Vector z,
+ SUNDIALS::realtype delta,
+ int lr,
+ void *user_data) -> int {
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
status = IDAInit(
ida_mem,
- [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tt,
-# else
- realtype tt,
-# endif
- N_Vector yy,
- N_Vector yp,
- N_Vector rr,
- void *user_data) -> int {
+ [](SUNDIALS::realtype tt,
+ N_Vector yy,
+ N_Vector yp,
+ N_Vector rr,
+ void *user_data) -> int {
IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
LS = SUNLinSolNewEmpty();
# else
- LS = SUNLinSolNewEmpty(ida_ctx);
+ LS = SUNLinSolNewEmpty(ida_ctx);
# endif
LS->content = this;
ExcFunctionNotProvided("solve_with_jacobian"));
LS->ops->solve = [](SUNLinearSolver LS,
SUNMatrix /*ignored*/,
- N_Vector x,
- N_Vector b,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tol
-# else
- realtype tol
-# endif
- ) -> int {
+ N_Vector x,
+ N_Vector b,
+ SUNDIALS::realtype tol) -> int {
IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
J = SUNMatNewEmpty();
# else
- J = SUNMatNewEmpty(ida_ctx);
+ J = SUNMatNewEmpty(ida_ctx);
# endif
J->content = this;
// calling IDASetLinearSolver
status = IDASetJacFn(
ida_mem,
- [](
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tt,
- sunrealtype cj,
-# else
- realtype tt,
- realtype cj,
-# endif
- N_Vector yy,
- N_Vector yp,
- N_Vector /* residual */,
- SUNMatrix /* ignored */,
- void *user_data,
- N_Vector /* tmp1 */,
- N_Vector /* tmp2 */,
- N_Vector /* tmp3 */) -> int {
+ [](SUNDIALS::realtype tt,
+ SUNDIALS::realtype cj,
+ N_Vector yy,
+ N_Vector yp,
+ N_Vector /* residual */,
+ SUNMatrix /* ignored */,
+ void *user_data,
+ N_Vector /* tmp1 */,
+ N_Vector /* tmp2 */,
+ N_Vector /* tmp3 */) -> int {
Assert(user_data != nullptr, ExcInternalError());
IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
LS->ops->solve = [](SUNLinearSolver LS,
SUNMatrix /*ignored*/,
- N_Vector x,
- N_Vector b,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tol
-# else
- realtype tol
-# endif
- ) -> int {
+ N_Vector x,
+ N_Vector b,
+ SUNDIALS::realtype tol) -> int {
// Receive the object that describes the linear solver and
// unpack the pointer to the KINSOL object from which we can then
// get the 'reinit' and 'solve' functions.
# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
J = SUNMatNewEmpty();
# else
- J = SUNMatNewEmpty(kinsol_ctx);
+ J = SUNMatNewEmpty(kinsol_ctx);
# endif
J->content = this;
# endif
# include <deal.II/sundials/n_vector.h>
+# include <deal.II/sundials/sundials_types.h>
DEAL_II_NAMESPACE_OPEN
int
arkode_linsol_solve(SUNLinearSolver LS,
SUNMatrix /*ignored*/,
- N_Vector x,
- N_Vector b,
-# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype tol
-# else
- realtype tol
-# endif
- )
+ N_Vector x,
+ N_Vector b,
+ SUNDIALS::realtype tol)
{
LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
};
- ode.jacobian_times_setup = [&](
-#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-#else
- realtype t,
-#endif
- const VectorType &y,
- const VectorType &fy) {
- J = 0;
- J(2, 2) = -1.0 / eps;
- };
+ ode.jacobian_times_setup =
+ [&](SUNDIALS::realtype t, const VectorType &y, const VectorType &fy) {
+ J = 0;
+ J(2, 2) = -1.0 / eps;
+ };
ode.jacobian_times_vector = [&](const VectorType &v,
VectorType &Jv,
};
- ode.jacobian_times_setup = [&](
-#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- sunrealtype t,
-#else
- realtype t,
-#endif
- const VectorType &y,
- const VectorType &fy) {
- J = 0;
- J(2, 2) = -1.0 / eps;
- };
+ ode.jacobian_times_setup =
+ [&](SUNDIALS::realtype t, const VectorType &y, const VectorType &fy) {
+ J = 0;
+ J(2, 2) = -1.0 / eps;
+ };
ode.jacobian_times_vector = [&](const VectorType &v,
VectorType &Jv,