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 <limits>
       9             : #include <optional>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : 
      13             : /// \cond
      14             : namespace PUP {
      15             : class er;
      16             : }  // namespace PUP
      17             : /// \endcond
      18             : 
      19             : namespace domain {
      20             : namespace CoordinateMaps {
      21             : 
      22             : /*!
      23             :  * \ingroup CoordinateMapsGroup
      24             :  *
      25             :  * \brief Redistributes gridpoints within the unit sphere.
      26             :  * \image html SpecialMobius.png "A sphere with a `mu` of 0.25."
      27             :  *
      28             :  * \details A special case of the conformal Mobius transformation that
      29             :  * maps the unit ball to itself. This map depends on a single
      30             :  * parameter, `mu` \f$ = \mu\f$, which is the x-coordinate of the preimage
      31             :  * of the origin under this map. This map has the fixed points \f$x=1\f$ and
      32             :  * \f$x=-1\f$. The map is singular for \f$\mu=1\f$ but we have found that this
      33             :  * map is accurate up to 12 decimal places for values of \f$\mu\f$ up to 0.96.
      34             :  *
      35             :  * We define the auxiliary variables
      36             :  * \f[ r := \sqrt{x^2 + y^2 +z^2}\f]
      37             :  * and
      38             :  * \f[ \lambda := \frac{1}{1 - 2 x \mu + \mu^2 r^2}\f]
      39             :  *
      40             :  * The map corresponding to this transformation in cartesian coordinates
      41             :  * is then given by:
      42             :  *
      43             :  * \f[\vec{x}'(x,y,z) =
      44             :  * \lambda\begin{bmatrix}
      45             :  * x(1+\mu^2) - \mu(1+r^2)\\
      46             :  * y(1-\mu^2)\\
      47             :  * z(1-\mu^2)\\
      48             :  * \end{bmatrix}\f]
      49             :  *
      50             :  * The inverse map is the same as the forward map with \f$\mu\f$
      51             :  * replaced by \f$-\mu\f$.
      52             :  *
      53             :  * This map is intended to be used only inside the unit sphere.  A
      54             :  * point inside the unit sphere maps to another point inside the unit
      55             :  * sphere. The map can have undesirable behavior at certain points
      56             :  * outside the unit sphere: The map is singular at
      57             :  * \f$(x,y,z) = (1/\mu, 0, 0)\f$ (which is outside the unit sphere
      58             :  * since \f$|\mu| < 1\f$). Moreover, a point on the \f$x\f$-axis
      59             :  * arbitrarily close to the singularity maps to an arbitrarily large
      60             :  * value on the \f$\pm x\f$-axis, where the sign depends on which side
      61             :  * of the singularity the point is on.
      62             :  *
      63             :  * A general Mobius transformation is a function on the complex plane, and
      64             :  * takes the form \f$ f(z) = \frac{az+b}{cz+d}\f$, where
      65             :  * \f$z, a, b, c, d \in \mathbb{C}\f$, and \f$ad-bc\neq 0\f$.
      66             :  *
      67             :  * The special case used in this map is the function
      68             :  * \f$ f(z) = \frac{z - \mu}{1 - z\mu}\f$. This has the desired properties:
      69             :  * - The unit disk in the complex plane is mapped to itself.
      70             :  *
      71             :  * - The x-axis is mapped to itself.
      72             :  *
      73             :  * - \f$f(\mu) = 0\f$.
      74             :  *
      75             :  * The three-dimensional version of this map is obtained by rotating the disk
      76             :  * in the plane about the x-axis.
      77             :  *
      78             :  * This map is useful for performing transformations along the x-axis
      79             :  * that preserve the unit disk. A concrete example of this is in the BBH
      80             :  * domain, where two BBHs with a center-of-mass at x=\f$\mu\f$ can be shifted
      81             :  * such that the new center of mass is now located at x=0. Additionally,
      82             :  * the spherical shape of the outer wave-zone is preserved and, as a mobius
      83             :  * map, the spherical coordinate shapes of the black holes is also preserved.
      84             :  */
      85           1 : class SpecialMobius {
      86             :  public:
      87           0 :   static constexpr size_t dim = 3;
      88           0 :   explicit SpecialMobius(double mu);
      89           0 :   SpecialMobius() = default;
      90           0 :   ~SpecialMobius() = default;
      91           0 :   SpecialMobius(SpecialMobius&&) = default;
      92           0 :   SpecialMobius(const SpecialMobius&) = default;
      93           0 :   SpecialMobius& operator=(const SpecialMobius&) = default;
      94           0 :   SpecialMobius& operator=(SpecialMobius&&) = default;
      95             : 
      96             :   template <typename T>
      97           0 :   std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
      98             : 
      99             :   /// Returns std::nullopt for target_coords outside the unit sphere.
     100             :   /// The inverse function is only callable with doubles because the inverse
     101             :   /// might fail if called for a point out of range, and it is unclear
     102             :   /// what should happen if the inverse were to succeed for some points in a
     103             :   /// DataVector but fail for other points.
     104           1 :   std::optional<std::array<double, 3>> inverse(
     105             :       const std::array<double, 3>& target_coords) const;
     106             : 
     107             :   template <typename T>
     108           0 :   tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
     109             :       const std::array<T, 3>& source_coords) const;
     110             : 
     111             :   template <typename T>
     112           0 :   tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
     113             :       const std::array<T, 3>& source_coords) const;
     114             : 
     115             :   // NOLINTNEXTLINE(google-runtime-references)
     116           0 :   void pup(PUP::er& p);
     117             : 
     118           0 :   bool is_identity() const { return is_identity_; }
     119             : 
     120           0 :   static constexpr bool supports_hessian{true};
     121             : 
     122             :  private:
     123             :   template <typename T>
     124           0 :   std::array<T, 3> mobius_distortion(const std::array<T, 3>& coords,
     125             :                                      double mu) const;
     126             :   template <typename T>
     127           0 :   tnsr::Ij<T, 3, Frame::NoFrame> mobius_distortion_jacobian(
     128             :       const std::array<T, 3>& coords, double mu) const;
     129           0 :   friend bool operator==(const SpecialMobius& lhs, const SpecialMobius& rhs);
     130             : 
     131           0 :   double mu_{std::numeric_limits<double>::signaling_NaN()};
     132           0 :   bool is_identity_{false};
     133             : };
     134           0 : bool operator!=(const SpecialMobius& lhs, const SpecialMobius& rhs);
     135             : }  // namespace CoordinateMaps
     136             : }  // namespace domain