SpECTRE Documentation Coverage Report
Current view: top level - Domain/CoordinateMaps - Wedge.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 3 37 8.1 %
Date: 2024-03-29 00:33:31
Legend: Lines: hit not hit

          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 <optional>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      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           1 : class Wedge {
     256             :  public:
     257           0 :   static constexpr size_t dim = Dim;
     258           0 :   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             :    * \param opening_angles Determines the angular size of the wedge. The
     300             :    * default value is pi/2, which corresponds to a wedge size of pi/2. For this
     301             :    * setting, four Wedges can be put together to cover 2pi in angle along a
     302             :    * great circle. This option is meant to be used with the equiangular map
     303             :    * option turned on.
     304             :    * \param with_adapted_equiangular_map Determines whether to adapt the
     305             :    * point distribution in the wedge to match its physical angular size. When
     306             :    * `true`, angular distances are proportional to logical distances. Note
     307             :    * that it is not possible to use adapted maps in every Wedge of a Sphere
     308             :    * unless each Wedge has the same size along both angular directions.
     309             :    */
     310           1 :   Wedge(double radius_inner, double radius_outer, double sphericity_inner,
     311             :         double sphericity_outer, OrientationMap<Dim> orientation_of_wedge,
     312             :         bool with_equiangular_map,
     313             :         WedgeHalves halves_to_use = WedgeHalves::Both,
     314             :         Distribution radial_distribution = Distribution::Linear,
     315             :         const std::array<double, Dim - 1>& opening_angles =
     316             :             make_array<Dim - 1>(M_PI_2),
     317             :         bool with_adapted_equiangular_map = true);
     318             : 
     319           0 :   Wedge() = default;
     320           0 :   ~Wedge() = default;
     321           0 :   Wedge(Wedge&&) = default;
     322           0 :   Wedge(const Wedge&) = default;
     323           0 :   Wedge& operator=(const Wedge&) = default;
     324           0 :   Wedge& operator=(Wedge&&) = default;
     325             : 
     326             :   template <typename T>
     327           0 :   std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
     328             :       const std::array<T, Dim>& source_coords) const;
     329             : 
     330             :   /// For a \f$+z\f$-oriented `Wedge`, returns invalid if \f$z<=0\f$
     331             :   /// or if \f$(x,y,z)\f$ is on or outside the cone defined
     332             :   /// by \f$(x^2/z^2 + y^2/z^2+1)^{1/2} = -S/F\f$, where
     333             :   /// \f$S = \frac{1}{2}(s_1 r_1 - s_0 r_0)\f$ and
     334             :   /// \f$F = \frac{1}{2\sqrt{3}}((1-s_1) r_1 - (1-s_0) r_0)\f$.
     335             :   /// Here \f$s_0,s_1\f$ and \f$r_0,r_1\f$ are the specified sphericities
     336             :   /// and radii of the inner and outer \f$z\f$ surfaces.  The map is singular on
     337             :   /// the cone and on the xy plane.
     338             :   /// The inverse function is only callable with doubles because the inverse
     339             :   /// might fail if called for a point out of range, and it is unclear
     340             :   /// what should happen if the inverse were to succeed for some points in a
     341             :   /// DataVector but fail for other points.
     342           1 :   std::optional<std::array<double, Dim>> inverse(
     343             :       const std::array<double, Dim>& target_coords) const;
     344             : 
     345             :   template <typename T>
     346           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
     347             :       const std::array<T, Dim>& source_coords) const;
     348             : 
     349             :   template <typename T>
     350           0 :   tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
     351             :       const std::array<T, Dim>& source_coords) const;
     352             : 
     353             :   // NOLINTNEXTLINE(google-runtime-references)
     354           0 :   void pup(PUP::er& p);
     355             : 
     356           0 :   static constexpr bool is_identity() { return false; }
     357             : 
     358             :  private:
     359             :   // maps between 2D and 3D choices for coordinate axis orientations
     360           0 :   static constexpr size_t radial_coord =
     361             :       detail::WedgeCoordOrientation<Dim>::radial_coord;
     362           0 :   static constexpr size_t polar_coord =
     363             :       detail::WedgeCoordOrientation<Dim>::polar_coord;
     364           0 :   static constexpr size_t azimuth_coord =
     365             :       detail::WedgeCoordOrientation<Dim>::azimuth_coord;
     366             : 
     367             :   // factors out calculation of z needed for mapping and jacobian
     368             :   template <typename T>
     369           0 :   tt::remove_cvref_wrap_t<T> default_physical_z(const T& zeta,
     370             :                                                 const T& one_over_rho) const;
     371             : 
     372             :   template <size_t LocalDim>
     373             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     374           0 :   friend bool operator==(const Wedge<LocalDim>& lhs,
     375             :                          const Wedge<LocalDim>& rhs);
     376             : 
     377           0 :   double radius_inner_{std::numeric_limits<double>::signaling_NaN()};
     378           0 :   double radius_outer_{std::numeric_limits<double>::signaling_NaN()};
     379           0 :   double sphericity_inner_{std::numeric_limits<double>::signaling_NaN()};
     380           0 :   double sphericity_outer_{std::numeric_limits<double>::signaling_NaN()};
     381           0 :   OrientationMap<Dim> orientation_of_wedge_{};
     382           0 :   bool with_equiangular_map_ = false;
     383           0 :   WedgeHalves halves_to_use_ = WedgeHalves::Both;
     384           0 :   Distribution radial_distribution_ = Distribution::Linear;
     385           0 :   double scaled_frustum_zero_{std::numeric_limits<double>::signaling_NaN()};
     386           0 :   double sphere_zero_{std::numeric_limits<double>::signaling_NaN()};
     387           0 :   double scaled_frustum_rate_{std::numeric_limits<double>::signaling_NaN()};
     388           0 :   double sphere_rate_{std::numeric_limits<double>::signaling_NaN()};
     389           0 :   std::array<double, Dim - 1> opening_angles_{
     390             :       make_array<Dim - 1>(std::numeric_limits<double>::signaling_NaN())};
     391           0 :   std::array<double, Dim - 1> opening_angles_distribution_{
     392             :       make_array<Dim - 1>(std::numeric_limits<double>::signaling_NaN())};
     393             : };
     394             : 
     395             : template <size_t Dim>
     396           0 : bool operator!=(const Wedge<Dim>& lhs, const Wedge<Dim>& rhs);
     397             : }  // namespace domain::CoordinateMaps

Generated by: LCOV version 1.14