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 <pup.h>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/DataBox/Prefixes.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "Elliptic/Systems/SelfForce/GeneralRelativity/Tags.hpp"
      15             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      16             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
      19             : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
      20             : #include "Utilities/Gsl.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : #include "Utilities/TaggedTuple.hpp"
      23             : 
      24           0 : namespace GrSelfForce::AnalyticData {
      25             : 
      26             : /*!
      27             :  * \brief Gravitational self-force of a point mass on a circular equatorial
      28             :  * orbit in Kerr.
      29             :  *
      30             :  * This class defines the gravitational self-force equations for a circular
      31             :  * orbit by setting the coefficients $\alpha$, $\beta$, and $\gamma$ (see
      32             :  * `GrSelfForce::FirstOrderSystem`). It also sets the effective source
      33             :  * $S_m^\mathrm{eff}$ and the singular field $\Psi_m^P$ in the regularized
      34             :  * region. The coefficients are computed using Mathematica-generated functions
      35             :  * (see CircularOrbitCoeffs.hpp) and the effective source is computed using the
      36             :  * GravitationalEffectiveSource code by Wardell et. al.
      37             :  * (https://github.com/barrywardell/GravitationalEffectiveSource) and then
      38             :  * transformed to our form of the equations with more Mathematica-generated
      39             :  * functions (see CircularOrbitConvertEffsource.hpp). The derivation of these
      40             :  * equations will be presented in a future publication. A very strong test of
      41             :  * the validity of these equations is evaluating them on the singular field and
      42             :  * the corresponding effective source provided by the external
      43             :  * GravitationalEffectiveSource code (see Test_CircularOrbit.cpp).
      44             :  */
      45           1 : class CircularOrbit : public elliptic::analytic_data::Background,
      46             :                       public elliptic::analytic_data::InitialGuess {
      47             :  public:
      48           0 :   struct BlackHoleMass {
      49           0 :     static constexpr Options::String help =
      50             :         "Kerr mass parameter 'M' of the black hole";
      51           0 :     using type = double;
      52             :   };
      53           0 :   struct BlackHoleSpin {
      54           0 :     static constexpr Options::String help =
      55             :         "Kerr dimensionless spin parameter 'chi' of the black hole";
      56           0 :     using type = double;
      57             :   };
      58           0 :   struct OrbitalRadius {
      59           0 :     static constexpr Options::String help =
      60             :         "Radius 'r_0' of the circular orbit";
      61           0 :     using type = double;
      62             :   };
      63           0 :   struct MModeNumber {
      64           0 :     static constexpr Options::String help =
      65             :         "Mode number 'm' of the m-mode decomposition";
      66           0 :     using type = int;
      67             :   };
      68           0 :   using options =
      69             :       tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius, MModeNumber>;
      70           0 :   static constexpr Options::String help =
      71             :       "Quasicircular orbit of a point mass in Kerr spacetime";
      72             : 
      73           0 :   CircularOrbit() = default;
      74           0 :   CircularOrbit(const CircularOrbit&) = default;
      75           0 :   CircularOrbit& operator=(const CircularOrbit&) = default;
      76           0 :   CircularOrbit(CircularOrbit&&) = default;
      77           0 :   CircularOrbit& operator=(CircularOrbit&&) = default;
      78           0 :   ~CircularOrbit() override = default;
      79             : 
      80           0 :   CircularOrbit(double black_hole_mass, double black_hole_spin,
      81             :                 double orbital_radius, int m_mode_number);
      82             : 
      83           0 :   explicit CircularOrbit(CkMigrateMessage* m);
      84             :   using PUP::able::register_constructor;
      85           0 :   WRAPPED_PUPable_decl_template(CircularOrbit);
      86             : 
      87           0 :   tnsr::I<double, 2> puncture_position() const;
      88           0 :   double black_hole_mass() const { return black_hole_mass_; }
      89           0 :   double black_hole_spin() const { return black_hole_spin_; }
      90           0 :   double orbital_radius() const { return orbital_radius_; }
      91           0 :   int m_mode_number() const { return m_mode_number_; }
      92             : 
      93           0 :   using background_tags =
      94             :       tmpl::list<Tags::Alpha, Tags::Beta, Tags::GammaRstar, Tags::GammaTheta>;
      95           0 :   using source_tags = tmpl::list<
      96             :       ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
      97             :       ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
      98             :       Tags::BoyerLindquistRadius>;
      99             : 
     100             :   // Background
     101           0 :   tuples::tagged_tuple_from_typelist<background_tags> variables(
     102             :       const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
     103             : 
     104             :   // Initial guess
     105           0 :   static tuples::TaggedTuple<Tags::MMode> variables(
     106             :       const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
     107             : 
     108             :   // Fixed sources
     109           0 :   tuples::tagged_tuple_from_typelist<source_tags> variables(
     110             :       const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
     111             : 
     112             :   template <typename... RequestedTags>
     113           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     114             :       const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
     115             :       const InverseJacobian<DataVector, 2, Frame::ElementLogical,
     116             :                             Frame::Inertial>& /*inv_jacobian*/,
     117             :       tmpl::list<RequestedTags...> /*meta*/) const {
     118             :     return variables(x, tmpl::list<RequestedTags...>{});
     119             :   }
     120             : 
     121             :   // NOLINTNEXTLINE
     122           0 :   void pup(PUP::er& p) override;
     123             : 
     124             :  private:
     125           0 :   friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
     126             : 
     127           0 :   double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
     128           0 :   double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
     129           0 :   double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
     130           0 :   int m_mode_number_{};
     131             : };
     132             : 
     133           0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
     134             : 
     135             : }  // namespace GrSelfForce::AnalyticData