Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/String.hpp"
12 : #include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
13 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
14 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
15 : #include "Utilities/TMPL.hpp"
16 :
17 : /// \cond
18 : namespace gsl {
19 : template <typename T>
20 : class not_null;
21 : } // namespace gsl
22 : /// \endcond
23 :
24 : namespace grmhd::AnalyticData::InitialMagneticFields {
25 :
26 : /*!
27 : * \brief %Poloidal magnetic field for GRMHD initial data.
28 : *
29 : * The vector potential has the form
30 : *
31 : * \f{align*}{
32 : * A_{\phi} = A_b \varpi^2 \max(p-p_{\mathrm{cut}}, 0)^{n_s} ,
33 : * \f}
34 : *
35 : * where \f$A_b\f$ controls the amplitude of the magnetic field,
36 : * \f$\varpi^2=x^2+y^2=r^2-z^2\f$ is the cylindrical radius,
37 : * \f$n_s\f$ controls the degree of differentiability, and
38 : * \f$p_{\mathrm{cut}}\f$ controls the pressure cutoff below which the magnetic
39 : * field is zero.
40 : *
41 : * In Cartesian coordinates the vector potential is:
42 : *
43 : * \f{align*}{
44 : * A_x & = -\frac{y}{\varpi^2}A_\phi
45 : * = -y A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
46 : * A_y & = \frac{x}{\varpi^2}A_\phi
47 : * = x A_b\max(p-p_{\mathrm{cut}}, 0)^{n_s}, \\
48 : * A_z & = 0,
49 : * \f}
50 : *
51 : * On the region where the field is non-zero, the magnetic field is given by:
52 : *
53 : * \f{align*}{
54 : * B^x & = -\frac{A_b n_s}{\sqrt{\gamma}}
55 : * (p-p_{\mathrm{cut}})^{n_s-1} x \partial_z p \\
56 : * B^x & = -\frac{A_b n_s}{\sqrt{\gamma}}
57 : * (p-p_{\mathrm{cut}})^{n_s-1} y \partial_z p \\
58 : * B^z & = \frac{A_b}{\sqrt{\gamma}}\left[
59 : * 2(p-p_{\mathrm{cut}})^{n_s}
60 : * + n_s (p-p_{\mathrm{cut}})^{n_s-1} (x \partial_x p + y \partial_y p)
61 : * \right]
62 : * \f}
63 : *
64 : * Taking the small-\f$r\f$ limit gives the magnetic field at the origin:
65 : *
66 : * \f{align*}{
67 : * B^x&=0, \\
68 : * B^y&=0, \\
69 : * B^z&=\frac{A_b}{\sqrt{\gamma}}
70 : * 2(p-p_{\mathrm{cut}})^{n_s}.
71 : * \f}
72 : *
73 : * Note that the coordinates are relative to the `Center` passed in, so the
74 : * field can be centered about any arbitrary point. The field is also zero
75 : * outside of `MaxDistanceFromCenter`, so that compact support can be imposed if
76 : * necessary.
77 : *
78 : * \warning This assumes the magnetic field is initialized, both in size and
79 : * value, before being passed into the `variables` function. This is so that
80 : * multiple magnetic fields can be superposed. Each magnetic field
81 : * configuration does a `+=` to make this possible.
82 : */
83 1 : class Poloidal : public InitialMagneticField {
84 : public:
85 0 : struct PressureExponent {
86 0 : using type = size_t;
87 0 : static constexpr Options::String help = {
88 : "The exponent n_s controlling the smoothness of the field"};
89 : };
90 :
91 0 : struct CutoffPressure {
92 0 : using type = double;
93 0 : static constexpr Options::String help = {
94 : "The pressure below which there is no magnetic field."};
95 0 : static type lower_bound() { return 0.0; }
96 : };
97 :
98 0 : struct VectorPotentialAmplitude {
99 0 : using type = double;
100 0 : static constexpr Options::String help = {
101 : "The amplitude A_b of the phi-component of the vector potential. This "
102 : "controls the magnetic field strength."};
103 0 : static type lower_bound() { return 0.0; }
104 : };
105 :
106 0 : struct Center {
107 0 : using type = std::array<double, 3>;
108 0 : static constexpr Options::String help = {
109 : "The center of the magnetic field."};
110 : };
111 :
112 0 : struct MaxDistanceFromCenter {
113 0 : using type = double;
114 0 : static constexpr Options::String help = {
115 : "The maximum distance from the center to compute the magnetic field. "
116 : "Everywhere outside the field is set to zero."};
117 0 : static type lower_bound() { return 0.0; }
118 : };
119 :
120 0 : using options =
121 : tmpl::list<PressureExponent, CutoffPressure, VectorPotentialAmplitude,
122 : Center, MaxDistanceFromCenter>;
123 :
124 0 : static constexpr Options::String help = {"Poloidal initial magnetic field"};
125 :
126 0 : Poloidal() = default;
127 0 : Poloidal(const Poloidal& /*rhs*/) = default;
128 0 : Poloidal& operator=(const Poloidal& /*rhs*/) = default;
129 0 : Poloidal(Poloidal&& /*rhs*/) = default;
130 0 : Poloidal& operator=(Poloidal&& /*rhs*/) = default;
131 0 : ~Poloidal() override = default;
132 :
133 0 : Poloidal(size_t pressure_exponent, double cutoff_pressure,
134 : double vector_potential_amplitude, std::array<double, 3> center,
135 : double max_distance_from_center);
136 :
137 0 : auto get_clone() const -> std::unique_ptr<InitialMagneticField> override;
138 :
139 : /// \cond
140 : explicit Poloidal(CkMigrateMessage* msg);
141 : using PUP::able::register_constructor;
142 : WRAPPED_PUPable_decl_template(Poloidal);
143 : /// \endcond
144 :
145 : // NOLINTNEXTLINE(google-runtime-references)
146 0 : void pup(PUP::er& p) override;
147 :
148 : /// Retrieve magnetic fields at `(x)`
149 1 : void variables(gsl::not_null<tnsr::I<DataVector, 3>*> result,
150 : const tnsr::I<DataVector, 3>& coords,
151 : const Scalar<DataVector>& pressure,
152 : const Scalar<DataVector>& sqrt_det_spatial_metric,
153 : const tnsr::i<DataVector, 3>& deriv_pressure) const override;
154 :
155 : /// Retrieve magnetic fields at `(x)`
156 1 : void variables(gsl::not_null<tnsr::I<double, 3>*> result,
157 : const tnsr::I<double, 3>& coords,
158 : const Scalar<double>& pressure,
159 : const Scalar<double>& sqrt_det_spatial_metric,
160 : const tnsr::i<double, 3>& deriv_pressure) const override;
161 :
162 0 : bool is_equal(const InitialMagneticField& rhs) const override;
163 :
164 : private:
165 : template <typename DataType>
166 0 : void variables_impl(gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
167 : const tnsr::I<DataType, 3>& coords,
168 : const Scalar<DataType>& pressure,
169 : const Scalar<DataType>& sqrt_det_spatial_metric,
170 : const tnsr::i<DataType, 3>& deriv_pressure) const;
171 :
172 0 : size_t pressure_exponent_ = std::numeric_limits<size_t>::max();
173 0 : double cutoff_pressure_ = std::numeric_limits<double>::signaling_NaN();
174 0 : double vector_potential_amplitude_ =
175 : std::numeric_limits<double>::signaling_NaN();
176 0 : std::array<double, 3> center_{{std::numeric_limits<double>::signaling_NaN(),
177 : std::numeric_limits<double>::signaling_NaN(),
178 : std::numeric_limits<double>::signaling_NaN()}};
179 0 : double max_distance_from_center_ =
180 : std::numeric_limits<double>::signaling_NaN();
181 :
182 0 : friend bool operator==(const Poloidal& lhs, const Poloidal& rhs);
183 0 : friend bool operator!=(const Poloidal& lhs, const Poloidal& rhs);
184 : };
185 :
186 : } // namespace grmhd::AnalyticData::InitialMagneticFields
|
