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 : // NOLINTNEXTLINE(google-runtime-references) 113 0 : void pup(PUP::er& p) override; 114 : 115 : /// Returns the quaternion at an arbitrary time `t`. 116 1 : std::array<DataVector, 1> func(const double t) const override { 117 : return quat_func(t); 118 : } 119 : 120 : /// Returns the quaternion and its first derivative at an arbitrary time `t`. 121 1 : std::array<DataVector, 2> func_and_deriv(const double t) const override { 122 : return quat_func_and_deriv(t); 123 : } 124 : 125 : /// Returns the quaternion and the first two derivatives at an arbitrary 126 : /// time `t`. 127 1 : std::array<DataVector, 3> func_and_2_derivs(const double t) const override { 128 : return quat_func_and_2_derivs(t); 129 : } 130 : 131 : /// Returns the quaternion at an arbitrary time `t`. 132 1 : std::array<DataVector, 1> quat_func(double t) const; 133 : 134 : /// Returns the quaternion and its first derivative at an arbitrary time `t`. 135 1 : std::array<DataVector, 2> quat_func_and_deriv(double t) const; 136 : 137 : /// Returns the quaternion and the first two derivatives at an arbitrary 138 : /// time `t`. 139 1 : std::array<DataVector, 3> quat_func_and_2_derivs(double t) const; 140 : 141 : /// Returns the quaternion and the first three derivatives at an arbitrary 142 : /// time `t`. 143 1 : std::array<DataVector, 4> quat_func_and_3_derivs(double t) const; 144 : 145 : /// Returns stored angle at an arbitrary time `t`. 146 1 : std::array<DataVector, 1> angle_func(const double t) const { 147 : return angle_f_of_t_.func(t); 148 : } 149 : 150 : /// Returns stored angle and its first derivative (omega) at an arbitrary time 151 : /// `t`. 152 1 : std::array<DataVector, 2> angle_func_and_deriv(const double t) const { 153 : return angle_f_of_t_.func_and_deriv(t); 154 : } 155 : 156 : /// Returns stored angle and the first two derivatives at an arbitrary 157 : /// time `t`. 158 1 : std::array<DataVector, 3> angle_func_and_2_derivs(const double t) const { 159 : return angle_f_of_t_.func_and_2_derivs(t); 160 : } 161 : 162 0 : std::vector<DataVector> angle_func_and_all_derivs(const double t) const { 163 : return angle_f_of_t_.func_and_all_derivs(t); 164 : } 165 : 166 : private: 167 : template <size_t LocalMaxDeriv> 168 0 : friend bool operator==( // NOLINT(readability-redundant-declaration) 169 : const QuaternionFunctionOfTime<LocalMaxDeriv>& lhs, 170 : const QuaternionFunctionOfTime<LocalMaxDeriv>& rhs); 171 : 172 : template <size_t LocalMaxDeriv> 173 0 : friend std::ostream& operator<<( // NOLINT(readability-redundant-declaration) 174 : std::ostream& os, 175 : const QuaternionFunctionOfTime<LocalMaxDeriv>& quaternion_f_of_t); 176 : 177 : FunctionOfTimeHelpers::ThreadsafeList<boost::math::quaternion<double>> 178 0 : stored_quaternions_and_times_{}; 179 0 : domain::FunctionsOfTime::PiecewisePolynomial<MaxDeriv> angle_f_of_t_{}; 180 0 : std::map<double, double> update_backlog_{}; 181 : 182 0 : void unpack_old_version(PUP::er& p, size_t version); 183 : 184 : /// Integrates the ODE \f$ \dot{q} = \frac{1}{2} q \times \omega \f$ from time 185 : /// `t0` to time `t`. On input, `quaternion_to_integrate` is the initial 186 : /// quaternion at time `t0` and on output, it stores the result at time `t` 187 1 : void solve_quaternion_ode( 188 : gsl::not_null<boost::math::quaternion<double>*> quaternion_to_integrate, 189 : double t0, double t) const; 190 : 191 : /// Does common operations to all the `func` functions such as updating stored 192 : /// info, solving the ODE, and returning the normalized quaternion as a boost 193 : /// quaternion for easy calculations 194 1 : boost::math::quaternion<double> setup_func(double t) const; 195 : }; 196 : 197 : template <size_t MaxDeriv> 198 0 : bool operator!=(const QuaternionFunctionOfTime<MaxDeriv>& lhs, 199 : const QuaternionFunctionOfTime<MaxDeriv>& rhs); 200 : 201 : /// \cond 202 : template <size_t MaxDeriv> 203 : PUP::able::PUP_ID QuaternionFunctionOfTime<MaxDeriv>::my_PUP_ID = 0; // NOLINT 204 : /// \endcond 205 : } // namespace domain::FunctionsOfTime