Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <boost/math/quaternion.hpp>
       8             : #include <cstddef>
       9             : #include <map>
      10             : #include <memory>
      11             : #include <ostream>
      12             : #include <pup.h>
      13             : #include <vector>
      14             : 
      15             : #include "DataStructures/DataVector.hpp"
      16             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      17             : #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
      18             : #include "Domain/FunctionsOfTime/ThreadsafeList.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : #include "Utilities/Serialization/CharmPupable.hpp"
      21             : 
      22             : namespace domain::FunctionsOfTime {
      23             : /*!
      24             :  * \ingroup ComputationalDomainGroup
      25             :  * \brief A FunctionOfTime that stores quaternions for the rotation map
      26             :  *
      27             :  * \details This FunctionOfTime stores quaternions that will be used in the
      28             :  * time-dependent rotation map as well as the orbital angular velocity that
      29             :  * will be controlled by the rotation control sytem. To get the quaternion, an
      30             :  * ODE is solved of the form \f$ \dot{q} = \frac{1}{2} q \times \omega \f$
      31             :  * where \f$ \omega \f$ is the orbital angular velocity which is stored
      32             :  * internally as the derivative of an angle `PiecewisePolynomial`, and
      33             :  * \f$ \times \f$ here is quaternion multiplication.
      34             :  *
      35             :  * Different from a `PiecewisePolynomial`, only the quaternion
      36             :  * itself is stored, not any of the derivatives because the derivatives must be
      37             :  * calculated from the solved ODE at every function call. Because
      38             :  * derivatives of the quaternion are not stored, the template parameter
      39             :  * `MaxDeriv` refers to both the max derivative of the stored angle
      40             :  * PiecewisePolynomial and the max derivative returned by the
      41             :  * QuaternionFunctionOfTime. The `update` function is then just a wrapper
      42             :  * around the internal `PiecewisePolynomial::update` function with the addition
      43             :  * that it then updates the stored quaternions as well.
      44             :  *
      45             :  * The angle PiecewisePolynomial is accessible through the `angle_func`,
      46             :  * `angle_func_and_deriv`, and `angle_func_and_2_derivs` functions which
      47             :  * correspond to the function calls of a normal PiecewisePolynomial except
      48             :  * without the `angle_` prefix.
      49             :  *
      50             :  * It is encouraged to use `quat_func` and `angle_func` when you want the
      51             :  * specific values of the functions to avoid ambiguity in what you are
      52             :  * calling. However, the original three `func` functions inherited from the
      53             :  * FunctionOfTime base class are necessary because the maps use the generic
      54             :  * `func` functions, thus they return the quaternion and its derivatives (which
      55             :  * are needed for the map). This is all to keep the symmetry of naming
      56             :  * `angle_func` and `quat_func` so that function calls won't be ambiguous.
      57             :  *
      58             :  * \note This class conforms to the requirements of the
      59             :  * `Parallel::GlobalCache` for objects held by mutable global cache tags.
      60             :  */
      61             : template <size_t MaxDeriv>
      62           1 : class QuaternionFunctionOfTime : public FunctionOfTime {
      63             :  public:
      64           0 :   QuaternionFunctionOfTime();
      65           0 :   QuaternionFunctionOfTime(QuaternionFunctionOfTime&&);
      66           0 :   QuaternionFunctionOfTime(const QuaternionFunctionOfTime&);
      67           0 :   QuaternionFunctionOfTime& operator=(QuaternionFunctionOfTime&&);
      68           0 :   QuaternionFunctionOfTime& operator=(const QuaternionFunctionOfTime&);
      69           0 :   ~QuaternionFunctionOfTime() override;
      70             : 
      71           0 :   QuaternionFunctionOfTime(
      72             :       double t, const std::array<DataVector, 1>& initial_quat_func,
      73             :       std::array<DataVector, MaxDeriv + 1> initial_angle_func,
      74             :       double expiration_time);
      75             : 
      76             :   // LCOV_EXCL_START
      77           0 :   explicit QuaternionFunctionOfTime(CkMigrateMessage* /*unused*/);
      78             :   // LCOV_EXCL_STOP
      79             : 
      80           0 :   auto get_clone() const -> std::unique_ptr<FunctionOfTime> override;
      81             : 
      82             :   // clang-tidy: google-runtime-references
      83             :   // clang-tidy: cppcoreguidelines-owning-memory,-warnings-as-errors
      84           0 :   WRAPPED_PUPable_decl_template(QuaternionFunctionOfTime<MaxDeriv>);  // NOLINT
      85             : 
      86             :   /// Returns domain of validity for the function of time
      87           1 :   std::array<double, 2> time_bounds() const override;
      88             : 
      89           1 :   double expiration_after(double time) const override;
      90             : 
      91             :   /// Updates the `MaxDeriv`th derivative of the angle piecewisepolynomial at
      92             :   /// the given time, then updates the stored quaternions.
      93             :   ///
      94             :   /// `updated_max_deriv` is a datavector of the `MaxDeriv`s for each component.
      95             :   /// `next_expiration_time` is the next expiration time.
      96             :   ///
      97             :   /// The \p time_of_update must be the same as the old expiration
      98             :   /// time.  It is passed as a check that the calling code is
      99             :   /// computing the other arguments with the correct value.
     100           1 :   void update(double time_of_update, DataVector updated_max_deriv,
     101             :               double next_expiration_time) override;
     102             : 
     103             :   // NOLINTNEXTLINE(google-runtime-references)
     104           0 :   void pup(PUP::er& p) override;
     105             : 
     106             :   /// Returns the quaternion at an arbitrary time `t`.
     107           1 :   std::array<DataVector, 1> func(const double t) const override {
     108             :     return quat_func(t);
     109             :   }
     110             : 
     111             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     112           1 :   std::array<DataVector, 2> func_and_deriv(const double t) const override {
     113             :     return quat_func_and_deriv(t);
     114             :   }
     115             : 
     116             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     117             :   /// time `t`.
     118           1 :   std::array<DataVector, 3> func_and_2_derivs(const double t) const override {
     119             :     return quat_func_and_2_derivs(t);
     120             :   }
     121             : 
     122             :   /// Returns the quaternion at an arbitrary time `t`.
     123           1 :   std::array<DataVector, 1> quat_func(double t) const;
     124             : 
     125             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     126           1 :   std::array<DataVector, 2> quat_func_and_deriv(double t) const;
     127             : 
     128             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     129             :   /// time `t`.
     130           1 :   std::array<DataVector, 3> quat_func_and_2_derivs(double t) const;
     131             : 
     132             :   /// Returns the quaternion and the first three derivatives at an arbitrary
     133             :   /// time `t`.
     134           1 :   std::array<DataVector, 4> quat_func_and_3_derivs(double t) const;
     135             : 
     136             :   /// Returns stored angle at an arbitrary time `t`.
     137           1 :   std::array<DataVector, 1> angle_func(const double t) const {
     138             :     return angle_f_of_t_.func(t);
     139             :   }
     140             : 
     141             :   /// Returns stored angle and its first derivative (omega) at an arbitrary time
     142             :   /// `t`.
     143           1 :   std::array<DataVector, 2> angle_func_and_deriv(const double t) const {
     144             :     return angle_f_of_t_.func_and_deriv(t);
     145             :   }
     146             : 
     147             :   /// Returns stored angle and the first two derivatives at an arbitrary
     148             :   /// time `t`.
     149           1 :   std::array<DataVector, 3> angle_func_and_2_derivs(const double t) const {
     150             :     return angle_f_of_t_.func_and_2_derivs(t);
     151             :   }
     152             : 
     153           0 :   std::vector<DataVector> angle_func_and_all_derivs(const double t) const {
     154             :     return angle_f_of_t_.func_and_all_derivs(t);
     155             :   }
     156             : 
     157             :  private:
     158             :   template <size_t LocalMaxDeriv>
     159           0 :   friend bool operator==(  // NOLINT(readability-redundant-declaration)
     160             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& lhs,
     161             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& rhs);
     162             : 
     163             :   template <size_t LocalMaxDeriv>
     164           0 :   friend std::ostream& operator<<(  // NOLINT(readability-redundant-declaration)
     165             :       std::ostream& os,
     166             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& quaternion_f_of_t);
     167             : 
     168             :   FunctionOfTimeHelpers::ThreadsafeList<boost::math::quaternion<double>>
     169           0 :       stored_quaternions_and_times_{};
     170           0 :   domain::FunctionsOfTime::PiecewisePolynomial<MaxDeriv> angle_f_of_t_{};
     171           0 :   std::map<double, double> update_backlog_{};
     172             : 
     173           0 :   void unpack_old_version(PUP::er& p, size_t version);
     174             : 
     175             :   /// Integrates the ODE \f$ \dot{q} = \frac{1}{2} q \times \omega \f$ from time
     176             :   /// `t0` to time `t`. On input, `quaternion_to_integrate` is the initial
     177             :   /// quaternion at time `t0` and on output, it stores the result at time `t`
     178           1 :   void solve_quaternion_ode(
     179             :       gsl::not_null<boost::math::quaternion<double>*> quaternion_to_integrate,
     180             :       double t0, double t) const;
     181             : 
     182             :   /// Does common operations to all the `func` functions such as updating stored
     183             :   /// info, solving the ODE, and returning the normalized quaternion as a boost
     184             :   /// quaternion for easy calculations
     185           1 :   boost::math::quaternion<double> setup_func(double t) const;
     186             : };
     187             : 
     188             : template <size_t MaxDeriv>
     189           0 : bool operator!=(const QuaternionFunctionOfTime<MaxDeriv>& lhs,
     190             :                 const QuaternionFunctionOfTime<MaxDeriv>& rhs);
     191             : 
     192             : /// \cond
     193             : template <size_t MaxDeriv>
     194             : PUP::able::PUP_ID QuaternionFunctionOfTime<MaxDeriv>::my_PUP_ID = 0;  // NOLINT
     195             : /// \endcond
     196             : }  // namespace domain::FunctionsOfTime