|
SpECTRE
v2023.10.11
|
Distorts a distribution of points radially according to a spherical harmonic expansion while preserving angles. More...
#include <Shape.hpp>
Public Types | |
| using | FunctionsOfTimeMap = std::unordered_map< std::string, std::unique_ptr< domain::FunctionsOfTime::FunctionOfTime > > |
Public Member Functions | |
| Shape (const std::array< double, 3 > ¢er, size_t l_max, size_t m_max, std::unique_ptr< ShapeMapTransitionFunctions::ShapeMapTransitionFunction > transition_func, std::string shape_function_of_time_name, std::optional< std::string > size_function_of_time_name=std::nullopt) | |
| Shape (Shape &&)=default | |
| Shape & | operator= (Shape &&)=default |
| Shape (const Shape &rhs) | |
| Shape & | operator= (const Shape &rhs) |
| template<typename T > | |
| std::array< tt::remove_cvref_wrap_t< T >, 3 > | operator() (const std::array< T, 3 > &source_coords, double time, const FunctionsOfTimeMap &functions_of_time) const |
| std::optional< std::array< double, 3 > > | inverse (const std::array< double, 3 > &target_coords, double time, const FunctionsOfTimeMap &functions_of_time) const |
| template<typename T > | |
| std::array< tt::remove_cvref_wrap_t< T >, 3 > | frame_velocity (const std::array< T, 3 > &source_coords, double time, const FunctionsOfTimeMap &functions_of_time) const |
| template<typename T > | |
| tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | jacobian (const std::array< T, 3 > &source_coords, double time, const FunctionsOfTimeMap &functions_of_time) const |
| template<typename T > | |
| tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame > | inv_jacobian (const std::array< T, 3 > &source_coords, double time, const FunctionsOfTimeMap &functions_of_time) const |
| void | pup (PUP::er &p) |
| const std::unordered_set< std::string > & | function_of_time_names () const |
Static Public Member Functions | |
| static bool | is_identity () |
Static Public Attributes | |
| static constexpr size_t | dim = 3 |
Friends | |
| bool | operator== (const Shape &lhs, const Shape &rhs) |
Distorts a distribution of points radially according to a spherical harmonic expansion while preserving angles.
The shape map distorts the distance \(r\) between a point and the center while leaving the angles \(\theta\), \(\phi\) between them preserved by applying a spherical harmonic expansion with time-dependent coefficients \(\lambda_{lm}(t)\). There are two ways to specify the time-dependent coefficients \(\lambda_{lm}(t)\):
ylm::Spherepack::spectral_size() number of components. These are in Spherepack order and should be the Spherepack coefficients, not the spherical harmonic coefficients. See the note below. To use this, set the size_function_of_time_name argument of the constructor to std::nullopt.size_function_of_time_name argument of the constructor to the name of the FunctionOfTime that's in the cache. This method is useful if we have control systems because we have a separate control system controlling a separate function of time for the \(l = 0\) coefficient than we do for the other coefficients.shape_function_of_time_name argument in the constructor that must always be specified) are not the complex spherical-harmonic coefficients \(\lambda_{lm}(t)\), but instead are the real-valued SPHEREPACK coefficients \(a_{lm}(t)\) and \(b_{lm}(t)\) used by Spherepack. This is the same for both methods of specifying FunctionOfTimes above. The relationship between these two sets of coefficients is \begin{align} a_{l0} & = \sqrt{\frac{2}{\pi}}\lambda_{l0}&\qquad l\geq 0,\\ a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(\lambda_{lm}) &\qquad l\geq 1, m\geq 1, \\ b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(\lambda_{lm}) &\qquad l\geq 1, m\geq 1. \end{align}
The "shape" FunctionOfTime stores coefficients only for non-negative \(m\); this is because the function we are expanding is real, so the coefficients for \(m<0\) can be obtained from \(m>0\) coefficients by complex conjugation. If thesize_function_of_time_name argument is given to the constructor, then it is asserted that the \(l=0\) coefficient of the "shape" function of time is exactly 0. The \(l=0\) coefficient is then controlled by the "size" FunctionOfTime. Unlike the "shape" FunctionOfTime, the quantity in the "size" FunctionOfTime is the "complex" spherical harmonic coefficient \(\lambda_{00}(t)\), and not the SPHEREPACK coefficient \(a_{00}(t)\) ("complex" is in quotes because all \(m=0\) coefficients are always real.) Here and below we write the equations in terms of \(\lambda_{lm}(t)\) instead of \(a_{lm}(t)\) and \(b_{lm}(t)\), regardless of which FunctionOfTime representation we are using, because the resulting expressions are much shorter.An additional dimensionless domain-dependent transition function \(f(r, \theta, \phi)\) ensures that the distortion falls off correctly to zero at a certain boundary (must be a block boundary). The function \(f(r, \theta, \phi)\) is restricted such that
\begin{equation} 0 \leq f(r, \theta, \phi) \leq 1 \end{equation}
Given a point with cartesian coordinates \(\xi^i\), let the polar coordinates \((r, \theta, \phi)\) with respect to a center \(x_c^i\) be defined in the usual way:
\begin{align} \xi^0 - x_c^0 &= r \sin(\theta) \cos(\phi)\\ \xi^1 - x_c^1 &= r \sin(\theta) \sin(\phi)\\ \xi^2 - x_c^2 &= r \cos(\theta) \end{align}
The shape map maps the unmapped coordinates \(\xi^i\) to coordinates \(x^i\):
\begin{equation} x^i = \xi^i - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi). \end{equation}
Or written another way
\begin{equation}\label{eq:map_center_plus_distortion} x^i = x_c^i + (\xi^i - x_c^i) \left(1 - \frac{f(r, \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right). \end{equation}
The form in Eq. \(\ref{eq:map_center_plus_distortion}\) makes two things clearer
The inverse map is given by:
\begin{equation} \xi^i = x_c^i + (x^i-x_c^i)*(r/\tilde{r}), \end{equation}
where \(\tilde{r}\) is the radius of \(\xi\), calculated by the transition map. In order to compute \(r/\tilde{r}\), the following equation must be solved
\begin{equation} \frac{r}{\tilde{r}} = \frac{1}{1-f(r,\theta,\phi)\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r} \end{equation}
For more details, see ShapeMapTransitionFunction::original_radius_over_radius .
The frame velocity \(v^i\ = dx^i / dt\) is calculated trivially:
\begin{equation} v^i = - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm} \dot{\lambda}_{lm}(t)Y_{lm}(\theta, \phi). \end{equation}
The Jacobian is given by:
\begin{align} \frac{\partial x^i}{\partial \xi^j} = \delta_j^i &\left( 1 - \frac{f(r, \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right) \nonumber \\ &- (\xi^i - x_c^i) \left[\left(\frac{1}{r}\frac{\partial}{\partial \xi^j} f(r, \theta, \phi) - \frac{\xi_j}{r^3} f(r, \theta, \phi)\right) \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi) + \frac{f(r, \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t) \frac{\partial}{\partial \xi^j} Y_{lm}(\theta, \phi) \right]. \end{align}
where \(\xi_j = \xi^j\).
The inverse Jacobian is computed by numerically inverting the Jacobian.
For future optimization, the interpolation_info objects calculated in all functions of this class could be cached. Since every element should evaluate the same grid coordinates most time steps, this might greatly decrease computation. Every element has their own clone of the shape map so the caching could be done with member variables. Care must be taken that jacobian currently calculates the interpolation_info with an order higher.