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