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 <cstddef>
8 : #include <memory>
9 : #include <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <unordered_set>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "PointwiseFunctions/MathFunctions/MathFunction.hpp"
16 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
17 :
18 : /// \cond
19 : namespace domain::FunctionsOfTime {
20 : class FunctionOfTime;
21 : } // namespace domain::FunctionsOfTime
22 : namespace PUP {
23 : class er;
24 : } // namespace PUP
25 : /// \endcond
26 :
27 : namespace domain::CoordinateMaps::TimeDependent {
28 : /*!
29 : * \ingroup CoordMapsTimeDependentGroup
30 : * \brief RotScaleTrans map which applies a combination of rotation, expansion,
31 : * and translation based on which maps are supplied.
32 : *
33 : * \details This map adds a rotation, expansion and translation based on what
34 : * types of maps are needed. Translation and expansion have piecewise functions
35 : * that map the coordinates $\vec{\xi}$ based on what region
36 : * $|\vec{\xi}|$ is in. Coordinates within the inner radius are translated by
37 : * the translation function of time $\vec{T}(t)$ and expanded by the inner
38 : * expansion function of time $E_{a}(t)$. Coordinates in between the inner
39 : * radius and outer radius have a linear radial falloff applied to them.
40 : * Coordinates at or beyond the outer radius have no translation applied and are
41 : * expanded by the outer expansion function of time $E_{b}(t)$. This map assumes
42 : * that the center of the map is at (0., 0., 0.). There is an enum class to set
43 : * which region your block is in. Specifying RotScaleTrans::BlockRegion::Inner
44 : * treats coordinates as if they're inside the inner radius, setting
45 : * RotScaleTrans::BlockRegion::Transition treats the points as if they're
46 : * between the inner and outer radius, and setting
47 : * RotScaleTrans::BlockRegion::Outer treats points as if they're on or beyond
48 : * the outer boundary.
49 : *
50 : * \note $E_{a}(t)$ and $E_{b}(t)$ are $a(t)$ and $b(t)$ from the CubicScale
51 : * documentation.
52 : *
53 : * ## Mapped Coordinates
54 : * The RotScaleTrans map takes the coordinates $\vec{\xi}$ to the target
55 : * coordinates $\vec{\bar{\xi}}$ through
56 : * \f{equation}{
57 : * \vec{\bar{\xi}} = \left\{\begin{array}{ll}E_{a}(t)R(t)\vec{\xi} + \vec{T}(t),
58 : * & |\vec{\xi}| \leq R_{in}, \\ (w_{E} + E_{a}(t))R(t)\vec{\xi} + (1 +
59 : * w_{T})\vec{T}(t), & R_{in} < |\vec{\xi}| \leq 0.5(R_{in} + R_{out}),
60 : * \\ (w_{E} + E_{b}(t))R(t)\vec{\xi} + w_{T}\vec{T}(t), & 0.5(R_{in} + R_{out})
61 : * < |\vec{\xi}| < R_{out}, \\ E_{b}(t)R(t)\vec{\xi}, & |\vec{\xi}| \geq R_{out}
62 : * \end{array}\right.
63 : * \f}
64 : *
65 : * Where $R_{in}$ is the inner radius, $R_{out}$ is the outer radius, and
66 : * $w_{T}$ is the translation falloff factor and $w_{E}$ is the expansion
67 : * falloff factor found through
68 : * \f{equation}{
69 : * w_{E} = \left\{\begin{array}{ll}\frac{R_{in}(R_{out} - |\vec{\xi}|)(E_{a}(t)
70 : * - E_{b}(t))}{|\vec{\xi}|(R_{out} - R_{in})}, & R_{in} < |\vec{\xi}| \leq
71 : * 0.5(R_{in} + R_{out}), \\ \frac{R_{out}(R_{in} - |\vec{\xi}|)(E_{a}(t)
72 : * - E_{b}(t))}{|\vec{\xi}|(R_{out} - R_{in})}, & 0.5(R_{in} + R_{out}) <
73 : * |\vec{\xi}| < R_{out} \end{array}\right.
74 : * \f}
75 : *
76 : * and
77 : *
78 : * \f{equation}{
79 : * w_{T} = \left\{\begin{array}{ll}\frac{R_{in} - |\vec{\xi}|}{R_{out} -
80 : * R_{in}}, & R_{in} < |\vec{\xi}| \leq 0.5(R_{in} + R_{out}), \\ \frac{R_{out}
81 : * - |\vec{\xi}|}{R_{out} - R_{in}}, & 0.5(R_{in} + R_{out}) < |\vec{\xi}| <
82 : * R_{out} \end{array}\right.
83 : * \f}
84 : *
85 : * $w_{E}$ and $w_{T}$ are calculated differently based on if you're closer
86 : * to the inner radius or outer radius to reduce roundoff error.
87 : *
88 : * ## Inverse
89 : * The inverse function maps the coordinates $\vec{\bar{\xi}}$ back to the
90 : * original coordinates $\vec{\xi}$ through different equations based on
91 : * which maps are supplied.
92 : *
93 : * If Rotation, Expansion, and Translation Maps are supplied then the
94 : * inverse is given by
95 : *
96 : * \f{equation}{
97 : * \label{eq:full_inverse}
98 : * \vec{\xi} = \left\{\begin{array}{ll}R^{T}(t)(\frac{(\vec{\bar{\xi}} -
99 : * \vec{T}(t))}{E_{a}(t)}), & |\vec{\bar{\xi}}| \leq R_{in}E_{a}(t),
100 : * \\ R^{T}(t)\frac{\vec{\bar{\xi}} - w_{T}\vec{T}(t)}{w_{E}},
101 : * & R_{in}E_{a}(t) < |\vec{\bar{\xi}}| \leq 0.5(R_{in}E_{a}(t) +
102 : * R_{out}E_{b}(t)), \\ R^{T}(t)\frac{\vec{\bar{\xi}} - (1.0 -
103 : * w_{T})\vec{T}(t)}{w_{E}}, & 0.5(R_{in}E_{a}(t) + R_{out}E_{b}(t)) <
104 : * |\vec{\bar{\xi}}| < R_{out}E_{b}(t),
105 : * \\ R^{T}(t)\frac{\vec{\bar{\xi}}}{E_{b}(t)}, & |\vec{\bar{\xi}}| \geq
106 : * R_{out}E_{b}(t) \end{array}\right.
107 : * \f}
108 : *
109 : * Where $w_{T}$ and $w_{E}$ are found through different quadratic solves.
110 : *
111 : * When closer to $R_{in}$ the quadratic has the form
112 : * \f{equation}{
113 : * w_{T}^2((E_{a}(t)R_{in} - E_{b}(t)R_{out})^2 - T(t)^2) +
114 : * 2w_{T}(E_{a}(t)R_{in}(E_{b}(t)R_{out} - E_{a}(t)R_{in}) + \vec{T}(t) \cdot
115 : * (\vec{T}(t) - \vec{\bar{\xi}})) + \vec{T}(t) \cdot \vec{\bar{\xi}})
116 : * + E_{a}(t)^2 R_{in}^2 - (\vec{\bar{\xi}} - \vec{T}(t))^2
117 : * \f}
118 : *
119 : * where $w_{E} = \frac{(1.0 - w_{T})R_{in}E_{a}(t) +
120 : * w_{T}R_{out}E_{b}(t)}{(1.0 - w_{T})R_{in} + w_{T}R_{out}}$
121 : *
122 : * When closer to $R_{out}$ the quadratic has the form
123 : * \f{equation}{
124 : * w_{T}^2((E_{a}(t)R_{in} - E_{b}(t)R_{out})^2 - T(t)^2) +
125 : * 2w_{T}(E_{b}(t)R_{out}(E_{a}(t)R_{in} - E_{b}(t)R_{out}) +
126 : * \vec{T}(t) \cdot \vec{\bar{\xi}})
127 : * + E_{b}(t)^2 R_{out}^2 - \vec{\bar{\xi}}^2
128 : * \f}
129 : *
130 : * where $w_{E} = \frac{w_{T}R_{in}E_{a}(t) + (1.0 -
131 : * w_{T})R_{out}E_{b}(t)}{w_{T}R_{in} + (1.0 - w_{T})R_{out}}$
132 : *
133 : * If Rotation and Expansion are supplied then the inverse is given by
134 : *
135 : * \f{equation}{
136 : * \vec{\xi} =
137 : * \left\{\begin{array}{ll}R^{T}(t)(\frac{\vec{\bar{\xi}}}{E_{a}(t)}), &
138 : * |\vec{\bar{\xi}}| \leq R_{in}E_{a}(t),
139 : * \\ R^{T}(t)\frac{\vec{\bar{\xi}}}{w_{E}}, & R_{in}E_{a}(t) <
140 : * |\vec{\bar{\xi}}| \leq 0.5(R_{in}E_{a}(t) + R_{out}E_{b}(t)),
141 : * \\ R^{T}(t)\frac{\vec{\bar{\xi}}}{w_{E}}, & 0.5(R_{in}E_{a}(t)
142 : * + R_{out}E_{b}(t)) < |\vec{\bar{\xi}}| < R_{out}E_{b}(t),
143 : * \\ R^{T}(t)\frac{\vec{\bar{\xi}}}{E_{b}(t)}, & |\vec{\bar{\xi}}| \geq
144 : * R_{out}E_{b}(t) \end{array}\right.
145 : * \f}
146 : *
147 : * Where $w_{E}$ is found through different quadratic solves.
148 : *
149 : * When closer to $R_{in}$ the quadratic has the form
150 : * \f{equation}{
151 : * w^2(E_{a}(t)R_{in} - E_{b}(t)R_{out})^2 +
152 : * 2wE_{a}(t)R_{in}(E_{b}(t)R_{out} - E_{a}(t)R_{in})
153 : * + (E_{a}(t) R_{in})^2 - \bar{\xi}^2
154 : * \f}
155 : * with $w_{E} = \frac{E_{a}(t)R_{in}(1.0 - w) + wE_{b}(t)R_{out}}{R_{in}(1.0 -
156 : * w) + wR_{out}}$
157 : *
158 : * When closer to $R_{out}$ the quadratic has the form
159 : * \f{equation}{
160 : * w^2(E_{a}(t)R_{in} - E_{b}(t)R_{out})^2 +
161 : * 2wE_{b}(t)R_{out}(E_{a}(t)R_{in} - E_{b}(t)R_{out})
162 : * + (E_{b}(t) R_{out})^2 - \bar{\xi}^2
163 : * \f}
164 : * with $w_{E} = \frac{wE_{a}(t)R_{in} + E_{b}(t)R_{out}(1.0 - w)}{wR_{in} +
165 : * R_{out}(1.0 - w)}$
166 : *
167 : * If Rotation and Translation are supplied, then the inverse is given by
168 : *
169 : * \f{equation}{
170 : * \vec{\xi} = \left\{\begin{array}{ll}R^{T}(t)(\vec{\bar{\xi}} -
171 : * \vec{T}(t)), & |\vec{\bar{\xi}}| \leq R_{in},
172 : * \\ R^{T}(t)(\vec{\bar{\xi}} - w_{T}\vec{T}(t)), & R_{in} < |\vec{\bar{\xi}}|
173 : * \leq 0.5(R_{in} + R_{out}), \\ R^{T}(t)(\vec{\bar{\xi}} - (1.0 -
174 : * w_{T})\vec{T}(t)), & 0.5(R_{in} + R_{out}) < |\vec{\bar{\xi}}| < R_{out},
175 : * \\ R^{T}(t)\vec{\bar{\xi}}, & |\vec{\bar{\xi}}| \geq R_{out}
176 : * \end{array}\right.
177 : * \f}
178 : *
179 : * Where $w_{T}$ is found through different quadratic solves.
180 : *
181 : * When closer to $R_{in}$ the quadratic has the form
182 : * \f{equation}{
183 : * w_{T}^2(T(t)^2 - (R_{out} - R_{in})^2) - 2w_{T}(\vec{T}(t) \cdot
184 : * (\vec{T}(t) - \vec{\bar{\xi}}) + R_{in}(R_{out} - R_{in})
185 : * + (\vec{T}(t) - \vec{\bar{\xi}})^2 - R_{in}^2
186 : * \f}
187 : *
188 : * When closer to $R_{out}$ the quadratic has the form
189 : * \f{equation}{
190 : * w_{T}^2(T(t)^2 - (R_{out} - R_{in})^2) +
191 : * 2w_{T}(R_{out}(R_{out} - R_{in}) - \vec{T}(t) \cdot \vec{\bar{\xi}})
192 : * + \vec{\bar{\xi}}^2 - R_{out}^2
193 : * \f}
194 : *
195 : * If Expansion and Translation are supplied, then the inverse is given by
196 : * Eq. $\ref{eq:full_inverse}$, with no transpose of rotation applied.
197 : *
198 : * \note For all the maps with rotation, the inverse of rotation is the
199 : * transpose of the original rotation. For maps with translation, the inverse
200 : * map also assumes that if $\vec{\bar{\xi}} - \vec{T}(t) \leq R_{in}$ then the
201 : * translated point originally came from within the inner radius so it'll be
202 : * translated back without a quadratic solve.
203 : *
204 : * ## Frame Velocity
205 : * The Frame Velocity is found through different equations based on which maps
206 : * are supplied.
207 : *
208 : * If Rotation, Expansion, and Translation are supplied then the frame
209 : * velocity is found through
210 : *
211 : * \f{equation}{
212 : * \vec{v} = \left\{\begin{array}{ll}(E_{a}(t)dR(t) +
213 : * dE_{a}(t)R(t))\vec{\xi} + d \vec{T}(t), & |\vec{\xi}| \leq R_{in},
214 : * \\ ((E_{a}(t) + w_{E})dR(t) + (dE_{a}(t) + dw_{E})R(t))\vec{\xi} + (1 +
215 : * w_{T})d \vec{T}(t), & R_{in} < |\vec{\xi}| \leq 0.5(R_{in} + R_{out}),
216 : * \\ ((E_{b}(t) + w_{E})dR(t) + (dE_{b}(t) + dw_{E})R(t))\vec{\xi} + w_{T}d
217 : * \vec{T}(t), & 0.5(R_{in} + R_{out}) < |\vec{\xi}| < R_{out},
218 : * \\ (E_{b}(t)dR(t) + dE_{b}(t)R(t))\vec{\xi}, & |\vec{\xi}| \geq R_{out}
219 : * \end{array}\right.
220 : * \f}
221 : *
222 : * where $dw_{E}$ is the derivative of the $w_{E}$ given by
223 : *
224 : * \f{equation}{
225 : * dw_{E} = \left\{\begin{array}{ll}\frac{R_{out}(R_{in} -
226 : * |\vec{\xi}|)(dE_{a}(t)
227 : * - dE_{b}(t))}{|\vec{\xi}|(R_{out} - R_{in})}, & R_{in} < |\vec{\xi}| \leq
228 : * 0.5(R_{in} + R_{out}), \\ \frac{R_{in}(R_{out} - |\vec{\xi}|)(dE_{a}(t) -
229 : * dE_{b}(t))}{|\vec{\xi}|(R_{out} - R_{in})}, & 0.5(R_{in} + R_{out}) <
230 : * |\vec{\xi}| < R_{out} \end{array}\right.
231 : * \f}
232 : *
233 : * ## Jacobian
234 : * The jacobian is also found through different equations based on which maps
235 : * are supplied.
236 : *
237 : * If Rotation, Expansion and Translation maps are supplied then the
238 : * jacobian is found through
239 : *
240 : * \f{equation}{
241 : * {J^{i}}_{j} = \left\{\begin{array}{ll}E_{a}(t){R^{i}}_{j}(t), & |\vec{\xi}|
242 : * \leq R_{in}, \\ {R^{i}}_{j}(t)E_{a}(t) + \frac{\alpha
243 : * {R^{i}}_{l}(t)\vec{\xi}^{l}\vec{\xi}_{j}(E_{A}(t) - E_{B}(t))}{|\vec{\xi}|} +
244 : * w_{E}{R^{i}}_{j}(t) + \frac{dw_{T}T^{i}\xi_{j}}{|\vec{\xi}|}, & R_{in} <
245 : * |\vec{\xi}| \leq 0.5(R_{in} + R_{out}), \\ {R^{i}}_{j}(t)E_{b}(t) +
246 : * \frac{\alpha {R^{i}}_{l}(t)\vec{\xi}^{l}\vec{\xi}_{j}(E_{A}(t) -
247 : * E_{B}(t))}{|\vec{\xi}|} + w_{E}{R^{i}}_{j}(t) +
248 : * \frac{dw_{T}T^{i}\xi_{j}}{|\vec{\xi}|}, & 0.5(R_{in} + R_{out}) < |\vec{\xi}|
249 : * < R_{out}, \\ E_{b}(t){R^{i}}_{j}(t), & |\vec{\xi}| \geq R_{out}
250 : * \end{array}\right.
251 : * \f}
252 : *
253 : * where $\alpha = \frac{R_{in}R_{out}}{\vec{\xi}^2(R_{in} - R_{out})}$ and
254 : * $dw_{T} = \frac{-1.0}{R_{out} - R_{in}}$
255 : *
256 : * \note For the translation map, the map returns the identity for all regions
257 : * except between $R_{in}$ and $R_{out}$
258 : *
259 : * ## Inverse Jacobian
260 : * The inverse jacobian is computed numerically by inverting the jacobian.
261 : */
262 : template <size_t Dim>
263 1 : class RotScaleTrans {
264 : public:
265 0 : enum class BlockRegion {
266 : /// Within the inner radius
267 : Inner,
268 : /// Between inner and outer radius
269 : Transition,
270 : /// At or beyond outer boundary
271 : Outer
272 : };
273 :
274 0 : static constexpr size_t dim = Dim;
275 :
276 0 : explicit RotScaleTrans(
277 : std::optional<std::pair<std::string, std::string>> scale_f_of_t_names,
278 : std::optional<std::string> rot_f_of_t_name,
279 : std::optional<std::string> trans_f_of_t_name, double inner_radius,
280 : double outer_radius, BlockRegion region);
281 :
282 0 : RotScaleTrans() = default;
283 0 : ~RotScaleTrans() = default;
284 0 : RotScaleTrans(const RotScaleTrans<Dim>& RotScaleTrans_Map) = default;
285 0 : RotScaleTrans(RotScaleTrans&&) = default;
286 0 : RotScaleTrans& operator=(RotScaleTrans&&) = default;
287 0 : RotScaleTrans& operator=(const RotScaleTrans& RotScaleTrans_Map) = default;
288 :
289 : template <typename T>
290 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
291 : const std::array<T, Dim>& source_coords, double time,
292 : const std::unordered_map<
293 : std::string,
294 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
295 : functions_of_time) const;
296 :
297 : /// The inverse function is only callable with doubles because the inverse
298 : /// might fail if called for a point out of range, and it is unclear
299 : /// what should happen if the inverse were to succeed for some points in a
300 : /// DataVector but fail for other points.
301 1 : std::optional<std::array<double, Dim>> inverse(
302 : const std::array<double, Dim>& target_coords, double time,
303 : const std::unordered_map<
304 : std::string,
305 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
306 : functions_of_time) const;
307 :
308 : template <typename T>
309 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> frame_velocity(
310 : const std::array<T, Dim>& source_coords, double time,
311 : const std::unordered_map<
312 : std::string,
313 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
314 : functions_of_time) const;
315 :
316 : template <typename T>
317 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
318 : const std::array<T, Dim>& source_coords, double time,
319 : const std::unordered_map<
320 : std::string,
321 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
322 : functions_of_time) const;
323 :
324 : template <typename T>
325 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
326 : const std::array<T, Dim>& source_coords, double time,
327 : const std::unordered_map<
328 : std::string,
329 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
330 : functions_of_time) const;
331 :
332 : // NOLINTNEXTLINE(google-runtime-references)
333 0 : void pup(PUP::er& p);
334 :
335 0 : static bool is_identity() { return false; }
336 :
337 0 : const std::unordered_set<std::string>& function_of_time_names() const {
338 : return f_of_t_names_;
339 : }
340 :
341 : private:
342 : template <size_t LocalDim>
343 0 : friend bool operator==( // NOLINT(readability-redundant-declaration)
344 : const RotScaleTrans<LocalDim>& lhs, const RotScaleTrans<LocalDim>& rhs);
345 :
346 : // The root helper returns the correct root of the roots found during the
347 : // quadratic solve in the inverse function.
348 0 : double root_helper(std::optional<std::array<double, 2>> roots) const;
349 :
350 0 : std::optional<std::string> scale_f_of_t_a_{};
351 0 : std::optional<std::string> scale_f_of_t_b_{};
352 0 : std::optional<std::string> rot_f_of_t_{};
353 0 : std::optional<std::string> trans_f_of_t_{};
354 0 : std::unordered_set<std::string> f_of_t_names_;
355 0 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
356 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
357 0 : BlockRegion region_ = BlockRegion::Inner;
358 : };
359 :
360 : template <size_t Dim>
361 0 : bool operator!=(const RotScaleTrans<Dim>& lhs, const RotScaleTrans<Dim>& rhs) {
362 : return not(lhs == rhs);
363 : }
364 :
365 : } // namespace domain::CoordinateMaps::TimeDependent
|
