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