Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
arkode.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_sundials_arkode_h
18 #define dealii_sundials_arkode_h
19 
20 #include <deal.II/base/config.h>
21 
22 
23 #ifdef DEAL_II_WITH_SUNDIALS
24 
26 # include <deal.II/base/exceptions.h>
27 # include <deal.II/base/logstream.h>
28 # include <deal.II/base/mpi_stub.h>
30 
31 # ifdef DEAL_II_WITH_PETSC
34 # endif
35 # include <deal.II/lac/vector.h>
37 
38 # include <arkode/arkode.h>
39 # include <nvector/nvector_serial.h>
40 # ifdef DEAL_II_WITH_MPI
41 # include <nvector/nvector_parallel.h>
42 # endif
44 
48 
49 # include <boost/signals2.hpp>
50 
51 # include <sundials/sundials_linearsolver.h>
52 # include <sundials/sundials_math.h>
53 
54 # include <exception>
55 # include <memory>
56 
57 
59 
60 
61 // Shorthand notation for ARKODE error codes.
62 # define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
63 
67 namespace SUNDIALS
68 {
329  template <typename VectorType = Vector<double>>
330  class ARKode
331  {
332  public:
337  {
338  public:
370  // Initial parameters
371  const double initial_time = 0.0,
372  const double final_time = 1.0,
373  const double initial_step_size = 1e-2,
374  const double output_period = 1e-1,
375  // Running parameters
376  const double minimum_step_size = 1e-6,
377  const unsigned int maximum_order = 5,
378  const unsigned int maximum_non_linear_iterations = 10,
379  const bool implicit_function_is_linear = false,
380  const bool implicit_function_is_time_independent = false,
381  const bool mass_is_time_independent = false,
382  const int anderson_acceleration_subspace = 3,
383  // Error parameters
384  const double absolute_tolerance = 1e-6,
385  const double relative_tolerance = 1e-5);
386 
401  void
403 
407  double initial_time;
408 
412  double final_time;
413 
418 
423 
428 
433 
437  unsigned int maximum_order;
438 
444 
450 
455 
461 
467 
473  };
474 
485 
493  ARKode(const AdditionalData &data, const MPI_Comm mpi_comm);
494 
498  ~ARKode();
499 
507  unsigned int
508  solve_ode(VectorType &solution);
509 
536  unsigned int
537  solve_ode_incrementally(VectorType &solution,
538  const double intermediate_time,
539  const bool reset_solver = false);
540 
555  void
556  reset(const double t, const double h, const VectorType &y);
557 
574  void *
575  get_arkode_memory() const;
576 
593  std::function<
594  void(const double t, const VectorType &y, VectorType &explicit_f)>
596 
613  std::function<void(const double t, const VectorType &y, VectorType &res)>
615 
633  std::function<void(const double t, const VectorType &v, VectorType &Mv)>
635 
670  std::function<void(const double t)> mass_times_setup;
671 
698  std::function<void(const VectorType &v,
699  VectorType &Jv,
700  const double t,
701  const VectorType &y,
702  const VectorType &fy)>
704 
739  std::function<
740  void(const double t, const VectorType &y, const VectorType &fy)>
742 
770 
794 
826  std::function<void(const double t,
827  const VectorType &y,
828  const VectorType &fy,
829  const VectorType &r,
830  VectorType &z,
831  const double gamma,
832  const double tol,
833  const int lr)>
835 
878  std::function<void(const double t,
879  const VectorType &y,
880  const VectorType &fy,
881  const int jok,
882  int &jcur,
883  const double gamma)>
885 
913  std::function<void(const double t,
914  const VectorType &r,
915  VectorType &z,
916  const double tol,
917  const int lr)>
919 
944  std::function<void(const double t)> mass_preconditioner_setup;
945 
960  std::function<void(const double t,
961  const VectorType &sol,
962  const unsigned int step_number)>
964 
982  std::function<bool(const double t, VectorType &sol)> solver_should_restart;
983 
990  std::function<VectorType &()> get_local_tolerances;
991 
1016  std::function<void(void *arkode_mem)> custom_setup;
1017 
1018  private:
1024  std::string,
1025  << "Please provide an implementation for the function \""
1026  << arg1 << "\"");
1027 
1031  int
1032  do_evolve_time(VectorType &solution,
1033  ::DiscreteTime &time,
1034  const bool do_reset);
1035 
1042  void
1043  setup_system_solver(const VectorType &solution);
1044 
1051  void
1052  setup_mass_solver(const VectorType &solution);
1053 
1059  void
1061 
1066 
1070  void *arkode_mem;
1071 
1072 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
1077 # endif
1078 
1084 
1089 
1090  std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
1091  std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
1092 
1098  mutable std::exception_ptr pending_exception;
1099 
1100 # ifdef DEAL_II_WITH_PETSC
1101 # ifdef PETSC_USE_COMPLEX
1102  static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
1103  "Sundials does not support complex scalar types, "
1104  "but PETSc is configured to use a complex scalar type!");
1105 
1106  static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
1107  "Sundials does not support complex scalar types, "
1108  "but PETSc is configured to use a complex scalar type!");
1109 # endif // PETSC_USE_COMPLEX
1110 # endif // DEAL_II_WITH_PETSC
1111  };
1112 
1113 
1117  DeclException1(ExcARKodeError,
1118  int,
1119  << "One of the SUNDIALS ARKode internal functions "
1120  << " returned a negative error code: " << arg1
1121  << ". Please consult SUNDIALS manual.");
1122 
1123 
1124  template <typename VectorType>
1126  // Initial parameters
1127  const double initial_time,
1128  const double final_time,
1129  const double initial_step_size,
1130  const double output_period,
1131  // Running parameters
1132  const double minimum_step_size,
1133  const unsigned int maximum_order,
1134  const unsigned int maximum_non_linear_iterations,
1135  const bool implicit_function_is_linear,
1136  const bool implicit_function_is_time_independent,
1137  const bool mass_is_time_independent,
1138  const int anderson_acceleration_subspace,
1139  // Error parameters
1140  const double absolute_tolerance,
1141  const double relative_tolerance)
1142  : initial_time(initial_time)
1143  , final_time(final_time)
1144  , initial_step_size(initial_step_size)
1145  , minimum_step_size(minimum_step_size)
1146  , absolute_tolerance(absolute_tolerance)
1147  , relative_tolerance(relative_tolerance)
1148  , maximum_order(maximum_order)
1149  , output_period(output_period)
1150  , maximum_non_linear_iterations(maximum_non_linear_iterations)
1151  , implicit_function_is_linear(implicit_function_is_linear)
1152  , implicit_function_is_time_independent(
1153  implicit_function_is_time_independent)
1154  , mass_is_time_independent(mass_is_time_independent)
1155  , anderson_acceleration_subspace(anderson_acceleration_subspace)
1156  {}
1157 
1158 
1159 
1160  template <typename VectorType>
1161  void
1163  {
1164  prm.add_parameter("Initial time", initial_time);
1165  prm.add_parameter("Final time", final_time);
1166  prm.add_parameter("Time interval between each output", output_period);
1167  prm.enter_subsection("Running parameters");
1168  prm.add_parameter("Initial step size", initial_step_size);
1169  prm.add_parameter("Minimum step size", minimum_step_size);
1170  prm.add_parameter("Maximum order of ARK", maximum_order);
1171  prm.add_parameter("Maximum number of nonlinear iterations",
1172  maximum_non_linear_iterations);
1173  prm.add_parameter("Implicit function is linear",
1174  implicit_function_is_linear);
1175  prm.add_parameter("Implicit function is time independent",
1176  implicit_function_is_time_independent);
1177  prm.add_parameter("Mass is time independent", mass_is_time_independent);
1178  prm.add_parameter("Anderson-acceleration subspace",
1179  anderson_acceleration_subspace);
1180  prm.leave_subsection();
1181  prm.enter_subsection("Error control");
1182  prm.add_parameter("Absolute error tolerance", absolute_tolerance);
1183  prm.add_parameter("Relative error tolerance", relative_tolerance);
1184  prm.leave_subsection();
1185  }
1186 
1187 } // namespace SUNDIALS
1188 
1190 #endif
1191 
1192 
1193 #endif
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection)
void add_parameters(ParameterHandler &prm)
Definition: arkode.h:1162
AdditionalData(const double initial_time=0.0, const double final_time=1.0, const double initial_step_size=1e-2, const double output_period=1e-1, const double minimum_step_size=1e-6, const unsigned int maximum_order=5, const unsigned int maximum_non_linear_iterations=10, const bool implicit_function_is_linear=false, const bool implicit_function_is_time_independent=false, const bool mass_is_time_independent=false, const int anderson_acceleration_subspace=3, const double absolute_tolerance=1e-6, const double relative_tolerance=1e-5)
Definition: arkode.h:1125
unsigned int maximum_non_linear_iterations
Definition: arkode.h:449
void setup_mass_solver(const VectorType &solution)
Definition: arkode.cc:578
std::function< void(const double t)> mass_preconditioner_setup
Definition: arkode.h:944
std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
Definition: arkode.h:595
void * arkode_mem
Definition: arkode.h:1070
std::exception_ptr pending_exception
Definition: arkode.h:1098
std::function< void(const double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, const double gamma, const double tol, const int lr)> jacobian_preconditioner_solve
Definition: arkode.h:834
ARKode(const AdditionalData &data=AdditionalData())
Definition: arkode.cc:54
void * get_arkode_memory() const
Definition: arkode.cc:727
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
Definition: arkode.h:963
std::function< VectorType &()> get_local_tolerances
Definition: arkode.h:990
AdditionalData data
Definition: arkode.h:1065
LinearSolveFunction< VectorType > solve_linearized_system
Definition: arkode.h:769
std::function< void(const double t, const VectorType &y, VectorType &res)> implicit_function
Definition: arkode.h:614
std::function< void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> mass_preconditioner_solve
Definition: arkode.h:918
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
Definition: arkode.cc:140
MPI_Comm mpi_communicator
Definition: arkode.h:1083
DeclException1(ExcFunctionNotProvided, std::string,<< "Please provide an implementation for the function \""<< arg1<< "\"")
void reset(const double t, const double h, const VectorType &y)
Definition: arkode.cc:235
SUNContext arkode_ctx
Definition: arkode.h:1076
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
Definition: arkode.h:1091
std::function< void(const double t, const VectorType &y, const VectorType &fy, const int jok, int &jcur, const double gamma)> jacobian_preconditioner_setup
Definition: arkode.h:884
std::function< bool(const double t, VectorType &sol)> solver_should_restart
Definition: arkode.h:982
std::function< void(const double t)> mass_times_setup
Definition: arkode.h:670
void set_functions_to_trigger_an_assert()
Definition: arkode.cc:716
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
Definition: arkode.h:1090
std::function< void(const double t, const VectorType &v, VectorType &Mv)> mass_times_vector
Definition: arkode.h:634
LinearSolveFunction< VectorType > solve_mass
Definition: arkode.h:793
std::function< void(const double t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
Definition: arkode.h:741
unsigned int solve_ode(VectorType &solution)
Definition: arkode.cc:108
void setup_system_solver(const VectorType &solution)
Definition: arkode.cc:369
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
Definition: arkode.cc:121
double last_end_time
Definition: arkode.h:1088
std::function< void(void *arkode_mem)> custom_setup
Definition: arkode.h:1016
std::function< void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
Definition: arkode.h:703
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
DeclException1(ExcARKodeError, int,<< "One of the SUNDIALS ARKode internal functions "<< " returned a negative error code: "<< arg1<< ". Please consult SUNDIALS manual.")
long double gamma(const unsigned int n)