|
SpECTRE
v2024.04.12
|
A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects portions of two spherical surfaces. More...
#include <FocallyLiftedEndcap.hpp>
Public Member Functions | |
| Endcap (const std::array< double, 3 > ¢er, double radius, double z_plane) | |
| Endcap (Endcap &&)=default | |
| Endcap (const Endcap &)=default | |
| Endcap & | operator= (const Endcap &)=default |
| Endcap & | operator= (Endcap &&)=default |
| template<typename T > | |
| void | forward_map (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > target_coords, const std::array< T, 3 > &source_coords) const |
| std::optional< std::array< double, 3 > > | inverse (const std::array< double, 3 > &target_coords, double sigma_in) const |
| The inverse function is only callable with doubles because the inverse might fail if called for a point out of range, and it is unclear what should happen if the inverse were to succeed for some points in a DataVector but fail for other points. | |
| template<typename T > | |
| void | jacobian (const gsl::not_null< tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > * > jacobian_out, const std::array< T, 3 > &source_coords) const |
| template<typename T > | |
| void | inv_jacobian (const gsl::not_null< tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > * > inv_jacobian_out, const std::array< T, 3 > &source_coords) const |
| template<typename T > | |
| void | sigma (const gsl::not_null< tt::remove_cvref_wrap_t< T > * > sigma_out, const std::array< T, 3 > &source_coords) const |
| template<typename T > | |
| void | deriv_sigma (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > deriv_sigma_out, const std::array< T, 3 > &source_coords) const |
| template<typename T > | |
| void | dxbar_dsigma (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > dxbar_dsigma_out, const std::array< T, 3 > &source_coords) const |
| std::optional< double > | lambda_tilde (const std::array< double, 3 > &parent_mapped_target_coords, const std::array< double, 3 > &projection_point, bool source_is_between_focus_and_target) const |
| template<typename T > | |
| void | deriv_lambda_tilde (const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > deriv_lambda_tilde_out, const std::array< T, 3 > &target_coords, const T &lambda_tilde, const std::array< double, 3 > &projection_point) const |
| void | pup (PUP::er &p) |
Static Public Member Functions | |
| static bool | is_identity () |
Friends | |
| bool | operator== (const Endcap &lhs, const Endcap &rhs) |
A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects portions of two spherical surfaces.
Because FocallyLiftedEndcap is a FocallyLiftedInnerMap, it is meant to be a template parameter of FocallyLiftedMap, and its member functions are meant to be used by FocallyLiftedMap. See FocallyLiftedMap for further documentation.
The domain of the map is a 3D unit right cylinder with coordinates \((\bar{x},\bar{y},\bar{z})\) such that \(-1\leq\bar{z}\leq 1\) and \(\bar{x}^2+\bar{y}^2 \leq 1\). The range of the map has coordinates \((x,y,z)\).
Consider a sphere with center \(C^i\) and radius \(R\) that is intersected by a plane normal to the \(z\) axis and located at \(z = z_\mathrm{P}\). In the figure above, every point \(\bar{x}^i\) in the blue region \(\sigma=0\) maps to a point \(x_0^i\) on a portion of the surface of the sphere.
Endcap provides the following functions:
forward_map maps \((\bar{x},\bar{y},\bar{z}=-1)\) to the portion of the sphere with \(z \geq z_\mathrm{P}\). The arguments to forward_map are \((\bar{x},\bar{y},\bar{z})\), but \(\bar{z}\) is ignored. forward_map returns \(x_0^i\), the 3D coordinates on that sphere, which are given by
\begin{align} x_0^0 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max}) \bar{x}}{\bar{\rho}} + C^0,\\ x_0^1 &= R \frac{\sin(\bar{\rho} \theta_\mathrm{max}) \bar{y}}{\bar{\rho}} + C^1,\\ x_0^2 &= R \cos(\bar{\rho} \theta_\mathrm{max}) + C^2. \end{align}
Here \(\bar{\rho}^2 \equiv (\bar{x}^2+\bar{y}^2)/\bar{R}^2\), where \(\bar{R}\) is the radius of the cylinder in barred coordinates, which is always unity, and where \(\theta_\mathrm{max}\) is defined by \(\cos(\theta_\mathrm{max}) = (z_\mathrm{P}-C^2)/R\). Note that when \(\bar{\rho}=0\), we must evaluate \(\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\) as \(\theta_\mathrm{max}\).
\(\sigma\) is a function that is zero at \(\bar{z}=-1\) (which maps onto the sphere \(x^i=x_0^i\)) and unity at \(\bar{z}=+1\) (corresponding to the upper surface of the FocallyLiftedMap). We define
\begin{align} \sigma &= \frac{\bar{z}+1}{2}. \end{align}
deriv_sigma returns
\begin{align} \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2). \end{align}
jacobian returns \(\partial x_0^k/\partial \bar{x}^j\). The arguments to jacobian are \((\bar{x},\bar{y},\bar{z})\), but \(\bar{z}\) is ignored.
Differentiating Eqs.(1–3) above yields
\begin{align*} \frac{\partial x_0^2}{\partial \bar{x}} &= - R \theta_\mathrm{max} \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{x}\\ \frac{\partial x_0^2}{\partial \bar{y}} &= - R \theta_\mathrm{max} \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\bar{y}\\ \frac{\partial x_0^0}{\partial \bar{x}} &= R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} + R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}^2\\ \frac{\partial x_0^0}{\partial \bar{y}} &= R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}\bar{y}\\ \frac{\partial x_0^1}{\partial \bar{x}} &= R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{x}\bar{y}\\ \frac{\partial x_0^1}{\partial \bar{y}} &= R \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}} + R \frac{1}{\bar{\rho}}\frac{d}{d\bar{\rho}} \left(\frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}\right) \bar{y}^2\\ \frac{\partial x_0^i}{\partial \bar{z}} &=0. \end{align*}
Note that \(\sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\) and its derivative appear in the above equations. We evaluate \(\sin(ax)/x\) in a straightforward way, except we are careful to evaluate \(\sin(ax)/x = a\) for the special case \(x=0\).
The derivative of the sync function is more complicated to evaluate because of roundoff. Note that we can expand
\begin{align*} \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &= \frac{a}{x^2}\left(1 - 1 - \frac{2 (ax)^2}{3!} + \frac{4(ax)^4}{5!} - \frac{5(ax)^6}{7!} + \cdots \right), \end{align*}
where we kept the "1 - 1" above as a reminder that when evaluating this function directly as \((a \cos(ax)/x^2 - \sin(ax)/x^3)\), there can be significant roundoff because of the "1" in each of the two terms that are subtracted. The relative error in direct evaluation is \(3\epsilon/(ax)^2\), where \(\epsilon\) is machine precision (This expression comes from replacing "1 - 1" with \(\epsilon\) and comparing the lowest-order contribution to the correct answer, i.e. the \(2(ax)^2/3!\) term, with \(\epsilon\), the error contribution). This means the error is 100% if \((ax)^2=3\epsilon\).
To avoid roundoff, we evaluate the series if \(ax\) is small enough. Suppose we keep terms up to and including the \((ax)^{2n}\) term in the series. Then we evaluate the series if the next term, the \((ax)^{2n+2}\) term, is roundoff, i.e. if \((2n+2)(ax)^{2n+2}/(2n+3)! < \epsilon\). In this case, the direct evaluation has the maximum error if \((2n+2)(ax)^{2n+2}/(2n+3)! = \epsilon\). We showed above that the relative error in direct evaluation is \(3\epsilon/(ax)^2\), which evaluates to \((\epsilon^n (2n+2)/(2n+3)!)^{1/(n+1)}\).
\begin{align*} n=1 \qquad& \mathrm{error}=3(\epsilon/30)^{1/2} &\qquad \sim \mathrm{5e-09}\\ n=2 \qquad& \mathrm{error}=3(\epsilon^2/840)^{1/3} &\qquad \sim \mathrm{7e-12}\\ n=3 \qquad& \mathrm{error}=3(\epsilon^3/45360)^{1/4} &\qquad \sim \mathrm{2e-13}\\ n=4 \qquad& \mathrm{error}=3(\epsilon^4/3991680)^{1/5} &\qquad \sim \mathrm{2e-14}\\ n=5 \qquad& \mathrm{error}=3(\epsilon^5/518918400)^{1/6} &\qquad \sim \mathrm{5e-15}\\ n=6 \qquad& \mathrm{error}=3(\epsilon^6/93405312000)^{1/7} &\qquad \sim \mathrm{1e-15} \end{align*}
We gain less and less with each order, so we choose \(n=3\). Then the series above can be rewritten to this order in the form
\begin{align*} \frac{1}{x}\frac{d}{dx}\left(\frac{\sin(ax)}{x}\right) &= -\frac{a^3}{3}\left(1 - \frac{3\cdot 4 (ax)^2}{5!} + \frac{3 \cdot 6(ax)^4}{7!}\right). \end{align*}
inverse takes \(x_0^i\) and \(\sigma\) as arguments, and returns \((\bar{x},\bar{y},\bar{z})\), or boost::none if \(x_0^i\) or \(\sigma\) are outside the range of the map. For example, if \(x_0^i\) does not lie on the sphere, we return boost::none.
The easiest to compute is \(\bar{z}\), which is given by inverting Eq. (4):
\begin{align} \bar{z} &= 2\sigma - 1. \end{align}
If \(\bar{z}\) is outside the range \([-1,1]\) then we return boost::none.
To get \(\bar{x}\) and \(\bar{y}\), we invert Eqs (1–3). If \(x_0^0=x_0^1=0\), then \(\bar{x}=\bar{y}=0\). Otherwise, we compute
\begin{align} \bar{\rho} = \theta_\mathrm{max}^{-1} \tan^{-1}\left(\frac{\rho}{x_0^2-C^2}\right), \end{align}
where \(\rho^2 = (x_0^0-C^0)^2+(x_0^1-C^1)^2\). Then
\begin{align} \bar{x} &= (x_0^0-C^0)\frac{\bar{\rho}}{\rho},\\ \bar{y} &= (x_0^1-C^1)\frac{\bar{\rho}}{\rho}. \end{align}
Note that if \(\bar{x}^2+\bar{y}^2 > 1\), the original point is outside the range of the map so we return boost::none.
lambda_tilde takes as arguments a point \(x^i\) and a projection point \(P^i\), and computes \(\tilde{\lambda}\), the solution to
\begin{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\end{align}
Since \(x_0^i\) must lie on the sphere, \(\tilde{\lambda}\) is the solution of the quadratic equation
\begin{align} |P^i + (x^i - P^i) \tilde{\lambda} - C^i |^2 - R^2 = 0. \end{align}
In solving the quadratic, we demand a root that is positive and less than or equal to unity, since \(x_0^i\) is always between the projection point and \(x^i\). If there are two suitable roots, this means that the entire sphere lies between \(P^i\) and \(x^i\); in this case if \(x^2 \geq z_\mathrm{P}\) we choose the larger root, otherwise we choose the smaller one: this gives us the root with \(x_0^2 \geq z_\mathrm{P}\), the portion of the sphere that is the range of Endcap. If there is no suitable root, this means that the point \(x^i\) is not in the range of the map so we return a default-constructed std::optional.
deriv_lambda_tilde takes as arguments \(x_0^i\), a projection point \(P^i\), and \(\tilde{\lambda}\), and returns \(\partial \tilde{\lambda}/\partial x^i\). By differentiating Eq. (11), we find
\begin{align} \frac{\partial\tilde{\lambda}}{\partial x^j} &= \tilde{\lambda}^2 \frac{C^j - x_0^j}{|x_0^i - P^i|^2 + (x_0^i - P^i)(P_i - C_{i})}. \end{align}
inv_jacobian returns \(\partial \bar{x}^i/\partial x_0^k\), where \(\sigma\) is held fixed. The arguments to inv_jacobian are \((\bar{x},\bar{y},\bar{z})\), but \(\bar{z}\) is ignored.
Note that \(\bar{x}\) and \(\bar{y}\) can be considered to depend only on \(x_0^0\) and \(x_0^1\) but not on \(x_0^2\), because the point \(x_0^i\) is constrained to lie on a sphere of radius \(R\). Note that there is an alternative way to compute Eqs. (8) and (9) using only \(x_0^0\) and \(x_0^1\). To do this, define
\begin{align} \upsilon \equiv \sin(\bar{\rho}\theta_\mathrm{max}) &= \sqrt{\frac{(x_0^0-C^0)^2+(x_0^1-C^1)^2}{R^2}}. \end{align}
Then we can write
\begin{align} \frac{1}{\bar{\rho}}\sin(\bar{\rho}\theta_\mathrm{max}) &= \frac{\theta_\mathrm{max}\upsilon}{\arcsin(\upsilon)}, \end{align}
so that
\begin{align} \bar{x} &= \frac{x_0^0-C^0}{R}\left(\frac{1}{\bar{\rho}} \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1} \\ \bar{y} &= \frac{x_0^1-C^1}{R}\left(\frac{1}{\bar{\rho}} \sin(\bar{\rho}\theta_\mathrm{max})\right)^{-1}. \end{align}
We will compute \(\partial \bar{x}^i/\partial x_0^k\) by differentiating Eqs. (15) and (16). Because those equations involve \(\bar{\rho}\), we first establish some relations involving derivatives of \(\bar{\rho}\). For ease of notation, we define
\begin{align} q \equiv \frac{\sin(\bar{\rho}\theta_\mathrm{max})}{\bar{\rho}}. \end{align}
First observe that
\begin{align} \frac{dq}{d\upsilon} = \frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1}, \end{align}
where \(\upsilon\) is the quantity defined by Eq. (13). Therefore
\begin{align} \frac{\partial q}{\partial x_0^0} &= \frac{\bar{x}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1},\\ \frac{\partial q}{\partial x_0^1} &= \frac{\bar{y}}{\bar{\rho}R}\frac{dq}{d\bar{\rho}} \left(\bar{\rho} \frac{dq}{d\bar{\rho}} + q\right)^{-1}, \end{align}
where we have differentiated Eq. (13), and where we have used Eqs. (15) and (16) to eliminate \(x_0^0\) and \(x_0^1\) in favor of \(\bar{x}\) and \(\bar{y}\) in the final result.
Let
\begin{align} \Sigma \equiv \frac{1}{\bar\rho} \frac{dq}{d\bar{\rho}}, \end{align}
since that combination will appear frequently in formulas below. Note that \(\Sigma\) has a finite limit as \(\bar{\rho}\to 0\), and it is evaluated according to the section on evaluating sinc functions above.
By differentiating Eqs. (15) and (16), and using Eqs. (19) and (20), we find
\begin{align} \frac{\partial \bar{x}}{\partial x_0^0} &= \frac{1}{R q} - \frac{\bar{x}^2 \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{x}}{\partial x_0^1} &= - \frac{\bar{x}\bar{y} \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{x}}{\partial x_0^2} &= 0,\\ \frac{\partial \bar{y}}{\partial x_0^0} &= \frac{\partial \bar{x}}{\partial x_0^1},\\ \frac{\partial \bar{y}}{\partial x_0^1} &= \frac{1}{R q} - \frac{\bar{y}^2 \Sigma}{R q} \left(\bar{\rho}^2 \Sigma + q\right)^{-1},\\ \frac{\partial \bar{y}}{\partial x_0^2} &= 0,\\ \frac{\partial \bar{z}}{\partial x_0^i} &= 0. \end{align}
Note that care must be taken to evaluate \(q = \sin(\bar{\rho}\theta_\mathrm{max})/\bar{\rho}\) and its derivative \(\Sigma\) near \(\bar{\rho}=0\); see the discussion above on evaluating sinc functions.
dxbar_dsigma returns \(\partial \bar{x}^i/\partial \sigma\), where \(x_0^i\) is held fixed.
From Eq. (6) we have
\begin{align} \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2). \end{align}