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 F^i = 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             :  * Note that we use $\cos\theta$ as angular coordinate but \cite Osburn:2022bby
      44             :  * uses $\theta$. We also multiply Eq. (2.9) by the factor $\Sigma^2 / (r^2 +
      45             :  * a^2)^2$ so we can easily write it in first-order flux form. The resulting
      46             :  * factors in the equation are:
      47             :  *
      48             :  * \begin{align}
      49             :  * &\alpha = \frac{\Delta}{(r^2 + a^2)^2} \sin^2\theta \\
      50             :  * &\beta = \left(-m^2\Omega^2 \Sigma^2 + 4a m^2 \Omega M r + \Delta \left[
      51             :  *   \frac{m^2}{\sin^2\theta} + \frac{2M}{r}(1-\frac{a^2}{Mr}) + \frac{2iam}{r}
      52             :  *   \right]\right) \frac{1}{(r^2 + a^2)^2} \\
      53             :  * &\gamma_{r_\star} = -\frac{2iam}{r^2+a^2} + \frac{2a^2}{r}
      54             :  *   \frac{\alpha}{\sin^2\theta} \\
      55             :  * &\gamma_{\cos\theta} = 0
      56             :  * \end{align}
      57             :  *
      58             :  * This class also provides the effective source $S_m^\mathrm{eff} = \Delta_m
      59             :  * \Psi_m^P$ and the singular field $\Psi_m^P$ in the regularized region (see
      60             :  * `ScalarSelfForce::FirstOrderSystem` and Sec. III in \cite Osburn:2022bby ).
      61             :  * The effective source is computed using the scalar EffectiveSource code by
      62             :  * Wardell et. al. (https://github.com/barrywardell/EffectiveSource and
      63             :  * \cite Wardell:2011gb ). It is then transformed to correspond to the m-mode
      64             :  * decomposition used in \cite Osburn:2022bby Eq. (2.8) as:
      65             :  * \begin{align}
      66             :  *   &\Psi_m^P = \frac{r}{2 \pi} e^{-i m \Delta\phi} \Phi_m^\mathrm{Wardell} \\
      67             :  *   &S_m^\mathrm{eff} = \frac{r}{2 \pi} e^{-i m \Delta\phi}
      68             :  *     \frac{\Delta\,(r^2 + a^2\cos^2\theta)}{(r^2 + a^2)^2}
      69             :  *     S_m^\mathrm{Wardell}
      70             :  *   \text{,}
      71             :  * \end{align}
      72             :  * where $\Delta\phi = \frac{a}{r_+ - r_-} \ln(\frac{r-r_+}{r-r_-})$ (Eq. (2.7)
      73             :  * in \cite Osburn:2022bby ).
      74             :  *
      75             :  * \par Hyperboloidal slicing
      76             :  * Transforming to a hyperboloidal time coordinate $s = t - h(r_*)$ can simplify
      77             :  * the solution a lot because the slice will only intersect a finite number of
      78             :  * wave fronts. In frequency domain this transformation means that partial
      79             :  * derivatives transform as:
      80             :  * \begin{align}
      81             :  *   &\partial_{r_*} \rightarrow \partial_{r_*} + i m\Omega H(r_*) \\
      82             :  *   &\partial_{r_*}^2 \rightarrow \partial_{r_*}^2 +
      83             :  *     2 i m\Omega H(r_*) \partial_{r_*} + i m\Omega H'(r_*)
      84             :  *     -m^2\Omega^2 H(r_*)^2
      85             :  * \end{align}
      86             :  * where $H(r_*) = h'(r_*)$ is the boost function that asymptotes to
      87             :  * $H(r_* \rightarrow \pm \infty) = \pm 1$.
      88             :  * This maps to the following additional terms in $\beta$ and $\gamma_{r_*}$:
      89             :  * \begin{align}
      90             :  *   &\beta \rightarrow \beta + i m\Omega \gamma_{r_*} H(r_*)
      91             :  *     - i m\Omega H'(r_*) + m^2\Omega^2 H(r_*)^2 \\
      92             :  *   &\gamma_{r_*} \rightarrow \gamma_{r_*} - 2 i m\Omega H(r_*)
      93             :  * \end{align}
      94             :  * For the boost function we choose here a simple sigmoid so that it is exactly
      95             :  * zero in a finite region around the puncture where we apply the effective
      96             :  * source, and transitions to -1 and 1 outside of this region. We choose
      97             :  * the $C^1$ continuous `smoothstep<1>` for this. Reasonable transition points
      98             :  * are from the boundaries of the effective source region to about $20M$
      99             :  * away from it, though this hasn't been tested very much yet. We may also want
     100             :  * to try jumping from $H=0$ to $H=\pm 1$ discontinuously at the boundary of the
     101             :  * effective source region, which is supported by the DG scheme (see
     102             :  * `fluxes_computer::is_discontinuous` in
     103             :  * `elliptic::protocols::FirstOrderSystem`), so (1) we don't have to resolve the
     104             :  * transition, and (2) we can more easily support the 2nd order self-force
     105             :  * source.
     106             :  */
     107           1 : class CircularOrbit : public elliptic::analytic_data::Background,
     108             :                       public elliptic::analytic_data::InitialGuess {
     109             :  public:
     110           0 :   struct BlackHoleMass {
     111           0 :     static constexpr Options::String help =
     112             :         "Kerr mass parameter 'M' of the black hole";
     113           0 :     using type = double;
     114             :   };
     115           0 :   struct BlackHoleSpin {
     116           0 :     static constexpr Options::String help =
     117             :         "Kerr dimensionless spin parameter 'chi' of the black hole";
     118           0 :     using type = double;
     119             :   };
     120           0 :   struct OrbitalRadius {
     121           0 :     static constexpr Options::String help =
     122             :         "Radius 'r_0' of the circular orbit";
     123           0 :     using type = double;
     124             :   };
     125           0 :   struct MModeNumber {
     126           0 :     static constexpr Options::String help =
     127             :         "Mode number 'm' of the scalar field";
     128           0 :     using type = int;
     129             :   };
     130           0 :   struct HyperboloidalSlicingTransitions {
     131           0 :     static constexpr Options::String help =
     132             :         "Enable hyperboloidal slicing by specifying the transition points for "
     133             :         "the boost function. The boost function transitions from -1 to zero "
     134             :         "between the first two points and from zero to 1 between the last "
     135             :         "two points. The effective source can only be evaluated where the "
     136             :         "boost function is zero, so the regularized region must be between "
     137             :         "the second and third points.";
     138           0 :     using type = Options::Auto<std::array<double, 4>, Options::AutoLabel::None>;
     139             :   };
     140           0 :   using options = tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius,
     141             :                              MModeNumber, HyperboloidalSlicingTransitions>;
     142           0 :   static constexpr Options::String help =
     143             :       "Quasicircular orbit of a scalar point charge in Kerr spacetime";
     144             : 
     145           0 :   CircularOrbit() = default;
     146           0 :   CircularOrbit(const CircularOrbit&) = default;
     147           0 :   CircularOrbit& operator=(const CircularOrbit&) = default;
     148           0 :   CircularOrbit(CircularOrbit&&) = default;
     149           0 :   CircularOrbit& operator=(CircularOrbit&&) = default;
     150           0 :   ~CircularOrbit() override = default;
     151             : 
     152           0 :   CircularOrbit(
     153             :       double black_hole_mass, double black_hole_spin, double orbital_radius,
     154             :       int m_mode_number,
     155             :       std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions);
     156             : 
     157           0 :   explicit CircularOrbit(CkMigrateMessage* m);
     158             :   using PUP::able::register_constructor;
     159           0 :   WRAPPED_PUPable_decl_template(CircularOrbit);
     160             : 
     161           0 :   tnsr::I<double, 2> puncture_position() const;
     162           0 :   double black_hole_mass() const { return black_hole_mass_; }
     163           0 :   double black_hole_spin() const { return black_hole_spin_; }
     164           0 :   double orbital_radius() const { return orbital_radius_; }
     165           0 :   int m_mode_number() const { return m_mode_number_; }
     166           0 :   std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions()
     167             :       const {
     168             :     return hyperboloidal_slicing_transitions_;
     169             :   }
     170             : 
     171           0 :   using background_tags = tmpl::list<Tags::Alpha, Tags::Beta, Tags::Gamma>;
     172           0 :   using source_tags = tmpl::list<
     173             :       ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
     174             :       ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
     175             :       Tags::BoyerLindquistRadius>;
     176             : 
     177             :   // Background
     178           0 :   tuples::tagged_tuple_from_typelist<background_tags> variables(
     179             :       const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
     180             : 
     181             :   // Initial guess
     182           0 :   static tuples::TaggedTuple<Tags::MMode> variables(
     183             :       const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
     184             : 
     185             :   // Fixed sources
     186           0 :   tuples::tagged_tuple_from_typelist<source_tags> variables(
     187             :       const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
     188             : 
     189             :   template <typename... RequestedTags>
     190           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     191             :       const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
     192             :       const InverseJacobian<DataVector, 2, Frame::ElementLogical,
     193             :                             Frame::Inertial>& /*inv_jacobian*/,
     194             :       tmpl::list<RequestedTags...> /*meta*/) const {
     195             :     return variables(x, tmpl::list<RequestedTags...>{});
     196             :   }
     197             : 
     198             :   // NOLINTNEXTLINE
     199           0 :   void pup(PUP::er& p) override;
     200             : 
     201             :  private:
     202           0 :   friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
     203             : 
     204           0 :   double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
     205           0 :   double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
     206           0 :   double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
     207           0 :   int m_mode_number_{};
     208           0 :   std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions_{};
     209             : };
     210             : 
     211           0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
     212             : 
     213             : }  // namespace ScalarSelfForce::AnalyticData