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