Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <boost/math/quaternion.hpp>
7 : #include <cstddef>
8 : #include <pup.h>
9 :
10 : #include "ControlSystem/ControlErrors/Expansion.hpp"
11 : #include "ControlSystem/ControlErrors/Rotation.hpp"
12 : #include "ControlSystem/Protocols/ControlError.hpp"
13 : #include "ControlSystem/Tags/QueueTags.hpp"
14 : #include "ControlSystem/Tags/SystemTags.hpp"
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/Matrix.hpp"
17 : #include "DataStructures/TaggedTuple.hpp"
18 : #include "Domain/CoordinateMaps/TimeDependent/RotationMatrixHelpers.hpp"
19 : #include "Domain/Creators/Tags/ObjectCenter.hpp"
20 : #include "Domain/FunctionsOfTime/QuaternionHelpers.hpp"
21 : #include "Domain/Structure/ObjectLabel.hpp"
22 : #include "Options/String.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Utilities/EqualWithinRoundoff.hpp"
25 : #include "Utilities/ErrorHandling/Assert.hpp"
26 : #include "Utilities/GetOutput.hpp"
27 : #include "Utilities/ProtocolHelpers.hpp"
28 : #include "Utilities/TMPL.hpp"
29 :
30 : /// \cond
31 : namespace domain::Tags {
32 : template <size_t Dim>
33 : struct Domain;
34 : struct FunctionsOfTime;
35 : } // namespace domain::Tags
36 : /// \endcond
37 :
38 : namespace control_system {
39 : namespace ControlErrors {
40 : /*!
41 : * \brief Control error in the 3D \link
42 : * domain::CoordinateMaps::TimeDependent::Translation Translation \endlink
43 : * coordinate map
44 : *
45 : * \details Computes the error in how much the system has translated. When there
46 : * are two excisions, it does this by using a modified version of Eq. (42) from
47 : * \cite Ossokine2013zga. The equation is
48 : *
49 : * \begin{equation}
50 : * \left(0, \delta\vec{T}\right) = a\mathbf{q}\left(\frac{1}{2}(\mathbf{x}_A
51 : * + \mathbf{x}_B - \frac{1}{2}(\mathbf{c}_A + \mathbf{c}_B)) - \mathbf{\delta
52 : * q}\wedge\frac{1}{2}(\mathbf{c}_A + \mathbf{c}_B) - \frac{\delta
53 : * a}{a}\frac{1}{2}(\mathbf{c}_A + \mathbf{c}_B) \right)\mathbf{q}^*
54 : * \end{equation}
55 : *
56 : * where object A is located on the positive x-axis in the grid frame, bold face
57 : * letters are quaternions, vectors are promoted to quaternions as \f$
58 : * \mathbf{v} = (0, \vec{v}) \f$, \f$ \mathbf{q} \f$ is the quaternion from the
59 : * \link domain::CoordinateMaps::TimeDependent::Rotation Rotation \endlink map,
60 : * \f$ a \f$ is the function \f$ a(t) \f$ from the \link
61 : * domain::CoordinateMaps::TimeDependent::CubicScale CubicScale \endlink map,
62 : * \f$ \mathbf{\delta q}\wedge\mathbf{c}_A \equiv (0, \delta\vec{q} \times
63 : * \vec{c}_A) \f$, \f$ \delta\vec{q} \f$ is the \link
64 : * control_system::ControlErrors::Rotation Rotation \endlink control error, and
65 : * \f$ \delta a\f$ is the \link control_system::ControlErrors::Expansion
66 : * Expansion \endlink control error.
67 : *
68 : * When there is only a single excision, the control error assumes that the
69 : * center of the excision is at the origin. Thus, the control error is just
70 : * taken to be the current center of the horizon mapped through the expansion
71 : * and rotation maps if there are any.
72 : *
73 : * Requirements:
74 : * - This control error requires that there be either one or two objects in the
75 : * simulation
76 : * - Currently this control error can only be used with the \link
77 : * control_system::Systems::Translation Translation \endlink control system
78 : * - There must exist an expansion map and a quaternion rotation map in the
79 : * coordinate map with names "Expansion" and "Rotation", respectively if there
80 : * are two object. If there is a single object, then the "Expansion" and
81 : * "Rotation" maps may exist, but don't need to.
82 : */
83 : template <size_t NumberOfObjects>
84 1 : struct Translation : tt::ConformsTo<protocols::ControlError> {
85 : static_assert(NumberOfObjects == 1 or NumberOfObjects == 2,
86 : "Translation control can only work with 1 or 2 objects.");
87 :
88 0 : using object_centers = tmpl::conditional_t<
89 : NumberOfObjects == 1, domain::object_list<domain::ObjectLabel::None>,
90 : domain::object_list<domain::ObjectLabel::A, domain::ObjectLabel::B>>;
91 :
92 0 : using options = tmpl::list<>;
93 0 : static constexpr Options::String help{
94 : "Computes the control error for translation control. This should not "
95 : "take any options."};
96 :
97 : // NOLINTNEXTLINE(readability-convert-member-functions-to-static)
98 0 : std::optional<double> get_suggested_timescale() const { return std::nullopt; }
99 :
100 0 : void reset() {}
101 :
102 0 : void pup(PUP::er& /*p*/) {}
103 :
104 : template <typename Metavariables, typename... TupleTags>
105 0 : DataVector operator()(const ::TimescaleTuner<true>& /*unused*/,
106 : const Parallel::GlobalCache<Metavariables>& cache,
107 : const double time,
108 : const std::string& /*function_of_time_name*/,
109 : const tuples::TaggedTuple<TupleTags...>& measurements) {
110 : const auto& functions_of_time = get<domain::Tags::FunctionsOfTime>(cache);
111 :
112 : if constexpr (NumberOfObjects == 2) {
113 : using quat = boost::math::quaternion<double>;
114 :
115 : const quat quaternion = datavector_to_quaternion(
116 : functions_of_time.at("Rotation")->func(time)[0]);
117 : const double expansion_factor =
118 : functions_of_time.at("Expansion")->func(time)[0][0];
119 :
120 : using center_A =
121 : control_system::QueueTags::Center<::domain::ObjectLabel::A,
122 : Frame::Grid>;
123 : using center_B =
124 : control_system::QueueTags::Center<::domain::ObjectLabel::B,
125 : Frame::Grid>;
126 :
127 : const tnsr::I<double, 3, Frame::Grid>& grid_position_of_A_tnsr =
128 : Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::A>>(
129 : cache);
130 : const DataVector grid_position_of_A{{grid_position_of_A_tnsr[0],
131 : grid_position_of_A_tnsr[1],
132 : grid_position_of_A_tnsr[2]}};
133 : const DataVector& current_position_of_A = get<center_A>(measurements);
134 :
135 : const tnsr::I<double, 3, Frame::Grid>& grid_position_of_B_tnsr =
136 : Parallel::get<domain::Tags::ObjectCenter<domain::ObjectLabel::B>>(
137 : cache);
138 : const DataVector grid_position_of_B{{grid_position_of_B_tnsr[0],
139 : grid_position_of_B_tnsr[1],
140 : grid_position_of_B_tnsr[2]}};
141 : const DataVector& current_position_of_B = get<center_B>(measurements);
142 :
143 : const DataVector grid_position_average =
144 : 0.5 * (grid_position_of_A + grid_position_of_B);
145 : const DataVector current_position_average =
146 : 0.5 * (current_position_of_A + current_position_of_B);
147 :
148 : const DataVector grid_separation =
149 : grid_position_of_A - grid_position_of_B;
150 : const DataVector current_separation =
151 : current_position_of_A - current_position_of_B;
152 :
153 : // These quantities come from the translation control implementation in
154 : // SpEC
155 : const double current_separation_dot_grid_separation =
156 : dot(current_separation, grid_separation);
157 : const double current_separation_dot_grid_average =
158 : dot(current_separation, grid_position_average);
159 : const double grid_separation_dot_grid_average =
160 : dot(grid_separation, grid_position_average);
161 : const double grid_separation_dot_grid_separation =
162 : dot(grid_separation, grid_separation);
163 :
164 : // From eq. 42 in 1304.3067 where the grid and current position are
165 : // swapped from only object A to the average grid and current positions of
166 : // both objects.
167 : const DataVector translation_control =
168 : expansion_factor *
169 : (grid_separation_dot_grid_separation * current_position_average -
170 : current_separation_dot_grid_separation * grid_position_average -
171 : grid_separation_dot_grid_average * current_separation +
172 : current_separation_dot_grid_average * grid_separation) /
173 : grid_separation_dot_grid_separation;
174 : const quat middle_expression =
175 : datavector_to_quaternion(translation_control);
176 :
177 : // Because we are converting from a quaternion to a DataVector, there will
178 : // be four components in the DataVector. However, translation control only
179 : // requires three (the latter three to be exact, because the first
180 : // component should be 0. We ASSERT this also.)
181 : const DataVector result_with_four_components = quaternion_to_datavector(
182 : quaternion * middle_expression * conj(quaternion));
183 : ASSERT(equal_within_roundoff(result_with_four_components[0], 0.0),
184 : "Error in computing translation control error. First component of "
185 : "resulting quaternion should be 0.0, but is " +
186 : get_output(result_with_four_components[0]) + " instead.");
187 :
188 : return {result_with_four_components[1], result_with_four_components[2],
189 : result_with_four_components[3]};
190 : } else {
191 : const DataVector& current_position =
192 : get<control_system::QueueTags::Center<::domain::ObjectLabel::None>>(
193 : measurements);
194 : DataVector result{3, 0.0};
195 :
196 : double expansion_factor = 1.0;
197 : if (functions_of_time.count("Expansion") == 1) {
198 : expansion_factor = functions_of_time.at("Expansion")->func(time)[0][0];
199 : }
200 :
201 : if (functions_of_time.count("Rotation") == 1) {
202 : const Matrix rot_matrix =
203 : rotation_matrix<3>(time, *functions_of_time.at("Rotation").get());
204 :
205 : for (size_t i = 0; i < 3; i++) {
206 : for (size_t j = 0; j < 3; j++) {
207 : result[i] += rot_matrix(i, j) * current_position[j];
208 : }
209 : }
210 : } else {
211 : result = current_position;
212 : }
213 :
214 : result *= expansion_factor;
215 :
216 : return result;
217 : }
218 : }
219 : };
220 : } // namespace ControlErrors
221 : } // namespace control_system
|
