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