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 <optional>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Utilities/Gsl.hpp"
13 :
14 : /// \cond
15 : namespace PUP {
16 : class er;
17 : } // namespace PUP
18 : /// \endcond
19 :
20 : /// Contains FocallyLiftedInnerMaps
21 1 : namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {
22 : /*!
23 : * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylinder
24 : * to a volume that connects portions of two spherical surfaces.
25 : *
26 : * Because FocallyLiftedEndcap is a FocallyLiftedInnerMap, it is meant
27 : * to be a template parameter of FocallyLiftedMap, and its member functions
28 : * are meant to be used by FocallyLiftedMap. See FocallyLiftedMap for further
29 : * documentation.
30 : *
31 : * \image html FocallyLiftedEndcap.svg "Focally Lifted Endcap."
32 : *
33 : * \details The domain of the map is a 3D unit right cylinder with
34 : * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
35 : * \f$-1\leq\bar{z}\leq 1\f$ and \f$\bar{x}^2+\bar{y}^2 \leq
36 : * 1\f$. The range of the map has coordinates \f$(x,y,z)\f$.
37 : *
38 : * Consider a sphere with center \f$C^i\f$ and radius \f$R\f$ that is
39 : * intersected by a plane normal to the \f$z\f$ axis and located at
40 : * \f$z = z_\mathrm{P}\f$. In the figure above, every point
41 : * \f$\bar{x}^i\f$ in the blue region \f$\sigma=0\f$ maps to a point
42 : * \f$x_0^i\f$ on a portion of the surface of the sphere.
43 : *
44 : * `Endcap` provides the following functions:
45 : *
46 : * ### forward_map
47 : * `forward_map` maps \f$(\bar{x},\bar{y},\bar{z}=-1)\f$ to the portion of
48 : * the sphere with \f$z \geq z_\mathrm{P}\f$. The arguments to `forward_map`
49 : * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
50 : * `forward_map` returns \f$x_0^i\f$,
51 : * the 3D coordinates on that sphere, which are given by
52 : *
53 : * \f{align}
54 : * x_0^0 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max})
55 : * \bar{x}}{\bar{\rho}} + C^0,\\
56 : * x_0^1 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max})
57 : * \bar{y}}{\bar{\rho}} + C^1,\\
58 : * x_0^2 &= R \cos(\bar{\rho} \theta_\mathrm{max}) + C^2.
59 : * \f}
60 : *
61 : * Here \f$\bar{\rho}^2 \equiv (\bar{x}^2+\bar{y}^2)/\bar{R}^2\f$, where
62 : * \f$\bar{R}\f$ is the radius of the cylinder in barred coordinates,
63 : * which is always unity,
64 : * and where
65 : * \f$\theta_\mathrm{max}\f$ is defined by
66 : * \f$\cos(\theta_\mathrm{max}) = (z_\mathrm{P}-C^2)/R\f$.
67 : * Note that when \f$\bar{\rho}=0\f$, we must evaluate
68 : * \f$\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\f$
69 : * as \f$\theta_\mathrm{max}\f$.
70 : *
71 : * ### sigma
72 : *
73 : * \f$\sigma\f$ is a function that is zero at \f$\bar{z}=-1\f$
74 : * (which maps onto the sphere \f$x^i=x_0^i\f$) and
75 : * unity at \f$\bar{z}=+1\f$ (corresponding to the
76 : * upper surface of the FocallyLiftedMap). We define
77 : *
78 : * \f{align}
79 : * \sigma &= \frac{\bar{z}+1}{2}.
80 : * \f}
81 : *
82 : * ### deriv_sigma
83 : *
84 : * `deriv_sigma` returns
85 : *
86 : * \f{align}
87 : * \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2).
88 : * \f}
89 : *
90 : * ### jacobian
91 : *
92 : * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
93 : * The arguments to `jacobian`
94 : * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
95 : *
96 : * Differentiating Eqs.(1--3) above yields
97 : *
98 : * \f{align*}
99 : * \frac{\partial x_0^2}{\partial \bar{x}} &=
100 : * - R \theta_\mathrm{max}
101 : * \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{x}\\
102 : * \frac{\partial x_0^2}{\partial \bar{y}} &=
103 : * - R \theta_\mathrm{max}
104 : * \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{y}\\
105 : * \frac{\partial x_0^0}{\partial \bar{x}} &=
106 : * R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} +
107 : * R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}}
108 : * \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right)
109 : * \bar{x}^2\\
110 : * \frac{\partial x_0^0}{\partial \bar{y}} &=
111 : * R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}}
112 : * \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right)
113 : * \bar{x}\bar{y}\\
114 : * \frac{\partial x_0^1}{\partial \bar{x}} &=
115 : * R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}}
116 : * \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right)
117 : * \bar{x}\bar{y}\\
118 : * \frac{\partial x_0^1}{\partial \bar{y}} &=
119 : * R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} +
120 : * R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}}
121 : * \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right)
122 : * \bar{y}^2\\
123 : * \frac{\partial x_0^i}{\partial \bar{z}} &=0.
124 : * \f}
125 : *
126 : * ### Evaluating sinc functions
127 : *
128 : * Note that \f$\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\f$ and
129 : * its derivative appear in the above equations. We evaluate
130 : * \f$\sin(ax)/x\f$ in a straightforward way, except we are careful
131 : * to evaluate \f$\sin(ax)/x = a\f$ for the special case \f$x=0\f$.
132 : *
133 : * The derivative of the sync function is more complicated to evaluate
134 : * because of roundoff. Note that we can expand
135 : *
136 : * \f{align*}
137 : * \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &=
138 : * \frac{a}{x^2}\left(1 - 1 - \frac{2 (ax)^2}{3!} +
139 : * \frac{4(ax)^4}{5!} - \frac{5(ax)^6}{7!} + \cdots \right),
140 : * \f}
141 : *
142 : * where we kept the "1 - 1" above as a reminder that when evaluating
143 : * this function directly as \f$(a \cos(ax)/x^2 - \sin(ax)/x^3)\f$, there
144 : * can be significant roundoff because of the "1" in each of the two
145 : * terms that are subtracted. The relative error in direct evaluation is
146 : * \f$3\epsilon/(ax)^2\f$, where \f$\epsilon\f$ is machine precision
147 : * (This expression comes from replacing "1 - 1" with \f$\epsilon\f$
148 : * and comparing the lowest-order contribution to the correct answer, i.e.
149 : * the \f$2(ax)^2/3!\f$ term, with \f$\epsilon\f$, the error contribution).
150 : * This means the error is 100% if \f$(ax)^2=3\epsilon\f$.
151 : *
152 : * To avoid roundoff, we evaluate the series if \f$ax\f$ is small
153 : * enough. Suppose we keep terms up to and including the \f$(ax)^{2n}\f$
154 : * term in the series. Then we evaluate the series if the
155 : * next term, the \f$(ax)^{2n+2}\f$ term, is roundoff,
156 : * i.e. if \f$(2n+2)(ax)^{2n+2}/(2n+3)! < \epsilon\f$. In this case,
157 : * the direct evaluation has the maximum error if
158 : * \f$(2n+2)(ax)^{2n+2}/(2n+3)! = \epsilon\f$. We showed above that the
159 : * relative error in direct evaluation is \f$3\epsilon/(ax)^2\f$,
160 : * which evaluates to \f$(\epsilon^n (2n+2)/(2n+3)!)^{1/(n+1)}\f$.
161 : *
162 : * \f{align*}
163 : * n=1 \qquad& \mathrm{error}=3(\epsilon/30)^{1/2} &\qquad
164 : * \sim \mathrm{5e-09}\\
165 : * n=2 \qquad& \mathrm{error}=3(\epsilon^2/840)^{1/3} &\qquad
166 : * \sim \mathrm{7e-12}\\
167 : * n=3 \qquad& \mathrm{error}=3(\epsilon^3/45360)^{1/4} &\qquad
168 : * \sim \mathrm{2e-13}\\
169 : * n=4 \qquad& \mathrm{error}=3(\epsilon^4/3991680)^{1/5} &\qquad
170 : * \sim \mathrm{2e-14}\\
171 : * n=5 \qquad& \mathrm{error}=3(\epsilon^5/518918400)^{1/6} &\qquad
172 : * \sim \mathrm{5e-15}\\
173 : * n=6 \qquad& \mathrm{error}=3(\epsilon^6/93405312000)^{1/7} &\qquad
174 : * \sim \mathrm{1e-15}
175 : * \f}
176 : * We gain less and less with each order, so we choose \f$n=3\f$.
177 : * Then the series above can be rewritten to this order in the form
178 : * \f{align*}
179 : * \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &=
180 : * -\frac{a^3}{3}\left(1 - \frac{3\cdot 4 (ax)^2}{5!} +
181 : * \frac{3 \cdot 6(ax)^4}{7!}\right).
182 : * \f}
183 : *
184 : * ### inverse
185 : *
186 : * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
187 : * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or boost::none if
188 : * \f$x_0^i\f$ or \f$\sigma\f$ are outside the range of the map.
189 : * For example, if \f$x_0^i\f$ does not lie on the sphere,
190 : * we return boost::none.
191 : *
192 : * The easiest to compute is \f$\bar{z}\f$, which is given by inverting
193 : * Eq. (4):
194 : *
195 : * \f{align}
196 : * \bar{z} &= 2\sigma - 1.
197 : * \f}
198 : *
199 : * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ then we return
200 : * boost::none.
201 : *
202 : * To get \f$\bar{x}\f$ and \f$\bar{y}\f$,
203 : * we invert
204 : * Eqs (1--3). If \f$x_0^0=x_0^1=0\f$, then \f$\bar{x}=\bar{y}=0\f$.
205 : * Otherwise, we compute
206 : *
207 : * \f{align}
208 : * \bar{\rho} = \theta_\mathrm{max}^{-1}
209 : * \tan^{-1}\left(\frac{\rho}{x_0^2-C^2}\right),
210 : * \f}
211 : *
212 : * where \f$\rho^2 = (x_0^0-C^0)^2+(x_0^1-C^1)^2\f$. Then
213 : *
214 : * \f{align}
215 : * \bar{x} &= (x_0^0-C^0)\frac{\bar{\rho}}{\rho},\\
216 : * \bar{y} &= (x_0^1-C^1)\frac{\bar{\rho}}{\rho}.
217 : * \f}
218 : *
219 : * Note that if \f$\bar{x}^2+\bar{y}^2 > 1\f$, the original point is outside
220 : * the range of the map so we return boost::none.
221 : *
222 : * ### lambda_tilde
223 : *
224 : * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
225 : * \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
226 : *
227 : * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
228 : *
229 : * Since \f$x_0^i\f$ must lie on the sphere, \f$\tilde{\lambda}\f$ is the
230 : * solution of the quadratic equation
231 : *
232 : * \f{align}
233 : * |P^i + (x^i - P^i) \tilde{\lambda} - C^i |^2 - R^2 = 0.
234 : * \f}
235 : *
236 : * In solving the quadratic, we demand a root that is positive and
237 : * less than or equal to unity, since \f$x_0^i\f$ is always between
238 : * the projection point and \f$x^i\f$. If there are two suitable
239 : * roots, this means that the entire sphere lies between \f$P^i\f$ and
240 : * \f$x^i\f$; in this case if \f$x^2 \geq z_\mathrm{P}\f$ we choose the
241 : * larger root, otherwise we choose the smaller one: this gives
242 : * us the root with \f$x_0^2 \geq z_\mathrm{P}\f$, the portion of the sphere
243 : * that is the range of `Endcap`. If there is no suitable root,
244 : * this means that the point \f$x^i\f$ is not in the range of the map
245 : * so we return a default-constructed std::optional.
246 : *
247 : * ### deriv_lambda_tilde
248 : *
249 : * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
250 : * \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
251 : * returns \f$\partial \tilde{\lambda}/\partial x^i\f$.
252 : * By differentiating Eq. (11), we find
253 : *
254 : * \f{align}
255 : * \frac{\partial\tilde{\lambda}}{\partial x^j} &=
256 : * \tilde{\lambda}^2 \frac{C^j - x_0^j}{|x_0^i - P^i|^2
257 : * + (x_0^i - P^i)(P_i - C_{i})}.
258 : * \f}
259 : *
260 : * ### inv_jacobian
261 : *
262 : * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
263 : * where \f$\sigma\f$ is held fixed.
264 : * The arguments to `inv_jacobian`
265 : * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
266 : *
267 : * Note that \f$\bar{x}\f$ and \f$\bar{y}\f$ can be considered to
268 : * depend only on \f$x_0^0\f$ and \f$x_0^1\f$ but not on \f$x_0^2\f$,
269 : * because the point \f$x_0^i\f$ is constrained to lie on a sphere of
270 : * radius \f$R\f$. Note that there is an alternative way to compute
271 : * Eqs. (8) and (9) using only \f$x_0^0\f$ and \f$x_0^1\f$. To do
272 : * this, define
273 : *
274 : * \f{align}
275 : * \upsilon \equiv \sin(\bar{\rho}\theta_\mathrm{max})
276 : * &= \sqrt{\frac{(x_0^0-C^0)^2+(x_0^1-C^1)^2}{R^2}}.
277 : * \f}
278 : *
279 : * Then we can write
280 : *
281 : * \f{align}
282 : * \frac{1}{\bar{\rho}}\sin(\bar{\rho}\theta_\mathrm{max})
283 : * &= \frac{\theta_\mathrm{max}\upsilon}{\arcsin(\upsilon)},
284 : * \f}
285 : *
286 : * so that
287 : *
288 : * \f{align}
289 : * \bar{x} &= \frac{x_0^0-C^0}{R}\left(\frac{1}{\bar{\rho}}
290 : * \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1} \\
291 : * \bar{y} &= \frac{x_0^1-C^1}{R}\left(\frac{1}{\bar{\rho}}
292 : * \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1}.
293 : * \f}
294 : *
295 : * We will compute \f$\partial \bar{x}^i/\partial
296 : * x_0^k\f$ by differentiating Eqs. (15) and (16). Because those equations
297 : * involve \f$\bar{\rho}\f$, we first establish some relations
298 : * involving derivatives of \f$\bar{\rho}\f$. For ease of notation, we define
299 : *
300 : * \f{align}
301 : * q \equiv \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}.
302 : * \f}
303 : *
304 : * First observe that
305 : * \f{align}
306 : * \frac{dq}{d\upsilon}
307 : * = \frac{dq}{d\bar{\rho}}
308 : * \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1},
309 : * \f}
310 : *
311 : * where \f$\upsilon\f$ is the quantity defined by Eq. (13). Therefore
312 : *
313 : * \f{align}
314 : * \frac{\partial q}{\partial x_0^0} &=
315 : * \frac{\bar{x}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}}
316 : * \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1},\\
317 : * \frac{\partial q}{\partial x_0^1} &=
318 : * \frac{\bar{y}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}}
319 : * \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1},
320 : * \f}
321 : *
322 : * where we have differentiated Eq. (13), and where we have
323 : * used Eqs. (15) and (16) to eliminate \f$x_0^0\f$ and
324 : * \f$x_0^1\f$ in favor of \f$\bar{x}\f$ and
325 : * \f$\bar{y}\f$ in the final result.
326 : *
327 : * Let
328 : * \f{align}
329 : * \Sigma \equiv \frac{1}{\bar\rho} \frac{dq}{d\bar{\rho}},
330 : * \f}
331 : * since that combination will appear frequently in formulas below. Note that
332 : * \f$\Sigma\f$ has a finite limit as \f$\bar{\rho}\to 0\f$, and it is evaluated
333 : * according to the section on evaluating sinc functions above.
334 : *
335 : * By differentiating Eqs. (15) and (16), and using Eqs. (19) and (20), we
336 : * find
337 : *
338 : * \f{align}
339 : * \frac{\partial \bar{x}}{\partial x_0^0} &=
340 : * \frac{1}{R q}
341 : * - \frac{\bar{x}^2 \Sigma}{R q}
342 : * \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\
343 : * \frac{\partial \bar{x}}{\partial x_0^1} &=
344 : * - \frac{\bar{x}\bar{y} \Sigma}{R q}
345 : * \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\
346 : * \frac{\partial \bar{x}}{\partial x_0^2} &= 0,\\
347 : * \frac{\partial \bar{y}}{\partial x_0^0} &=
348 : * \frac{\partial \bar{x}}{\partial x_0^1},\\
349 : * \frac{\partial \bar{y}}{\partial x_0^1} &=
350 : * \frac{1}{R q}
351 : * - \frac{\bar{y}^2 \Sigma}{R q}
352 : * \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\
353 : * \frac{\partial \bar{y}}{\partial x_0^2} &= 0,\\
354 : * \frac{\partial \bar{z}}{\partial x_0^i} &= 0.
355 : * \f}
356 : * Note that care must be taken to evaluate
357 : * \f$q = \sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\f$ and its
358 : * derivative \f$\Sigma\f$ near \f$\bar{\rho}=0\f$; see the discussion above on
359 : * evaluating sinc functions.
360 : *
361 : * ### dxbar_dsigma
362 : *
363 : * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
364 : * where \f$x_0^i\f$ is held fixed.
365 : *
366 : * From Eq. (6) we have
367 : *
368 : * \f{align}
369 : * \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2).
370 : * \f}
371 : *
372 : */
373 1 : class Endcap {
374 : public:
375 0 : Endcap(const std::array<double, 3>& center, double radius, double z_plane);
376 :
377 0 : Endcap() = default;
378 0 : ~Endcap() = default;
379 0 : Endcap(Endcap&&) = default;
380 0 : Endcap(const Endcap&) = default;
381 0 : Endcap& operator=(const Endcap&) = default;
382 0 : Endcap& operator=(Endcap&&) = default;
383 :
384 : template <typename T>
385 0 : void forward_map(gsl::not_null<std::array<T, 3>*> target_coords,
386 : const std::array<T, 3>& source_coords) const;
387 :
388 : /// The inverse function is only callable with doubles because the inverse
389 : /// might fail if called for a point out of range, and it is unclear
390 : /// what should happen if the inverse were to succeed for some points in a
391 : /// DataVector but fail for other points.
392 1 : std::optional<std::array<double, 3>> inverse(
393 : const std::array<double, 3>& target_coords, double sigma_in) const;
394 :
395 : template <typename T>
396 0 : void jacobian(gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> jacobian_out,
397 : const std::array<T, 3>& source_coords) const;
398 :
399 : template <typename T>
400 0 : void inv_jacobian(
401 : gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> inv_jacobian_out,
402 : const std::array<T, 3>& source_coords) const;
403 :
404 : template <typename T>
405 0 : void sigma(gsl::not_null<T*> sigma_out,
406 : const std::array<T, 3>& source_coords) const;
407 :
408 : template <typename T>
409 0 : void deriv_sigma(gsl::not_null<std::array<T, 3>*> deriv_sigma_out,
410 : const std::array<T, 3>& source_coords) const;
411 :
412 : template <typename T>
413 0 : void dxbar_dsigma(gsl::not_null<std::array<T, 3>*> dxbar_dsigma_out,
414 : const std::array<T, 3>& source_coords) const;
415 :
416 0 : std::optional<double> lambda_tilde(
417 : const std::array<double, 3>& parent_mapped_target_coords,
418 : const std::array<double, 3>& projection_point,
419 : bool source_is_between_focus_and_target) const;
420 :
421 : template <typename T>
422 0 : void deriv_lambda_tilde(
423 : gsl::not_null<std::array<T, 3>*> deriv_lambda_tilde_out,
424 : const std::array<T, 3>& target_coords, const T& lambda_tilde,
425 : const std::array<double, 3>& projection_point) const;
426 :
427 : // NOLINTNEXTLINE(google-runtime-references)
428 0 : void pup(PUP::er& p);
429 :
430 0 : static bool is_identity() { return false; }
431 :
432 0 : static constexpr bool supports_hessian{false};
433 :
434 : private:
435 0 : friend bool operator==(const Endcap& lhs, const Endcap& rhs);
436 0 : std::array<double, 3> center_{};
437 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
438 0 : double theta_max_{std::numeric_limits<double>::signaling_NaN()};
439 : };
440 0 : bool operator!=(const Endcap& lhs, const Endcap& rhs);
441 : } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
|
