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 : * An additional dimensionless domain-dependent transition function
90 : * \f$f(r, \theta, \phi)\f$ ensures that the distortion falls off correctly to
91 : * zero at a certain boundary (must be a block boundary). The function \f$f(r,
92 : * \theta, \phi)\f$ is restricted such that
93 : *
94 : * \f{equation}{
95 : * 0 \leq f(r, \theta, \phi) \leq 1
96 : * \f}
97 : *
98 : * ### Mapped coordinates
99 : *
100 : * Given a point with cartesian coordinates \f$\xi^i\f$, let the polar
101 : * coordinates \f$(r, \theta, \phi)\f$ with respect to a center \f$x_c^i\f$ be
102 : * defined in the usual way:
103 : *
104 : * \f{align}{
105 : * \xi^0 - x_c^0 &= r \sin(\theta) \cos(\phi)\\
106 : * \xi^1 - x_c^1 &= r \sin(\theta) \sin(\phi)\\
107 : * \xi^2 - x_c^2 &= r \cos(\theta)
108 : * \f}
109 : *
110 : * The shape map maps the unmapped
111 : * coordinates \f$\xi^i\f$ to coordinates \f$x^i\f$:
112 : *
113 : * \f{equation}{
114 : * x^i = \xi^i - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm}
115 : * \lambda_{lm}(t)Y_{lm}(\theta, \phi).
116 : * \f}
117 : *
118 : * Or written another way
119 : *
120 : * \f{equation}{\label{eq:map_center_plus_distortion}
121 : * x^i = x_c^i + (\xi^i - x_c^i) \left(1 - \frac{f(r, \theta, \phi)}{r}
122 : * \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right).
123 : * \f}
124 : *
125 : * The form in Eq. \f$\ref{eq:map_center_plus_distortion}\f$ makes two things
126 : * clearer
127 : *
128 : * 1. This shape map is just a radial distortion about \f$x_c^i\f$
129 : * 2. The coefficients \f$\lambda_{lm}\f$ have units of distance because
130 : * \f$\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r\f$ must be dimensionless
131 : * (because \f$f\f$ is dimensionless).
132 : *
133 : * ### Inverse map
134 : *
135 : * The inverse map is given by:
136 : * \f{equation}{
137 : * \xi^i = x_c^i + (x^i-x_c^i)*(r/\tilde{r}),
138 : * \f}
139 : * where \f$\tilde{r}\f$ is the radius of \f$\xi\f$, calculated by the
140 : * transition map. In order to compute \f$r/\tilde{r}\f$, the following equation
141 : * must be solved
142 : *
143 : * \f{equation}{
144 : * \frac{r}{\tilde{r}} =
145 : * \frac{1}{1-f(r,\theta,\phi)\sum\lambda_{lm}(t)Y_{lm}(\theta,\phi) / r}
146 : * \f}
147 : *
148 : * For more details, see
149 : * \link domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
150 : * ShapeMapTransitionFunction::original_radius_over_radius \endlink.
151 : *
152 : * ### Frame velocity
153 : *
154 : * The frame velocity \f$v^i\ = dx^i / dt\f$ is calculated trivially:
155 : * \f{equation}{
156 : * v^i = - (\xi^i - x_c^i) \frac{f(r, \theta, \phi)}{r} \sum_{lm}
157 : * \dot{\lambda}_{lm}(t)Y_{lm}(\theta, \phi).
158 : * \f}
159 : *
160 : * ### Jacobian
161 : *
162 : * The Jacobian is given by:
163 : * \f{align}{
164 : * \frac{\partial x^i}{\partial \xi^j} = \delta_j^i &\left( 1 - \frac{f(r,
165 : * \theta, \phi)}{r} \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)\right)
166 : * \nonumber \\
167 : * &- (\xi^i - x_c^i)
168 : * \left[\left(\frac{1}{r}\frac{\partial}{\partial \xi^j} f(r, \theta, \phi) -
169 : * \frac{\xi_j}{r^3} f(r, \theta, \phi)\right) \sum_{lm}
170 : * \lambda_{lm}(t)Y_{lm}(\theta, \phi) + \frac{f(r, \theta, \phi)}{r}
171 : * \sum_{lm} \lambda_{lm}(t) \frac{\partial}{\partial \xi^j} Y_{lm}(\theta,
172 : * \phi) \right].
173 : * \f}
174 : *
175 : * where \f$\xi_j = \xi^j\f$.
176 : *
177 : * ### Inverse Jacobian
178 : *
179 : * The inverse Jacobian is computed by numerically inverting the Jacobian.
180 : *
181 : * For future optimization, the `interpolation_info` objects calculated in all
182 : * functions of this class could be cached. Since every element should evaluate
183 : * the same grid coordinates most time steps, this might greatly decrease
184 : * computation. Every element has their own clone of the shape map so the
185 : * caching could be done with member variables. Care must be taken that
186 : * `jacobian` currently calculates the `interpolation_info` with an order
187 : * higher.
188 : */
189 1 : class Shape {
190 : public:
191 0 : using FunctionsOfTimeMap = std::unordered_map<
192 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
193 :
194 0 : explicit Shape(
195 : const std::array<double, 3>& center, size_t l_max, size_t m_max,
196 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
197 : transition_func,
198 : std::string shape_function_of_time_name,
199 : std::optional<std::string> size_function_of_time_name = std::nullopt);
200 :
201 0 : Shape() = default;
202 0 : ~Shape() = default;
203 0 : Shape(Shape&&) = default;
204 0 : Shape& operator=(Shape&&) = default;
205 0 : Shape(const Shape& rhs);
206 0 : Shape& operator=(const Shape& rhs);
207 :
208 : template <typename T>
209 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
210 : const std::array<T, 3>& source_coords, double time,
211 : const FunctionsOfTimeMap& functions_of_time) const;
212 :
213 0 : std::optional<std::array<double, 3>> inverse(
214 : const std::array<double, 3>& target_coords, double time,
215 : const FunctionsOfTimeMap& functions_of_time) const;
216 :
217 : template <typename T>
218 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> frame_velocity(
219 : const std::array<T, 3>& source_coords, double time,
220 : const FunctionsOfTimeMap& functions_of_time) const;
221 :
222 : template <typename T>
223 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
224 : const std::array<T, 3>& source_coords, double time,
225 : const FunctionsOfTimeMap& functions_of_time) const;
226 :
227 : template <typename T>
228 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
229 : const std::array<T, 3>& source_coords, double time,
230 : const FunctionsOfTimeMap& functions_of_time) const;
231 :
232 : // NOLINTNEXTLINE(google-runtime-references)
233 0 : void pup(PUP::er& p);
234 0 : static bool is_identity() { return false; }
235 0 : static constexpr size_t dim = 3;
236 :
237 0 : const std::unordered_set<std::string>& function_of_time_names() const {
238 : return f_of_t_names_;
239 : }
240 :
241 : private:
242 0 : std::string shape_f_of_t_name_;
243 0 : std::optional<std::string> size_f_of_t_name_;
244 0 : std::unordered_set<std::string> f_of_t_names_;
245 0 : std::array<double, 3> center_{};
246 0 : size_t l_max_ = 2;
247 0 : size_t m_max_ = 2;
248 0 : ylm::Spherepack ylm_{2, 2};
249 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
250 0 : transition_func_;
251 :
252 : template <typename T>
253 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> center_coordinates(
254 : const std::array<T, 3>& coords) const {
255 : return {coords[0] - center_[0], coords[1] - center_[1],
256 : coords[2] - center_[2]};
257 : }
258 :
259 0 : void check_size(const gsl::not_null<DataVector*>& coefs,
260 : const FunctionsOfTimeMap& functions_of_time, double time,
261 : bool use_deriv) const;
262 :
263 : // Checks that the vector of coefficients has the right size and that the
264 : // monopole and dipole coefficients are zero.
265 0 : void check_coefficients(const DataVector& coefs) const;
266 :
267 : template <typename T>
268 0 : T check_and_compute_one_over_radius(
269 : const std::array<T, 3>& centered_coords) const;
270 :
271 0 : friend bool operator==(const Shape& lhs, const Shape& rhs);
272 : };
273 0 : bool operator!=(const Shape& lhs, const Shape& rhs);
274 :
275 : } // namespace domain::CoordinateMaps::TimeDependent
|