Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <memory>
       9             : #include <optional>
      10             : #include <string>
      11             : #include <unordered_map>
      12             : #include <unordered_set>
      13             : 
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      16             : 
      17             : /// \cond
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : /// \endcond
      22             : 
      23             : namespace domain::CoordinateMaps::TimeDependent {
      24             : /*!
      25             :  * \ingroup CoordMapsTimeDependentGroup
      26             :  * \brief %Time dependent coordinate map that keeps the $y$ and $z$ coordinates
      27             :  * unchanged, but will distort (or skew) the $x$ coordinate based on functions
      28             :  * of time and radial distance to the origin.
      29             :  *
      30             :  * \details This coordinate map is only available in 3 dimensions and is
      31             :  * intended to be used in the BinaryCompactObject domain.
      32             :  *
      33             :  * ### Mapped Coordinates
      34             :  * The Skew coordinate map is given by the mapping
      35             :  *
      36             :  * \begin{align}
      37             :  *   \label{eq:map}
      38             :  *   \bar{x} &= x - W(\vec{x})\left(\tan(F_y(t))(y-y_C) +
      39             :  *       \tan(F_z(t))(z-z_C)\right) \\
      40             :  *   \bar{y} &= y \\
      41             :  *   \bar{z} &= z
      42             :  * \end{align}
      43             :  *
      44             :  * where $\vec{x}_C = (x_C, y_C, z_C)$ is the \p center of the skew map (which
      45             :  * is different than the origin of the coordinate system), and $F_y(t)$ and
      46             :  * $F_z(t)$ are the angles within the $(x,y)$ plane between the undistorted
      47             :  * $x$-axis and the skewed $\bar{y}$-axis at the origin, represented by
      48             :  * `domain::FunctionsOfTime::FunctionOfTime`s. The actual function of time
      49             :  * should have two components; the first corresponds to $y$ and the second
      50             :  * corresponds to $z$.
      51             :  *
      52             :  * $W(\vec{x})$ is a spatial function that should be 1 at $\vec{x}_C$ (i.e.
      53             :  * maximally skewed between the two objects) and fall off to 0 at the \p
      54             :  * outer_radius, $R$. Typically the \p outer_radius should be set to the
      55             :  * envelope radius for the `domain::creators::BinaryCompactObject`. The reason
      56             :  * that $W(\vec{x})$ *should* be 1 at $\vec{x}_C$, and doesn't *need* to be 1 at
      57             :  * $\vec{x}_C$ is because having $W(\vec{x})$ centered at the origin is much
      58             :  * better for the smoothness of higher derivatives of $W(\vec{x})$ than if it
      59             :  * were centered at $\vec{x}_C$. The important part is that the map is left
      60             :  * invariant along the $x=x_C$ line and that $W(\vec{x}_C)\approx 1$ for the
      61             :  * skew control system, both of which are satisfied if $W(\vec{x})$ is centered
      62             :  * at the origin and $|\vec{x}_C|\ll R$.
      63             :  *
      64             :  * With that in mind, $W$ is chosen to be
      65             :  *
      66             :  * \begin{equation}
      67             :  *   \label{eq:W}
      68             :  *   W(\vec{x}) = \frac{1}{2}\left(1 + \cos(\pi\lambda(\vec{x}))\right).
      69             :  * \end{equation}
      70             :  *
      71             :  * When the $\cos$ term is $-1$, then $W = 0$, and when the $\cos$ term is 1,
      72             :  * then $W = 1$. Therefore, the function $\lambda(\vec{x})$ must go from 0 at
      73             :  * the origin, to 1 at $R$. We choose $\lambda(\vec{x})$ to quadratically go
      74             :  * from 0 at the origin to 1 at $R$ with the form
      75             :  *
      76             :  * \begin{equation}
      77             :  *   \label{eq:lambda}
      78             :  *   \lambda(\vec{x}) = \frac{|\vec{x}|^2}{R^2}.
      79             :  * \end{equation}
      80             :  *
      81             :  * A quadratic form was chosen for $\lambda$ rather than higher powers of 2 to
      82             :  * avoid needing too much resolution to resolve a steep falloff in the function.
      83             :  *
      84             :  * \note If the quantity $S = \tan(F_y(t))(y-y_C) + \tan(F_z(t))(z-z_C)$ becomes
      85             :  * too large, then the map becomes singular because multiple $x$ will be mapped
      86             :  * to the same $\bar{x}$. This is due to the fact we are adding a linear term to
      87             :  * a cosine term with different weights. If $S$ is too large, the cosine term
      88             :  * will overpower the linear and the map will become singular.
      89             :  *
      90             :  * ### Inverse
      91             :  * To find the inverse, we need to solve a 1D root find for the $x$ component of
      92             :  * the coordinate. The inverse of the $\bar{y}$ and $\bar{z}$ coordinates are
      93             :  * trivial because there was no mapping. The equation we need to find the root
      94             :  * of is
      95             :  *
      96             :  * \begin{equation}
      97             :  *   0 = x - W(\vec{x})\left(\tan(F_y(t))(y-y_C) + \tan(F_z(t))(z-z_C)\right) -
      98             :  *   \bar{x}
      99             :  * \end{equation}
     100             :  *
     101             :  * We can bound the root by noticing that in $\ref{eq:map}$, if we substitute
     102             :  * the extremal values of $W = 0$ and $W = 1$, we get bounds on $\bar{x}$ which
     103             :  * we can turn into bounds on $x$:
     104             :  *
     105             :  * \begin{align}
     106             :  *   x &<=& \bar{x} &<=& x - (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\
     107             :  *   0 &<=& \bar{x} - x &<=& -(\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\
     108             :  *   -\bar{x} &<=& - x &<=& -\bar{x} - (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\
     109             :  *   \bar{x} &>=& x &>=& \bar{x} + (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C))
     110             :  * \end{align}
     111             :  *
     112             :  * where on each line we just made simple arithmetic operations. We pad each
     113             :  * bound by $10^{-14}$ just to avoid roundoff issues. If either of the bounds is
     114             :  * within roundoff of zero, the map is the identity at that point and we forgo
     115             :  * the root find. The root that is found is the original $x$ coordinate.
     116             :  *
     117             :  * ### Frame Velocity
     118             :  * Taking the time derivative of $\ref{eq:map}$, the frame velocity is
     119             :  *
     120             :  * \begin{align}
     121             :  *   \label{eq:frame_vel}
     122             :  *   \dot{\bar{x}} &= -W(\vec{x})\left(\dot{F}_y(t)(1 + \tan^2(F_y(t)))(y-y_C) +
     123             :  *     \dot{F}_z(t)(1 + \tan^2(F_z(t))(z-z_C))\right) \\
     124             :  *   \dot{\bar{y}} &= 0 \\
     125             :  *   \dot{\bar{z}} &= 0
     126             :  * \end{align}
     127             :  *
     128             :  * ### Jacobian and Inverse Jacobian
     129             :  * Considering the first terms in each equation of $\ref{eq:map}$, part of the
     130             :  * jacobian will be the identity matrix. The rest will come from only the $x$
     131             :  * equation. Therefore we can express the jacobian as
     132             :  *
     133             :  * \begin{equation}
     134             :  *   \frac{\partial\bar{x}^i}{\partial x^j} = \delta^i_j + {W^i}_j
     135             :  * \end{equation}
     136             :  *
     137             :  * where all components of ${W^i}_j$ are zero except the following
     138             :  *
     139             :  * \begin{align}
     140             :  *   {W^0}_0 &= \frac{\partial(\bar{x}-x)}{\partial x} &= -\frac{\partial
     141             :  *     W(\vec{x})}{\partial x}&\left(\tan(F_y(t))(y-y_C) +
     142             :  *     \tan(F_z(t))(z-z_C)\right), \\
     143             :  *   {W^0}_1 &= \frac{\partial(\bar{x}-x)}{\partial y} &= -\frac{\partial
     144             :  *     W(\vec{x})}{\partial y}&\left(\tan(F_y(t))(y-y_C) +
     145             :  *     \tan(F_z(t))(z-z_C)\right) -
     146             :  *     W\tan(F_y(t)), \\
     147             :  *   {W^0}_2 &= \frac{\partial(\bar{x}-x)}{\partial z} &= -\frac{\partial
     148             :  *     W(\vec{x})}{\partial z}&\left(\tan(F_y(t))(y-y_C) +
     149             :  *     \tan(F_z(t))(z-z_C)\right) - W\tan(F_z(t)).
     150             :  * \end{align}
     151             :  *
     152             :  * The gradient of $W(\vec{x})$ (Eq. $\ref{eq:W}$) is given by
     153             :  *
     154             :  * \begin{equation}
     155             :  *   \frac{\partial W(\vec{x})}{\partial x^i} =
     156             :  *   -\frac{\pi}{2}\frac{\partial\lambda(\vec{x})}{\partial
     157             :  *   x^i}\sin(\pi\lambda(\vec{x})).
     158             :  * \end{equation}
     159             :  *
     160             :  * The gradient of $\lambda(\vec{x})$ (Eq. $\ref{eq:lambda}$) is given by
     161             :  *
     162             :  * \begin{equation}
     163             :  *   \frac{\partial \lambda(\vec{x})}{\partial x^i} = \frac{2x^i}{R^2}
     164             :  * \end{equation}
     165             :  *
     166             :  * The inverse jacobian is computed by numerically inverting the jacobian.
     167             :  */
     168           1 : class Skew {
     169             :  public:
     170           0 :   static constexpr size_t dim = 3;
     171             : 
     172           0 :   Skew(std::string function_of_time_name, const std::array<double, 3>& center,
     173             :        double outer_radius);
     174           0 :   Skew() = default;
     175             : 
     176             :   template <typename T>
     177           0 :   std::array<T, 3> operator()(
     178             :       const std::array<T, 3>& source_coords, double time,
     179             :       const domain::FunctionsOfTimeMap& functions_of_time) const;
     180             : 
     181             :   /// The inverse function is only callable with doubles because the inverse
     182             :   /// might fail if called for a point out of range, and it is unclear
     183             :   /// what should happen if the inverse were to succeed for some points in a
     184             :   /// DataVector but fail for other points.
     185           1 :   std::optional<std::array<double, 3>> inverse(
     186             :       const std::array<double, 3>& target_coords, double time,
     187             :       const domain::FunctionsOfTimeMap& functions_of_time) const;
     188             : 
     189             :   template <typename T>
     190           0 :   std::array<T, 3> frame_velocity(
     191             :       const std::array<T, 3>& source_coords, double time,
     192             :       const domain::FunctionsOfTimeMap& functions_of_time) const;
     193             : 
     194             :   template <typename T>
     195           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     196             :       const std::array<T, 3>& source_coords, double time,
     197             :       const domain::FunctionsOfTimeMap& functions_of_time) const;
     198             : 
     199             :   template <typename T>
     200           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     201             :       const std::array<T, 3>& source_coords, double time,
     202             :       const domain::FunctionsOfTimeMap& functions_of_time) const;
     203             : 
     204             :   // NOLINTNEXTLINE(google-runtime-references)
     205           0 :   void pup(PUP::er& p);
     206             : 
     207           0 :   static bool is_identity() { return false; }
     208             : 
     209           0 :   static constexpr bool supports_hessian{false};
     210             : 
     211           0 :   const std::unordered_set<std::string>& function_of_time_names() const {
     212             :     return f_of_t_names_;
     213             :   }
     214             : 
     215             :  private:
     216             :   template <typename T>
     217           0 :   T get_width(const std::array<T, 3>& source_coords,
     218             :               bool ignore_error = false) const;
     219             : 
     220             :   template <typename T>
     221           0 :   std::array<T, 3> get_width_deriv(const std::array<T, 3>& source_coords) const;
     222             : 
     223             :   template <typename T>
     224           0 :   std::array<T, 3> map_and_velocity_helper(
     225             :       const std::array<T, 3>& source_coords, double time,
     226             :       const domain::FunctionsOfTimeMap& functions_of_time,
     227             :       bool return_velocity) const;
     228             : 
     229             :   template <typename T>
     230           0 :   void check_for_singular_map(
     231             :       const std::array<T, 3>& source_coords, double time,
     232             :       const FunctionsOfTimeMap& functions_of_time) const;
     233             : 
     234             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     235           0 :   friend bool operator==(const Skew& lhs, const Skew& rhs);
     236           0 :   std::string f_of_t_name_;
     237           0 :   std::array<double, 3> center_{};
     238           0 :   double outer_radius_{};
     239           0 :   double one_over_outer_radius_squared_{};
     240           0 :   std::unordered_set<std::string> f_of_t_names_;
     241             : };
     242             : 
     243           0 : bool operator!=(const Skew& lhs, const Skew& rhs);
     244             : 
     245             : }  // namespace domain::CoordinateMaps::TimeDependent