Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <limits>
       7             : #include <memory>
       8             : #include <pup.h>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Evolution/Systems/ForceFree/Tags.hpp"
      12             : #include "Options/Options.hpp"
      13             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      14             : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp"
      15             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : namespace PUP {
      21             : class er;
      22             : }  // namespace PUP
      23             : /// \endcond
      24             : 
      25             : namespace ForceFree::AnalyticData {
      26             : /*!
      27             :  * \brief The magnetospheric Wald problem proposed in \cite Komissarov2004
      28             :  *
      29             :  * This is an initial value problem that evolves the magnetosphere of a rotating
      30             :  * black hole. The initial condition is given as same as the exact Wald solution
      31             :  * \cite Wald1974 (see also documentation of ForceFree::Solutions::ExactWald)
      32             :  *
      33             :  * \begin{equation}
      34             :  *  A_\mu = \frac{B_0}{2}(\phi_\mu + 2a t_\mu) ,
      35             :  * \end{equation}
      36             :  *
      37             :  * but electric field is set to zero at $t=0$.
      38             :  *
      39             :  * In the cartesian projection of the spherical Kerr-Schild coordinates
      40             :  * (which we use in the code for representing tensors), initial magnetic fields
      41             :  * is given as
      42             :  *
      43             :  * \begin{align}
      44             :  * \tilde{B}^{x} &= a B_0 z \left[ (ax-ry) \left\{
      45             :  *      \frac{1}{r^4} + \frac{2M r (r^2-a^2)}{(r^4+a^2z^2)^2} \right\}
      46             :  *      + a M r x \left\{ \frac{r^2-z^2}{r^4(r^4+a^2z^2)}
      47             :  *      - \frac{4(r^2+z^2)}{(r^4+a^2z^2)^2} \right\} \right] \\
      48             :  * \tilde{B}^{y} &= a B_0 z \left[ (rx+ay) \left\{
      49             :  *      \frac{1}{r^4} + \frac{2M r (r^2-a^2)}{(r^4+a^2z^2)^2} \right\}
      50             :  *      + a M r y \left\{ \frac{r^2-z^2}{r^4(r^4+a^2z^2)}
      51             :  *      - \frac{4(r^2+z^2)}{(r^4+a^2z^2)^2} \right\} \right] \\
      52             :  * \tilde{B}^{z} &= B_0 \left[
      53             :  *      1 + \frac{a^2z^2}{r^4} + \frac{M a^2}{r^3}\left\{
      54             :  *      1 - \frac{z^2(a^2+z^2)(5r^4+a^2z^2)}{(r^4+a^2z^2)^2} \right\} \right] .
      55             :  * \end{align}
      56             :  *
      57             :  * where $M$ and $a$ are mass and (dimensionless) spin of the Kerr black hole,
      58             :  * $B_0$ is the amplitude of magnetic field, and $r$ is the radial coordinate
      59             :  * defined in the spherical Kerr-Schild coordinate system (see the documentation
      60             :  * of gr::Solutions::SphericalKerrSchild). All other variables are set to zero
      61             :  * at $t=0$.
      62             :  *
      63             :  * There is no known exact solution to this problem, but numerical simulations
      64             :  * \cite Komissarov2004 \cite Paschalidis2013 \cite Etienne2017 report that the
      65             :  * system converges to a steady state with an equatorial current sheet inside
      66             :  * the ergosphere.
      67             :  *
      68             :  * \note We set $M=1$ and $B_0=1$ in the initial data to fix scales.
      69             :  *
      70             :  */
      71           1 : class MagnetosphericWald : public evolution::initial_data::InitialData,
      72             :                            public MarkAsAnalyticData {
      73             :  public:
      74           0 :   struct Spin {
      75           0 :     using type = double;
      76           0 :     static constexpr Options::String help = {
      77             :         "The dimensionless spin of the Kerr BH"};
      78           0 :     static type upper_bound() { return 1.0; }
      79           0 :     static type lower_bound() { return -1.0; }
      80             :   };
      81             : 
      82           0 :   using options = tmpl::list<Spin>;
      83           0 :   static constexpr Options::String help{
      84             :       "Magnetospheric Wald initial value problem"};
      85             : 
      86           0 :   MagnetosphericWald() = default;
      87           0 :   MagnetosphericWald(const MagnetosphericWald&) = default;
      88           0 :   MagnetosphericWald& operator=(const MagnetosphericWald&) = default;
      89           0 :   MagnetosphericWald(MagnetosphericWald&&) = default;
      90           0 :   MagnetosphericWald& operator=(MagnetosphericWald&&) = default;
      91           0 :   ~MagnetosphericWald() override = default;
      92             : 
      93           0 :   explicit MagnetosphericWald(double spin,
      94             :                               const Options::Context& context = {});
      95             : 
      96           0 :   auto get_clone() const
      97             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
      98             : 
      99             :   /// \cond
     100             :   explicit MagnetosphericWald(CkMigrateMessage* msg);
     101             :   using PUP::able::register_constructor;
     102             :   WRAPPED_PUPable_decl_template(MagnetosphericWald);
     103             :   /// \endcond
     104             : 
     105             :   // NOLINTNEXTLINE(google-runtime-references)
     106           0 :   void pup(PUP::er& p) override;
     107             : 
     108             :   /// @{
     109             :   /// Retrieve the EM variables at (x,t).
     110           1 :   static auto variables(const tnsr::I<DataVector, 3>& x,
     111             :                         tmpl::list<Tags::TildeE> /*meta*/)
     112             :       -> tuples::TaggedTuple<Tags::TildeE>;
     113             : 
     114           1 :   auto variables(const tnsr::I<DataVector, 3>& x,
     115             :                  tmpl::list<Tags::TildeB> /*meta*/) const
     116             :       -> tuples::TaggedTuple<Tags::TildeB>;
     117             : 
     118           1 :   static auto variables(const tnsr::I<DataVector, 3>& x,
     119             :                         tmpl::list<Tags::TildePsi> /*meta*/)
     120             :       -> tuples::TaggedTuple<Tags::TildePsi>;
     121             : 
     122           1 :   static auto variables(const tnsr::I<DataVector, 3>& x,
     123             :                         tmpl::list<Tags::TildePhi> /*meta*/)
     124             :       -> tuples::TaggedTuple<Tags::TildePhi>;
     125             : 
     126           1 :   static auto variables(const tnsr::I<DataVector, 3>& x,
     127             :                         tmpl::list<Tags::TildeQ> /*meta*/)
     128             :       -> tuples::TaggedTuple<Tags::TildeQ>;
     129             :   /// @}
     130             : 
     131             :   /// Retrieve a collection of EM variables at position x
     132             :   template <typename... Tags>
     133           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataVector, 3>& x,
     134             :                                          tmpl::list<Tags...> /*meta*/) const {
     135             :     static_assert(sizeof...(Tags) > 1,
     136             :                   "The generic template will recurse infinitely if only one "
     137             :                   "tag is being retrieved.");
     138             :     return {get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     139             :   }
     140             : 
     141             :   /// Retrieve the metric variables
     142             :   template <typename Tag>
     143           1 :   tuples::TaggedTuple<Tag> variables(const tnsr::I<DataVector, 3>& x,
     144             :                                      tmpl::list<Tag> /*meta*/) const {
     145             :     constexpr double dummy_time = 0.0;
     146             :     return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
     147             :   }
     148             : 
     149             :  private:
     150           0 :   double spin_ = std::numeric_limits<double>::signaling_NaN();
     151           0 :   gr::Solutions::SphericalKerrSchild background_spacetime_{};
     152             : 
     153           0 :   friend bool operator==(const MagnetosphericWald& lhs,
     154             :                          const MagnetosphericWald& rhs);
     155             : };
     156             : 
     157           0 : bool operator!=(const MagnetosphericWald& lhs, const MagnetosphericWald& rhs);
     158             : 
     159             : }  // namespace ForceFree::AnalyticData