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             :  * \note The initial rotation angles passed to the angle PiecewisePolynomial
      46             :  * don't matter as we never actually use the angles themselves. We only use
      47             :  * their derivatives (angular velocity) to determine map parameters. In theory
      48             :  * we could determine each initial angle from the input axis-angle
      49             :  * representation, but we don't need to.
      50             :  *
      51             :  * The angle PiecewisePolynomial is accessible through the `angle_func`,
      52             :  * `angle_func_and_deriv`, and `angle_func_and_2_derivs` functions which
      53             :  * correspond to the function calls of a normal PiecewisePolynomial except
      54             :  * without the `angle_` prefix.
      55             :  *
      56             :  * It is encouraged to use `quat_func` and `angle_func` when you want the
      57             :  * specific values of the functions to avoid ambiguity in what you are
      58             :  * calling. However, the original three `func` functions inherited from the
      59             :  * FunctionOfTime base class are necessary because the maps use the generic
      60             :  * `func` functions, thus they return the quaternion and its derivatives (which
      61             :  * are needed for the map). This is all to keep the symmetry of naming
      62             :  * `angle_func` and `quat_func` so that function calls won't be ambiguous.
      63             :  *
      64             :  * \note This class conforms to the requirements of the
      65             :  * `Parallel::GlobalCache` for objects held by mutable global cache tags.
      66             :  */
      67             : template <size_t MaxDeriv>
      68           1 : class QuaternionFunctionOfTime : public FunctionOfTime {
      69             :  public:
      70           0 :   QuaternionFunctionOfTime();
      71           0 :   QuaternionFunctionOfTime(QuaternionFunctionOfTime&&);
      72           0 :   QuaternionFunctionOfTime(const QuaternionFunctionOfTime&);
      73           0 :   QuaternionFunctionOfTime& operator=(QuaternionFunctionOfTime&&);
      74           0 :   QuaternionFunctionOfTime& operator=(const QuaternionFunctionOfTime&);
      75           0 :   ~QuaternionFunctionOfTime() override;
      76             : 
      77           0 :   QuaternionFunctionOfTime(
      78             :       double t, const std::array<DataVector, 1>& initial_quat_func,
      79             :       std::array<DataVector, MaxDeriv + 1> initial_angle_func,
      80             :       double expiration_time);
      81             : 
      82             :   // LCOV_EXCL_START
      83           0 :   explicit QuaternionFunctionOfTime(CkMigrateMessage* /*unused*/);
      84             :   // LCOV_EXCL_STOP
      85             : 
      86           0 :   auto get_clone() const -> std::unique_ptr<FunctionOfTime> override;
      87             : 
      88           1 :   std::unique_ptr<FunctionOfTime> create_at_time(
      89             :       double t, double expiration_time) const override;
      90             : 
      91             :   // clang-tidy: google-runtime-references
      92             :   // clang-tidy: cppcoreguidelines-owning-memory,-warnings-as-errors
      93           0 :   WRAPPED_PUPable_decl_template(QuaternionFunctionOfTime<MaxDeriv>);  // NOLINT
      94             : 
      95             :   /// Returns domain of validity for the function of time
      96           1 :   std::array<double, 2> time_bounds() const override;
      97             : 
      98           1 :   double expiration_after(double time) const override;
      99             : 
     100             :   /// Updates the `MaxDeriv`th derivative of the angle piecewisepolynomial at
     101             :   /// the given time, then updates the stored quaternions.
     102             :   ///
     103             :   /// `updated_max_deriv` is a datavector of the `MaxDeriv`s for each component.
     104             :   /// `next_expiration_time` is the next expiration time.
     105             :   ///
     106             :   /// The \p time_of_update must be the same as the old expiration
     107             :   /// time.  It is passed as a check that the calling code is
     108             :   /// computing the other arguments with the correct value.
     109           1 :   void update(double time_of_update, DataVector updated_max_deriv,
     110             :               double next_expiration_time) override;
     111             : 
     112           1 :   void truncate_at_time(double time) override;
     113             : 
     114             :   // NOLINTNEXTLINE(google-runtime-references)
     115           0 :   void pup(PUP::er& p) override;
     116             : 
     117             :   /// Returns the quaternion at an arbitrary time `t`.
     118           1 :   std::array<DataVector, 1> func(const double t) const override {
     119             :     return quat_func(t);
     120             :   }
     121             : 
     122             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     123           1 :   std::array<DataVector, 2> func_and_deriv(const double t) const override {
     124             :     return quat_func_and_deriv(t);
     125             :   }
     126             : 
     127             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     128             :   /// time `t`.
     129           1 :   std::array<DataVector, 3> func_and_2_derivs(const double t) const override {
     130             :     return quat_func_and_2_derivs(t);
     131             :   }
     132             : 
     133             :   /// Returns the quaternion at an arbitrary time `t`.
     134           1 :   std::array<DataVector, 1> quat_func(double t) const;
     135             : 
     136             :   /// Returns the quaternion and its first derivative at an arbitrary time `t`.
     137           1 :   std::array<DataVector, 2> quat_func_and_deriv(double t) const;
     138             : 
     139             :   /// Returns the quaternion and the first two derivatives at an arbitrary
     140             :   /// time `t`.
     141           1 :   std::array<DataVector, 3> quat_func_and_2_derivs(double t) const;
     142             : 
     143             :   /// Returns the quaternion and the first three derivatives at an arbitrary
     144             :   /// time `t`.
     145           1 :   std::array<DataVector, 4> quat_func_and_3_derivs(double t) const;
     146             : 
     147             :   /// Returns stored angle at an arbitrary time `t`.
     148           1 :   std::array<DataVector, 1> angle_func(const double t) const {
     149             :     return angle_f_of_t_.func(t);
     150             :   }
     151             : 
     152             :   /// Returns stored angle and its first derivative (omega) at an arbitrary time
     153             :   /// `t`.
     154           1 :   std::array<DataVector, 2> angle_func_and_deriv(const double t) const {
     155             :     return angle_f_of_t_.func_and_deriv(t);
     156             :   }
     157             : 
     158             :   /// Returns stored angle and the first two derivatives at an arbitrary
     159             :   /// time `t`.
     160           1 :   std::array<DataVector, 3> angle_func_and_2_derivs(const double t) const {
     161             :     return angle_f_of_t_.func_and_2_derivs(t);
     162             :   }
     163             : 
     164           0 :   std::vector<DataVector> angle_func_and_all_derivs(const double t) const {
     165             :     return angle_f_of_t_.func_and_all_derivs(t);
     166             :   }
     167             : 
     168             :  private:
     169             :   template <size_t LocalMaxDeriv>
     170           0 :   friend bool operator==(  // NOLINT(readability-redundant-declaration)
     171             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& lhs,
     172             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& rhs);
     173             : 
     174             :   template <size_t LocalMaxDeriv>
     175           0 :   friend std::ostream& operator<<(  // NOLINT(readability-redundant-declaration)
     176             :       std::ostream& os,
     177             :       const QuaternionFunctionOfTime<LocalMaxDeriv>& quaternion_f_of_t);
     178             : 
     179             :   FunctionOfTimeHelpers::ThreadsafeList<boost::math::quaternion<double>>
     180           0 :       stored_quaternions_and_times_{};
     181           0 :   domain::FunctionsOfTime::PiecewisePolynomial<MaxDeriv> angle_f_of_t_{};
     182           0 :   std::map<double, double> update_backlog_{};
     183             : 
     184           0 :   void unpack_old_version(PUP::er& p, size_t version);
     185             : 
     186             :   /// Integrates the ODE \f$ \dot{q} = \frac{1}{2} q \times \omega \f$ from time
     187             :   /// `t0` to time `t`. On input, `quaternion_to_integrate` is the initial
     188             :   /// quaternion at time `t0` and on output, it stores the result at time `t`
     189           1 :   void solve_quaternion_ode(
     190             :       gsl::not_null<boost::math::quaternion<double>*> quaternion_to_integrate,
     191             :       double t0, double t) const;
     192             : 
     193             :   /// Does common operations to all the `func` functions such as updating stored
     194             :   /// info, solving the ODE, and returning the normalized quaternion as a boost
     195             :   /// quaternion for easy calculations
     196           1 :   boost::math::quaternion<double> setup_func(double t) const;
     197             : };
     198             : 
     199             : template <size_t MaxDeriv>
     200           0 : bool operator!=(const QuaternionFunctionOfTime<MaxDeriv>& lhs,
     201             :                 const QuaternionFunctionOfTime<MaxDeriv>& rhs);
     202             : 
     203             : /// \cond
     204             : template <size_t MaxDeriv>
     205             : PUP::able::PUP_ID QuaternionFunctionOfTime<MaxDeriv>::my_PUP_ID = 0;  // NOLINT
     206             : /// \endcond
     207             : }  // namespace domain::FunctionsOfTime