SpECTRE  v2025.03.17
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
domain::CoordinateMaps::TimeDependent::Skew Class Reference

Time dependent coordinate map that keeps the y and z coordinates unchanged, but will distort (or skew) the x coordinate based on functions of time and radial distance to the origin. More...

#include <Skew.hpp>

Public Member Functions

 Skew (std::string function_of_time_name, const std::array< double, 3 > &center, double outer_radius)
 
template<typename T >
std::array< tt::remove_cvref_wrap_t< T >, 3 > operator() (const std::array< T, 3 > &source_coords, double time, const domain::FunctionsOfTimeMap &functions_of_time) const
 
std::optional< std::array< double, 3 > > inverse (const std::array< double, 3 > &target_coords, double time, const domain::FunctionsOfTimeMap &functions_of_time) 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 >
std::array< tt::remove_cvref_wrap_t< T >, 3 > frame_velocity (const std::array< T, 3 > &source_coords, double time, const domain::FunctionsOfTimeMap &functions_of_time) const
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFramejacobian (const std::array< T, 3 > &source_coords, double time, const domain::FunctionsOfTimeMap &functions_of_time) const
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrameinv_jacobian (const std::array< T, 3 > &source_coords, double time, const domain::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 Skew &lhs, const Skew &rhs)
 

Detailed Description

Time dependent coordinate map that keeps the y and z coordinates unchanged, but will distort (or skew) the x coordinate based on functions of time and radial distance to the origin.

Details

This coordinate map is only available in 3 dimensions and is intended to be used in the BinaryCompactObject domain.

Mapped Coordinates

The Skew coordinate map is given by the mapping

(1)x¯=xW(x)(tan(Fy(t))(yyC)+tan(Fz(t))(zzC))(2)y¯=y(3)z¯=z

where xC=(xC,yC,zC) is the center of the skew map (which is different than the origin of the coordinate system), and Fy(t) and Fz(t) are the angles within the (x,y) plane between the undistorted x-axis and the skewed y¯-axis at the origin, represented by domain::FunctionsOfTime::FunctionOfTimes. The actual function of time should have two components; the first corresponds to y and the second corresponds to z.

W(x) is a spatial function that should be 1 at xC (i.e. maximally skewed between the two objects) and fall off to 0 at the outer_radius, R. Typically the outer_radius should be set to the envelope radius for the domain::creators::BinaryCompactObject. The reason that W(x) should be 1 at xC, and doesn't need to be 1 at xC is because having W(x) centered at the origin is much better for the smoothness of higher derivatives of W(x) than if it were centered at xC. The important part is that the map is left invariant along the x=xC line and that W(xC)1 for the skew control system, both of which are satisfied if W(x) is centered at the origin and |xC|R.

With that in mind, W is chosen to be

(4)W(x)=12(1+cos(πλ(x))).

When the cos term is 1, then W=0, and when the cos term is 1, then W=1. Therefore, the function λ(x) must go from 0 at the origin, to 1 at R. We choose λ(x) to quadratically go from 0 at the origin to 1 at R with the form

(5)λ(x)=|x|2R2.

A quadratic form was chosen for λ rather than higher powers of 2 to avoid needing too much resolution to resolve a steep falloff in the function.

Note
If the quantity S=tan(Fy(t))(yyC)+tan(Fz(t))(zzC) becomes too large, then the map becomes singular because multiple x will be mapped to the same x¯. This is due to the fact we are adding a linear term to a cosine term with different weights. If S is too large, the cosine term will overpower the linear and the map will become singular.

Inverse

To find the inverse, we need to solve a 1D root find for the x component of the coordinate. The inverse of the y¯ and z¯ coordinates are trivial because there was no mapping. The equation we need to find the root of is

(6)0=xW(x)(tan(Fy(t))(yyC)+tan(Fz(t))(zzC))x¯

We can bound the root by noticing that in 1, if we substitute the extremal values of W=0 and W=1, we get bounds on x¯ which we can turn into bounds on x:

(7)x<=x¯<=x(tan(Fy)(yyC)+tan(Fz)(zzC))(8)0<=x¯x<=(tan(Fy)(yyC)+tan(Fz)(zzC))(9)x¯<=x<=x¯(tan(Fy)(yyC)+tan(Fz)(zzC))(10)x¯>=x>=x¯+(tan(Fy)(yyC)+tan(Fz)(zzC))

where on each line we just made simple arithmetic operations. We pad each bound by 1014 just to avoid roundoff issues. If either of the bounds is within roundoff of zero, the map is the identity at that point and we forgo the root find. The root that is found is the original x coordinate.

Frame Velocity

Taking the time derivative of 1, the frame velocity is

(11)x¯˙=W(x)(F˙y(t)(1+tan2(Fy(t)))(yyC)+F˙z(t)(1+tan2(Fz(t))(zzC)))(12)y¯˙=0(13)z¯˙=0

Jacobian and Inverse Jacobian

Considering the first terms in each equation of 1, part of the jacobian will be the identity matrix. The rest will come from only the x equation. Therefore we can express the jacobian as

(14)x¯ixj=δji+Wij

where all components of Wij are zero except the following

(15)W00=(x¯x)x=W(x)x(tan(Fy(t))(yyC)+tan(Fz(t))(zzC)),(16)W01=(x¯x)y=W(x)y(tan(Fy(t))(yyC)+tan(Fz(t))(zzC))Wtan(Fy(t)),(17)W02=(x¯x)z=W(x)z(tan(Fy(t))(yyC)+tan(Fz(t))(zzC))Wtan(Fz(t)).

The gradient of W(x) (Eq. 4) is given by

(18)W(x)xi=π2λ(x)xisin(πλ(x)).

The gradient of λ(x) (Eq. 5) is given by

(19)λ(x)xi=2xiR2

The inverse jacobian is computed by numerically inverting the jacobian.


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