OdeIntegration.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <array>
8 #include <complex>
9 #include <cstddef>
10 
11 #include "DataStructures/BoostMultiArray.hpp"
12 #include "DataStructures/ComplexDataVector.hpp"
13 #include "DataStructures/ComplexModalVector.hpp"
14 #include "DataStructures/DataVector.hpp"
15 #include "DataStructures/ModalVector.hpp"
16 #include "DataStructures/VectorImpl.hpp"
17 
18 // Boost MultiArray is used internally in odeint, so odeint must be included
19 // later
20 #include <boost/numeric/odeint.hpp>
21 
22 namespace detail {
23 struct vector_impl_algebra : boost::numeric::odeint::vector_space_algebra {
24  template <typename VectorType>
25  static double norm_inf(const VectorType& vector) noexcept {
26  return max(abs(vector));
27  }
28 };
29 
30 struct vector_impl_array_algebra : boost::numeric::odeint::array_algebra {
31  template <typename VectorType, size_t Size>
32  static double norm_inf(
33  const std::array<VectorType, Size>& vector_array) noexcept {
34  return std::accumulate(
35  vector_array.begin(), vector_array.end(), 0.0,
36  [](const double maximum_so_far, const VectorType& element) noexcept {
37  return std::max(maximum_so_far, max(abs(element)));
38  });
39  }
40 };
41 } // namespace detail
42 
43 namespace boost {
44 namespace numeric {
45 namespace odeint {
46 
47 namespace detail {
48 template <typename V>
49 struct set_unit_value_impl<ComplexDataVector, V, void> {
50  static void set_value(ComplexDataVector& t, const V& v) noexcept { // NOLINT
51  // this ensures that the blaze expression template resolve the result to a
52  // complex vector
53  t = std::complex<double>(1.0, 0.0) * v;
54  }
55 };
56 template <typename V>
57 struct set_unit_value_impl<ComplexModalVector, V, void> {
58  static void set_value(ComplexModalVector& t, const V& v) noexcept { // NOLINT
59  // this ensures that the blaze expression template resolve the result to a
60  // complex vector
61  t = std::complex<double>(1.0, 0.0) * v;
62  }
63 };
64 } // namespace detail
65 
66 template <>
67 struct algebra_dispatcher<DataVector> {
68  using algebra_type = ::detail::vector_impl_algebra;
69 };
70 
71 template <>
72 struct algebra_dispatcher<ComplexDataVector> {
73  using algebra_type = ::detail::vector_impl_algebra;
74 };
75 
76 template <>
77 struct algebra_dispatcher<ModalVector> {
78  using algebra_type = ::detail::vector_impl_algebra;
79 };
80 
81 template <>
82 struct algebra_dispatcher<ComplexModalVector> {
83  using algebra_type = ::detail::vector_impl_algebra;
84 };
85 
86 template <size_t Size>
87 struct algebra_dispatcher<std::array<DataVector, Size>> {
88  using algebra_type = ::detail::vector_impl_array_algebra;
89 };
90 
91 template <size_t Size>
92 struct algebra_dispatcher<std::array<ComplexDataVector, Size>> {
93  using algebra_type = ::detail::vector_impl_array_algebra;
94 };
95 
96 template <size_t Size>
97 struct algebra_dispatcher<std::array<ModalVector, Size>> {
98  using algebra_type = ::detail::vector_impl_array_algebra;
99 };
100 
101 template <size_t Size>
102 struct algebra_dispatcher<std::array<ComplexModalVector, Size>> {
103  using algebra_type = ::detail::vector_impl_array_algebra;
104 };
105 } // namespace odeint
106 } // namespace numeric
107 } // namespace boost
108 
109 /*!
110  * \ingroup NumericalAlgorithmsGroup
111  * \brief For ODE integration, we suggest using the boost libraries whenever
112  * possible.
113  *
114  * \details Here we describe briefly the suggested setup to use a boost ODE
115  * integrator in SpECTRE. The boost utilities have a number of more elaborate
116  * features, including features that may result in simpler function calls than
117  * the specific recipes below. For full documentation, see the Boost odeint
118  * documentation:
119  * https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/
120  *
121  * The integration methods can be used largely as in the boost documentation. We
122  * summarize the salient points necessary to get an example running. The main
123  * steps in using a boost integrator are:
124  * - define the system
125  * - construct the stepper
126  * - initialize (might be performed implicitly, depending on the stepper)
127  * - perform steps
128  * - (compute dense output, for dense-output steppers)
129  *
130  * For most cases the Dormand-Prince fifth order controlled, dense-output
131  * stepper is recommended. That stepper is used in the section "SpECTRE vectors
132  * or `std::array` thereof in boost integrators" below.
133  *
134  * #### Fundamental types or `std::array` thereof in fixed-step integrators
135  *
136  * Let us consider a simple oscillator system, which we'll declare as a lambda:
137  *
138  * \snippet Test_OdeIntegration.cpp explicit_fundamental_array_system
139  *
140  * Note that the first argument is a `const` lvalue reference to the state type,
141  * `std::array<double, 2>`, the second argument is an lvalue reference for the
142  * time derivatives that is written to by the system function, and the final
143  * argument is the current time.
144  *
145  * The construction and initialization of the stepper is simple:
146  *
147  * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_construction
148  *
149  * Finally, we can perform the steps and examine the output,
150  *
151  * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_use
152  *
153  * #### Fundamental types or `std::array` thereof in Dense-output integrators
154  *
155  * The dense-output and controlled-step-size ODE integrators in boost comply
156  * with a somewhat different interface, as significantly more state information
157  * must be managed. The result is a somewhat simpler, but more opaque, interface
158  * to the user code.
159  *
160  * Once again, we start by constructing the system we'd like to integrate,
161  *
162  * \snippet Test_OdeIntegration.cpp dense_output_fundamental_system
163  *
164  * The constructor for dense steppers takes optional arguments for tolerance,
165  * and dense steppers need to be initialized.
166  *
167  * \snippet Test_OdeIntegration.cpp dense_output_fundamental_construction
168  *
169  * It is important to take note of the tolerance arguments, as the defaults are
170  * `1.0e-6`, which is often looser than we want for calculations in SpECTRE.
171  *
172  * We then perform the step supplying the system function, and can retrieve the
173  * dense output state with `calc_state`, which returns by reference by the
174  * second argument.
175  *
176  * \snippet Test_OdeIntegration.cpp dense_output_fundamental_stepper_use
177  *
178  * #### SpECTRE vectors or `std::array` thereof in boost integrators
179  *
180  * The additional template specializations present in the `OdeIntegration.hpp`
181  * file ensure that the boost ODE methods work transparently for SpECTRE vectors
182  * as well. We'll run through a brief example to emphasize the functionality.
183  *
184  * \snippet Test_OdeIntegration.cpp dense_output_vector_stepper
185  */
186 namespace OdeIntegration {}
algorithm
ModalVector
A class for storing spectral coefficients on a spectral grid.
Definition: ModalVector.hpp:30
cstddef
OdeIntegration
For ODE integration, we suggest using the boost libraries whenever possible.
Definition: OdeIntegration.hpp:186
array
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
cpp2b::accumulate
constexpr T accumulate(InputIt first, InputIt last, T init)
Definition: Numeric.hpp:33
ComplexModalVector
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
complex
std::max
T max(T... args)
ComplexDataVector
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47