SpECTRE
v2024.05.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 timedependent coefficients \(\lambda_{lm}(t)\). There are two ways to specify the timedependent 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 sphericalharmonic coefficients \(\lambda_{lm}(t)\), but instead are the realvalued 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 nonnegative \(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.Strahlkorper
. This is because you would typically represent the expansion of a strahlkorper as \(S(r) = +\sum S_{lm} Y_{lm}\). However, in equation \(\ref{eq:map_form_2}\) there is a minus sign on the \(\sum \lambda_{lm} Y_{lm}\), not a plus sign. Therefore, \(\lambda_{lm}(t)\) picks up an extra factor of \(1\). This is purely a choice of convention.An additional dimensionless domaindependent 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}\label{eq:map_form_1} 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_form_2} 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_form_2}\) makes two things clearer
The inverse map is given by:
\begin{equation} \xi^i = x_c^i + (x^ix_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}{1f(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.