Wedge.hpp
1 // 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 <optional>
10 
12 #include "Domain/CoordinateMaps/Distribution.hpp"
13 #include "Domain/Structure/OrientationMap.hpp"
14 #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
15 
16 /// \cond
17 namespace PUP {
18 class er;
19 } // namespace PUP
20 /// \endcond
21 
22 namespace domain::CoordinateMaps {
23 
24 namespace detail {
25 // This mapping can be deleted once the 2D and 3D wedges are oriented the same
26 // (see issue https://github.com/sxs-collaboration/spectre/issues/2988)
27 template <size_t Dim>
28 struct WedgeCoordOrientation;
29 template <>
30 struct WedgeCoordOrientation<2> {
31  static constexpr size_t radial_coord = 0;
32  static constexpr size_t polar_coord = 1;
33  static constexpr size_t azimuth_coord = 2; // unused
34 };
35 template <>
36 struct WedgeCoordOrientation<3> {
37  static constexpr size_t radial_coord = 2;
38  static constexpr size_t polar_coord = 0;
39  static constexpr size_t azimuth_coord = 1;
40 };
41 } // namespace detail
42 
43 /*!
44  * \ingroup CoordinateMapsGroup
45  *
46  * \brief Map from a square or cube to a wedge.
47  * \image html Shell.png "A shell can be constructed out of six wedges."
48  *
49  * \details The mapping that goes from a reference cube (in 3D) or square (in
50  * 2D) to a wedge centered on a coordinate axis covering a volume between an
51  * inner surface and outer surface. Each surface can be given a curvature
52  * between flat (a sphericity of 0) or spherical (a sphericity of 1).
53  *
54  * In 2D, the first logical coordinate corresponds to the radial coordinate,
55  * and the second logical coordinates correspond to the angular coordinate. In
56  * 3D, the first two logical coordinates correspond to the two angular
57  * coordinates, and the third to the radial coordinate. This difference
58  * originates from separate implementations for the 2D and 3D map that were
59  * merged. The 3D implementation can be changed to use the first logical
60  * coordinate as radial direction, but this requires propagating the change
61  * through the rest of the domain code (see issue
62  * https://github.com/sxs-collaboration/spectre/issues/2988).
63  *
64  * The following documentation is for the 3D map. The 2D map is obtained by
65  * setting either of the two angular coordinates to zero (and using \f$\xi\f$
66  * as radial coordinate).
67  *
68  * The Wedge map is constructed by linearly interpolating between a bulged
69  * face of radius `radius_of_inner_surface` to a bulged face of
70  * radius `radius_of_outer_surface`, where the radius of each bulged face
71  * is defined to be the radius of the sphere circumscribing the bulge.
72  *
73  * We make a choice here as to whether we wish to use the logical coordinates
74  * parameterizing these surface as they are, in which case we have the
75  * equidistant choice of coordinates, or whether to apply a tangent map to them
76  * which leads us to the equiangular choice of coordinates. In terms of the
77  * logical coordinates, the equiangular coordinates are:
78  *
79  * \f[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\f]
80  *
81  * \f[\textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
82  *
83  * With derivatives:
84  *
85  * \f[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\f]
86  *
87  * \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
88  *
89  * The equidistant coordinates are:
90  *
91  * \f[ \textrm{equidistant xi} : \Xi = \xi\f]
92  *
93  * \f[ \textrm{equidistant eta} : \mathrm{H} = \eta\f]
94  *
95  * with derivatives:
96  *
97  * <center>\f$\Xi'(\xi) = 1\f$, and \f$\mathrm{H}'(\eta) = 1\f$</center>
98  *
99  * We also define the variable \f$\rho\f$, given by:
100  *
101  * \f[\textrm{rho} : \rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\f]
102  *
103  * ### The Spherical Face Map
104  * The surface map for the spherical face of radius \f$R\f$ lying in the
105  * \f$+z\f$
106  * direction in either choice of coordinates is then given by:
107  *
108  * \f[\vec{\sigma}_{spherical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
109  * Where
110  * \f[
111  * \vec{x}(\xi,\eta) =
112  * \begin{bmatrix}
113  * x(\xi,\eta)\\
114  * y(\xi,\eta)\\
115  * z(\xi,\eta)\\
116  * \end{bmatrix} = \frac{R}{\rho}
117  * \begin{bmatrix}
118  * \Xi\\
119  * \mathrm{H}\\
120  * 1\\
121  * \end{bmatrix}\f]
122  *
123  * ### The Bulged Face Map
124  * The bulged surface is itself constructed by linearly interpolating between
125  * a cubical face and a spherical face. The surface map for the cubical face
126  * of side length \f$2L\f$ lying in the \f$+z\f$ direction is given by:
127  *
128  * \f[\vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
129  * Where
130  * \f[
131  * \vec{x}(\xi,\eta) =
132  * \begin{bmatrix}
133  * x(\xi,\eta)\\
134  * y(\xi,\eta)\\
135  * L\\
136  * \end{bmatrix} = L
137  * \begin{bmatrix}
138  * \Xi\\
139  * \mathrm{H}\\
140  * 1\\
141  * \end{bmatrix}\f]
142  *
143  * To construct the bulged map we interpolate between this cubical face map
144  * and a spherical face map of radius \f$R\f$, with the
145  * interpolation parameter being \f$s\f$. The surface map for the bulged face
146  * lying in the \f$+z\f$ direction is then given by:
147  *
148  * \f[\vec{\sigma}_{bulged}(\xi,\eta) = {(1-s)L + \frac{sR}{\rho}}
149  * \begin{bmatrix}
150  * \Xi\\
151  * \mathrm{H}\\
152  * 1\\
153  * \end{bmatrix}\f]
154  *
155  * We constrain L by demanding that the spherical face circumscribe the cube.
156  * With this condition, we have \f$L = R/\sqrt3\f$.
157  * \note This differs from the choice in SpEC where it is demanded that the
158  * surfaces touch at the center, which leads to \f$L = R\f$.
159  *
160  * ### The Full Volume Map
161  * The final map for the wedge which lies along the \f$+z\f$ is obtained by
162  * interpolating between the two surfaces with the
163  * interpolation parameter being the logical coordinate \f$\zeta\f$. This
164  * results in:
165  *
166  * \f[\vec{x}(\xi,\eta,\zeta) =
167  * \frac{1}{2}\left\{(1-\zeta)\Big[(1-s_{inner})\frac{R_{inner}}{\sqrt 3}
168  * + s_{inner}\frac{R_{inner}}{\rho}\Big] +
169  * (1+\zeta)\Big[(1-s_{outer})\frac{R_{outer}}{\sqrt 3} +s_{outer}
170  * \frac{R_{outer}}{\rho}\Big] \right\}\begin{bmatrix}
171  * \Xi\\
172  * \mathrm{H}\\
173  * 1\\
174  * \end{bmatrix}\f]
175  *
176  * We will define the variables \f$F(\zeta)\f$ and \f$S(\zeta)\f$, the frustum
177  * and sphere factors: \f[F(\zeta) = F_0 + F_1\zeta\f] \f[S(\zeta) = S_0 +
178  * S_1\zeta\f]
179  * Where \f{align*}F_0 &= \frac{1}{2} \big\{ (1-s_{outer})R_{outer} +
180  * (1-s_{inner})R_{inner}\big\}\\
181  * F_1 &= \partial_{\zeta} F = \frac{1}{2} \big\{ (1-s_{outer})R_{outer} -
182  * (1-s_{inner})R_{inner}\big\}\\
183  * S_0 &= \frac{1}{2} \big\{ s_{outer}R_{outer} + s_{inner}R_{inner}\big\}\\
184  * S_1 &= \partial_{\zeta} S = \frac{1}{2} \big\{ s_{outer}R_{outer} -
185  * s_{inner}R_{inner}\big\}\f}
186  *
187  * The map can then be rewritten as:
188  * \f[\vec{x}(\xi,\eta,\zeta) = \left\{\frac{F(\zeta)}{\sqrt 3} +
189  * \frac{S(\zeta)}{\rho}\right\}\begin{bmatrix}
190  * \Xi\\
191  * \mathrm{H}\\
192  * 1\\
193  * \end{bmatrix}\f]
194  *
195  * We provide some common derivatives:
196  * \f[\partial_{\xi}z = \frac{-S(\zeta)\Xi\Xi'}{\rho^3}\f]
197  * \f[\partial_{\eta}z = \frac{-S(\zeta)\mathrm{H}\mathrm{H}'}{\rho^3}\f]
198  * \f[\partial_{\zeta}z = \frac{F'}{\sqrt 3} + \frac{S'}{\rho}\f]
199  * The Jacobian then is: \f[J =
200  * \begin{bmatrix}
201  * \Xi'z + \Xi\partial_{\xi}z & \Xi\partial_{\eta}z & \Xi\partial_{\zeta}z \\
202  * \mathrm{H}\partial_{\xi}z & \mathrm{H}'z +
203  * \mathrm{H}\partial_{\eta}z & \mathrm{H}\partial_{\zeta}z\\
204  * \partial_{\xi}z&\partial_{\eta}z &\partial_{\zeta}z \\
205  * \end{bmatrix}
206  * \f]
207  *
208  * A common factor that shows up in the inverse jacobian is:
209  * \f[ T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3}\f]
210  *
211  * The inverse Jacobian then is: \f[J^{-1} =
212  * \frac{1}{z}\begin{bmatrix}
213  * \Xi'^{-1} & 0 & -\Xi\Xi'^{-1}\\
214  * 0 & \mathrm{H}'^{-1} & -\mathrm{H}\mathrm{H}'^{-1}\\
215  * T\Xi &
216  * T\mathrm{H} &
217  * T + F(\partial_{\zeta}z)^{-1}/\sqrt 3\\
218  * \end{bmatrix}
219  * \f]
220  *
221  * ### Changing the radial distribution of the gridpoints
222  * By default, Wedge linearly distributes its gridpoints in the radial
223  * direction. An exponential distribution of gridpoints can be obtained by
224  * linearly interpolating in the logarithm of the radius, in order to obtain
225  * a relatively higher resolution at smaller radii. Since this is a radial
226  * rescaling of Wedge, this option is only supported for fully spherical
227  * wedges with `sphericity_inner` = `sphericity_outer` = 1.
228  *
229  * The linear interpolation done is:
230  * \f[
231  * \log r = \frac{1-\zeta}{2}\log R_{inner} +
232  * \frac{1+\zeta}{2}\log R_{outer}
233  * \f]
234  *
235  * The map then is:
236  * \f[\vec{x}(\xi,\eta,\zeta) =
237  * \frac{\sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}}}{\rho}\begin{bmatrix}
238  * \Xi\\
239  * \mathrm{H}\\
240  * 1\\
241  * \end{bmatrix}\f]
242  *
243  * The jacobian simplifies similarly.
244  *
245  * Alternatively, an inverse radial distribution can be chosen where the linear
246  * interpolation is:
247  *
248  * \f[
249  * \frac{1}{r} = \frac{R_\mathrm{inner} + R_\mathrm{outer}}{2 R_\mathrm{inner}
250  * R_\mathrm{outer}} + \frac{R_\mathrm{inner} - R_\mathrm{outer}}{2
251  * R_\mathrm{inner} R_\mathrm{outer}} \zeta
252  * \f]
253  */
254 template <size_t Dim>
255 class Wedge {
256  public:
257  static constexpr size_t dim = Dim;
258  enum class WedgeHalves {
259  /// Use the entire wedge
260  Both,
261  /// Use only the upper logical half
262  UpperOnly,
263  /// Use only the lower logical half
264  LowerOnly
265  };
266 
267  /*!
268  * Constructs a 3D wedge.
269  * \param radius_inner Distance from the origin to one of the
270  * corners which lie on the inner surface.
271  * \param radius_outer Distance from the origin to one of the
272  * corners which lie on the outer surface.
273  * \param orientation_of_wedge The orientation of the desired wedge relative
274  * to the orientation of the default wedge which is a wedge that has its
275  * curved surfaces pierced by the upper-z axis. The logical xi and eta
276  * coordinates point in the cartesian x and y directions, respectively.
277  * \param sphericity_inner Value between 0 and 1 which determines
278  * whether the inner surface is flat (value of 0), spherical (value of 1) or
279  * somewhere in between
280  * \param sphericity_outer Value between 0 and 1 which determines
281  * whether the outer surface is flat (value of 0), spherical (value of 1) or
282  * somewhere in between
283  * \param with_equiangular_map Determines whether to apply a tangent function
284  * mapping to the logical coordinates (for `true`) or not (for `false`).
285  * \param halves_to_use Determines whether to construct a full wedge or only
286  * half a wedge. If constructing only half a wedge, the resulting shape has a
287  * face normal to the x direction (assuming default OrientationMap). If
288  * constructing half a wedge, an intermediate affine map is applied to the
289  * logical xi coordinate such that the interval [-1,1] is mapped to the
290  * corresponding logical half of the wedge. For example, if `UpperOnly` is
291  * specified, [-1,1] is mapped to [0,1], and if `LowerOnly` is specified,
292  * [-1,1] is mapped to [-1,0]. The case of `Both` means a full wedge, with no
293  * intermediate map applied. In all cases, the logical points returned by the
294  * inverse map will lie in the range [-1,1] in each dimension. Half wedges are
295  * currently only useful in constructing domains for binary systems.
296  * \param radial_distribution Determines how to distribute grid points along
297  * the radial direction. For wedges that are not exactly spherical, only
298  * `Distribution::Linear` is currently supported.
299  */
300  Wedge(double radius_inner, double radius_outer, double sphericity_inner,
301  double sphericity_outer, OrientationMap<Dim> orientation_of_wedge,
302  bool with_equiangular_map,
303  WedgeHalves halves_to_use = WedgeHalves::Both,
304  Distribution radial_distribution = Distribution::Linear) noexcept;
305 
306  Wedge() = default;
307  ~Wedge() = default;
308  Wedge(Wedge&&) = default;
309  Wedge(const Wedge&) = default;
310  Wedge& operator=(const Wedge&) = default;
311  Wedge& operator=(Wedge&&) = default;
312 
313  template <typename T>
314  std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
315  const std::array<T, Dim>& source_coords) const noexcept;
316 
317  /// For a \f$+z\f$-oriented `Wedge`, returns invalid if \f$z<=0\f$
318  /// or if \f$(x,y,z)\f$ is on or outside the cone defined
319  /// by \f$(x^2/z^2 + y^2/z^2+1)^{1/2} = -S/F\f$, where
320  /// \f$S = \frac{1}{2}(s_1 r_1 - s_0 r_0)\f$ and
321  /// \f$F = \frac{1}{2\sqrt{3}}((1-s_1) r_1 - (1-s_0) r_0)\f$.
322  /// Here \f$s_0,s_1\f$ and \f$r_0,r_1\f$ are the specified sphericities
323  /// and radii of the inner and outer \f$z\f$ surfaces. The map is singular on
324  /// the cone and on the xy plane.
325  /// The inverse function is only callable with doubles because the inverse
326  /// might fail if called for a point out of range, and it is unclear
327  /// what should happen if the inverse were to succeed for some points in a
328  /// DataVector but fail for other points.
330  const std::array<double, Dim>& target_coords) const noexcept;
331 
332  template <typename T>
333  tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
334  const std::array<T, Dim>& source_coords) const noexcept;
335 
336  template <typename T>
337  tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
338  const std::array<T, Dim>& source_coords) const noexcept;
339 
340  // clang-tidy: google runtime references
341  void pup(PUP::er& p) noexcept; // NOLINT
342 
343  static constexpr bool is_identity() noexcept { return false; }
344 
345  private:
346  // maps between 2D and 3D choices for coordinate axis orientations
347  static constexpr size_t radial_coord =
348  detail::WedgeCoordOrientation<Dim>::radial_coord;
349  static constexpr size_t polar_coord =
350  detail::WedgeCoordOrientation<Dim>::polar_coord;
351  static constexpr size_t azimuth_coord =
352  detail::WedgeCoordOrientation<Dim>::azimuth_coord;
353 
354  // factors out calculation of z needed for mapping and jacobian
355  template <typename T>
356  tt::remove_cvref_wrap_t<T> default_physical_z(
357  const T& zeta, const T& one_over_rho) const noexcept;
358 
359  template <size_t LocalDim>
360  // NOLINTNEXTLINE(readability-redundant-declaration)
361  friend bool operator==(const Wedge<LocalDim>& lhs,
362  const Wedge<LocalDim>& rhs) noexcept;
363 
364  double radius_inner_{std::numeric_limits<double>::signaling_NaN()};
365  double radius_outer_{std::numeric_limits<double>::signaling_NaN()};
366  double sphericity_inner_{std::numeric_limits<double>::signaling_NaN()};
367  double sphericity_outer_{std::numeric_limits<double>::signaling_NaN()};
368  OrientationMap<Dim> orientation_of_wedge_{};
369  bool with_equiangular_map_ = false;
370  WedgeHalves halves_to_use_ = WedgeHalves::Both;
371  Distribution radial_distribution_ = Distribution::Linear;
372  double scaled_frustum_zero_{std::numeric_limits<double>::signaling_NaN()};
373  double sphere_zero_{std::numeric_limits<double>::signaling_NaN()};
374  double scaled_frustum_rate_{std::numeric_limits<double>::signaling_NaN()};
375  double sphere_rate_{std::numeric_limits<double>::signaling_NaN()};
376 };
377 
378 template <size_t Dim>
379 bool operator!=(const Wedge<Dim>& lhs, const Wedge<Dim>& rhs) noexcept;
380 } // namespace domain::CoordinateMaps
std::rel_ops::operator!=
T operator!=(T... args)
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
OrientationMap< Dim >
domain::CoordinateMaps::Wedge::inverse
std::optional< std::array< double, Dim > > inverse(const std::array< double, Dim > &target_coords) const noexcept
For a -oriented Wedge, returns invalid if or if is on or outside the cone defined by ,...
domain::CoordinateMaps::Distribution
Distribution
Distribution of grid points in one dimension.
Definition: Distribution.hpp:25
domain::CoordinateMaps::Wedge::WedgeHalves
WedgeHalves
Definition: Wedge.hpp:258
cstddef
domain::CoordinateMaps::Wedge::WedgeHalves::Both
@ Both
Use the entire wedge.
array
domain::CoordinateMaps::Wedge::WedgeHalves::UpperOnly
@ UpperOnly
Use only the upper logical half.
domain::CoordinateMaps
Contains all coordinate maps.
Definition: Affine.hpp:23
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
limits
TypeAliases.hpp
optional
domain::CoordinateMaps::Wedge::WedgeHalves::LowerOnly
@ LowerOnly
Use only the lower logical half.
domain::CoordinateMaps::Wedge
Map from a square or cube to a wedge.
Definition: Wedge.hpp:255