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
|
