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 "PointwiseFunctions/MathFunctions/MathFunction.hpp"
16 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
17 :
18 : /// \cond
19 : namespace domain::FunctionsOfTime {
20 : class FunctionOfTime;
21 : } // namespace domain::FunctionsOfTime
22 : namespace PUP {
23 : class er;
24 : } // namespace PUP
25 : /// \endcond
26 :
27 : namespace domain::CoordinateMaps::TimeDependent {
28 : /*!
29 : * \ingroup CoordMapsTimeDependentGroup
30 : * \brief Translation map defined by \f$\vec{x} = \vec{\xi}+F(r)\vec{T}(t)\f$
31 : * where $F(r)$ takes on different forms based on which constructor is used.
32 : *
33 : * \details The map adds a translation to the coordinates $\vec{\xi}$ based on
34 : * what type of translation is needed. For the piecewise translation, a
35 : * translation $F(r)\vec{T}(t)$ is added to $\vec{\xi}$ based on what region
36 : * $|\vec{\xi}|$ is in. For coordinates within the inner radius, $F(r) = 1$
37 : * causing a uniform translation. Coordinates in between the inner and outer
38 : * radius have a linear radial falloff applied to them. Coordinates beyond the
39 : * outer radius have no translation applied to them $F(r) = 0$. The piecewise
40 : * translation assumes that the center of your map is at (0., 0., 0.). For the
41 : * radial MathFunction translation, a radial translation \f$F(r)\vec{T}(t)\f$ is
42 : * added to the coordinates \f$\vec{\xi}\f$, where \f$\vec{T}(t)\f$ is a
43 : * FunctionOfTime and\f$F(r)\f$ is a 1D radial MathFunction. The radius of each
44 : * point is found by subtracting the center map argument from the coordinates
45 : * \f$\vec{\xi}\f$ or the target coordinates \f$\vec{\bar{\xi}}\f$. The
46 : * Translation Map class is overloaded so that the user can choose between a
47 : * piecewise translation, radial translation or a uniform translation based on
48 : * their problem. If a radial dependence is not specified, this sets \f$F(r) =
49 : * 1\f$.
50 : *
51 : * ### Mapped Coordinates
52 : * The piecewise translation translates the coordinates $\vec{\xi}$
53 : * to the target coordinates $\vec{\bar{\xi}}$ based on the region $\vec{\xi}$
54 : * is in.
55 : * \f{equation}{
56 : * \vec{\bar{\xi}} = \left\{\begin{array}{ll}\vec{\xi} + \vec{T}(t), &
57 : * |\vec{\xi}| \leq R_{in}, \\ \vec{\xi} + wT(t), & R_{in} < |\vec{\xi}| <
58 : * R_{out}, \\ \vec{\xi}, & |\vec{\xi}| \geq R_{out} \end{array}\right.
59 : * \f}
60 : *
61 : * Where $R_{in}$ is the inner radius, $R_{out}$ is the outer radius, and $w$ is
62 : * the radial falloff factor found through
63 : * \f{equation}{
64 : * w = \frac{R_{out} - |\vec{\xi}|}{R_{out} - R_{in}}
65 : * \f}
66 : *
67 : * The radial MathFunction translation translates the coordinates
68 : * \f$\vec{\xi}\f$ to the target coordinates \f{equation}{\vec{\bar{\xi}} =
69 : * \vec{\xi} + F(r)\vec{T}(t) \f}
70 : *
71 : * If you only supply a FunctionOfTime to the constructor of this class, the
72 : * radial function will be set to 1.0 causing a uniform translation for your
73 : * coordinates. If a FunctionOfTime, MathFunction, and map center are passed in,
74 : * the radius will be found through
75 : * \f{equation}{
76 : * r = |\vec{\xi} - \vec{c}|
77 : * \f}
78 : * where r is the radius and \f$\vec{c}\f$ is the center argument.
79 : *
80 : * ### Inverse Translation
81 : * The piecewise inverse translates the coordinates
82 : * \f$\vec{\bar{\xi}}\f$ to the original coordinates based on what region
83 : * $\vec{\bar{\xi}}$ is in.
84 : * \f{equation}{
85 : * \vec{\xi} = \left\{\begin{array}{ll}\vec{\bar{\xi}} -
86 : * \vec{T}(t), & |\vec{\bar{\xi}}| \leq R_{in}, or, |\vec{\bar{\xi}} - T(t)|
87 : * \leq R_{in}, \\
88 : * \vec{\bar{\xi}} - wT(t), & R_{in} < |\vec{\bar{\xi}}| < R_{out}, \\
89 : * \vec{\bar{\xi}}, & |\vec{\bar{\xi}}| \geq R_{out}\end{array}\right.
90 : * \f}
91 : * Where $w$ is the radial falloff factor found through a quadratic solve of the
92 : * form
93 : * \f{equation}{
94 : * w^2(\vec{T}(t)^2 - (R_{out} - R_{in})^2) - 2w(\vec{T}(t)\vec{\bar{\xi}} -
95 : * R_{out}(R_{out} - R_{in})) + \vec{\bar{\xi}}^2 - R_{out}^2
96 : * \f}
97 : * The inverse map also assumes that if $\vec{\bar{\xi}}
98 : * - \vec{T}(t) \leq R_{in}$ then the translated point originally came from
99 : * within the inner radius so it'll be translated back without a quadratic
100 : * solve.
101 : *
102 : * The radial MathFunction inverse translates the coordinates
103 : * \f$\vec{\bar{\xi}}\f$ to the original coordinates using
104 : * \f{equation}{
105 : * \vec{\xi} = \vec{\bar{\xi}} - F(r)\vec{T}(t)
106 : * \f}
107 : * where \f$r^2\f$ is found as the root of
108 : * \f{equation}{
109 : * r^2 = \Big(\vec{\bar{\xi}} - \vec{c} - F(r) \vec{T}(t)\Big)^2.
110 : * \f}
111 : *
112 : * ### Frame Velocity
113 : * For the piecewise translation, the frame velocity is found through
114 : * \f{equation}{
115 : * \vec{v} = \left\{\begin{array}{ll}\frac{\vec{dT}(t)}{dt}, & |\vec{\xi}| \leq
116 : * R_{in}, \\ w\frac{\vec{dT}(t)}{dt}, & R_{in} < |\vec{\xi}| < R_{out}, \\ 0,
117 : * & |\vec{\xi}| \geq R_{out} \end{array}\right.
118 : * \f}
119 : *
120 : * For the radial MathFunction translation, the frame velocity is found through
121 : * \f{equation}{
122 : * \vec{v} = \frac{\vec{dT}(t)}{dt} F(r)
123 : * \f}
124 : * where \f$\frac{\vec{dT}(t)}{dt}\f$ is the first derivative of the
125 : * FunctionOfTime.
126 : *
127 : * ### Jacobian
128 : * For the piecewise translation, the jacobian is computed based on what region
129 : * the coordinates $\vec{\xi}$ is in.
130 : * \f{equation}{
131 : * {J^{i}}_{j} = \frac{dw}{dr} T(t)^i \frac{\xi_j}{r}, R_{in} <
132 : * |\vec{\bar{\xi}}| < R_{out}
133 : * \f}
134 : * otherwise, it will return the identity matrix.
135 : *
136 : * For the radial MathFunction translation, the jacobian is computed through the
137 : * first derivative when the radius is bigger than 1.e-13:
138 : * \f{equation}{
139 : * {J^{i}}_{j} = \frac{dF(r)}{dr} T(t)^i \frac{(\xi_j - c_j)}{r}
140 : * \f}
141 : * Where \f$\frac{dF(r)}{dr}\f$ is the first derivative of the MathFunction,
142 : * \f$\vec{\xi_j}\f$ is the source coordinates, \f$\vec{c}\f$ is the center of
143 : * your map, and r is the radius.
144 : *
145 : * At a radius smaller than 1e-13, we ASSERT that the radial MathFunction is
146 : * smooth $\frac{dF(r)}{dr} \approx 0$, so return the identity matrix.
147 : *
148 : *
149 : * ### Inverse Jacobian
150 : * The inverse jacobian is computed numerically by inverting the jacobian.
151 : */
152 : template <size_t Dim>
153 1 : class Translation {
154 : public:
155 0 : static constexpr size_t dim = Dim;
156 :
157 0 : Translation() = default;
158 0 : explicit Translation(std::string function_of_time_name);
159 :
160 0 : explicit Translation(std::string function_of_time_name, double inner_radius,
161 : double outer_radius);
162 :
163 0 : explicit Translation(
164 : std::string function_of_time_name,
165 : std::unique_ptr<MathFunction<1, Frame::Inertial>> radial_function,
166 : std::array<double, Dim>& center);
167 :
168 0 : Translation(const Translation<Dim>& Translation_Map);
169 :
170 0 : ~Translation() = default;
171 0 : Translation(Translation&&) = default;
172 0 : Translation& operator=(Translation&&) = default;
173 0 : Translation& operator=(const Translation& Translation_Map);
174 :
175 : template <typename T>
176 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
177 : const std::array<T, Dim>& source_coords, double time,
178 : const std::unordered_map<
179 : std::string,
180 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
181 : functions_of_time) const;
182 :
183 : /// The inverse function is only callable with doubles because the inverse
184 : /// might fail if called for a point out of range, and it is unclear
185 : /// what should happen if the inverse were to succeed for some points in a
186 : /// DataVector but fail for other points.
187 1 : std::optional<std::array<double, Dim>> inverse(
188 : const std::array<double, Dim>& target_coords, double time,
189 : const std::unordered_map<
190 : std::string,
191 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
192 : functions_of_time) const;
193 :
194 : template <typename T>
195 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> frame_velocity(
196 : const std::array<T, Dim>& source_coords, double time,
197 : const std::unordered_map<
198 : std::string,
199 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
200 : functions_of_time) const;
201 :
202 : template <typename T>
203 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
204 : const std::array<T, Dim>& source_coords, double time,
205 : const std::unordered_map<
206 : std::string,
207 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
208 : functions_of_time) const;
209 :
210 : template <typename T>
211 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
212 : const std::array<T, Dim>& source_coords, double time,
213 : const std::unordered_map<
214 : std::string,
215 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
216 : functions_of_time) const;
217 :
218 : // NOLINTNEXTLINE(google-runtime-references)
219 0 : void pup(PUP::er& p);
220 :
221 0 : static bool is_identity() { return false; }
222 :
223 0 : const std::unordered_set<std::string>& function_of_time_names() const {
224 : return f_of_t_names_;
225 : }
226 :
227 : private:
228 : template <size_t LocalDim>
229 0 : friend bool operator==( // NOLINT(readability-redundant-declaration)
230 : const Translation<LocalDim>& lhs, const Translation<LocalDim>& rhs);
231 :
232 : // These 2 helper functions compute the translated coordinates or frame
233 : // velocity based on the option passed in, 0 for translated coordinates, and
234 : // frame velocity for any other number.
235 : template <typename T>
236 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> math_function_helper(
237 : const std::array<T, Dim>& source_coords, double time,
238 : const std::unordered_map<
239 : std::string,
240 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
241 : functions_of_time,
242 : size_t function_or_deriv_index) const;
243 :
244 : template <typename T>
245 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> piecewise_helper(
246 : const std::array<T, Dim>& source_coords, double time,
247 : const std::unordered_map<
248 : std::string,
249 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
250 : functions_of_time,
251 : size_t function_or_deriv_index) const;
252 :
253 0 : double root_finder(const std::array<double, Dim>& distance_to_center,
254 : const DataVector& function_of_time) const;
255 :
256 0 : std::string f_of_t_name_{};
257 0 : std::unordered_set<std::string> f_of_t_names_;
258 0 : std::optional<double> inner_radius_;
259 0 : std::optional<double> outer_radius_;
260 0 : std::unique_ptr<MathFunction<1, Frame::Inertial>> f_of_r_{};
261 0 : std::array<double, Dim> center_{};
262 : };
263 :
264 : template <size_t Dim>
265 0 : inline bool operator!=(const Translation<Dim>& lhs,
266 : const Translation<Dim>& rhs) {
267 : return not(lhs == rhs);
268 : }
269 :
270 : } // namespace domain::CoordinateMaps::TimeDependent
|