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 "Parallel/FifoCache.hpp"
19 : #include "Utilities/Gsl.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 : * \warning The Shape map uses a mutable spherepack cache that can be read and
217 : * mutated concurrently. It is up to the user to guarantee that
218 : * the cache is in a valid state while threads are calling the map methods.
219 : * Special care needs to be taken when using move and copy assignment.
220 : */
221 1 : class Shape {
222 : public:
223 0 : using FunctionsOfTimeMap = std::unordered_map<
224 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
225 :
226 0 : explicit Shape(
227 : const std::array<double, 3>& center, double truncation_limit,
228 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
229 : transition_func,
230 : std::string shape_function_of_time_name,
231 : std::optional<std::string> size_function_of_time_name = std::nullopt);
232 :
233 0 : Shape() = default;
234 0 : ~Shape() = default;
235 0 : Shape(const Shape& rhs);
236 0 : Shape& operator=(const Shape& rhs);
237 0 : Shape(Shape&& rhs);
238 0 : Shape& operator=(Shape&& rhs);
239 :
240 : template <typename T>
241 0 : std::array<T, 3> operator()(
242 : const std::array<T, 3>& source_coords, double time,
243 : const FunctionsOfTimeMap& functions_of_time) const;
244 :
245 0 : std::optional<std::array<double, 3>> inverse(
246 : const std::array<double, 3>& target_coords, double time,
247 : const FunctionsOfTimeMap& functions_of_time) const;
248 :
249 : template <typename T>
250 0 : std::array<T, 3> frame_velocity(
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<T, 3, Frame::NoFrame> jacobian(
256 : const std::array<T, 3>& source_coords, double time,
257 : const FunctionsOfTimeMap& functions_of_time) const;
258 :
259 : template <typename T>
260 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
261 : const std::array<T, 3>& source_coords, double time,
262 : const FunctionsOfTimeMap& functions_of_time) const;
263 :
264 : /*!
265 : * \brief An optimized call that computes the target coordinates, frame
266 : * velocity and jacobian at once to avoid duplicate calculations.
267 : *
268 : * \details The first argument `source_and_target_coords` should contain
269 : * the source coordinates and will be overwritten in place with the target
270 : * coordinates.
271 : */
272 1 : void coords_frame_velocity_jacobian(
273 : gsl::not_null<std::array<DataVector, 3>*> source_and_target_coords,
274 : gsl::not_null<std::array<DataVector, 3>*> frame_vel,
275 : gsl::not_null<tnsr::Ij<DataVector, 3, Frame::NoFrame>*> jac, double time,
276 : const FunctionsOfTimeMap& functions_of_time) const;
277 :
278 : // NOLINTNEXTLINE(google-runtime-references)
279 0 : void pup(PUP::er& p);
280 0 : static bool is_identity() { return false; }
281 0 : static constexpr bool supports_hessian{false};
282 0 : static constexpr size_t dim = 3;
283 :
284 0 : const std::unordered_set<std::string>& function_of_time_names() const {
285 : return f_of_t_names_;
286 : }
287 :
288 : private:
289 0 : std::string shape_f_of_t_name_;
290 0 : std::optional<std::string> size_f_of_t_name_;
291 0 : std::unordered_set<std::string> f_of_t_names_;
292 0 : std::array<double, 3> center_{};
293 0 : double truncation_limit_{0.};
294 : std::unique_ptr<ShapeMapTransitionFunctions::ShapeMapTransitionFunction>
295 0 : transition_func_;
296 :
297 0 : using SpherepackEntry = std::pair<ylm::SpherepackIterator, ylm::Spherepack>;
298 0 : using SpherepackCache = Parallel::FifoCache<std::unique_ptr<SpherepackEntry>>;
299 0 : using CachedItem = SpherepackCache::Cached;
300 0 : static constexpr size_t cache_capacity_ = 5;
301 :
302 : // NOLINTNEXTLINE(spectre-mutable)
303 0 : mutable SpherepackCache spherepack_cache_{cache_capacity_};
304 :
305 : template <typename T>
306 0 : std::array<T, 3> center_coordinates(const std::array<T, 3>& coords) const {
307 : return {coords[0] - center_[0], coords[1] - center_[1],
308 : coords[2] - center_[2]};
309 : }
310 :
311 : template <typename T>
312 0 : void center_coordinates(gsl::not_null<std::array<T, 3>*> result,
313 : const std::array<T, 3>& coords) const {
314 : for (size_t i = 0; i < 3; ++i) {
315 : gsl::at(*result, i) = gsl::at(coords, i) - gsl::at(center_, i);
316 : }
317 : }
318 :
319 : template <typename T>
320 0 : void jacobian_helper(
321 : gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> result,
322 : const ylm::Spherepack::InterpolationInfo<T>& interpolation_info,
323 : const DataVector& extended_coefs, const std::array<T, 3>& centered_coords,
324 : const T& radial_distortion, const T& transition_func,
325 : const ylm::Spherepack& ylm) const;
326 :
327 0 : void check_size(const gsl::not_null<DataVector*>& coefs,
328 : const FunctionsOfTimeMap& functions_of_time, double time,
329 : bool use_deriv) const;
330 :
331 0 : size_t find_truncated_l_max(const DataVector& coefs,
332 : const DataVector& coef_derivs,
333 : const DataVector& coef_dderivs,
334 : ylm::SpherepackIterator iterator) const;
335 :
336 0 : CachedItem get_spherepack_cache_entry(size_t l_max) const;
337 :
338 0 : friend bool operator==(const Shape& lhs, const Shape& rhs);
339 : };
340 0 : bool operator!=(const Shape& lhs, const Shape& rhs);
341 :
342 : } // namespace domain::CoordinateMaps::TimeDependent
|
