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