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 <limits>
9 : #include <memory>
10 : #include <optional>
11 : #include <string>
12 : #include <unordered_map>
13 : #include <unordered_set>
14 :
15 : #include "DataStructures/Tensor/TypeAliases.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 Time-dependent compression of a finite 3D spherical volume.
31 : *
32 : * \details Let \f$\xi^i\f$ be the unmapped coordinates, and let \f$\rho\f$ be
33 : * the Euclidean radius corresponding to these coordinates with respect to
34 : * some center \f$C^i\f$. The transformation implemented by this map is
35 : * equivalent to the following transformation: at each point, the mapped
36 : * coordinates are the same as the unmapped coordinates, except in
37 : * a spherical region \f$\rho \leq \rho_{\rm max}\f$, where instead coordinates
38 : * are mapped using a compression that is spherically symmetric about the center
39 : * \f$C^i\f$. The amount of compression decreases linearly from a maximum at
40 : * \f$\rho = \rho_{\rm min}\f$ to zero at \f$\rho = \rho_{\rm max}\f$. A
41 : * scalar domain::FunctionsOfTime::FunctionOfTime \f$\lambda_{00}(t)\f$ controls
42 : * the amount of compression.
43 : *
44 : * The mapped coordinates are a continuous function of the unmapped
45 : * coordinates, but the Jacobians are not continuous at \f$\rho_{\rm min}\f$
46 : * and \f$\rho_{\rm max}\f$. Therefore, \f$\rho_{\rm min}\f$ and \f$\rho_{\rm
47 : * max}\f$ should both be surfaces corresponding to block boundaries. Therefore,
48 : * this class implements the transformation described above as follows: the
49 : * if the template parameter `InteriorMap` is true, the map is the one
50 : * appropriate for \f$\rho < \rho_{\rm min}\f$, while if `InteriorMap` is false,
51 : * the map is the one appropriate for \f$\rho_{\rm min} \leq \rho \leq \rho_{\rm
52 : * max}\f$. To use this map, add it to the blocks where the transformation is
53 : * not the identity, using the appropriate template parameter, depending on
54 : * which region the block is in.
55 : *
56 : * \note This map performs a only a spherical compression. A
57 : * generalization of this map that changes the region's shape as well as
58 : * its size, by including more terms than the spherically symmetric
59 : * term included here, can be found in the
60 : * domain::CoordinateMaps::TimeDependent::Shape map.
61 : *
62 : * \note The quantity stored in the FunctionOfTime is really
63 : * the spherical-harmonic coefficient \f$\lambda_{00}(t)\f$. This is
64 : * different from the Shape map, which stores ylm::Spherepack coefficients
65 : * \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$ instead of \f$\lambda_{lm}(t)\f$.
66 : * See domain::CoordinateMaps::TimeDependent::Shape for more details.
67 : *
68 : * ### Mapped coordinates
69 : *
70 : * The mapped coordinates
71 : * \f$x^i\f$ are related to the unmapped coordinates \f$\xi^i\f$
72 : * as follows:
73 : * \f{align}{
74 : * x^i &= \left\{\begin{array}{ll}\xi^i - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
75 : * \frac{\rho^i}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
76 : * \xi^i - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
77 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}} \rho^i, &
78 : * \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
79 : * \xi^i, & \rho_{\rm max} < \rho,\end{array}\right.
80 : * \f}
81 : * where \f$\rho^i = \xi^i - C^i\f$ is the Euclidean radial position vector in
82 : * the unmapped coordinates with respect to the center \f$C^i\f$, \f$\rho =
83 : * \sqrt{\delta_{kl}\left(\xi^k - C^l\right)\left(\xi^l - C^l\right)}\f$ is the
84 : * Euclidean magnitude of \f$\rho^i\f$, and \f$\rho_j = \delta_{ij} \rho^i\f$.
85 : *
86 : * ### Frame velocity
87 : *
88 : * The frame velocity \f$v^i \equiv dx^i/dt\f$ is then
89 : * \f{align}{
90 : * v^i &= \left\{\begin{array}{ll} - \frac{\lambda_{00}^\prime(t)}{\sqrt{4\pi}}
91 : * \frac{\rho^i}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
92 : * - \frac{\lambda_{00}^\prime(t)}{\sqrt{4\pi}}
93 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}} \rho^i,
94 : * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
95 : * 0, & \rho_{\rm max} < \rho,\end{array}\right.
96 : * \f} where \f$\lambda_{00}^\prime(t) \equiv d\lambda_{00}/dt\f$.
97 : *
98 : * ### Jacobian
99 : *
100 : * Differentiating the equations for \f$x^i\f$ gives the Jacobian
101 : * \f$\partial x^i / \partial \xi^j\f$. Using the result
102 : * \f{align}{
103 : * \frac{\partial \rho^i}{\partial \xi^j} &= \frac{\partial}{\partial \xi^j}
104 : * \left(\xi^i - C^i\right) = \frac{\partial \xi^i}{\partial \xi^j}
105 : * = \delta^i_{j}
106 : * \f}
107 : * and taking the derivatives yields
108 : * \f{align}{
109 : * \frac{\partial x^i}{\partial \xi^j} &= \left\{\begin{array}{ll}
110 : * \delta^i_j \left(1
111 : * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm min}}\right),
112 : * & \rho < \rho_{\rm min},\\
113 : * \delta^i_j
114 : * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
115 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
116 : * - \rho^i \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
117 : * \frac{\partial}{\partial \xi^j}\left(
118 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right),
119 : * & \rho_{\rm min} \leq \rho < \rho_{\rm max},\\
120 : * \delta^i_j, & \rho_{\rm max} < \rho.\end{array}\right.
121 : * \f}
122 : * Inserting
123 : * \f{align}{
124 : * \frac{\partial}{\partial \xi^j}\left(
125 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
126 : * &= \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}
127 : * \frac{\partial}{\partial \xi^j}\left(\frac{1}{\rho}\right)
128 : * = - \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}} \frac{1}{\rho^2}
129 : * \frac{\partial \rho}{\partial \xi^j}
130 : * \f}
131 : * and
132 : * \f{align}{
133 : * \frac{\partial \rho}{\partial \xi^j} &= \frac{\rho_j}{\rho}.
134 : * \f}
135 : * into the Jacobian yields
136 : * \f{align}{
137 : * \frac{\partial x^i}{\partial \xi^j} &= \left\{\begin{array}{ll}
138 : * \delta^i_j \left(1
139 : * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm min}}\right),
140 : * & \rho < \rho_{\rm min},\\
141 : * \delta^i_j
142 : * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
143 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
144 : * + \rho^i \rho_j \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
145 : * \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}\frac{1}{\rho^3},
146 : * & \rho_{\rm min} \leq \rho < \rho_{\rm max},\\
147 : * \delta^i_j, & \rho_{\rm max} < \rho.\end{array}\right.
148 : * \f}
149 : *
150 : * ### Inverse Jacobian
151 : *
152 : * This map finds the inverse Jacobian by first finding the Jacobian and then
153 : * numerically inverting it.
154 : *
155 : * ### Inverse map
156 : *
157 : * For \f$\lambda_{00}(t)\f$ that satisfy
158 : * \f{align}{
159 : * \rho_{\rm min} - \rho_{\rm max} < \lambda_{00}(t) / \sqrt{4\pi} <
160 : * \rho_{\rm min},
161 : * \f}
162 : * the map will be invertible and nonsingular. For simplicity, here we
163 : * enforce this condition, even though perhaps the map might be generalized to
164 : * handle cases that are still invertible but violate this condition. This
165 : * avoids the need to specially handle the cases
166 : * \f$\lambda_{00}(t) / \sqrt{4\pi} = \rho_{\rm min} - \rho_{\rm max}\f$
167 : * and \f$\lambda_{00}(t) / \sqrt{4\pi} = \rho_{\rm min}\f$, both of which
168 : * yield a singular map, and it also avoids cases where the map behaves
169 : * in undesirable ways (such as a larger \f$\lambda_{00}(t)\f$ leading to
170 : * an expansion and a coordinate inversion instead of a compression).
171 : *
172 : * After
173 : * requiring the above inequality to be satisfied, however, the inverse mapping
174 : * can be derived as follows. Let \f$r^i \equiv x^i - C^i\f$. In terms of
175 : * \f$r^i\f$, the map is \f{align}{ r^i &= \left\{\begin{array}{ll}\rho^i
176 : * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
177 : * \frac{1}{\rho_{\rm min}}\right), & \rho < \rho_{\rm min}, \\
178 : * \rho^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
179 : * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right),
180 : * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
181 : * \rho^i, & \rho_{\rm max} < \rho.\end{array}\right.
182 : * \f}
183 : *
184 : * Taking the Euclidean magnitude of both sides and simplifying yields
185 : * \f{align}{
186 : * \frac{r}{\rho} &= \left\{\begin{array}{ll}
187 : * 1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
188 : * \frac{1}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
189 : * 1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
190 : * \frac{\rho_{\rm max}/\rho - 1}{\rho_{\rm max} - \rho_{\rm min}},
191 : * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
192 : * 1, & \rho_{\rm max} < \rho,\end{array}\right.
193 : * \f}
194 : * which implies
195 : * \f{align}{
196 : * r^i = \rho^i \frac{r}{\rho} \Rightarrow \rho^i = r^i \frac{\rho}{r}.
197 : * \f}
198 : *
199 : * Inserting \f$\rho_{\rm min}\f$ or \f$\rho_{\rm max}\f$ then gives the
200 : * corresponding bounds in the mapped coordinates: \f{align}{
201 : * r_{\rm min} &= \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
202 : * r_{\rm max} &= \rho_{\rm max}.
203 : * \f}
204 : *
205 : * In the regime \f$\rho_{\rm min} \leq \rho < \rho_{\rm max}\f$, rearranging
206 : * yields a linear relationship between \f$\rho\f$ and \f$r\f$, which
207 : * can then be solved for \f$\rho(r)\f$:
208 : * \f{align}{
209 : * r &= \rho - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
210 : * \frac{\rho_{\rm max} - \rho}{\rho_{\rm max} - \rho_{\rm min}}\\
211 : * \Rightarrow r &= \rho \left(1 + \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
212 : * \frac{1}{\rho_{\rm max} - \rho_{\rm min}}\right)
213 : * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
214 : * \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}.
215 : * \f}
216 : * Solving this linear equation for \f$\rho\f$ yields
217 : * \f{align}{
218 : * \rho &= \left(r+\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
219 : * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)
220 : * \left(1 + \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
221 : * \frac{1}{\rho_{\rm max} - \rho_{\rm min}}\right)^{-1}.
222 : * \f}
223 : *
224 : * Inserting the expressions for \f$\rho\f$ into the equation
225 : * \f{align}{
226 : * \rho^i = r^i \frac{\rho}{r}
227 : * \f}
228 : * then gives
229 : * \f{align}{
230 : * \rho^i &= \left\{\begin{array}{ll}
231 : * r^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
232 : * \frac{1}{\rho_{\rm min}}\right)^{-1},
233 : * & r < \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
234 : * r^i
235 : * \left(1+\frac{1}{r}\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
236 : * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)\left(1 +
237 : * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm max} - \rho_{\rm
238 : * min}}\right)^{-1}, & \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
239 : * \leq r
240 : * \leq \rho_{\rm max},\\
241 : * r^i, & \rho_{\rm max} < r.\end{array}\right.
242 : * \f}
243 : * Finally, inserting \f$\rho^i = \xi^i - C^i\f$ yields the inverse map:
244 : * \f{align}{
245 : * \xi^i &= \left\{\begin{array}{ll}
246 : * r^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
247 : * \frac{1}{\rho_{\rm min}}\right)^{-1} + C^i,
248 : * & r < \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
249 : * r^i
250 : * \left(1+\frac{1}{r}\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
251 : * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)\left(1 +
252 : * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm max} - \rho_{\rm
253 : * min}}\right)^{-1} + C^i, & \rho_{\rm min} -
254 : * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \leq r
255 : * \leq \rho_{\rm max},\\
256 : * r^i + C^i = x^i, & \rho_{\rm max} < r.\end{array}\right.
257 : * \f}
258 : *
259 : */
260 : template <bool InteriorMap>
261 1 : class SphericalCompression {
262 : public:
263 0 : static constexpr size_t dim = 3;
264 :
265 0 : explicit SphericalCompression(std::string function_of_time_name,
266 : double min_radius, double max_radius,
267 : const std::array<double, 3>& center);
268 0 : SphericalCompression() = default;
269 :
270 : template <typename T>
271 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
272 : const std::array<T, 3>& source_coords, double time,
273 : const std::unordered_map<
274 : std::string,
275 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
276 : functions_of_time) const;
277 :
278 : /// The inverse function is only callable with doubles because the inverse
279 : /// might fail if called for a point out of range, and it is unclear
280 : /// what should happen if the inverse were to succeed for some points in a
281 : /// DataVector but fail for other points.
282 1 : std::optional<std::array<double, 3>> inverse(
283 : const std::array<double, 3>& target_coords, double time,
284 : const std::unordered_map<
285 : std::string,
286 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
287 : functions_of_time) const;
288 :
289 : template <typename T>
290 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> frame_velocity(
291 : const std::array<T, 3>& 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 : template <typename T>
298 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
299 : const std::array<T, 3>& source_coords, double time,
300 : const std::unordered_map<
301 : std::string,
302 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
303 : functions_of_time) const;
304 :
305 : template <typename T>
306 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
307 : const std::array<T, 3>& source_coords, double time,
308 : const std::unordered_map<
309 : std::string,
310 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
311 : functions_of_time) const;
312 :
313 : // NOLINTNEXTLINE(google-runtime-references)
314 0 : void pup(PUP::er& p);
315 :
316 0 : static bool is_identity() { return false; }
317 :
318 0 : const std::unordered_set<std::string>& function_of_time_names() const {
319 : return f_of_t_names_;
320 : }
321 :
322 : private:
323 0 : friend bool operator==(const SphericalCompression& lhs,
324 : const SphericalCompression& rhs) {
325 : return lhs.f_of_t_name_ == rhs.f_of_t_name_ and
326 : lhs.min_radius_ == rhs.min_radius_ and
327 : lhs.max_radius_ == rhs.max_radius_ and lhs.center_ == rhs.center_;
328 : }
329 0 : std::string f_of_t_name_;
330 0 : std::unordered_set<std::string> f_of_t_names_;
331 0 : double min_radius_ = std::numeric_limits<double>::signaling_NaN();
332 0 : double max_radius_ = std::numeric_limits<double>::signaling_NaN();
333 0 : std::array<double, 3> center_;
334 : };
335 :
336 : template <bool InteriorMap>
337 0 : bool operator!=(const SphericalCompression<InteriorMap>& lhs,
338 : const SphericalCompression<InteriorMap>& rhs) {
339 : return not(lhs == rhs);
340 : }
341 : } // namespace domain::CoordinateMaps::TimeDependent
|
