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
|
