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