SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/OdeIntegration - OdeIntegration.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 1 13 7.7 %
Date: 2024-03-29 00:33:31
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14