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/Minkowski.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             : /*!
      28             :  * \brief The magnetosphere of an isolated rotating star with dipolar initial
      29             :  * magnetic field in the flat spacetime. This is a toy model of a pulsar
      30             :  * magnetosphere.
      31             :  *
      32             :  * \note Coordinate radius of the star is rescaled to 1.0 in code units.
      33             :  *
      34             :  * The vector potential of the initial magnetic field has the form
      35             :  * \cite Most2022
      36             :  *
      37             :  * \begin{equation}
      38             :  *  A_\phi = \frac{A_0 \varpi_0 (x^2+y^2)}{(r^2 + \delta^2)^{3/2}}
      39             :  * \end{equation}
      40             :  *
      41             :  * where $A_0$ is the vector potential amplitude, $\varpi_0$ is a constant with
      42             :  * the unit of length, $r^2 = x^2 + y^2 + z^2$, and $\delta$ is a small number
      43             :  * for regularization of the dipole magnetic field at the origin ($r=0$).
      44             :  *
      45             :  * In the Cartesian coordinates, components of densitized magnetic fields are
      46             :  * given as
      47             :  *
      48             :  * \begin{align}
      49             :  *  \tilde{B}^x & = A_0 \varpi_0 \frac{3xz}{(r^2 + \delta^2)^{5/2}} , \\
      50             :  *  \tilde{B}^y & = A_0 \varpi_0 \frac{3yz}{(r^2 + \delta^2)^{5/2}} , \\
      51             :  *  \tilde{B}^z & = A_0 \varpi_0 \frac{3z^2 - r^2 + 2\delta^2}{(r^2 +
      52             :  *                    \delta^2)^{5/2}} .
      53             :  * \end{align}
      54             :  *
      55             :  * Rotation of the star is switched on at $t=0$ with the angular velocity
      56             :  * specified in the input file. The grid points inside the star ($x^2 + y^2 +
      57             :  * z^2 < 1.0$) are identified as the interior and masked by `interior_mask()`
      58             :  * member function. In the masked region we impose the MHD condition
      59             :  * $\mathbf{E} + \mathbf{v} \times \mathbf{B} = 0$, where the velocity field is
      60             :  * given by $\mathbf{v} \equiv \Omega \hat{z} \times \mathbf{r}$. By this means,
      61             :  * initial dipolar magnetic field is effectively "anchored" inside the star
      62             :  * while electromagnetic fields are evolved self-consistently outside the star.
      63             :  *
      64             :  * \note We impose the MHD condition stated above and $q=0$ inside the masked
      65             :  * interior region during the evolution phase. While this initial data class
      66             :  * sets both electric field and charge density to zero at $t=0$, those variables
      67             :  * are immediately overwritten with proper values ($\mathbf{E} =
      68             :  * -\mathbf{v}\times\mathbf{B}$, $q=0$) once the simulation begins.
      69             :  *
      70             :  * When the system reaches a stationary state, magnetic field lines far from the
      71             :  * star are opening up while the field lines close to the star are corotating.
      72             :  * The light cylinder and the Y-point, which marks the boundary between these
      73             :  * two regions with different magnetic field topology, are expected to be formed
      74             :  * at $r_\text{LC} = c/\Omega$.
      75             :  *
      76             :  * The option `TiltAngle` controls the angle $\alpha$ between the rotation axis
      77             :  * ($z$) and the magnetic axis of initial magnetic field on the $x-z$ plane. An
      78             :  * aligned rotator ($\alpha = 0$) is a common test problem for FFE codes,
      79             :  * whereas an oblique rotator ($\alpha \neq 0$) is a more realistic model of
      80             :  * pulsars \cite Spitkovsky2006.
      81             :  *
      82             :  */
      83           1 : class RotatingDipole : public evolution::initial_data::InitialData,
      84             :                        public MarkAsAnalyticData {
      85             :  public:
      86           0 :   struct VectorPotentialAmplitude {
      87           0 :     using type = double;
      88           0 :     static constexpr Options::String help = {
      89             :         "The vector potential amplitude A_0"};
      90             :   };
      91             : 
      92           0 :   struct Varpi0 {
      93           0 :     using type = double;
      94           0 :     static constexpr Options::String help = {"The length constant varpi_0"};
      95           0 :     static type lower_bound() { return 0.0; }
      96             :   };
      97             : 
      98           0 :   struct Delta {
      99           0 :     using type = double;
     100           0 :     static constexpr Options::String help = {
     101             :         "A small value used to regularize magnetic fields at r=0."};
     102           0 :     static type lower_bound() { return 0.0; }
     103             :   };
     104             : 
     105           0 :   struct AngularVelocity {
     106           0 :     using type = double;
     107           0 :     static constexpr Options::String help = {
     108             :         "Rotation angular velocity of the star."};
     109           0 :     static type upper_bound() { return 1.0; }
     110           0 :     static type lower_bound() { return -1.0; }
     111             :   };
     112             : 
     113           0 :   struct TiltAngle {
     114           0 :     using type = double;
     115           0 :     static constexpr Options::String help = {
     116             :         "Angle between the rotation axis (z) and magnetic axis at t = 0."};
     117           0 :     static type upper_bound() { return M_PI; }
     118           0 :     static type lower_bound() { return 0.0; }
     119             :   };
     120             : 
     121           0 :   using options = tmpl::list<VectorPotentialAmplitude, Varpi0, Delta,
     122             :                              AngularVelocity, TiltAngle>;
     123           0 :   static constexpr Options::String help{
     124             :       "Magnetosphere of an isolated rotating star with dipole magnetic field."};
     125             : 
     126           0 :   RotatingDipole() = default;
     127           0 :   RotatingDipole(const RotatingDipole&) = default;
     128           0 :   RotatingDipole& operator=(const RotatingDipole&) = default;
     129           0 :   RotatingDipole(RotatingDipole&&) = default;
     130           0 :   RotatingDipole& operator=(RotatingDipole&&) = default;
     131           0 :   ~RotatingDipole() override = default;
     132             : 
     133           0 :   RotatingDipole(double vector_potential_amplitude, double varpi0, double delta,
     134             :                  double angular_velocity, double tilt_angle,
     135             :                  const Options::Context& context = {});
     136             : 
     137           0 :   auto get_clone() const
     138             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     139             : 
     140             :   /// \cond
     141             :   explicit RotatingDipole(CkMigrateMessage* msg);
     142             :   using PUP::able::register_constructor;
     143             :   WRAPPED_PUPable_decl_template(RotatingDipole);
     144             :   /// \endcond
     145             : 
     146             :   // NOLINTNEXTLINE(google-runtime-references)
     147           0 :   void pup(PUP::er& p) override;
     148             : 
     149             :   /// @{
     150             :   /// Retrieve the EM variables.
     151           1 :   static auto variables(const tnsr::I<DataVector, 3>& coords,
     152             :                         tmpl::list<Tags::TildeE> /*meta*/)
     153             :       -> tuples::TaggedTuple<Tags::TildeE>;
     154             : 
     155           1 :   auto variables(const tnsr::I<DataVector, 3>& coords,
     156             :                  tmpl::list<Tags::TildeB> /*meta*/) const
     157             :       -> tuples::TaggedTuple<Tags::TildeB>;
     158             : 
     159           1 :   static auto variables(const tnsr::I<DataVector, 3>& coords,
     160             :                         tmpl::list<Tags::TildePsi> /*meta*/)
     161             :       -> tuples::TaggedTuple<Tags::TildePsi>;
     162             : 
     163           1 :   static auto variables(const tnsr::I<DataVector, 3>& coords,
     164             :                         tmpl::list<Tags::TildePhi> /*meta*/)
     165             :       -> tuples::TaggedTuple<Tags::TildePhi>;
     166             : 
     167           1 :   static auto variables(const tnsr::I<DataVector, 3>& coords,
     168             :                         tmpl::list<Tags::TildeQ> /*meta*/)
     169             :       -> tuples::TaggedTuple<Tags::TildeQ>;
     170             :   /// @}
     171             : 
     172             :   /// Retrieve a collection of EM variables at position x
     173             :   template <typename... Tags>
     174           1 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataVector, 3>& x,
     175             :                                          tmpl::list<Tags...> /*meta*/) const {
     176             :     static_assert(sizeof...(Tags) > 1,
     177             :                   "The generic template will recurse infinitely if only one "
     178             :                   "tag is being retrieved.");
     179             :     return {get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     180             :   }
     181             : 
     182             :   /// Retrieve the metric variables
     183             :   template <typename Tag>
     184           1 :   tuples::TaggedTuple<Tag> variables(const tnsr::I<DataVector, 3>& x,
     185             :                                      tmpl::list<Tag> /*meta*/) const {
     186             :     constexpr double dummy_time = 0.0;
     187             :     return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
     188             :   }
     189             : 
     190             :   // Returns the value of NS interior mask
     191           0 :   static std::optional<Scalar<DataVector>> interior_mask(
     192             :       const tnsr::I<DataVector, 3, Frame::Inertial>& x);
     193             : 
     194             :   // Returns the value of angular velocity.
     195           0 :   double angular_velocity() const { return angular_velocity_; };
     196             : 
     197             :  private:
     198           0 :   double vector_potential_amplitude_ =
     199             :       std::numeric_limits<double>::signaling_NaN();
     200           0 :   double varpi0_ = std::numeric_limits<double>::signaling_NaN();
     201           0 :   double delta_ = std::numeric_limits<double>::signaling_NaN();
     202           0 :   double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
     203           0 :   double tilt_angle_ = std::numeric_limits<double>::signaling_NaN();
     204           0 :   gr::Solutions::Minkowski<3> background_spacetime_{};
     205             : 
     206           0 :   friend bool operator==(const RotatingDipole& lhs, const RotatingDipole& rhs);
     207             : };
     208             : 
     209           0 : bool operator!=(const RotatingDipole& lhs, const RotatingDipole& rhs);
     210             : 
     211             : }  // namespace ForceFree::AnalyticData