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 : : boost::true_type {}; 84 : 85 : template <> 86 0 : struct algebra_dispatcher<DataVector> { 87 0 : using algebra_type = ::detail::vector_impl_algebra; 88 : }; 89 : 90 : template <> 91 0 : struct algebra_dispatcher<ComplexDataVector> { 92 0 : using algebra_type = ::detail::vector_impl_algebra; 93 : }; 94 : 95 : template <> 96 0 : struct algebra_dispatcher<ModalVector> { 97 0 : using algebra_type = ::detail::vector_impl_algebra; 98 : }; 99 : 100 : template <> 101 0 : struct algebra_dispatcher<ComplexModalVector> { 102 0 : using algebra_type = ::detail::vector_impl_algebra; 103 : }; 104 : 105 : template <size_t Size> 106 : struct algebra_dispatcher<std::array<DataVector, Size>> { 107 : using algebra_type = ::detail::vector_impl_array_algebra; 108 : }; 109 : 110 : template <size_t Size> 111 : struct algebra_dispatcher<std::array<ComplexDataVector, Size>> { 112 : using algebra_type = ::detail::vector_impl_array_algebra; 113 : }; 114 : 115 : template <size_t Size> 116 : struct algebra_dispatcher<std::array<ModalVector, Size>> { 117 : using algebra_type = ::detail::vector_impl_array_algebra; 118 : }; 119 : 120 : template <size_t Size> 121 : struct algebra_dispatcher<std::array<ComplexModalVector, Size>> { 122 : using algebra_type = ::detail::vector_impl_array_algebra; 123 : }; 124 : } // namespace odeint 125 : } // namespace numeric 126 : } // namespace boost 127 : 128 : /*! 129 : * \ingroup NumericalAlgorithmsGroup 130 : * \brief For ODE integration, we suggest using the boost libraries whenever 131 : * possible. 132 : * 133 : * \details Here we describe briefly the suggested setup to use a boost ODE 134 : * integrator in SpECTRE. The boost utilities have a number of more elaborate 135 : * features, including features that may result in simpler function calls than 136 : * the specific recipes below. For full documentation, see the Boost odeint 137 : * documentation: 138 : * https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/ 139 : * 140 : * The integration methods can be used largely as in the boost documentation. We 141 : * summarize the salient points necessary to get an example running. The main 142 : * steps in using a boost integrator are: 143 : * - define the system 144 : * - construct the stepper 145 : * - initialize (might be performed implicitly, depending on the stepper) 146 : * - perform steps 147 : * - (compute dense output, for dense-output steppers) 148 : * 149 : * For most cases the Dormand-Prince fifth order controlled, dense-output 150 : * stepper is recommended. That stepper is used in the section "SpECTRE vectors 151 : * or `std::array` thereof in boost integrators" below. 152 : * 153 : * #### Fundamental types or `std::array` thereof in fixed-step integrators 154 : * 155 : * Let us consider a simple oscillator system, which we'll declare as a lambda: 156 : * 157 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_array_system 158 : * 159 : * Note that the first argument is a `const` lvalue reference to the state type, 160 : * `std::array<double, 2>`, the second argument is an lvalue reference for the 161 : * time derivatives that is written to by the system function, and the final 162 : * argument is the current time. 163 : * 164 : * The construction and initialization of the stepper is simple: 165 : * 166 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_construction 167 : * 168 : * Finally, we can perform the steps and examine the output, 169 : * 170 : * \snippet Test_OdeIntegration.cpp explicit_fundamental_stepper_use 171 : * 172 : * #### Fundamental types or `std::array` thereof in Dense-output integrators 173 : * 174 : * The dense-output and controlled-step-size ODE integrators in boost comply 175 : * with a somewhat different interface, as significantly more state information 176 : * must be managed. The result is a somewhat simpler, but more opaque, interface 177 : * to the user code. 178 : * 179 : * Once again, we start by constructing the system we'd like to integrate, 180 : * 181 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_system 182 : * 183 : * The constructor for dense steppers takes optional arguments for tolerance, 184 : * and dense steppers need to be initialized. 185 : * 186 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_construction 187 : * 188 : * It is important to take note of the tolerance arguments, as the defaults are 189 : * `1.0e-6`, which is often looser than we want for calculations in SpECTRE. 190 : * 191 : * We then perform the step supplying the system function, and can retrieve the 192 : * dense output state with `calc_state`, which returns by reference by the 193 : * second argument. 194 : * 195 : * \snippet Test_OdeIntegration.cpp dense_output_fundamental_stepper_use 196 : * 197 : * #### SpECTRE vectors or `std::array` thereof in boost integrators 198 : * 199 : * The additional template specializations present in the `OdeIntegration.hpp` 200 : * file ensure that the boost ODE methods work transparently for SpECTRE vectors 201 : * as well. We'll run through a brief example to emphasize the functionality. 202 : * 203 : * \snippet Test_OdeIntegration.cpp dense_output_vector_stepper 204 : */ 205 1 : namespace OdeIntegration {}