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