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