domain::CoordinateMaps::TimeDependent::Shape Class Reference

Distorts a distribution of points radially according to a spherical harmonic expansion while preserving angles. More...

#include <Shape.hpp>

Public Types
using	FunctionsOfTimeMap

Public Member Functions
	Shape (const std::array< double, 3 > &center, double truncation_limit, 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 (const Shape &rhs)
Shape &	operator= (const Shape &rhs)
	Shape (Shape &&rhs)
Shape &	operator= (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	coords_frame_velocity_jacobian (gsl::not_null< std::array< DataVector, 3 > * > source_and_target_coords, gsl::not_null< std::array< DataVector, 3 > * > frame_vel, gsl::not_null< tnsr::Ij< DataVector, 3, Frame::NoFrame > * > jac, double time, const FunctionsOfTimeMap &functions_of_time) const
	An optimized call that computes the target coordinates, frame velocity and jacobian at once to avoid duplicate calculations.
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 bool	supports_hessian {false}
static constexpr size_t	dim = 3

Friends
bool	operator== (const Shape &lhs, const Shape &rhs)

Detailed Description

Distorts a distribution of points radially according to a spherical harmonic expansion while preserving angles.

Details

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)\):

1. A single FunctionOfTime which specifies all coefficients. This FunctionOfTime should have 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.
2. Two different FunctionOfTimes. The first is similar to 1.) in that it should have the same number of components, be in Spherepack order, and be the Spherepack coefficients. The only difference is that the \(l = 0\) coefficient should be identically 0. The second FunctionOfTime should have a single component which will be the \(l = 0\) coefficient. This component should be stored as the spherical harmonic coefficient and not a Spherepack coefficient. See the note below. To use this method, set the 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.

Note

The quantities stored in the "shape" FunctionOfTime (the 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 the size_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.

Note

Also note that the FunctionOfTime coefficients \(\lambda_{lm}(t)\) are stored as negative of the coefficients you'd retrieve from a 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 domain-dependent transition function

\begin{equation} G(r,\theta,\phi) = \frac{f(r,\theta,\phi)}{r} \end{equation}

ensures that the distortion falls off correctly to zero at a certain boundary (must be a block boundary). The dimensionless function \(f(r, \theta, \phi)\) is restricted such that

\begin{equation}0 \leq f(r, \theta, \phi) \leq 1 \end{equation}

Mapped coordinates

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) G(r,\theta,\phi) \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 - G(r,\theta,\phi) \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

1. This shape map is just a radial distortion about \(x_c^i\)
2. The coefficients \(\lambda_{lm}\) have units of distance because \(\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r\) must be dimensionless (because \(f\) is dimensionless).

Inverse map

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 \(\vec{x}\), 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-G(r,\theta,\phi)\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi)} \end{equation}

For more details, see ShapeMapTransitionFunction::original_radius_over_radius .

Frame velocity

The frame velocity \(v^i\ = dx^i / dt\) is calculated trivially:

\begin{equation}v^i = - (\xi^i - x_c^i) G(r, \theta, \phi) \sum_{lm} \dot{\lambda}_{lm}(t)Y_{lm}(\theta, \phi). \end{equation}

Jacobian

The Jacobian is given by:

\begin{align}\frac{\partial x^i}{\partial \xi^j} = \delta_j^i &\left( 1 - G(r,\theta,\phi) \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right) \nonumber \\ &- (\xi^i - x_c^i) \left[\frac{\partial G(r,\theta,\phi)}{\partial\xi^j} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi) + G(r, \theta, \phi) \sum_{lm} \lambda_{lm}(t) \frac{\partial}{\partial \xi^j} Y_{lm}(\theta, \phi) \right]. \end{align}

where \(\xi_j = \xi^j\). It should be noted that there is an additional factor of \(1/r\) hidden in the \(\partial/\partial\xi^j Y_{lm}(\theta, \phi)\) term, so the transition function \(G(r,\theta,\phi)\) must have a functional form to avoid division by zero if \(r=0\).

Inverse Jacobian

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.

Warning

The Shape map uses a mutable spherepack cache that can be read and mutated concurrently. It is up to the user to guarantee that the cache is in a valid state while threads are calling the map methods. Special care needs to be taken when using move and copy assignment.

Member Typedef Documentation

FunctionsOfTimeMap

using domain::CoordinateMaps::TimeDependent::Shape::FunctionsOfTimeMap

Initial value:

 std::unordered_map<
      std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>

Member Function Documentation

coords_frame_velocity_jacobian()

void domain::CoordinateMaps::TimeDependent::Shape::coords_frame_velocity_jacobian	(	gsl::not_null< std::array< DataVector, 3 > * >	source_and_target_coords,
		gsl::not_null< std::array< DataVector, 3 > * >	frame_vel,
		gsl::not_null< tnsr::Ij< DataVector, 3, Frame::NoFrame > * >	jac,
		double	time,
		const FunctionsOfTimeMap &	functions_of_time ) const

An optimized call that computes the target coordinates, frame velocity and jacobian at once to avoid duplicate calculations.

Details

The first argument source_and_target_coords should contain the source coordinates and will be overwritten in place with the target coordinates.

---

The documentation for this class was generated from the following file:

• src/Domain/CoordinateMaps/TimeDependent/Shape.hpp