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