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
|
