SphericalCompression.hpp
1 // 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 
15 #include "Utilities/TypeTraits/RemoveReferenceWrapper.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 
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 Currently, this map only performs a compression. A generalization of
56  * this map could also change the region's shape as well as its size, by
57  * including more terms than the spherically symmetric term included here.
58  *
59  * ### Mapped coordinates
60  *
61  * The mapped coordinates
62  * \f$x^i\f$ are related to the unmapped coordinates \f$\xi^i\f$
63  * as follows:
64  * \f{align}{
65  * x^i &= \left\{\begin{array}{ll}\xi^i - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
66  * \frac{\rho^i}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
67  * \xi^i - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
68  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}} \rho^i, &
69  * \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
70  * \xi^i, & \rho_{\rm max} < \rho,\end{array}\right.
71  * \f}
72  * where \f$\rho^i = \xi^i - C^i\f$ is the Euclidean radial position vector in
73  * the unmapped coordinates with respect to the center \f$C^i\f$, \f$\rho =
74  * \sqrt{\delta_{kl}\left(\xi^k - C^l\right)\left(\xi^l - C^l\right)}\f$ is the
75  * Euclidean magnitude of \f$\rho^i\f$, and \f$\rho_j = \delta_{ij} \rho^i\f$.
76  *
77  * ### Frame velocity
78  *
79  * The frame velocity \f$v^i \equiv dx^i/dt\f$ is then
80  * \f{align}{
81  * v^i &= \left\{\begin{array}{ll} - \frac{\lambda_{00}^\prime(t)}{\sqrt{4\pi}}
82  * \frac{\rho^i}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
83  * - \frac{\lambda_{00}^\prime(t)}{\sqrt{4\pi}}
84  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}} \rho^i,
85  * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
86  * 0, & \rho_{\rm max} < \rho,\end{array}\right.
87  * \f} where \f$\lambda_{00}^\prime(t) \equiv d\lambda_{00}/dt\f$.
88  *
89  * ### Jacobian
90  *
91  * Differentiating the equations for \f$x^i\f$ gives the Jacobian
92  * \f$\partial x^i / \partial \xi^j\f$. Using the result
93  * \f{align}{
94  * \frac{\partial \rho^i}{\partial \xi^j} &= \frac{\partial}{\partial \xi^j}
95  * \left(\xi^i - C^i\right) = \frac{\partial \xi^i}{\partial \xi^j}
96  * = \delta^i_{j}
97  * \f}
98  * and taking the derivatives yields
99  * \f{align}{
100  * \frac{\partial x^i}{\partial \xi^j} &= \left\{\begin{array}{ll}
101  * \delta^i_j \left(1
102  * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm min}}\right),
103  * & \rho < \rho_{\rm min},\\
104  * \delta^i_j
105  * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
106  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
107  * - \rho^i \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
108  * \frac{\partial}{\partial \xi^j}\left(
109  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right),
110  * & \rho_{\rm min} \leq \rho < \rho_{\rm max},\\
111  * \delta^i_j, & \rho_{\rm max} < \rho.\end{array}\right.
112  * \f}
113  * Inserting
114  * \f{align}{
115  * \frac{\partial}{\partial \xi^j}\left(
116  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
117  * &= \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}
118  * \frac{\partial}{\partial \xi^j}\left(\frac{1}{\rho}\right)
119  * = - \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}} \frac{1}{\rho^2}
120  * \frac{\partial \rho}{\partial \xi^j}
121  * \f}
122  * and
123  * \f{align}{
124  * \frac{\partial \rho}{\partial \xi^j} &= \frac{\rho_j}{\rho}.
125  * \f}
126  * into the Jacobian yields
127  * \f{align}{
128  * \frac{\partial x^i}{\partial \xi^j} &= \left\{\begin{array}{ll}
129  * \delta^i_j \left(1
130  * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm min}}\right),
131  * & \rho < \rho_{\rm min},\\
132  * \delta^i_j
133  * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
134  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right)
135  * + \rho^i \rho_j \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
136  * \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}\frac{1}{\rho^3},
137  * & \rho_{\rm min} \leq \rho < \rho_{\rm max},\\
138  * \delta^i_j, & \rho_{\rm max} < \rho.\end{array}\right.
139  * \f}
140  *
141  * ### Inverse Jacobian
142  *
143  * This map finds the inverse Jacobian by first finding the Jacobian and then
144  * numerically inverting it.
145  *
146  * ### Inverse map
147  *
148  * For \f$\lambda_{00}(t)\f$ that satisfy
149  * \f{align}{
150  * \rho_{\rm min} - \rho_{\rm max} < \lambda_{00}(t) / \sqrt{4\pi} <
151  * \rho_{\rm min},
152  * \f}
153  * the map will be invertible and nonsingular. For simplicity, here we
154  * enforce this condition, even though perhaps the map might be generalized to
155  * handle cases that are still invertible but violate this condition. This
156  * avoids the need to specially handle the cases
157  * \f$\lambda_{00}(t) / \sqrt{4\pi} = \rho_{\rm min} - \rho_{\rm max}\f$
158  * and \f$\lambda_{00}(t) / \sqrt{4\pi} = \rho_{\rm min}\f$, both of which
159  * yield a singular map, and it also avoids cases where the map behaves
160  * in undesirable ways (such as a larger \f$\lambda_{00}(t)\f$ leading to
161  * an expansion and a coordinate inversion instead of a compression).
162  *
163  * After
164  * requiring the above inequality to be satisfied, however, the inverse mapping
165  * can be derived as follows. Let \f$r^i \equiv x^i - C^i\f$. In terms of
166  * \f$r^i\f$, the map is \f{align}{ r^i &= \left\{\begin{array}{ll}\rho^i
167  * \left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
168  * \frac{1}{\rho_{\rm min}}\right), & \rho < \rho_{\rm min}, \\
169  * \rho^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
170  * \frac{\rho_{\rm max} / \rho - 1}{\rho_{\rm max} - \rho_{\rm min}}\right),
171  * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
172  * \rho^i, & \rho_{\rm max} < \rho.\end{array}\right.
173  * \f}
174  *
175  * Taking the Euclidean magnitude of both sides and simplifying yields
176  * \f{align}{
177  * \frac{r}{\rho} &= \left\{\begin{array}{ll}
178  * 1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
179  * \frac{1}{\rho_{\rm min}}, & \rho < \rho_{\rm min}, \\
180  * 1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
181  * \frac{\rho_{\rm max}/\rho - 1}{\rho_{\rm max} - \rho_{\rm min}},
182  * & \rho_{\rm min} \leq \rho \leq \rho_{\rm max}, \\
183  * 1, & \rho_{\rm max} < \rho,\end{array}\right.
184  * \f}
185  * which implies
186  * \f{align}{
187  * r^i = \rho^i \frac{r}{\rho} \Rightarrow \rho^i = r^i \frac{\rho}{r}.
188  * \f}
189  *
190  * Inserting \f$\rho_{\rm min}\f$ or \f$\rho_{\rm max}\f$ then gives the
191  * corresponding bounds in the mapped coordinates: \f{align}{
192  * r_{\rm min} &= \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
193  * r_{\rm max} &= \rho_{\rm max}.
194  * \f}
195  *
196  * In the regime \f$\rho_{\rm min} \leq \rho < \rho_{\rm max}\f$, rearranging
197  * yields a linear relationship between \f$\rho\f$ and \f$r\f$, which
198  * can then be solved for \f$\rho(r)\f$:
199  * \f{align}{
200  * r &= \rho - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
201  * \frac{\rho_{\rm max} - \rho}{\rho_{\rm max} - \rho_{\rm min}}\\
202  * \Rightarrow r &= \rho \left(1 + \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
203  * \frac{1}{\rho_{\rm max} - \rho_{\rm min}}\right)
204  * - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
205  * \frac{\rho_{\rm max}}{\rho_{\rm max} - \rho_{\rm min}}.
206  * \f}
207  * Solving this linear equation for \f$\rho\f$ yields
208  * \f{align}{
209  * \rho &= \left(r+\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
210  * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)
211  * \left(1 + \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
212  * \frac{1}{\rho_{\rm max} - \rho_{\rm min}}\right)^{-1}.
213  * \f}
214  *
215  * Inserting the expressions for \f$\rho\f$ into the equation
216  * \f{align}{
217  * \rho^i = r^i \frac{\rho}{r}
218  * \f}
219  * then gives
220  * \f{align}{
221  * \rho^i &= \left\{\begin{array}{ll}
222  * r^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
223  * \frac{1}{\rho_{\rm min}}\right)^{-1},
224  * & r < \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
225  * r^i
226  * \left(1+\frac{1}{r}\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
227  * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)\left(1 +
228  * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm max} - \rho_{\rm
229  * min}}\right)^{-1}, & \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
230  * \leq r
231  * \leq \rho_{\rm max},\\
232  * r^i, & \rho_{\rm max} < r.\end{array}\right.
233  * \f}
234  * Finally, inserting \f$\rho^i = \xi^i - C^i\f$ yields the inverse map:
235  * \f{align}{
236  * \xi^i &= \left\{\begin{array}{ll}
237  * r^i\left(1 - \frac{\lambda_{00}(t)}{\sqrt{4\pi}}
238  * \frac{1}{\rho_{\rm min}}\right)^{-1} + C^i,
239  * & r < \rho_{\rm min} - \frac{\lambda_{00}(t)}{\sqrt{4\pi}},\\
240  * r^i
241  * \left(1+\frac{1}{r}\frac{\lambda_{00}(t)}{\sqrt{4\pi}}\frac{\rho_{\rm
242  * max}}{\rho_{\rm max}-\rho_{\rm min}}\right)\left(1 +
243  * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \frac{1}{\rho_{\rm max} - \rho_{\rm
244  * min}}\right)^{-1} + C^i, & \rho_{\rm min} -
245  * \frac{\lambda_{00}(t)}{\sqrt{4\pi}} \leq r
246  * \leq \rho_{\rm max},\\
247  * r^i + C^i = x^i, & \rho_{\rm max} < r.\end{array}\right.
248  * \f}
249  *
250  */
251 template <bool InteriorMap>
253  public:
254  static constexpr size_t dim = 3;
255 
256  explicit SphericalCompression(std::string function_of_time_name,
257  double min_radius, double max_radius,
258  const std::array<double, 3>& center) noexcept;
259  SphericalCompression() = default;
260 
261  template <typename T>
263  const std::array<T, 3>& source_coords, double time,
264  const std::unordered_map<
265  std::string,
267  functions_of_time) const noexcept;
268 
269  /// The inverse function is only callable with doubles because the inverse
270  /// might fail if called for a point out of range, and it is unclear
271  /// what should happen if the inverse were to succeed for some points in a
272  /// DataVector but fail for other points.
274  const std::array<double, 3>& target_coords, double time,
275  const std::unordered_map<
276  std::string,
278  functions_of_time) const noexcept;
279 
280  template <typename T>
281  std::array<tt::remove_cvref_wrap_t<T>, 3> frame_velocity(
282  const std::array<T, 3>& source_coords, double time,
283  const std::unordered_map<
284  std::string,
286  functions_of_time) const noexcept;
287 
288  template <typename T>
289  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
290  const std::array<T, 3>& source_coords, double time,
291  const std::unordered_map<
292  std::string,
294  functions_of_time) const noexcept;
295 
296  template <typename T>
297  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
298  const std::array<T, 3>& source_coords, double time,
299  const std::unordered_map<
300  std::string,
302  functions_of_time) const noexcept;
303 
304  void pup(PUP::er& p) noexcept; // NOLINT
305 
306  static bool is_identity() noexcept { return false; }
307 
308  private:
309  friend bool operator==(const SphericalCompression& lhs,
310  const SphericalCompression& rhs) noexcept {
311  return lhs.f_of_t_name_ == rhs.f_of_t_name_ and
312  lhs.min_radius_ == rhs.min_radius_ and
313  lhs.max_radius_ == rhs.max_radius_ and lhs.center_ == rhs.center_;
314  }
315  std::string f_of_t_name_;
316  double min_radius_ = std::numeric_limits<double>::signaling_NaN();
317  double max_radius_ = std::numeric_limits<double>::signaling_NaN();
318  std::array<double, 3> center_;
319 };
320 
321 template <bool InteriorMap>
322 bool operator!=(const SphericalCompression<InteriorMap>& lhs,
323  const SphericalCompression<InteriorMap>& rhs) noexcept {
324  return not(lhs == rhs);
325 }
326 } // namespace domain::CoordinateMaps::TimeDependent
std::string
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
domain::FunctionsOfTime
Contains functions of time to support the dual frame system.
Definition: FixedSpeedCubic.hpp:17
domain::CoordinateMaps::TimeDependent::SphericalCompression::inverse
std::optional< std::array< double, 3 > > inverse(const std::array< double, 3 > &target_coords, double time, const std::unordered_map< std::string, std::unique_ptr< domain::FunctionsOfTime::FunctionOfTime >> &functions_of_time) const noexcept
The inverse function is only callable with doubles because the inverse might fail if called for a poi...
cstddef
array
domain::CoordinateMaps::TimeDependent::SphericalCompression
Time-dependent compression of a finite 3D spherical volume.
Definition: SphericalCompression.hpp:252
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
memory
limits
TypeAliases.hpp
optional
domain::CoordinateMaps::TimeDependent
Contains the time-dependent coordinate maps.
Definition: CoordinateMap.hpp:35
std::unique_ptr
unordered_map
string