Line data Source code
1 0 : // 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) {
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(const std::array<VectorType, Size>& vector_array) {
33 : return std::accumulate(
34 : vector_array.begin(), vector_array.end(), 0.0,
35 : [](const double maximum_so_far, const VectorType& element) {
36 : return std::max(maximum_so_far, max(abs(element)));
37 : });
38 : }
39 : };
40 : } // namespace detail
41 :
42 : namespace boost {
43 : namespace numeric {
44 : namespace odeint {
45 :
46 : namespace detail {
47 : template <typename V>
48 : struct set_unit_value_impl<ComplexDataVector, V, void> {
49 : static void set_value(ComplexDataVector& t, const V& v) { // NOLINT
50 : // this ensures that the blaze expression template resolve the result to a
51 : // complex vector
52 : t = std::complex<double>(1.0, 0.0) * v;
53 : }
54 : };
55 : template <typename V>
56 : struct set_unit_value_impl<ComplexModalVector, V, void> {
57 : static void set_value(ComplexModalVector& t, const V& v) { // NOLINT
58 : // this ensures that the blaze expression template resolve the result to a
59 : // complex vector
60 : t = std::complex<double>(1.0, 0.0) * v;
61 : }
62 : };
63 : } // namespace detail
64 :
65 : // In some integration contexts, boost requires the ability to resize the
66 : // integration arguments in preparation for writing to output buffers passed by
67 : // reference. These specializations make the resize work correctly with spectre
68 : // vector types.
69 : template <class VectorType1, class VectorType2>
70 0 : struct resize_impl_sfinae<VectorType1, VectorType2,
71 : typename boost::enable_if_c<
72 : is_derived_of_vector_impl_v<VectorType1> and
73 : is_derived_of_vector_impl_v<VectorType2>>::type> {
74 0 : static void resize(VectorType1& x1, const VectorType2& x2) {
75 : x1.destructive_resize(x2.size());
76 : }
77 : };
78 :
79 : template <class VectorType>
80 0 : struct is_resizeable_sfinae<
81 : VectorType,
82 : typename boost::enable_if_c<is_derived_of_vector_impl_v<VectorType>>::type>
83 : :
84 : #if BOOST_VERSION >= 108600
85 : std::true_type
86 : #else
87 : boost::true_type
88 : #endif
89 : {
90 : };
91 :
92 : template <>
93 0 : struct algebra_dispatcher<DataVector> {
94 0 : using algebra_type = ::detail::vector_impl_algebra;
95 : };
96 :
97 : template <>
98 0 : struct algebra_dispatcher<ComplexDataVector> {
99 0 : using algebra_type = ::detail::vector_impl_algebra;
100 : };
101 :
102 : template <>
103 0 : struct algebra_dispatcher<ModalVector> {
104 0 : using algebra_type = ::detail::vector_impl_algebra;
105 : };
106 :
107 : template <>
108 0 : struct algebra_dispatcher<ComplexModalVector> {
109 0 : using algebra_type = ::detail::vector_impl_algebra;
110 : };
111 :
112 : template <size_t Size>
113 : struct algebra_dispatcher<std::array<DataVector, Size>> {
114 : using algebra_type = ::detail::vector_impl_array_algebra;
115 : };
116 :
117 : template <size_t Size>
118 : struct algebra_dispatcher<std::array<ComplexDataVector, Size>> {
119 : using algebra_type = ::detail::vector_impl_array_algebra;
120 : };
121 :
122 : template <size_t Size>
123 : struct algebra_dispatcher<std::array<ModalVector, Size>> {
124 : using algebra_type = ::detail::vector_impl_array_algebra;
125 : };
126 :
127 : template <size_t Size>
128 : struct algebra_dispatcher<std::array<ComplexModalVector, Size>> {
129 : using algebra_type = ::detail::vector_impl_array_algebra;
130 : };
131 : } // namespace odeint
132 : } // namespace numeric
133 : } // namespace boost
134 :
135 : /*!
136 : * \ingroup NumericalAlgorithmsGroup
137 : * \brief For ODE integration, we suggest using the boost libraries whenever
138 : * possible.
139 : *
140 : * \details Here we describe briefly the suggested setup to use a boost ODE
141 : * integrator in SpECTRE. The boost utilities have a number of more elaborate
142 : * features, including features that may result in simpler function calls than
143 : * the specific recipes below. For full documentation, see the Boost odeint
144 : * documentation:
145 : * https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/
146 : *
147 : * The integration methods can be used largely as in the boost documentation. We
148 : * summarize the salient points necessary to get an example running. The main
149 : * steps in using a boost integrator are:
150 : * - define the system
151 : * - construct the stepper
152 : * - initialize (might be performed implicitly, depending on the stepper)
153 : * - perform steps
154 : * - (compute dense output, for dense-output steppers)
155 : *
156 : * For most cases the Dormand-Prince fifth order controlled, dense-output
157 : * stepper is recommended. That stepper is used in the section "SpECTRE vectors
158 : * or `std::array` thereof in boost integrators" below.
159 : *
160 : * #### Fundamental types or `std::array` thereof in fixed-step integrators
161 : *
162 : * Let us consider a simple oscillator system, which we'll declare as a lambda:
163 : *
164 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_array_system
165 : *
166 : * Note that the first argument is a `const` lvalue reference to the state type,
167 : * `std::array<double, 2>`, the second argument is an lvalue reference for the
168 : * time derivatives that is written to by the system function, and the final
169 : * argument is the current time.
170 : *
171 : * The construction and initialization of the stepper is simple:
172 : *
173 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_construction
174 : *
175 : * Finally, we can perform the steps and examine the output,
176 : *
177 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_use
178 : *
179 : * #### Fundamental types or `std::array` thereof in Dense-output integrators
180 : *
181 : * The dense-output and controlled-step-size ODE integrators in boost comply
182 : * with a somewhat different interface, as significantly more state information
183 : * must be managed. The result is a somewhat simpler, but more opaque, interface
184 : * to the user code.
185 : *
186 : * Once again, we start by constructing the system we'd like to integrate,
187 : *
188 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_system
189 : *
190 : * The constructor for dense steppers takes optional arguments for tolerance,
191 : * and dense steppers need to be initialized.
192 : *
193 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_construction
194 : *
195 : * It is important to take note of the tolerance arguments, as the defaults are
196 : * `1.0e-6`, which is often looser than we want for calculations in SpECTRE.
197 : *
198 : * We then perform the step supplying the system function, and can retrieve the
199 : * dense output state with `calc_state`, which returns by reference by the
200 : * second argument.
201 : *
202 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_stepper_use
203 : *
204 : * #### SpECTRE vectors or `std::array` thereof in boost integrators
205 : *
206 : * The additional template specializations present in the `OdeIntegration.hpp`
207 : * file ensure that the boost ODE methods work transparently for SpECTRE vectors
208 : * as well. We'll run through a brief example to emphasize the functionality.
209 : *
210 : * \snippet Test_OdeIntegration.cpp dense_output_vector_stepper
211 : */
212 1 : namespace OdeIntegration {}
|
