Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2749-g12bc158a88 2025-03-04 14:20:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
62namespace SUNDIALS
63{
324 template <typename VectorType = Vector<double>>
325 class ARKode
326 {
327 public:
332 {
333 public:
365 // Initial parameters
366 const double initial_time = 0.0,
367 const double final_time = 1.0,
368 const double initial_step_size = 1e-2,
369 const double output_period = 1e-1,
370 // Running parameters
371 const double minimum_step_size = 1e-6,
372 const unsigned int maximum_order = 5,
373 const unsigned int maximum_non_linear_iterations = 10,
374 const bool implicit_function_is_linear = false,
375 const bool implicit_function_is_time_independent = false,
376 const bool mass_is_time_independent = false,
377 const int anderson_acceleration_subspace = 3,
378 // Error parameters
379 const double absolute_tolerance = 1e-6,
380 const double relative_tolerance = 1e-5);
381
396 void
398
403
408
413
418
423
428
432 unsigned int maximum_order;
433
439
445
450
456
462
468 };
469
480
488 ARKode(const AdditionalData &data, const MPI_Comm mpi_comm);
489
493 ~ARKode();
494
502 unsigned int
503 solve_ode(VectorType &solution);
504
531 unsigned int
532 solve_ode_incrementally(VectorType &solution,
533 const double intermediate_time,
534 const bool reset_solver = false);
535
550 void
551 reset(const double t, const double h, const VectorType &y);
552
569 void *
570 get_arkode_memory() const;
571
588 std::function<
589 void(const double t, const VectorType &y, VectorType &explicit_f)>
591
608 std::function<void(const double t, const VectorType &y, VectorType &res)>
610
628 std::function<void(const double t, const VectorType &v, VectorType &Mv)>
630
665 std::function<void(const double t)> mass_times_setup;
666
693 std::function<void(const VectorType &v,
694 VectorType &Jv,
695 const double t,
696 const VectorType &y,
697 const VectorType &fy)>
699
734 std::function<
735 void(const double t, const VectorType &y, const VectorType &fy)>
737
765
789
821 std::function<void(const double t,
822 const VectorType &y,
823 const VectorType &fy,
824 const VectorType &r,
825 VectorType &z,
826 const double gamma,
827 const double tol,
828 const int lr)>
830
873 std::function<void(const double t,
874 const VectorType &y,
875 const VectorType &fy,
876 const int jok,
877 int &jcur,
878 const double gamma)>
880
908 std::function<void(const double t,
909 const VectorType &r,
910 VectorType &z,
911 const double tol,
912 const int lr)>
914
939 std::function<void(const double t)> mass_preconditioner_setup;
940
955 std::function<void(const double t,
956 const VectorType &sol,
957 const unsigned int step_number)>
959
977 std::function<bool(const double t, VectorType &sol)> solver_should_restart;
978
985 std::function<VectorType &()> get_local_tolerances;
986
1011 std::function<void(void *arkode_mem)> custom_setup;
1012
1013 private:
1019 std::string,
1020 << "Please provide an implementation for the function \""
1021 << arg1 << "\"");
1022
1026 unsigned int
1027 do_evolve_time(VectorType &solution,
1028 ::DiscreteTime &time,
1029 const bool do_reset);
1030
1037 void
1038 setup_system_solver(const VectorType &solution);
1039
1046 void
1047 setup_mass_solver(const VectorType &solution);
1048
1054 void
1056
1061
1066
1067# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
1072# endif
1073
1079
1084
1085 std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
1086 std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
1087
1093 mutable std::exception_ptr pending_exception;
1094
1095# ifdef DEAL_II_WITH_PETSC
1096# ifdef PETSC_USE_COMPLEX
1097 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
1098 "Sundials does not support complex scalar types, "
1099 "but PETSc is configured to use a complex scalar type!");
1100
1101 static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
1102 "Sundials does not support complex scalar types, "
1103 "but PETSc is configured to use a complex scalar type!");
1104# endif // PETSC_USE_COMPLEX
1105# endif // DEAL_II_WITH_PETSC
1106 };
1107
1108
1113 int,
1114 << "One of the SUNDIALS ARKode internal functions "
1115 << " returned a negative error code: " << arg1
1116 << ". Please consult SUNDIALS manual.");
1117
1118
1119 template <typename VectorType>
1121 // Initial parameters
1122 const double initial_time,
1123 const double final_time,
1124 const double initial_step_size,
1125 const double output_period,
1126 // Running parameters
1127 const double minimum_step_size,
1128 const unsigned int maximum_order,
1129 const unsigned int maximum_non_linear_iterations,
1130 const bool implicit_function_is_linear,
1131 const bool implicit_function_is_time_independent,
1132 const bool mass_is_time_independent,
1133 const int anderson_acceleration_subspace,
1134 // Error parameters
1135 const double absolute_tolerance,
1136 const double relative_tolerance)
1137 : initial_time(initial_time)
1138 , final_time(final_time)
1139 , initial_step_size(initial_step_size)
1140 , minimum_step_size(minimum_step_size)
1141 , absolute_tolerance(absolute_tolerance)
1142 , relative_tolerance(relative_tolerance)
1143 , maximum_order(maximum_order)
1144 , output_period(output_period)
1145 , maximum_non_linear_iterations(maximum_non_linear_iterations)
1146 , implicit_function_is_linear(implicit_function_is_linear)
1147 , implicit_function_is_time_independent(
1148 implicit_function_is_time_independent)
1149 , mass_is_time_independent(mass_is_time_independent)
1150 , anderson_acceleration_subspace(anderson_acceleration_subspace)
1151 {}
1152
1153
1154
1155 template <typename VectorType>
1156 void
1158 {
1159 prm.add_parameter("Initial time", initial_time);
1160 prm.add_parameter("Final time", final_time);
1161 prm.add_parameter("Time interval between each output", output_period);
1162 prm.enter_subsection("Running parameters");
1163 prm.add_parameter("Initial step size", initial_step_size);
1164 prm.add_parameter("Minimum step size", minimum_step_size);
1165 prm.add_parameter("Maximum order of ARK", maximum_order);
1166 prm.add_parameter("Maximum number of nonlinear iterations",
1167 maximum_non_linear_iterations);
1168 prm.add_parameter("Implicit function is linear",
1169 implicit_function_is_linear);
1170 prm.add_parameter("Implicit function is time independent",
1171 implicit_function_is_time_independent);
1172 prm.add_parameter("Mass is time independent", mass_is_time_independent);
1173 prm.add_parameter("Anderson-acceleration subspace",
1174 anderson_acceleration_subspace);
1175 prm.leave_subsection();
1176 prm.enter_subsection("Error control");
1177 prm.add_parameter("Absolute error tolerance", absolute_tolerance);
1178 prm.add_parameter("Relative error tolerance", relative_tolerance);
1179 prm.leave_subsection();
1180 }
1181
1182} // namespace SUNDIALS
1183
1185#endif
1186
1187
1188#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:1157
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:1120
unsigned int maximum_non_linear_iterations
Definition arkode.h:444
void setup_mass_solver(const VectorType &solution)
Definition arkode.cc:602
std::function< void(const double t)> mass_preconditioner_setup
Definition arkode.h:939
std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
Definition arkode.h:590
void * arkode_mem
Definition arkode.h:1065
std::exception_ptr pending_exception
Definition arkode.h:1093
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:829
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:958
std::function< VectorType &()> get_local_tolerances
Definition arkode.h:985
AdditionalData data
Definition arkode.h:1060
LinearSolveFunction< VectorType > solve_linearized_system
Definition arkode.h:764
std::function< void(const double t, const VectorType &y, VectorType &res)> implicit_function
Definition arkode.h:609
std::function< void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> mass_preconditioner_solve
Definition arkode.h:913
MPI_Comm mpi_communicator
Definition arkode.h:1078
void reset(const double t, const double h, const VectorType &y)
Definition arkode.cc:247
SUNContext arkode_ctx
Definition arkode.h:1071
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
Definition arkode.h:1086
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:879
std::function< bool(const double t, VectorType &sol)> solver_should_restart
Definition arkode.h:977
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:665
void set_functions_to_trigger_an_assert()
Definition arkode.cc:742
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
Definition arkode.h:1085
std::function< void(const double t, const VectorType &v, VectorType &Mv)> mass_times_vector
Definition arkode.h:629
LinearSolveFunction< VectorType > solve_mass
Definition arkode.h:788
std::function< void(const double t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
Definition arkode.h:736
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:1083
std::function< void(void *arkode_mem)> custom_setup
Definition arkode.h:1011
std::function< void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
Definition arkode.h:698
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
static ::ExceptionBase & ExcARKodeError(int arg1)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
std::function< void(SundialsOperator< VectorType > &op, SundialsPreconditioner< VectorType > &prec, VectorType &x, const VectorType &b, double tol)> LinearSolveFunction