Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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 - 2020 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 #ifndef dealii_sundials_arkode_h
17 #define dealii_sundials_arkode_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/mpi.h>
22 
23 #ifdef DEAL_II_WITH_SUNDIALS
24 
26 # include <deal.II/base/exceptions.h>
27 # include <deal.II/base/logstream.h>
29 # ifdef DEAL_II_WITH_PETSC
32 # endif
33 # include <deal.II/lac/vector.h>
35 
36 # include <arkode/arkode.h>
37 # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
38 # include <arkode/arkode_impl.h>
39 # endif
40 # include <nvector/nvector_serial.h>
41 # ifdef DEAL_II_WITH_MPI
42 # include <nvector/nvector_parallel.h>
43 # endif
45 
48 
49 # include <boost/signals2.hpp>
50 
51 # include <sundials/sundials_linearsolver.h>
52 # include <sundials/sundials_math.h>
53 # include <sundials/sundials_types.h>
54 
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 {
332  template <typename VectorType = Vector<double>>
333  class ARKode
334  {
335  public:
340  {
341  public:
373  // Initial parameters
374  const double initial_time = 0.0,
375  const double final_time = 1.0,
376  const double initial_step_size = 1e-2,
377  const double output_period = 1e-1,
378  // Running parameters
379  const double minimum_step_size = 1e-6,
380  const unsigned int maximum_order = 5,
381  const unsigned int maximum_non_linear_iterations = 10,
382  const bool implicit_function_is_linear = false,
383  const bool implicit_function_is_time_independent = false,
384  const bool mass_is_time_independent = false,
385  const int anderson_acceleration_subspace = 3,
386  // Error parameters
387  const double absolute_tolerance = 1e-6,
388  const double relative_tolerance = 1e-5)
403  {}
404 
419  void
421  {
422  prm.add_parameter("Initial time", initial_time);
423  prm.add_parameter("Final time", final_time);
424  prm.add_parameter("Time interval between each output", output_period);
425 
426  prm.enter_subsection("Running parameters");
427  prm.add_parameter("Initial step size", initial_step_size);
428  prm.add_parameter("Minimum step size", minimum_step_size);
429  prm.add_parameter("Maximum order of ARK", maximum_order);
430  prm.add_parameter("Maximum number of nonlinear iterations",
432  prm.add_parameter("Implicit function is linear",
434  prm.add_parameter("Implicit function is time independent",
436  prm.add_parameter("Mass is time independent", mass_is_time_independent);
437  prm.add_parameter("Anderson-acceleration subspace",
439  prm.leave_subsection();
440 
441  prm.enter_subsection("Error control");
442  prm.add_parameter("Absolute error tolerance", absolute_tolerance);
443  prm.add_parameter("Relative error tolerance", relative_tolerance);
444  prm.leave_subsection();
445  }
446 
450  double initial_time;
451 
455  double final_time;
456 
461 
466 
471 
476 
480  unsigned int maximum_order;
481 
487 
493 
498 
504 
510 
516  };
517 
530  const MPI_Comm & mpi_comm = MPI_COMM_WORLD);
531 
535  ~ARKode();
536 
544  unsigned int
545  solve_ode(VectorType &solution);
546 
573  unsigned int
575  const double intermediate_time,
576  const bool reset_solver = false);
577 
592  void
593  reset(const double t, const double h, const VectorType &y);
594 
611  void *
612  get_arkode_memory() const;
613 
623  std::function<void(VectorType &)> reinit_vector;
624 
641  std::function<
642  int(const double t, const VectorType &y, VectorType &explicit_f)>
644 
661  std::function<int(const double t, const VectorType &y, VectorType &res)>
663 
664 # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
665 
749  std::function<int(const int convfail,
750  const double t,
751  const double gamma,
752  const VectorType &ypred,
753  const VectorType &fpred,
754  bool & j_is_current)>
756 
801  std::function<int(const double t,
802  const double gamma,
803  const VectorType &ycur,
804  const VectorType &fcur,
805  const VectorType &rhs,
806  VectorType & dst)>
808 
809 
843  std::function<int(const double t)> setup_mass;
844 
863  std::function<int(const VectorType &rhs, VectorType &dst)>
865 # else
866 
884  std::function<int(double t, const VectorType &v, VectorType &Mv)>
886 
921  std::function<int(const double t)> mass_times_setup;
922 
949  std::function<int(const VectorType &v,
950  VectorType & Jv,
951  double t,
952  const VectorType &y,
953  const VectorType &fy)>
955 
990  std::function<int(realtype t, const VectorType &y, const VectorType &fy)>
992 
1013 
1030 
1031 
1063  std::function<int(double t,
1064  const VectorType &y,
1065  const VectorType &fy,
1066  const VectorType &r,
1067  VectorType & z,
1068  double gamma,
1069  double tol,
1070  int lr)>
1072 
1115  std::function<int(double t,
1116  const VectorType &y,
1117  const VectorType &fy,
1118  int jok,
1119  int & jcur,
1120  double gamma)>
1122 
1150  std::function<
1151  int(double t, const VectorType &r, VectorType &z, double tol, int lr)>
1153 
1178  std::function<int(double t)> mass_preconditioner_setup;
1179 # endif
1180 
1195  std::function<void(const double t,
1196  const VectorType & sol,
1197  const unsigned int step_number)>
1199 
1217  std::function<bool(const double t, VectorType &sol)> solver_should_restart;
1218 
1225  std::function<VectorType &()> get_local_tolerances;
1226 
1251  std::function<void(void *arkode_mem)> custom_setup;
1252 
1253  private:
1258  DeclException1(ExcFunctionNotProvided,
1259  std::string,
1260  << "Please provide an implementation for the function \""
1261  << arg1 << "\"");
1262 
1266  int
1267  do_evolve_time(VectorType & solution,
1268  ::DiscreteTime &time,
1269  const bool do_reset);
1270 
1271 # if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
1272 
1279  void
1280  setup_system_solver(const VectorType &solution);
1281 
1288  void
1289  setup_mass_solver(const VectorType &solution);
1290 
1291 # endif
1292 
1298  void
1300 
1305 
1309  void *arkode_mem;
1310 
1317 
1322 
1323 # if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
1324  std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
1325  std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
1326 # endif
1327 
1328 # ifdef DEAL_II_WITH_PETSC
1329 # ifdef PETSC_USE_COMPLEX
1330  static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
1331  "Sundials does not support complex scalar types, "
1332  "but PETSc is configured to use a complex scalar type!");
1333 
1334  static_assert(
1335  !std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::value,
1336  "Sundials does not support complex scalar types, "
1337  "but PETSc is configured to use a complex scalar type!");
1338 # endif // PETSC_USE_COMPLEX
1339 # endif // DEAL_II_WITH_PETSC
1340  };
1341 
1342 
1346  DeclException1(ExcARKodeError,
1347  int,
1348  << "One of the SUNDIALS ARKode internal functions "
1349  << " returned a negative error code: " << arg1
1350  << ". Please consult SUNDIALS manual.");
1351 
1352 } // namespace SUNDIALS
1353 
1355 #endif
1356 
1357 
1358 #endif
void setup_mass_solver(const VectorType &solution)
Definition: arkode.cc:761
void reset(const double t, const double h, const VectorType &y)
Definition: arkode.cc:505
void * arkode_mem
Definition: arkode.h:1309
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
LinearSolveFunction< VectorType > solve_linearized_system
Definition: arkode.h:1012
double last_end_time
Definition: arkode.h:1321
void * get_arkode_memory() const
Definition: arkode.cc:833
std::function< int(const VectorType &rhs, VectorType &dst)> solve_mass_system
Definition: arkode.h:864
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)
DEAL_II_DEPRECATED std::function< void(VectorType &)> reinit_vector
Definition: arkode.h:623
MPI_Comm communicator
Definition: arkode.h:1316
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
Definition: arkode.h:1198
std::function< VectorType &()> get_local_tolerances
Definition: arkode.h:1225
unsigned int maximum_non_linear_iterations
Definition: arkode.h:492
std::function< int(const double t)> setup_mass
Definition: arkode.h:843
std::function< bool(const double t, VectorType &sol)> solver_should_restart
Definition: arkode.h:1217
void enter_subsection(const std::string &subsection)
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
Definition: arkode.h:1325
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:372
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm &mpi_comm=MPI_COMM_WORLD)
Definition: arkode.cc:380
std::function< int(const double t, const double gamma, const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system
Definition: arkode.h:807
std::function< int(const double t)> mass_times_setup
Definition: arkode.h:921
AdditionalData data
Definition: arkode.h:1304
std::function< int(double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, double gamma, double tol, int lr)> jacobian_preconditioner_solve
Definition: arkode.h:1071
void setup_system_solver(const VectorType &solution)
Definition: arkode.cc:684
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
std::function< int(const double t, const VectorType &y, VectorType &res)> implicit_function
Definition: arkode.h:662
std::function< int(realtype t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
Definition: arkode.h:991
LinearSolveFunction< VectorType > solve_mass
Definition: arkode.h:1029
std::function< int(double t, const VectorType &v, VectorType &Mv)> mass_times_vector
Definition: arkode.h:885
void add_parameters(ParameterHandler &prm)
Definition: arkode.h:420
std::function< int(double t, const VectorType &y, const VectorType &fy, int jok, int &jcur, double gamma)> jacobian_preconditioner_setup
Definition: arkode.h:1121
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
Definition: arkode.h:643
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
Definition: arkode.cc:448
std::function< int(const VectorType &v, VectorType &Jv, double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
Definition: arkode.h:954
DeclException1(ExcFunctionNotProvided, std::string,<< "Please provide an implementation for the function << arg1<< "\")
std::function< int(double t)> mass_preconditioner_setup
Definition: arkode.h:1178
std::function< void(void *arkode_mem)> custom_setup
Definition: arkode.h:1251
unsigned int solve_ode(VectorType &solution)
Definition: arkode.cc:416
std::function< int(const int convfail, const double t, const double gamma, const VectorType &ypred, const VectorType &fpred, bool &j_is_current)> setup_jacobian
Definition: arkode.h:755
std::function< int(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction
std::function< int(double t, const VectorType &r, VectorType &z, double tol, int lr)> mass_preconditioner_solve
Definition: arkode.h:1152
long double gamma(const unsigned int n)
void set_functions_to_trigger_an_assert()
Definition: arkode.cc:818
#define DEAL_II_DEPRECATED
Definition: config.h:158
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
Definition: arkode.cc:429
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
Definition: arkode.h:1324