Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <optional>
8 : #include <pup.h>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "Utilities/Serialization/CharmPupable.hpp"
12 :
13 : namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
14 :
15 : /*!
16 : * \brief Abstract base class for the transition function $G(r,\theta,\phi)$
17 : * used by the domain::CoordinateMaps::TimeDependent::Shape map.
18 : *
19 : * \details This base class defines the required methods of a transition
20 : * function that is defined as
21 : *
22 : * \begin{equation}
23 : * G(r,\theta,\phi) = \frac{f(r,\theta,\phi)}{r}
24 : * \end{equation}
25 : *
26 : * where $0 <= f(r,\theta,\phi) <= 1$. Different domains may require the shape
27 : * map to fall off towards the boundary in different ways. This behavior is
28 : * controlled by the transition function. It is also needed to find the inverse
29 : * of the shape map. Since the shape map preserves angles, the problem of
30 : * finding its inverse reduces to the 1-dimensional problem of finding the
31 : * original radius from the mapped radius. The mapped radius $\tilde{r}$ is
32 : * related to the original $r$ radius by:
33 : * \begin{equation}
34 : * \label{eq:shape_map_radius}
35 : * \tilde{r} = r \left(1 - G(r,\theta,\phi)\sum_{lm}
36 : * \lambda_{lm}(t)Y_{lm}(\theta,\phi)\right),
37 : * \end{equation}
38 : * Depending on the format of the transition function, it should be possible to
39 : * analytically derive this map's inverse because it preserves angles and shifts
40 : * only the radius of each point. Otherwise the inverse has to be computed
41 : * numerically.
42 : *
43 : * The transition function must also be able to compute the gradient and the
44 : * value of the function divided by the radius. Care must be taken that this
45 : * does not divide by zero.
46 : *
47 : * All member functions with the exception of `original_radius_over_radius`
48 : * exist as overloads for types `double` and `DataVector` so that they work with
49 : * the templated shape map methods calling them. To avoid code duplication these
50 : * can be forwarded to templated implementation methods held by the derived
51 : * classes only.
52 : *
53 : * For an example, see SphereTransition.
54 : *
55 : * #### Design Decisions:
56 : *
57 : * It was decided to make the ShapeMapTransitionFunction an abstract base class
58 : * with overloads for types `double` and `DataVector` corresponding to the
59 : * template parameter `T` of the shape map's methods. The shape map holds a
60 : * `unique_ptr` to this abstract class using a common dynamic dispatch design
61 : * pattern. This approach avoids templating the shape map all together.
62 : *
63 : * An alternative approach would be to directly template the transition
64 : * functions onto the shape map so that no abstract base class is necessary.
65 : * These approaches can also be combined by making the transition function an
66 : * abstract base class but also templating it onto the shape map. In this way
67 : * the shape map does not need to hold a `unique_ptr` but can hold the
68 : * transition function directly as a member.
69 : */
70 1 : class ShapeMapTransitionFunction : public PUP::able {
71 : public:
72 0 : ShapeMapTransitionFunction() = default;
73 :
74 : /// @{
75 : /*!
76 : * Evaluate the transition function $G(r,\theta,\phi)$ at the
77 : * Cartesian coordinates `source_coords` and possibly divide by the radius
78 : * $r$.
79 : *
80 : * \details If \p one_over_radius_power has a value, then divide $G$ by $r$ to
81 : * that power. This is done here because the transition function knows how to
82 : * handle points in various regions of the map.
83 : */
84 1 : virtual double operator()(
85 : const std::array<double, 3>& source_coords,
86 : const std::optional<size_t>& one_over_radius_power) const = 0;
87 1 : virtual DataVector operator()(
88 : const std::array<DataVector, 3>& source_coords,
89 : const std::optional<size_t>& one_over_radius_power) const = 0;
90 : /// @}
91 :
92 : /*!
93 : * \brief The inverse of the transition function
94 : *
95 : * This method returns $r/\tilde{r}$ given the mapped coordinates
96 : * $\tilde{x}^i$ (`target_coords`) and the spherical harmonic expansion
97 : * $\Sigma(t, \theta, \phi) = \sum_{lm} \lambda_{lm}(t)Y_{lm}(\theta, \phi)$
98 : * (`radial_distortion`). See domain::CoordinateMaps::TimeDependent::Shape for
99 : * details on how this quantity is used to compute the inverse of the Shape
100 : * map.
101 : *
102 : * To derive the expression for this inverse, solve Eq.
103 : * ($\ref{eq:shape_map_radius}$) for $r$ after substituting
104 : * $G(r,\theta,\phi)$.
105 : *
106 : * \param target_coords The mapped Cartesian coordinates $\tilde{x}^i$.
107 : * \param radial_distortion The spherical harmonic expansion
108 : * $\Sigma(t, \theta, \phi)$.
109 : * \return The quantity $r/\tilde{r}$.
110 : */
111 1 : virtual std::optional<double> original_radius_over_radius(
112 : const std::array<double, 3>& target_coords,
113 : double radial_distortion) const = 0;
114 :
115 : /*!
116 : * Evaluate the gradient of the transition function with respect to the
117 : * Cartesian coordinates x, y and z at the Cartesian coordinates
118 : * `source_coords`.
119 : */
120 : /// @{
121 1 : virtual std::array<double, 3> gradient(
122 : const std::array<double, 3>& source_coords) const = 0;
123 1 : virtual std::array<DataVector, 3> gradient(
124 : const std::array<DataVector, 3>& source_coords) const = 0;
125 : /// @}
126 :
127 0 : virtual std::unique_ptr<ShapeMapTransitionFunction> get_clone() const = 0;
128 :
129 0 : virtual bool operator==(const ShapeMapTransitionFunction& other) const = 0;
130 0 : virtual bool operator!=(const ShapeMapTransitionFunction& other) const = 0;
131 :
132 0 : WRAPPED_PUPable_abstract(ShapeMapTransitionFunction);
133 0 : explicit ShapeMapTransitionFunction(CkMigrateMessage* m) : PUP::able(m) {}
134 : };
135 : } // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions
|
