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 <cmath>
       8             : #include <cstddef>
       9             : #include <limits>
      10             : #include <optional>
      11             : 
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "DataStructures/Tensor/TypeAliases.hpp"
      14             : #include "Options/Context.hpp"
      15             : #include "Options/String.hpp"
      16             : 
      17             : /// \cond
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : namespace gsl {
      22             : template <typename T>
      23             : class not_null;
      24             : }  // namespace gsl
      25             : /// \endcond
      26             : 
      27             : namespace grmhd::AnalyticData {
      28             : /*!
      29             :  * \brief Torus made by removing two polar cones from a spherical shell
      30             :  *
      31             :  * Maps source coordinates \f$(\xi, \eta, \zeta)\f$ to
      32             :  * \f{align}
      33             :  * \vec{x}(\xi, \eta, \zeta) =
      34             :  * \begin{bmatrix}
      35             :  *  r \sin\theta\cos\phi \\
      36             :  *  r \sin\theta\sin\phi \\
      37             :  *  r \cos\theta
      38             :  * \end{bmatrix}
      39             :  * \f}
      40             :  *
      41             :  * where
      42             :  * \f{align}
      43             :  *  r & = r_\mathrm{min}\frac{1-\xi}{2} + r_\mathrm{max}\frac{1+\xi}{2}, \\
      44             :  *  \theta & = \pi/2 - (\pi/2 - \theta_\mathrm{min}) \eta, \\
      45             :  *  \phi   & = f_\mathrm{torus} \pi \zeta.
      46             :  * \f}
      47             :  *
      48             :  *  - $r_\mathrm{min}$ and $r_\mathrm{max}$ are inner and outer radius of torus.
      49             :  *  - $\theta_\mathrm{min}\in(0,\pi/2)$ is the minimum polar angle (measured
      50             :  *    from +z axis) of torus, which is equal to the half of the apex angle of
      51             :  *    the removed polar cones.
      52             :  *  - $f_\mathrm{torus}\in(0, 1)$ is azimuthal fraction that the torus covers.
      53             :  *
      54             :  */
      55           1 : class SphericalTorus {
      56             :  public:
      57           0 :   static constexpr size_t dim = 3;
      58             : 
      59           0 :   struct RadialRange {
      60           0 :     using type = std::array<double, 2>;
      61           0 :     static constexpr Options::String help =
      62             :         "Radial extent of the torus, "
      63             :         "[min_radius, max_radius] ";
      64             :   };
      65             : 
      66           0 :   struct MinPolarAngle {
      67           0 :     using type = double;
      68           0 :     static constexpr Options::String help =
      69             :         "Half of the apex angle of excised polar cones. "
      70             :         "Polar angle (measured from +z axis) of torus has range "
      71             :         "[MinPolarAngle, pi - MinPolarAngle]";
      72           0 :     static type lower_bound() { return 0.0; }
      73           0 :     static type upper_bound() { return 0.5 * M_PI; }
      74             :   };
      75             : 
      76           0 :   struct FractionOfTorus {
      77           0 :     using type = double;
      78           0 :     static constexpr Options::String help =
      79             :         "Fraction of (azimuthal) orbit covered. Azimuthal angle has range "
      80             :         "[- pi * FractionOfTorus, pi * FractionOfTorus].";
      81           0 :     static type lower_bound() { return 0.0; }
      82           0 :     static type upper_bound() { return 1.0; }
      83             :   };
      84             : 
      85           0 :   static constexpr Options::String help =
      86             :       "Torus made by removing polar cones from a spherical shell";
      87             : 
      88           0 :   using options = tmpl::list<RadialRange, MinPolarAngle, FractionOfTorus>;
      89             : 
      90           0 :   SphericalTorus(const std::array<double, 2>& radial_range,
      91             :                  const double min_polar_angle, const double fraction_of_torus,
      92             :                  const Options::Context& context = {});
      93             : 
      94           0 :   SphericalTorus(double r_min, double r_max, double min_polar_angle,
      95             :                  double fraction_of_torus = 1.0,
      96             :                  const Options::Context& context = {});
      97             : 
      98           0 :   SphericalTorus() = default;
      99             : 
     100             :   template <typename T>
     101           0 :   tnsr::I<T, 3> operator()(const tnsr::I<T, 3>& source_coords) const;
     102             : 
     103           0 :   tnsr::I<double, 3> inverse(const tnsr::I<double, 3>& target_coords) const;
     104             : 
     105             :   template <typename T>
     106           0 :   Jacobian<T, 3, Frame::BlockLogical, Frame::Inertial> jacobian(
     107             :       const tnsr::I<T, 3>& source_coords) const;
     108             : 
     109             :   template <typename T>
     110           0 :   InverseJacobian<T, 3, Frame::BlockLogical, Frame::Inertial> inv_jacobian(
     111             :       const tnsr::I<T, 3>& source_coords) const;
     112             : 
     113             :   template <typename T>
     114           0 :   tnsr::Ijj<T, 3, Frame::NoFrame> hessian(
     115             :       const tnsr::I<T, 3>& source_coords) const;
     116             : 
     117             :   template <typename T>
     118           0 :   tnsr::Ijk<T, 3, Frame::NoFrame> derivative_of_inv_jacobian(
     119             :       const tnsr::I<T, 3>& source_coords) const;
     120             : 
     121             :   // NOLINTNEXTLINE(google-runtime-references)
     122           0 :   void pup(PUP::er& p);
     123             : 
     124           0 :   bool is_identity() const { return false; }
     125             : 
     126             :  private:
     127             :   template <typename T>
     128           0 :   T radius(const T& x) const;
     129             : 
     130             :   template <typename T>
     131           0 :   void radius(gsl::not_null<T*> r, const T& x) const;
     132             : 
     133             :   template <typename T>
     134           0 :   T radius_inverse(const T& r) const;
     135             : 
     136           0 :   friend bool operator==(const SphericalTorus& lhs, const SphericalTorus& rhs);
     137             : 
     138           0 :   double r_min_ = std::numeric_limits<double>::signaling_NaN();
     139           0 :   double r_max_ = std::numeric_limits<double>::signaling_NaN();
     140           0 :   double pi_over_2_minus_theta_min_ =
     141             :       std::numeric_limits<double>::signaling_NaN();
     142           0 :   double fraction_of_torus_ = std::numeric_limits<double>::signaling_NaN();
     143             : };
     144             : 
     145           0 : bool operator!=(const SphericalTorus& lhs, const SphericalTorus& rhs);
     146             : }  // namespace grmhd::AnalyticData