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
|
