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             : #include <pup.h>
      11             : #include <vector>
      12             : 
      13             : #include "DataStructures/DataBox/Prefixes.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "Elliptic/Systems/SelfForce/Scalar/Tags.hpp"
      16             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      17             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      18             : #include "Options/Auto.hpp"
      19             : #include "Options/String.hpp"
      20             : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
      21             : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
      22             : #include "Utilities/Gsl.hpp"
      23             : #include "Utilities/TMPL.hpp"
      24             : #include "Utilities/TaggedTuple.hpp"
      25             : 
      26           0 : namespace ScalarSelfForce::AnalyticData {
      27             : 
      28             : /*!
      29             :  * \brief Scalar self force for a scalar charge on a circular equatorial orbit.
      30             :  *
      31             :  * This class implements Eq. (2.9) in \cite Osburn:2022bby . It does so by
      32             :  * defining the background fields $\alpha$, $\beta$, and $\gamma_i$ in the
      33             :  * general form of the equations
      34             :  * \begin{equation}
      35             :  * -\partial_i F^i + \beta \Psi_m + \gamma_i \partial_i \Psi_m  = S_m
      36             :  * \text{.}
      37             :  * \end{equation}
      38             :  * with the flux
      39             :  * \begin{equation}
      40             :  * F^i = \{\partial_{r_\star}, \alpha \partial_{\cos\theta}\} \Psi_m
      41             :  * \text{.}
      42             :  * \end{equation}
      43             :  * We make the following changes compared to \cite Osburn:2022bby :
      44             :  *
      45             :  * - Multiply by the factor $\Sigma^2 / (r^2 + a^2)^2$ so that we can easily
      46             :  *   write the equations in first-order flux form.
      47             :  * - Use $\cos\theta$ as angular coordinate instead of $\theta$. This avoids
      48             :  *   the $\cot\theta$ term by rewriting the angular derivatives as:
      49             :  *   \f[
      50             :  *   \partial_\theta^2+\cot\theta\partial_\theta =
      51             :  *   \partial_{\cos\theta}\sin^2\theta\partial_{\cos\theta}
      52             :  *   \f]
      53             :  * - Decompose $\Psi_m = \sin(\theta)^m u_m(r_\star, \theta)$. This avoids the
      54             :  *   $m^2/\sin^2\theta$ term by factoring out the singular behavior at the
      55             :  *   poles. The equations transform as:
      56             :  *   \f[
      57             :  *   -\partial_{\cos\theta}\sin^2\theta\partial_{\cos\theta} \Psi_m +
      58             :  *   \frac{m^2}{\sin^2\theta}\Psi_m = \sin(\theta)^m \left( m(m+1)
      59             :  *   + 2m \cos\theta \partial_{\cos\theta}
      60             :  *   - \partial_{\cos\theta}\sin^2\theta\partial_{\cos\theta} \right) u_m
      61             :  *   \f]
      62             :  *   We divide by $\sin(\theta)^m$ to get the equations for $u_m$.
      63             :  *
      64             :  * Written this way, the equations are regular at the poles and converge
      65             :  * exponentially. We also don't have to apply angular boundary conditions
      66             :  * because regularity at the poles is automatically enforced by the
      67             :  * $\sin^2\theta$ factor in the flux.
      68             :  *
      69             :  * The resulting factors in the equation are:
      70             :  * \begin{align}
      71             :  * &\alpha = \frac{\Delta}{(r^2 + a^2)^2} \sin^2\theta \\
      72             :  * &\beta = \left(-m^2\Omega^2 \Sigma^2 + 4a m^2 \Omega M r + \Delta \left[
      73             :  *   m (m + 1) + \frac{2M}{r}(1-\frac{a^2}{Mr}) + \frac{2iam}{r}
      74             :  *   \right]\right) \frac{1}{(r^2 + a^2)^2} \\
      75             :  * &\gamma_{r_\star} = -\frac{2iam}{r^2+a^2} + \frac{2a^2}{r}
      76             :  *   \frac{\alpha}{\sin^2\theta} \\
      77             :  * &\gamma_{\cos\theta} = 2 m \cos(\theta) \frac{\Delta}{(r^2 + a^2)^2}
      78             :  * \end{align}
      79             :  *
      80             :  * This class also provides the effective source $S_m^\mathrm{eff} = \Delta_m
      81             :  * \Psi_m^P$ and the singular field $\Psi_m^P$ in the regularized region (see
      82             :  * `ScalarSelfForce::FirstOrderSystem` and Sec. III in \cite Osburn:2022bby ).
      83             :  * The effective source is computed using the scalar EffectiveSource code by
      84             :  * Wardell et. al. (https://github.com/barrywardell/EffectiveSource and
      85             :  * \cite Wardell:2011gb ). It is then transformed to correspond to the m-mode
      86             :  * decomposition used in \cite Osburn:2022bby Eq. (2.8) as:
      87             :  * \begin{align}
      88             :  *   &\Psi_m^P = \frac{r}{2 \pi} e^{-i m \Delta\phi} \Phi_m^\mathrm{Wardell} \\
      89             :  *   &S_m^\mathrm{eff} = \frac{r}{2 \pi} e^{-i m \Delta\phi}
      90             :  *     \frac{\Delta\,(r^2 + a^2\cos^2\theta)}{(r^2 + a^2)^2}
      91             :  *     S_m^\mathrm{Wardell}
      92             :  *   \text{,}
      93             :  * \end{align}
      94             :  * where $\Delta\phi = \frac{a}{r_+ - r_-} \ln(\frac{r-r_+}{r-r_-})$ (Eq. (2.7)
      95             :  * in \cite Osburn:2022bby ).
      96             :  * We also divide by $\sin(\theta)^m$ to account for the change of variable from
      97             :  * $\Psi_m$ to $u_m$ described above.
      98             :  *
      99             :  * \par Impose equatorial symmetry
     100             :  * Since this work is restricted to equatorial orbits, we can enforce equatorial
     101             :  * symmetry by reformulating the equations in terms of the coordinate
     102             :  * $z^2 = \cos^2\theta$ instead of $z = \cos\theta$. This transform the factors
     103             :  * of first-order flux formulation as:
     104             :  * \begin{equation}
     105             :  * F^{\cos^2\theta}  = 4 \cos^2\theta F^{\cos\theta}
     106             :  * \text{.}
     107             :  * \end{equation}
     108             :  * \begin{equation}
     109             :  * \gamma_{\cos^2\theta} = 2 \cos\theta \gamma_{\cos\theta} + 2 \alpha
     110             :  * \text{.}
     111             :  * \end{equation}
     112             :  *
     113             :  * \par Hyperboloidal slicing
     114             :  * Transforming to a hyperboloidal time coordinate $s = t - h(r_*)$ can simplify
     115             :  * the solution a lot because the slice will only intersect a finite number of
     116             :  * wave fronts. In frequency domain this transformation means that partial
     117             :  * derivatives transform as:
     118             :  * \begin{align}
     119             :  *   &\partial_{r_*} \rightarrow \partial_{r_*} + i m\Omega H(r_*) \\
     120             :  *   &\partial_{r_*}^2 \rightarrow \partial_{r_*}^2 +
     121             :  *     2 i m\Omega H(r_*) \partial_{r_*} + i m\Omega H'(r_*)
     122             :  *     -m^2\Omega^2 H(r_*)^2
     123             :  * \end{align}
     124             :  * where $H(r_*) = h'(r_*)$ is the boost function that asymptotes to
     125             :  * $H(r_* \rightarrow \pm \infty) = \pm 1$.
     126             :  * This maps to the following additional terms in $\beta$ and $\gamma_{r_*}$:
     127             :  * \begin{align}
     128             :  *   &\beta \rightarrow \beta + i m\Omega \gamma_{r_*} H(r_*)
     129             :  *     - i m\Omega H'(r_*) + m^2\Omega^2 H(r_*)^2 \\
     130             :  *   &\gamma_{r_*} \rightarrow \gamma_{r_*} - 2 i m\Omega H(r_*)
     131             :  * \end{align}
     132             :  * For the boost function we choose here a simple sigmoid so that it is exactly
     133             :  * zero in a finite region around the puncture where we apply the effective
     134             :  * source, and transitions to -1 and 1 outside of this region. We choose
     135             :  * the $C^1$ continuous `smoothstep<1>` for this. Reasonable transition points
     136             :  * are from the boundaries of the effective source region to about $20M$
     137             :  * away from it, though this hasn't been tested very much yet. We may also want
     138             :  * to try jumping from $H=0$ to $H=\pm 1$ discontinuously at the boundary of the
     139             :  * effective source region, which is supported by the DG scheme (see
     140             :  * `fluxes_computer::is_discontinuous` in
     141             :  * `elliptic::protocols::FirstOrderSystem`), so (1) we don't have to resolve the
     142             :  * transition, and (2) we can more easily support the 2nd order self-force
     143             :  * source.
     144             :  */
     145           1 : class CircularOrbit : public elliptic::analytic_data::Background,
     146             :                       public elliptic::analytic_data::InitialGuess {
     147             :  public:
     148           0 :   struct BlackHoleMass {
     149           0 :     static constexpr Options::String help =
     150             :         "Kerr mass parameter 'M' of the black hole";
     151           0 :     using type = double;
     152             :   };
     153           0 :   struct BlackHoleSpin {
     154           0 :     static constexpr Options::String help =
     155             :         "Kerr dimensionless spin parameter 'chi' of the black hole";
     156           0 :     using type = double;
     157             :   };
     158           0 :   struct OrbitalRadius {
     159           0 :     static constexpr Options::String help =
     160             :         "Radius 'r_0' of the circular orbit";
     161           0 :     using type = double;
     162             :   };
     163           0 :   struct MModeNumber {
     164           0 :     static constexpr Options::String help =
     165             :         "Mode number 'm' of the scalar field";
     166           0 :     using type = int;
     167             :   };
     168           0 :   struct HyperboloidalSlicingTransitions {
     169           0 :     static constexpr Options::String help =
     170             :         "Enable hyperboloidal slicing by specifying the transition points for "
     171             :         "the boost function. The boost function transitions from -1 to zero "
     172             :         "between the first two points and from zero to 1 between the last "
     173             :         "two points. The effective source can only be evaluated where the "
     174             :         "boost function is zero, so the regularized region must be between "
     175             :         "the second and third points.";
     176           0 :     using type = Options::Auto<std::array<double, 4>, Options::AutoLabel::None>;
     177             :   };
     178           0 :   struct ImposeEquatorialSymmetry {
     179           0 :     static constexpr Options::String help =
     180             :         "Impose symmetry across the equatorial plane by using cos(theta)^2 "
     181             :         "as the angular coordinate instead of cos(theta). This means the "
     182             :         "domain should span [0, 1] instead of [-1, 1].";
     183           0 :     using type = bool;
     184             :   };
     185           0 :   using options =
     186             :       tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius, MModeNumber,
     187             :                  HyperboloidalSlicingTransitions, ImposeEquatorialSymmetry>;
     188           0 :   static constexpr Options::String help =
     189             :       "Quasicircular orbit of a scalar point charge in Kerr spacetime";
     190             : 
     191           0 :   CircularOrbit() = default;
     192           0 :   CircularOrbit(const CircularOrbit&) = default;
     193           0 :   CircularOrbit& operator=(const CircularOrbit&) = default;
     194           0 :   CircularOrbit(CircularOrbit&&) = default;
     195           0 :   CircularOrbit& operator=(CircularOrbit&&) = default;
     196           0 :   ~CircularOrbit() override = default;
     197             : 
     198           0 :   CircularOrbit(
     199             :       double black_hole_mass, double black_hole_spin, double orbital_radius,
     200             :       int m_mode_number,
     201             :       std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions,
     202             :       bool impose_equatorial_symmetry);
     203             : 
     204           0 :   explicit CircularOrbit(CkMigrateMessage* m);
     205             :   using PUP::able::register_constructor;
     206           0 :   WRAPPED_PUPable_decl_template(CircularOrbit);
     207             : 
     208           0 :   tnsr::I<double, 2> puncture_position() const;
     209           0 :   double black_hole_mass() const { return black_hole_mass_; }
     210           0 :   double black_hole_spin() const { return black_hole_spin_; }
     211           0 :   double orbital_radius() const { return orbital_radius_; }
     212           0 :   int m_mode_number() const { return m_mode_number_; }
     213           0 :   std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions()
     214             :       const {
     215             :     return hyperboloidal_slicing_transitions_;
     216             :   }
     217           0 :   bool impose_equatorial_symmetry() const {
     218             :     return impose_equatorial_symmetry_;
     219             :   }
     220             : 
     221           0 :   using background_tags = tmpl::list<Tags::Alpha, Tags::Beta, Tags::Gamma>;
     222           0 :   using source_tags = tmpl::list<
     223             :       ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
     224             :       ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
     225             :       Tags::BoyerLindquistRadius>;
     226             : 
     227             :   // Background
     228           0 :   tuples::tagged_tuple_from_typelist<background_tags> variables(
     229             :       const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
     230             : 
     231             :   // Initial guess
     232           0 :   static tuples::TaggedTuple<Tags::MMode> variables(
     233             :       const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
     234             : 
     235             :   // Fixed sources
     236           0 :   tuples::tagged_tuple_from_typelist<source_tags> variables(
     237             :       const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
     238             : 
     239             :   template <typename... RequestedTags>
     240           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     241             :       const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
     242             :       const InverseJacobian<DataVector, 2, Frame::ElementLogical,
     243             :                             Frame::Inertial>& /*inv_jacobian*/,
     244             :       tmpl::list<RequestedTags...> /*meta*/) const {
     245             :     return variables(x, tmpl::list<RequestedTags...>{});
     246             :   }
     247             : 
     248             :   // NOLINTNEXTLINE
     249           0 :   void pup(PUP::er& p) override;
     250             : 
     251             :  private:
     252           0 :   friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
     253             : 
     254           0 :   double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
     255           0 :   double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
     256           0 :   double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
     257           0 :   int m_mode_number_{};
     258           0 :   std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions_{};
     259           0 :   bool impose_equatorial_symmetry_{false};
     260             : };
     261             : 
     262           0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
     263             : 
     264             : }  // namespace ScalarSelfForce::AnalyticData