Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
arkode.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
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
22#ifdef DEAL_II_WITH_SUNDIALS
23
28
29# ifdef DEAL_II_WITH_PETSC
32# endif
33# include <deal.II/lac/vector.h>
35
36# include <arkode/arkode.h>
37# include <nvector/nvector_serial.h>
38# ifdef DEAL_II_WITH_MPI
39# include <nvector/nvector_parallel.h>
40# endif
42
46
47# include <boost/signals2.hpp>
48
49# include <sundials/sundials_linearsolver.h>
50# include <sundials/sundials_math.h>
51
52# include <exception>
53# include <memory>
54
55
57
58
59// Shorthand notation for ARKODE error codes.
60# define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
61
65namespace SUNDIALS
66{
327 template <typename VectorType = Vector<double>>
328 class ARKode
329 {
330 public:
335 {
336 public:
368 // Initial parameters
369 const double initial_time = 0.0,
370 const double final_time = 1.0,
371 const double initial_step_size = 1e-2,
372 const double output_period = 1e-1,
373 // Running parameters
374 const double minimum_step_size = 1e-6,
375 const unsigned int maximum_order = 5,
376 const unsigned int maximum_non_linear_iterations = 10,
377 const bool implicit_function_is_linear = false,
378 const bool implicit_function_is_time_independent = false,
379 const bool mass_is_time_independent = false,
380 const int anderson_acceleration_subspace = 3,
381 // Error parameters
382 const double absolute_tolerance = 1e-6,
383 const double relative_tolerance = 1e-5);
384
399 void
401
406
411
416
421
426
431
435 unsigned int maximum_order;
436
442
448
453
459
465
471 };
472
483
491 ARKode(const AdditionalData &data, const MPI_Comm mpi_comm);
492
496 ~ARKode();
497
505 unsigned int
506 solve_ode(VectorType &solution);
507
534 unsigned int
535 solve_ode_incrementally(VectorType &solution,
536 const double intermediate_time,
537 const bool reset_solver = false);
538
553 void
554 reset(const double t, const double h, const VectorType &y);
555
572 void *
573 get_arkode_memory() const;
574
591 std::function<
592 void(const double t, const VectorType &y, VectorType &explicit_f)>
594
611 std::function<void(const double t, const VectorType &y, VectorType &res)>
613
631 std::function<void(const double t, const VectorType &v, VectorType &Mv)>
633
668 std::function<void(const double t)> mass_times_setup;
669
696 std::function<void(const VectorType &v,
697 VectorType &Jv,
698 const double t,
699 const VectorType &y,
700 const VectorType &fy)>
702
737 std::function<
738 void(const double t, const VectorType &y, const VectorType &fy)>
740
768
792
824 std::function<void(const double t,
825 const VectorType &y,
826 const VectorType &fy,
827 const VectorType &r,
828 VectorType &z,
829 const double gamma,
830 const double tol,
831 const int lr)>
833
876 std::function<void(const double t,
877 const VectorType &y,
878 const VectorType &fy,
879 const int jok,
880 int &jcur,
881 const double gamma)>
883
911 std::function<void(const double t,
912 const VectorType &r,
913 VectorType &z,
914 const double tol,
915 const int lr)>
917
942 std::function<void(const double t)> mass_preconditioner_setup;
943
958 std::function<void(const double t,
959 const VectorType &sol,
960 const unsigned int step_number)>
962
980 std::function<bool(const double t, VectorType &sol)> solver_should_restart;
981
988 std::function<VectorType &()> get_local_tolerances;
989
1014 std::function<void(void *arkode_mem)> custom_setup;
1015
1016 private:
1022 std::string,
1023 << "Please provide an implementation for the function \""
1024 << arg1 << "\"");
1025
1029 unsigned int
1030 do_evolve_time(VectorType &solution,
1031 ::DiscreteTime &time,
1032 const bool do_reset);
1033
1040 void
1041 setup_system_solver(const VectorType &solution);
1042
1049 void
1050 setup_mass_solver(const VectorType &solution);
1051
1057 void
1059
1064
1069
1070# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
1075# endif
1076
1082
1087
1088 std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
1089 std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
1090
1096 mutable std::exception_ptr pending_exception;
1097
1098# ifdef DEAL_II_WITH_PETSC
1099# ifdef PETSC_USE_COMPLEX
1100 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
1101 "Sundials does not support complex scalar types, "
1102 "but PETSc is configured to use a complex scalar type!");
1103
1104 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
1105 "Sundials does not support complex scalar types, "
1106 "but PETSc is configured to use a complex scalar type!");
1107# endif // PETSC_USE_COMPLEX
1108# endif // DEAL_II_WITH_PETSC
1109 };
1110
1111
1116 int,
1117 << "One of the SUNDIALS ARKode internal functions "
1118 << " returned a negative error code: " << arg1
1119 << ". Please consult SUNDIALS manual.");
1120
1121
1122 template <typename VectorType>
1124 // Initial parameters
1125 const double initial_time,
1126 const double final_time,
1127 const double initial_step_size,
1128 const double output_period,
1129 // Running parameters
1130 const double minimum_step_size,
1131 const unsigned int maximum_order,
1132 const unsigned int maximum_non_linear_iterations,
1133 const bool implicit_function_is_linear,
1134 const bool implicit_function_is_time_independent,
1135 const bool mass_is_time_independent,
1136 const int anderson_acceleration_subspace,
1137 // Error parameters
1138 const double absolute_tolerance,
1139 const double relative_tolerance)
1140 : initial_time(initial_time)
1141 , final_time(final_time)
1142 , initial_step_size(initial_step_size)
1143 , minimum_step_size(minimum_step_size)
1144 , absolute_tolerance(absolute_tolerance)
1145 , relative_tolerance(relative_tolerance)
1146 , maximum_order(maximum_order)
1147 , output_period(output_period)
1148 , maximum_non_linear_iterations(maximum_non_linear_iterations)
1149 , implicit_function_is_linear(implicit_function_is_linear)
1150 , implicit_function_is_time_independent(
1151 implicit_function_is_time_independent)
1152 , mass_is_time_independent(mass_is_time_independent)
1153 , anderson_acceleration_subspace(anderson_acceleration_subspace)
1154 {}
1155
1156
1157
1158 template <typename VectorType>
1159 void
1161 {
1162 prm.add_parameter("Initial time", initial_time);
1163 prm.add_parameter("Final time", final_time);
1164 prm.add_parameter("Time interval between each output", output_period);
1165 prm.enter_subsection("Running parameters");
1166 prm.add_parameter("Initial step size", initial_step_size);
1167 prm.add_parameter("Minimum step size", minimum_step_size);
1168 prm.add_parameter("Maximum order of ARK", maximum_order);
1169 prm.add_parameter("Maximum number of nonlinear iterations",
1170 maximum_non_linear_iterations);
1171 prm.add_parameter("Implicit function is linear",
1172 implicit_function_is_linear);
1173 prm.add_parameter("Implicit function is time independent",
1174 implicit_function_is_time_independent);
1175 prm.add_parameter("Mass is time independent", mass_is_time_independent);
1176 prm.add_parameter("Anderson-acceleration subspace",
1177 anderson_acceleration_subspace);
1178 prm.leave_subsection();
1179 prm.enter_subsection("Error control");
1180 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
1181 prm.add_parameter("Relative error tolerance", relative_tolerance);
1182 prm.leave_subsection();
1183 }
1184
1185} // namespace SUNDIALS
1186
1188#endif
1189
1190
1191#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, const bool create_path_if_needed=true)
void add_parameters(ParameterHandler &prm)
Definition arkode.h:1160
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:1123
unsigned int maximum_non_linear_iterations
Definition arkode.h:447
void setup_mass_solver(const VectorType &solution)
Definition arkode.cc:602
std::function< void(const double t)> mass_preconditioner_setup
Definition arkode.h:942
std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
Definition arkode.h:593
void * arkode_mem
Definition arkode.h:1068
std::exception_ptr pending_exception
Definition arkode.h:1096
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:832
ARKode(const AdditionalData &data=AdditionalData())
Definition arkode.cc:53
void * get_arkode_memory() const
Definition arkode.cc:753
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
Definition arkode.h:961
std::function< VectorType &()> get_local_tolerances
Definition arkode.h:988
AdditionalData data
Definition arkode.h:1063
LinearSolveFunction< VectorType > solve_linearized_system
Definition arkode.h:767
std::function< void(const double t, const VectorType &y, VectorType &res)> implicit_function
Definition arkode.h:612
std::function< void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> mass_preconditioner_solve
Definition arkode.h:916
MPI_Comm mpi_communicator
Definition arkode.h:1081
void reset(const double t, const double h, const VectorType &y)
Definition arkode.cc:247
SUNContext arkode_ctx
Definition arkode.h:1074
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
Definition arkode.h:1089
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:882
std::function< bool(const double t, VectorType &sol)> solver_should_restart
Definition arkode.h:980
unsigned int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
Definition arkode.cc:146
std::function< void(const double t)> mass_times_setup
Definition arkode.h:668
void set_functions_to_trigger_an_assert()
Definition arkode.cc:742
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
Definition arkode.h:1088
std::function< void(const double t, const VectorType &v, VectorType &Mv)> mass_times_vector
Definition arkode.h:632
LinearSolveFunction< VectorType > solve_mass
Definition arkode.h:791
std::function< void(const double t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
Definition arkode.h:739
unsigned int solve_ode(VectorType &solution)
Definition arkode.cc:114
void setup_system_solver(const VectorType &solution)
Definition arkode.cc:392
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
Definition arkode.cc:127
double last_end_time
Definition arkode.h:1086
std::function< void(void *arkode_mem)> custom_setup
Definition arkode.h:1014
std::function< void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
Definition arkode.h:701
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcARKodeError(int arg1)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction