Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
arkode.h
Go to the documentation of this file.
1//-----------------------------------------------------------
2//
3// Copyright (C) 2017 - 2021 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
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
67namespace 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
451
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
574 solve_ode_incrementally(VectorType & solution,
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)
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
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 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:420
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
unsigned int maximum_non_linear_iterations
Definition: arkode.h:492
void setup_mass_solver(const VectorType &solution)
Definition: arkode.cc:761
std::function< int(const double t, const VectorType &y, VectorType &res)> implicit_function
Definition: arkode.h:662
void * arkode_mem
Definition: arkode.h:1309
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(double t, const VectorType &r, VectorType &z, double tol, int lr)> mass_preconditioner_solve
Definition: arkode.h:1152
std::function< int(const VectorType &v, VectorType &Jv, double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
Definition: arkode.h:954
std::function< int(double t)> mass_preconditioner_setup
Definition: arkode.h:1178
void * get_arkode_memory() const
Definition: arkode.cc:833
std::function< VectorType &()> get_local_tolerances
Definition: arkode.h:1225
AdditionalData data
Definition: arkode.h:1304
std::function< int(realtype t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
Definition: arkode.h:991
LinearSolveFunction< VectorType > solve_linearized_system
Definition: arkode.h:1012
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
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm &mpi_comm=MPI_COMM_WORLD)
Definition: arkode.cc:380
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
Definition: arkode.cc:448
std::function< int(const double t)> mass_times_setup
Definition: arkode.h:921
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:505
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
Definition: arkode.h:1325
std::function< int(double t, const VectorType &y, const VectorType &fy, int jok, int &jcur, double gamma)> jacobian_preconditioner_setup
Definition: arkode.h:1121
std::function< int(const VectorType &rhs, VectorType &dst)> solve_mass_system
Definition: arkode.h:864
std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
Definition: arkode.h:643
MPI_Comm communicator
Definition: arkode.h:1316
std::function< bool(const double t, VectorType &sol)> solver_should_restart
Definition: arkode.h:1217
std::function< int(double t, const VectorType &v, VectorType &Mv)> mass_times_vector
Definition: arkode.h:885
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
Definition: arkode.h:1198
void set_functions_to_trigger_an_assert()
Definition: arkode.cc:818
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
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
Definition: arkode.h:1324
LinearSolveFunction< VectorType > solve_mass
Definition: arkode.h:1029
std::function< int(const double t)> setup_mass
Definition: arkode.h:843
unsigned int solve_ode(VectorType &solution)
Definition: arkode.cc:416
void setup_system_solver(const VectorType &solution)
Definition: arkode.cc:684
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
Definition: arkode.cc:429
double last_end_time
Definition: arkode.h:1321
std::function< void(void *arkode_mem)> custom_setup
Definition: arkode.h:1251
DEAL_II_DEPRECATED std::function< void(VectorType &)> reinit_vector
Definition: arkode.h:623
#define DEAL_II_DEPRECATED
Definition: config.h:162
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::function< int(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)