Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <optional>
9 : #include <string>
10 : #include <unordered_set>
11 :
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
14 : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
15 : #include "Utilities/Gsl.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 : /*!
30 : * \ingroup CoordMapsTimeDependentGroup
31 : * \brief Distorts a distribution of points radially according to a spherical
32 : * harmonic expansion while preserving angles.
33 : *
34 : * \details The shape map distorts the distance \f$r\f$ between a point and
35 : * the center while leaving the angles \f$\theta\f$, \f$\phi\f$ between them
36 : * preserved by applying a spherical harmonic expansion with time-dependent
37 : * coefficients \f$\lambda_{lm}(t)\f$. There are two ways to specify the
38 : * time-dependent coefficients \f$\lambda_{lm}(t)\f$:
39 : *
40 : * 1. A single FunctionOfTime which specifies all coefficients. This
41 : * FunctionOfTime should have `ylm::Spherepack::spectral_size()` number of
42 : * components. These are in Spherepack order and should be the Spherepack
43 : * coefficients, *not* the spherical harmonic coefficients. See the note
44 : * below. To use this, set the `size_function_of_time_name` argument of the
45 : * constructor to `std::nullopt`.
46 : * 2. Two different FunctionOfTime%s. The first is similar to 1.) in that it
47 : * should have the same number of components, be in Spherepack order, and be
48 : * the Spherepack coefficients. The only difference is that the \f$l = 0\f$
49 : * coefficient should be identically 0. The second FunctionOfTime should have
50 : * a single component which will be the \f$l = 0\f$ coefficient. This
51 : * component should be stored as the spherical harmonic coefficient and *not*
52 : * a Spherepack coefficient. See the note below. To use this method, set the
53 : * `size_function_of_time_name` argument of the constructor to the name of
54 : * the FunctionOfTime that's in the cache. This method is useful if we have
55 : * control systems because we have a separate control system controlling a
56 : * separate function of time for the \f$l = 0\f$ coefficient than we do for
57 : * the other coefficients.
58 : *
59 : * \note The quantities stored in the "shape" FunctionOfTime (the
60 : * `shape_function_of_time_name` argument in the constructor that must always be
61 : * specified) are ***not*** the complex spherical-harmonic coefficients
62 : * \f$\lambda_{lm}(t)\f$, but instead are the real-valued SPHEREPACK
63 : * coefficients \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$ used by Spherepack. This
64 : * is the same for both methods of specifying FunctionOfTime%s above. The
65 : * relationship between these two sets of coefficients is
66 : * \f{align}
67 : * a_{l0} & = \sqrt{\frac{2}{\pi}}\lambda_{l0}&\qquad l\geq 0,\\
68 : * a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(\lambda_{lm})
69 : * &\qquad l\geq 1, m\geq 1, \\
70 : * b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(\lambda_{lm})
71 : * &\qquad l\geq 1, m\geq 1.
72 : * \f}
73 : * The "shape" FunctionOfTime stores coefficients only for non-negative \f$m\f$;
74 : * this is because the function we are expanding is real, so the
75 : * coefficients for \f$m<0\f$ can be obtained from \f$m>0\f$ coefficients by
76 : * complex conjugation.
77 : * If the `size_function_of_time_name` argument is given to the constructor,
78 : * then it is asserted that the \f$l=0\f$ coefficient of the "shape" function of
79 : * time is exactly 0. The \f$l=0\f$ coefficient is then controlled by the "size"
80 : * FunctionOfTime. Unlike the "shape" FunctionOfTime, the quantity in the
81 : * "size" FunctionOfTime ***is*** the "complex" spherical harmonic coefficient
82 : * \f$\lambda_{00}(t)\f$, and not the SPHEREPACK coefficient \f$a_{00}(t)\f$
83 : * ("complex" is in quotes because all \f$m=0\f$ coefficients are always real.)
84 : * Here and below we write the equations in terms of \f$\lambda_{lm}(t)\f$
85 : * instead of \f$a_{lm}(t)\f$ and \f$b_{lm}(t)\f$, regardless of which
86 : * FunctionOfTime representation we are using, because the resulting expressions
87 : * are much shorter.
88 : *
89 : * \parblock
90 : *
91 : * \note Also note that the FunctionOfTime coefficients $\lambda_{lm}(t)$ are
92 : * stored as *negative* of the coefficients you'd retrieve from a
93 : * `Strahlkorper`. This is because you would typically represent the expansion
94 : * of a strahlkorper as $S(r) = +\sum S_{lm} Y_{lm}$. However, in equation
95 : * $\ref{eq:map_form_2}$ there is a minus sign on the $\sum \lambda_{lm}
96 : * Y_{lm}$, not a plus sign. Therefore, $\lambda_{lm}(t)$ picks up an extra
97 : * factor of $-1$. This is purely a choice of convention.
98 : *
99 : * \endparblock
100 : *
101 : * An additional domain-dependent transition function
102 : *
103 : * \begin{equation}
104 : * G(r,\theta,\phi) = \frac{f(r,\theta,\phi)}{r}
105 : * \end{equation}
106 : *
107 : * ensures that the distortion falls off correctly to zero at a certain boundary
108 : * (must be a block boundary). The dimensionless function \f$f(r, \theta,
109 : * \phi)\f$ is restricted such that
110 : *
111 : * \f{equation}{
112 : * 0 \leq f(r, \theta, \phi) \leq 1
113 : * \f}
114 : *
115 : * ### Mapped coordinates
116 : *
117 : * Given a point with cartesian coordinates \f$\xi^i\f$, let the polar
118 : * coordinates \f$(r, \theta, \phi)\f$ with respect to a center \f$x_c^i\f$ be
119 : * defined in the usual way:
120 : *
121 : * \f{align}{
122 : * \xi^0 - x_c^0 &= r \sin(\theta) \cos(\phi)\\
123 : * \xi^1 - x_c^1 &= r \sin(\theta) \sin(\phi)\\
124 : * \xi^2 - x_c^2 &= r \cos(\theta)
125 : * \f}
126 : *
127 : * The shape map maps the unmapped
128 : * coordinates \f$\xi^i\f$ to coordinates \f$x^i\f$:
129 : *
130 : * \f{equation}{\label{eq:map_form_1}
131 : * x^i = \xi^i - (\xi^i - x_c^i) G(r,\theta,\phi) \sum_{lm}
132 : * \lambda_{lm}(t)Y_{lm}(\theta, \phi).
133 : * \f}
134 : *
135 : * Or written another way
136 : *
137 : * \f{equation}{\label{eq:map_form_2}
138 : * x^i = x_c^i + (\xi^i - x_c^i) \left(1 - G(r,\theta,\phi)
139 : * \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right).
140 : * \f}
141 : *
142 : * The form in Eq. \f$\ref{eq:map_form_2}\f$ makes two things
143 : * clearer
144 : *
145 : * 1. This shape map is just a radial distortion about \f$x_c^i\f$
146 : * 2. The coefficients \f$\lambda_{lm}\f$ have units of distance because
147 : * \f$\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r\f$ must be dimensionless
148 : * (because \f$f\f$ is dimensionless).
149 : *
150 : * ### Inverse map
151 : *
152 : * The inverse map is given by:
153 : * \f{equation}{
154 : * \xi^i = x_c^i + (x^i-x_c^i)*(r/\tilde{r}),
155 : * \f}
156 : * where \f$\tilde{r}\f$ is the radius of $\vec{x}$, calculated by the
157 : * transition map. In order to compute $r/\tilde{r}$, the following equation
158 : * must be solved
159 : *
160 : * \f{equation}{
161 : * \frac{r}{\tilde{r}} =
162 : * \frac{1}{1-G(r,\theta,\phi)\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi)}
163 : * \f}
164 : *
165 : * For more details, see
166 : * \link domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
167 : * ShapeMapTransitionFunction::original_radius_over_radius \endlink.
168 : *
169 : * ### Frame velocity
170 : *
171 : * The frame velocity \f$v^i\ = dx^i / dt\f$ is calculated trivially:
172 : * \f{equation}{
173 : * v^i = - (\xi^i - x_c^i) G(r, \theta, \phi) \sum_{lm}
174 : * \dot{\lambda}_{lm}(t)Y_{lm}(\theta, \phi).
175 : * \f}
176 : *
177 : * ### Jacobian
178 : *
179 : * The Jacobian is given by:
180 : * \f{align}{
181 : * \frac{\partial x^i}{\partial \xi^j} = \delta_j^i &\left( 1 - G(r,\theta,\phi)
182 : * \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right) \nonumber \\
183 : * &- (\xi^i - x_c^i)
184 : * \left[\frac{\partial G(r,\theta,\phi)}{\partial\xi^j} \sum_{lm}
185 : * \lambda_{lm}(t)Y_{lm}(\theta, \phi) + G(r, \theta, \phi)
186 : * \sum_{lm} \lambda_{lm}(t) \frac{\partial}{\partial \xi^j} Y_{lm}(\theta,
187 : * \phi) \right].
188 : * \f}
189 : *
190 : * where \f$\xi_j = \xi^j\f$. It should be noted that there is an additional
191 : * factor of $1/r$ hidden in the $\partial/\partial\xi^j Y_{lm}(\theta, \phi)$
192 : * term, so the transition function $G(r,\theta,\phi)$ must have a functional
193 : * form to avoid division by zero if $r=0$.
194 : *
195 : * ### Inverse Jacobian
196 : *
197 : * The inverse Jacobian is computed by numerically inverting the Jacobian.
198 : *
199 : * For future optimization, the `interpolation_info` objects calculated in all
200 : * functions of this class could be cached. Since every element should evaluate
201 : * the same grid coordinates most time steps, this might greatly decrease
202 : * computation. Every element has their own clone of the shape map so the
203 : * caching could be done with member variables. Care must be taken that
204 : * `jacobian` currently calculates the `interpolation_info` with an order
205 : * higher.
206 : */
207 1 : class Shape {
208 : public:
209 0 : using FunctionsOfTimeMap = std::unordered_map<
210 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
211 :
212 0 : explicit Shape(
213 : const std::array<double, 3>& center, size_t l_max, size_t m_max,
214 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
215 : transition_func,
216 : std::string shape_function_of_time_name,
217 : std::optional<std::string> size_function_of_time_name = std::nullopt);
218 :
219 0 : Shape() = default;
220 0 : ~Shape() = default;
221 0 : Shape(Shape&&) = default;
222 0 : Shape& operator=(Shape&&) = default;
223 0 : Shape(const Shape& rhs);
224 0 : Shape& operator=(const Shape& rhs);
225 :
226 : template <typename T>
227 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
228 : const std::array<T, 3>& source_coords, double time,
229 : const FunctionsOfTimeMap& functions_of_time) const;
230 :
231 0 : std::optional<std::array<double, 3>> inverse(
232 : const std::array<double, 3>& target_coords, double time,
233 : const FunctionsOfTimeMap& functions_of_time) const;
234 :
235 : template <typename T>
236 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> frame_velocity(
237 : const std::array<T, 3>& source_coords, double time,
238 : const FunctionsOfTimeMap& functions_of_time) const;
239 :
240 : template <typename T>
241 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
242 : const std::array<T, 3>& source_coords, double time,
243 : const FunctionsOfTimeMap& functions_of_time) const;
244 :
245 : template <typename T>
246 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
247 : const std::array<T, 3>& source_coords, double time,
248 : const FunctionsOfTimeMap& functions_of_time) const;
249 :
250 : /*!
251 : * \brief An optimized call that computes the target coordinates, frame
252 : * velocity and jacobian at once to avoid duplicate calculations.
253 : *
254 : * \details The first argument `source_and_target_coords` should contain
255 : * the source coordinates and will be overwritten in place with the target
256 : * coordinates.
257 : */
258 1 : void coords_frame_velocity_jacobian(
259 : gsl::not_null<std::array<DataVector, 3>*> source_and_target_coords,
260 : gsl::not_null<std::array<DataVector, 3>*> frame_vel,
261 : gsl::not_null<tnsr::Ij<DataVector, 3, Frame::NoFrame>*> jac, double time,
262 : const FunctionsOfTimeMap& functions_of_time) const;
263 :
264 : // NOLINTNEXTLINE(google-runtime-references)
265 0 : void pup(PUP::er& p);
266 0 : static bool is_identity() { return false; }
267 0 : static constexpr size_t dim = 3;
268 :
269 0 : const std::unordered_set<std::string>& function_of_time_names() const {
270 : return f_of_t_names_;
271 : }
272 :
273 : private:
274 0 : std::string shape_f_of_t_name_;
275 0 : std::optional<std::string> size_f_of_t_name_;
276 0 : std::unordered_set<std::string> f_of_t_names_;
277 0 : std::array<double, 3> center_{};
278 0 : size_t l_max_ = 2;
279 0 : size_t m_max_ = 2;
280 0 : ylm::Spherepack ylm_{2, 2};
281 0 : ylm::Spherepack extended_ylm_{3, 3};
282 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
283 0 : transition_func_;
284 :
285 : template <typename T>
286 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> center_coordinates(
287 : const std::array<T, 3>& coords) const {
288 : return {coords[0] - center_[0], coords[1] - center_[1],
289 : coords[2] - center_[2]};
290 : }
291 :
292 : template <typename T>
293 0 : void center_coordinates(
294 : gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*> result,
295 : const std::array<T, 3>& coords) const {
296 : for (size_t i = 0; i < 3; ++i) {
297 : gsl::at(*result, i) = gsl::at(coords, i) - gsl::at(center_, i);
298 : }
299 : }
300 :
301 : template <typename T>
302 0 : void jacobian_helper(
303 : gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> result,
304 : const ylm::Spherepack::InterpolationInfo<T>& interpolation_info,
305 : const DataVector& extended_coefs, const std::array<T, 3>& centered_coords,
306 : const T& radial_distortion, const T& transition_func) const;
307 :
308 0 : void check_size(const gsl::not_null<DataVector*>& coefs,
309 : const FunctionsOfTimeMap& functions_of_time, double time,
310 : bool use_deriv) const;
311 :
312 : // Checks that the vector of coefficients has the right size and that the
313 : // monopole and dipole coefficients are zero.
314 0 : void check_coefficients(const DataVector& coefs) const;
315 :
316 0 : friend bool operator==(const Shape& lhs, const Shape& rhs);
317 : };
318 0 : bool operator!=(const Shape& lhs, const Shape& rhs);
319 :
320 : } // namespace domain::CoordinateMaps::TimeDependent
|
