MagnetizedTovStar.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <cstddef>
7 #include <limits>
8
10 #include "Options/Options.hpp"
11 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
12 #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
13 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Tov.hpp"
14 #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp"
15 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
16 #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 #include "Utilities/TMPL.hpp"
18 #include "Utilities/TaggedTuple.hpp"
19
20 /// \cond
21 namespace PUP {
22 class er; // IWYU pragma: keep
23 } // namespace PUP
24 /// \endcond
25
26 namespace grmhd::AnalyticData {
27 /*!
28  * \brief Magnetized TOV star initial data, where metric terms only account for
29  * the hydrodynamics not the magnetic fields.
30  *
31  * Superposes a poloidal magnetic field on top of a TOV solution where the
32  * vector potential has the form
33  *
34  * \f{align*}{
35  * A_{\phi} = A_b \varpi^2 \max(p-p_{\mathrm{cut}}, 0)^{n_s}
36  * \f}
37  *
38  * where \f$A_b\f$ controls the amplitude of the magnetic field,
39  * \f$\varpi^2=x^2+y^2=r^2-z^2\f$ is the cylindrical radius,
40  * \f$n_s\f$ controls the degree of differentiability, and
41  * \f$p_{\mathrm{cut}}\f$ controls the pressure cutoff below which the magnetic
42  * field is zero.
43  *
44  * In Cartesian coordinates the vector potential is:
45  *
46  * \f{align*}{
47  * A_x&=-\frac{y}{\varpi^2}A_\phi, \\
48  * A_y&=\frac{x}{\varpi^2}A_\phi, \\
49  * A_z&=0,
50  * \f}
51  *
52  * For the poloidal field this means
53  *
54  * \f{align*}{
55  * A_x &= -y A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
56  * A_y &= x A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}.
57  * \f}
58  *
59  * The magnetic field is computed from the vector potential as
60  *
61  * \f{align*}{
62  * B^i&=n_a\epsilon^{aijk}\partial_jA_k \\
63  * &=-\frac{1}{\sqrt{\gamma}}[ijk]\partial_j A_k,
64  * \f}
65  *
66  * where \f$[ijk]\f$ is the total anti-symmetric symbol. This means that
67  *
68  * \f{align*}{
69  * B^x&=\frac{1}{\sqrt{\gamma}} (\partial_z A_y-\partial_y A_z), \\
70  * B^y&=\frac{1}{\sqrt{\gamma}} (\partial_x A_z - \partial_z A_x), \\
71  * B^z&=\frac{1}{\sqrt{\gamma}} (\partial_y A_x - \partial_x A_y).
72  * \f}
73  *
74  * Focusing on the region where the field is non-zero we have:
75  *
76  * \f{align*}{
77  * \partial_x A_y
78  * &= A_b(p-p_{\mathrm{cut}})^{n_s}+
79  * \frac{x^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\
80  * \partial_y A_x
81  * &= -A_b(p-p_{\mathrm{cut}})^{n_s}-
82  * \frac{y^2}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p \\
83  * \partial_z A_x
84  * &= -\frac{yz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
85  * \partial_z A_y
86  * &= \frac{xz}{r}A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp
87  * \f}
88  *
89  * The magnetic field is given by:
90  *
91  * \f{align*}{
92  * B^x&=\frac{1}{\sqrt{\gamma}}\frac{xz}{r}
93  * A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
94  * B^y&=\frac{1}{\sqrt{\gamma}}\frac{yz}{r}
95  * A_bn_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_rp \\
96  * B^z&=-\frac{A_b}{\sqrt{\gamma}}\left[
97  * 2(p-p_{\mathrm{cut}})^{n_s} \phantom{\frac{a}{b}}\right. \\
98  * &\left.+\frac{x^2+y^2}{r}
99  * n_s(p-p_{\mathrm{cut}})^{n_s-1}\partial_r p
100  * \right]
101  * \f}
102  *
103  * Taking the small-\f$r\f$ limit gives the \f$r=0\f$ magnetic field:
104  *
105  * \f{align*}{
106  * B^x&=0, \\
107  * B^y&=0, \\
108  * B^z&=-\frac{A_b}{\sqrt{\gamma}}
109  * 2(p-p_{\mathrm{cut}})^{n_s}.
110  * \f}
111  *
112  * While the amplitude \f$A_b\f$ is specified in the code, it is more natural
113  * to work with the magnetic field strength, which is given by \f$\sqrt{b^2}\f$
114  * (where \f$b^a\f$ is the comoving magnetic field), and in CGS units is
115  *
116  * \f{align*}{
117  * |B_{\mathrm{CGS}}|&= \sqrt{4 \pi b^2}
118  * \left(\frac{c^2}{G M_\odot}\right) \left(\frac{c}{\sqrt{4 \pi \epsilon_0
119  * G}}\right) \\
120  * &= \sqrt{b^2} \times 8.352\times10^{19}\mathrm{G} \,.
121  * \f}
122  *
123  * We now give values used for standard tests of magnetized stars.
124  * - \f$\rho_c(0)=1.28\times10^{-3}\f$
125  * - \f$K=100\f$
126  * - \f$\Gamma=2\f$
127  * - %Domain \f$[-20,20]^3\f$
128  * - Units \f$M=M_\odot\f$
129  * - A target final time 20ms means \f$20\times10^{-3}/(5\times10^{-6})=4000M\f$
130  * - The mass of the star is \f$1.4M_{\odot}\f$
131  *
132  * Parameters for desired magnetic field strength:
133  * - For \f$n_s=0\f$ and \f$p_{\mathrm{cut}}=0.04p_{\max}\f$ setting
134  * \f$A_b=6\times10^{-5}\f$ yields a maximum mganetic field strength of
135  * \f$1.002\times10^{16}G\f$.
136  * - For \f$n_s=1\f$ and \f$p_{\mathrm{cut}}=0.04p_{\max}\f$ setting
137  * \f$A_b=0.4\f$ yields a maximum mganetic field strength of
138  * \f$1.05\times10^{16}G\f$.
139  * - For \f$n_s=2\f$ and \f$p_{\mathrm{cut}}=0.04p_{\max}\f$ setting
140  * \f$A_b=2500\f$ yields a maximum mganetic field strength of
141  * \f$1.03\times10^{16}G\f$.
142  * - For \f$n_s=3\f$ and \f$p_{\mathrm{cut}}=0.04p_{\max}\f$ setting
143  * \f$A_b=1.65\times10^{7}\f$ yields a maximum mganetic field strength of
144  * \f$1.07\times10^{16}G\f$.
145  *
146  * Note that the magnetic field strength goes as \f$A_b\f$ so any desired value
147  * can be achieved by a linear scaling.
148  */
149 class MagnetizedTovStar : public MarkAsAnalyticData,
151  gr::Solutions::TovSolution> {
152  private:
153  using tov_star =
155
156  public:
158  using type = size_t;
159  static constexpr Options::String help = {
160  "The exponent n_s controlling the smoothness of the field"};
161  };
162
164  using type = double;
165  static constexpr Options::String help = {
166  "The amplitude A_b of the phi-component of the vector potential. This "
167  "controls the magnetic field strength."};
168  static type lower_bound() { return 0.0; }
169  };
170
172  using type = double;
173  static constexpr Options::String help = {
174  "The fraction of the central pressure below which there is no magnetic "
175  "field. p_cut = Fraction * p_max."};
176  static type lower_bound() { return 0.0; }
177  static type upper_bound() { return 1.0; }
178  };
179
180  using options =
181  tmpl::push_back<tov_star::options, PressureExponent,
183
184  static constexpr Options::String help = {
185  "Magnetized TOV star in areal coordinates."};
186
187  static constexpr size_t volume_dim = 3_st;
188
189  template <typename DataType>
190  using tags = typename tov_star::template tags<DataType>;
191
192  MagnetizedTovStar() = default;
193  MagnetizedTovStar(const MagnetizedTovStar& /*rhs*/) = delete;
194  MagnetizedTovStar& operator=(const MagnetizedTovStar& /*rhs*/) = delete;
195  MagnetizedTovStar(MagnetizedTovStar&& /*rhs*/) noexcept = default;
196  MagnetizedTovStar& operator=(MagnetizedTovStar&& /*rhs*/) noexcept = default;
197  ~MagnetizedTovStar() = default;
198
199  MagnetizedTovStar(double central_rest_mass_density,
200  double polytropic_constant, double polytropic_exponent,
201  size_t pressure_exponent, double cutoff_pressure_fraction,
202  double vector_potential_amplitude) noexcept;
203
204  // Overload the variables function from the base class.
205  using tov_star::equation_of_state;
206  using tov_star::equation_of_state_type;
207  using tov_star::variables;
208
209  /// Retrieve a collection of variables at (x)
210  template <typename DataType, typename... Tags>
211  tuples::TaggedTuple<Tags...> variables(
212  const tnsr::I<DataType, 3>& x,
213  tmpl::list<Tags...> /*meta*/) const noexcept {
217  }
218
219  void pup(PUP::er& p) noexcept; // NOLINT
220
221  private:
222  friend bool operator==(const MagnetizedTovStar& lhs,
223  const MagnetizedTovStar& rhs) noexcept;
224
225  protected:
226  template <typename DataType>
228  variables(const tnsr::I<DataType, 3>& coords,
229  tmpl::list<hydro::Tags::MagneticField<DataType, 3,
230  Frame::Inertial>> /*meta*/,
232
233  size_t pressure_exponent_ = std::numeric_limits<size_t>::max();
234  double cutoff_pressure_ = std::numeric_limits<double>::signaling_NaN();
235  double vector_potential_amplitude_ =
237 };
238
239 bool operator!=(const MagnetizedTovStar& lhs,
240  const MagnetizedTovStar& rhs) noexcept;
241 } // namespace grmhd::AnalyticData
grmhd::AnalyticData
Holds classes implementing analytic data for the GrMhd system.
Definition: AnalyticData.hpp:33
Frame::Inertial
Definition: IndexType.hpp:44
Options.hpp
RelativisticEuler::Solutions::TovStar
A static spherically symmetric star.
Definition: TovStar.hpp:47
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
cstddef
grmhd::AnalyticData::MagnetizedTovStar::PressureExponent
Definition: MagnetizedTovStar.hpp:157
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
grmhd::AnalyticData::MagnetizedTovStar::variables
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, 3 > &x, tmpl::list< Tags... >) const noexcept
Retrieve a collection of variables at (x)
Definition: MagnetizedTovStar.hpp:211
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
limits
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
hydro::Tags::MagneticField
The magnetic field measured by an Eulerian observer, where is the normal to the spatial hypersurfac...
Definition: Tags.hpp:80
grmhd::AnalyticData::MagnetizedTovStar::VectorPotentialAmplitude
Definition: MagnetizedTovStar.hpp:163
std::numeric_limits::max
T max(T... args)
grmhd::AnalyticData::MagnetizedTovStar::CutoffPressureFraction
Definition: MagnetizedTovStar.hpp:171
TMPL.hpp
grmhd::AnalyticData::MagnetizedTovStar
Magnetized TOV star initial data, where metric terms only account for the hydrodynamics not the magne...
Definition: MagnetizedTovStar.hpp:149