Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <limits>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/Context.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
14 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
15 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
16 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
17 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
18 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::AnalyticData {
29 :
30 : /*!
31 : * \brief Analytic initial data for a magnetic rotor.
32 : *
33 : * This is a test first described in \cite Balsara1999 for classical MHD and
34 : * later generalised to relativistic MHD in \cite DelZanna2002rv
35 : *
36 : * This effectively 2D test initially consists of an infinitely long cylinder of
37 : * radius `RotorRadius` rotating about the z-axis with the given
38 : * `AngularVelocity`. The rest mass density of the fluid inside the rotor,
39 : * `RotorDensity`, is higher than the `BackgroundDensity` outside of the rotor.
40 : * The fluid is at a constant `Pressure`. The rotor is embedded in a constant
41 : * `MagneticField` (usually taken to be along the x-axis). The fluid is an
42 : * ideal fluid with the given `AdiabaticIndex`. Evolving the initial data,
43 : * magnetic braking will slow down the rotor, while dragging the magnetic field
44 : * lines.
45 : *
46 : * The standard test setup is done on a unit cube \f$[-0.5,0.5]^3\f$ with the
47 : * following values given for the options:
48 : * - RotorRadius: 0.1
49 : * - RotorDensity: 10.0
50 : * - BackgroundDensity: 1.0
51 : * - Pressure: 1.0
52 : * - AngularVelocity: 9.95
53 : * - MagneticField: [1.0, 0.0, 0.0]
54 : * - AdiabaticIndex: 1.66666666666666667
55 : *
56 : * Note that \cite Zanotti2016efficient uses different parameters,
57 : * - RotorRadius: 0.1
58 : * - RotorDensity: 10.0
59 : * - BackgroundDensity: 1.0
60 : * - Pressure: 1.0
61 : * - AngularVelocity: 9.3
62 : * - MagneticField: [1.0, 0.0, 0.0]
63 : * - AdiabaticIndex: 1.333333333333333
64 : *
65 : * The magnetic field in the disk should rotate by about 90 degrees by t = 0.4.
66 : */
67 1 : class MagneticRotor : public evolution::initial_data::InitialData,
68 : public MarkAsAnalyticData,
69 : public AnalyticDataBase,
70 : public hydro::TemperatureInitialization<MagneticRotor> {
71 : public:
72 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
73 :
74 : /// Radius of the rotor.
75 1 : struct RotorRadius {
76 0 : using type = double;
77 0 : static constexpr Options::String help = {
78 : "The initial radius of the rotor."};
79 0 : static type lower_bound() { return 0.0; }
80 : };
81 : /// Density inside the rotor.
82 1 : struct RotorDensity {
83 0 : using type = double;
84 0 : static constexpr Options::String help = {"Density inside RotorRadius."};
85 0 : static type lower_bound() { return 0.0; }
86 : };
87 : /// Density outside the rotor.
88 1 : struct BackgroundDensity {
89 0 : using type = double;
90 0 : static constexpr Options::String help = {"Density outside RotorRadius."};
91 0 : static type lower_bound() { return 0.0; }
92 : };
93 : /// Uniform pressure inside and outside the rotor.
94 1 : struct Pressure {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {"Pressure."};
97 0 : static type lower_bound() { return 0.0; }
98 : };
99 : /// Angular velocity inside the rotor.
100 1 : struct AngularVelocity {
101 0 : using type = double;
102 0 : static constexpr Options::String help = {
103 : "Angular velocity of matter inside RotorRadius"};
104 : };
105 : /// The x,y,z components of the uniform magnetic field threading the matter.
106 1 : struct MagneticField {
107 0 : using type = std::array<double, 3>;
108 0 : static constexpr Options::String help = {
109 : "The x,y,z components of the uniform magnetic field."};
110 : };
111 : /// The adiabatic index of the ideal fluid.
112 1 : struct AdiabaticIndex {
113 0 : using type = double;
114 0 : static constexpr Options::String help = {
115 : "The adiabatic index of the ideal fluid."};
116 0 : static type lower_bound() { return 1.0; }
117 : };
118 :
119 0 : using options =
120 : tmpl::list<RotorRadius, RotorDensity, BackgroundDensity, Pressure,
121 : AngularVelocity, MagneticField, AdiabaticIndex>;
122 :
123 0 : static constexpr Options::String help = {
124 : "Magnetic rotor analytic initial data."};
125 :
126 0 : MagneticRotor() = default;
127 0 : MagneticRotor(const MagneticRotor& /*rhs*/) = default;
128 0 : MagneticRotor& operator=(const MagneticRotor& /*rhs*/) = default;
129 0 : MagneticRotor(MagneticRotor&& /*rhs*/) = default;
130 0 : MagneticRotor& operator=(MagneticRotor&& /*rhs*/) = default;
131 0 : ~MagneticRotor() override = default;
132 :
133 0 : MagneticRotor(double rotor_radius, double rotor_density,
134 : double background_density, double pressure,
135 : double angular_velocity,
136 : const std::array<double, 3>& magnetic_field,
137 : double adiabatic_index, const Options::Context& context = {});
138 :
139 0 : auto get_clone() const
140 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
141 :
142 : /// \cond
143 : explicit MagneticRotor(CkMigrateMessage* msg);
144 : using PUP::able::register_constructor;
145 : WRAPPED_PUPable_decl_template(MagneticRotor);
146 : /// \endcond
147 :
148 : /// @{
149 : /// Retrieve the GRMHD variables at a given position.
150 : template <typename DataType>
151 1 : auto variables(const tnsr::I<DataType, 3>& x,
152 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
153 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
154 :
155 : template <typename DataType>
156 1 : auto variables(const tnsr::I<DataType, 3>& x,
157 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
158 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
159 :
160 : template <typename DataType>
161 1 : auto variables(
162 : const tnsr::I<DataType, 3>& x,
163 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
164 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
165 :
166 : template <typename DataType>
167 1 : auto variables(const tnsr::I<DataType, 3>& x,
168 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
169 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
170 :
171 : template <typename DataType>
172 1 : auto variables(const tnsr::I<DataType, 3>& x,
173 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
174 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
175 :
176 : template <typename DataType>
177 1 : auto variables(const tnsr::I<DataType, 3>& x,
178 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
179 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
180 :
181 : template <typename DataType>
182 1 : auto variables(
183 : const tnsr::I<DataType, 3>& x,
184 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
185 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
186 :
187 : template <typename DataType>
188 1 : auto variables(const tnsr::I<DataType, 3>& x,
189 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
190 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
191 :
192 : template <typename DataType>
193 1 : auto variables(const tnsr::I<DataType, 3>& x,
194 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
195 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
196 :
197 : template <typename DataType>
198 1 : auto variables(const tnsr::I<DataType, 3>& x,
199 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
200 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
201 : return TemperatureInitialization::variables(
202 : x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
203 : }
204 : /// @}
205 :
206 : /// Retrieve a collection of hydrodynamic variables at position x
207 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
208 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
209 : const tnsr::I<DataType, 3>& x,
210 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
211 : return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
212 : tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
213 : tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
214 : }
215 :
216 : /// Retrieve the metric variables
217 : template <typename DataType, typename Tag,
218 : Requires<tmpl::list_contains_v<
219 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
220 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
221 : tmpl::list<Tag> /*meta*/) const {
222 : constexpr double dummy_time = 0.0;
223 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
224 : }
225 :
226 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
227 : return equation_of_state_;
228 : }
229 :
230 : // NOLINTNEXTLINE(google-runtime-references)
231 0 : void pup(PUP::er& /*p*/) override;
232 :
233 : private:
234 0 : double rotor_radius_ = std::numeric_limits<double>::signaling_NaN();
235 0 : double rotor_density_ = std::numeric_limits<double>::signaling_NaN();
236 0 : double background_density_ = std::numeric_limits<double>::signaling_NaN();
237 0 : double pressure_ = std::numeric_limits<double>::signaling_NaN();
238 0 : double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
239 0 : std::array<double, 3> magnetic_field_{
240 : {std::numeric_limits<double>::signaling_NaN(),
241 : std::numeric_limits<double>::signaling_NaN(),
242 : std::numeric_limits<double>::signaling_NaN()}};
243 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
244 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
245 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
246 :
247 0 : friend bool operator==(const MagneticRotor& lhs, const MagneticRotor& rhs);
248 :
249 0 : friend bool operator!=(const MagneticRotor& lhs, const MagneticRotor& rhs);
250 : };
251 :
252 : } // namespace grmhd::AnalyticData
|
