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 {}