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