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 <cstddef>
8 : #include <limits>
9 : #include <memory>
10 : #include <optional>
11 : #include <pup.h>
12 :
13 : #include "DataStructures/DataVector.hpp"
14 : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
15 :
16 : namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
17 : /*!
18 : * \brief A transition function that falls off linearly from an inner surface of
19 : * a wedge to an outer surface of a wedge. Meant to be used in
20 : * `domain::CoordinateMaps::Wedge` blocks.
21 : *
22 : * \details The functional form of this transition is
23 : *
24 : * \begin{equation}
25 : * f(r, \theta, \phi) = \frac{D_{\text{out}}(r, \theta, \phi) -
26 : * r}{D_{\text{out}}(r, \theta, \phi) - D_{\text{in}}(r, \theta, \phi)},
27 : * \label{eq:transition_func}
28 : * \end{equation}
29 : *
30 : * where
31 : *
32 : * \begin{equation}
33 : * D(r, \theta, \phi) = R\left(\frac{1 -
34 : * s}{\sqrt{3}\max(|\sin{\theta}\cos{\phi}|,|\sin{\theta}\sin{\phi}|,
35 : * |\cos{\theta}|)} + s\right),
36 : * \label{eq:distance}
37 : * \end{equation}
38 : *
39 : * where $s$ is the sphericity of the surface which goes from 0 (flat) to 1
40 : * (spherical), $R$ is the radius of the spherical surface, $\text{out}$ is the
41 : * outer surface, and $\text{in}$ is the inner surface. If the sphericity is 1,
42 : * then the surface is a sphere so $D = R$. If the sphericity is 0, then the
43 : * surface is a cube. This cube is circumscribed by a sphere of radius $R$. See
44 : * `domain::CoordinateMaps::Wedge` for more of an explanation of these boundary
45 : * surfaces and their sphericities.
46 : *
47 : * \note Because the shape map distorts only radii and does not affect angles,
48 : * $D$ is not a function of $r$ so we have that $D(r,\theta,\phi) =
49 : * D(\theta, \phi)$.
50 : *
51 : * There are several assumptions made for this mapping:
52 : *
53 : * - The `domain::CoordinateMaps::Wedge`s have an opening angle of
54 : * $\frac{\pi}{2}$ in each direction.
55 : * - The coordinates $r, \theta, \phi$ are assumed to be from the center of the
56 : * wedge, not the center of the computational domain.
57 : * - The wedges are concentric. (see the constructor)
58 : * - The $\max$ in the denominator of $\ref{eq:distance}$ can be simplified a
59 : * bit to $\max(|x|, |y|, |z|)/r$. It was written the other way in
60 : * $\ref{eq:distance}$ to emphasize that $D$ has no radial dependence.
61 : *
62 : * ## Gradient
63 : *
64 : * The cartesian gradient of the transition function is
65 : *
66 : * \begin{equation}
67 : * \frac{\partial f}{\partial x_i} = \frac{\frac{\partial
68 : * D_{\text{out}}}{\partial x_i} - \frac{x_i}{r}}{D_{\text{out}} -
69 : * D_{\text{in}}} - \frac{\left(D_{\text{out}} - r\right)\left(\frac{\partial
70 : * D_{\text{out}}}{\partial x_i} - \frac{\partial
71 : * D_{\text{in}}}{\partial x_i} \right)}{\left(D_{\text{out}} -
72 : * D_{\text{in}}\right)^2}.
73 : * \end{equation}
74 : *
75 : * The computation of the gradient of $D$ depends on the result of the $\max$ in
76 : * the denominator of $\ref{eq:distance}$. To simplify the expression, let $j
77 : * \in \{0,1,2\}$ correspond to the $\max(|x|, |y|, |z|)$ such that $j=0$ if
78 : * $|x|$ is the $\max$, and so on. Then we can write the gradient of $D$ as
79 : *
80 : * \f{align}{
81 : * \frac{\partial D}{\partial x_j} &= -\text{sgn}(x_j)\frac{R(1-s)\left(r^2 -
82 : * x_j^2\right)}{r x_j^2 \sqrt{3}} \\
83 : * \frac{\partial D}{\partial x_{j+1}} &= \frac{R(1-s)x_{j+1}}{r |x_j| \sqrt{3}}
84 : * \\
85 : * \frac{\partial D}{\partial x_{j+2}} &= \frac{R(1-s)x_{j+2}}{r |x_j| \sqrt{3}}
86 : * \f}
87 : *
88 : * where $j+1$ and $j+2$ are understood to be $\mod 3$. Also since we don't
89 : * allow points at the origin of these wedges, we can be assured that for
90 : * whichever $j$ is max, that $x_j$ won't be zero.
91 : *
92 : * ## Original radius divided by mapped radius
93 : *
94 : * Given an already mapped point and the distorted radius $\Sigma(\theta, \phi)
95 : * = \sum_{\ell,m}\lambda_{\ell m}Y_{\ell m}(\theta, \phi)$, we can figure out
96 : * the ratio of the original radius $r$ to the mapped $\tilde{r}$ by solving
97 : *
98 : * \begin{equation}
99 : * \frac{r}{\tilde{r}} =
100 : * \frac{1}{1-\frac{f(r,\theta,\phi)}{r}\Sigma(\theta,\phi)}
101 : * \end{equation}
102 : *
103 : * After plugging in the transition and solving, we get
104 : *
105 : * \begin{equation}
106 : * \frac{r}{\tilde{r}} = \frac{1 + \frac{D_{\text{out}}\Sigma(\theta,
107 : * \phi)}{\tilde{r}(D_{\text{out}} - D_{\text{in}})}}{1 + \frac{\Sigma(\theta,
108 : * \phi)}{D_{\text{out}} - D_{\text{in}}}}
109 : * \label{eq:r_over_rtil}
110 : * \end{equation}
111 : *
112 : * \note Since $D$ is not a function of the radius, we can compute the angles
113 : * using $\tilde{r}$ in Eq. $\ref{eq:r_over_rtil}$.
114 : *
115 : * ## Special cases
116 : *
117 : * The equations above become simpler if the inner boundary, outer boundary, or
118 : * both are spherical ($s = 1$). If a boundary is spherical, then $D = R$ and
119 : * $\frac{\partial D}{\partial x_i} = 0$. This simplifies the expanded form of
120 : * the gradient significantly.
121 : */
122 1 : class Wedge final : public ShapeMapTransitionFunction {
123 0 : struct Surface {
124 0 : double radius{};
125 0 : double sphericity{};
126 :
127 : // This is the distance from the center (assumed to be 0,0,0) to this
128 : // surface in the same direction as coords. The calculation is cheaper if
129 : // you know the axis ahead of time
130 : template <typename T>
131 0 : T distance(const std::array<T, 3>& coords,
132 : const std::optional<size_t>& axis = std::nullopt) const;
133 :
134 0 : void pup(PUP::er& p);
135 :
136 0 : bool operator==(const Surface& other) const;
137 0 : bool operator!=(const Surface& other) const;
138 : };
139 :
140 : public:
141 0 : explicit Wedge() = default;
142 :
143 : /*!
144 : * \brief Construct a Wedge transition for a wedge block in a given direction.
145 : *
146 : * \details Many concentric wedges can be part of the same falloff from 1 at
147 : * the inner boundary of the innermost wedge, to 0 at the outer boundary of
148 : * the outermost wedge.
149 : *
150 : * \param inner_radius Inner radius of innermost wedge
151 : * \param outer_radius Outermost radius of outermost wedge
152 : * \param inner_sphericity Sphericity of innermost surface of innermost wedge
153 : * \param outer_sphericity Sphericity of outermost surface of outermost wedge
154 : * \param axis The direction that this wedge is in. Both the positive and
155 : * negative direction get the same axis.
156 : */
157 1 : Wedge(double inner_radius, double outer_radius, double inner_sphericity,
158 : double outer_sphericity, size_t axis);
159 :
160 1 : double operator()(const std::array<double, 3>& source_coords) const override;
161 1 : DataVector operator()(
162 : const std::array<DataVector, 3>& source_coords) const override;
163 :
164 1 : std::optional<double> original_radius_over_radius(
165 : const std::array<double, 3>& target_coords,
166 : double distorted_radius) const override;
167 :
168 1 : std::array<double, 3> gradient(
169 : const std::array<double, 3>& source_coords) const override;
170 1 : std::array<DataVector, 3> gradient(
171 : const std::array<DataVector, 3>& source_coords) const override;
172 :
173 0 : WRAPPED_PUPable_decl_template(Wedge);
174 0 : explicit Wedge(CkMigrateMessage* msg);
175 0 : void pup(PUP::er& p) override;
176 :
177 1 : std::unique_ptr<ShapeMapTransitionFunction> get_clone() const override {
178 : return std::make_unique<Wedge>(*this);
179 : }
180 :
181 0 : bool operator==(const ShapeMapTransitionFunction& other) const override;
182 0 : bool operator!=(const ShapeMapTransitionFunction& other) const override;
183 :
184 : private:
185 : template <typename T>
186 0 : T call_impl(const std::array<T, 3>& source_coords) const;
187 :
188 : template <typename T>
189 0 : std::array<T, 3> gradient_impl(const std::array<T, 3>& source_coords) const;
190 :
191 : template <typename T>
192 0 : void check_distances(const std::array<T, 3>& coords) const;
193 :
194 0 : Surface inner_surface_{};
195 0 : Surface outer_surface_{};
196 0 : size_t axis_{};
197 0 : static constexpr double eps_ = std::numeric_limits<double>::epsilon() * 100;
198 : };
199 : } // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions
|
