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 : #include <string>
9 :
10 : #include "DataStructures/TaggedTuple.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
15 : #include "PointwiseFunctions/AnalyticSolutions/GrMhd/Solutions.hpp"
16 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
17 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
18 : #include "PointwiseFunctions/Hydro/Temperature.hpp"
19 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
20 : #include "Utilities/Serialization/CharmPupable.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : /// \cond
24 : namespace PUP {
25 : class er;
26 : } // namespace PUP
27 : /// \endcond
28 :
29 1 : namespace grmhd::Solutions {
30 :
31 : /*!
32 : * \brief Circularly polarized Alfvén wave solution in Minkowski
33 : * spacetime travelling along a background magnetic field.
34 : *
35 : * An analytic solution to the 3-D GRMHD system. The user specifies the
36 : * wavenumber \f$k\f$ of the Alfvén wave, the constant pressure
37 : * throughout the fluid \f$P\f$, the constant rest mass density throughout the
38 : * fluid \f$\rho_0\f$, the adiabatic index for the ideal fluid equation of
39 : * state \f$\gamma\f$, the magnetic field parallel to the wavevector
40 : * \f$\vec{B}_0\f$, and the transverse magnetic field vector \f$\vec{B}_1\f$ at
41 : * \f$x=y=z=t=0\f$.
42 : *
43 : * We define the auxiliary velocities:
44 : * \f[v^2_{B0} = \frac{B_0^2}{\rho_0 h + B_0^2 + B_1^2}\f]
45 : * \f[v^2_{B1} = \frac{B_1^2}{\rho_0 h + B_0^2 + B_1^2}\f]
46 : *
47 : * The Alfvén wave phase speed that solves the GRMHD equations, even for
48 : * finite amplitudes \cite DelZanna2007pk, is given by:
49 : *
50 : * \f[v_A^2 = \frac{2v^2_{B0}}{1 + \sqrt{1 - 4 v^2_{B0}v^2_{B1}}}\f]
51 : *
52 : * The amplitude of the fluid velocity is given by:
53 : *
54 : * \f[v_f^2 = \frac{2v^2_{B1}}{1 + \sqrt{1 - 4 v^2_{B0}v^2_{B1}}}\f]
55 : *
56 : * The electromagnetic field vectors define a set of basis vectors:
57 : *
58 : * \f{align*}{
59 : * \hat{b}_0 &= \vec{B_0}/B_0 \\
60 : * \hat{b}_1 &= \vec{B_1}/B_1 \\
61 : * \hat{e} &= \hat{b}_1 \times \hat{b}_0
62 : * \f}
63 : *
64 : * We also define the auxiliary variable for the phase \f$\phi\f$:
65 : * \f[\phi = k(\vec{x}\cdot\hat{b}_0 - v_A t)\f]
66 : * In Cartesian coordinates \f$(x, y, z)\f$, and using
67 : * dimensionless units, the primitive quantities at a given time \f$t\f$ are
68 : * then
69 : *
70 : * \f{align*}
71 : * \rho(\vec{x},t) &= \rho_0 \\
72 : * \vec{v}(\vec{x},t) &= v_f(-\hat{b}_1\cos\phi
73 : * +\hat{e}\sin\phi)\\
74 : * P(\vec{x},t) &= P, \\
75 : * \epsilon(\vec{x}, t) &= \frac{P}{(\gamma - 1)\rho_0}\\
76 : * \vec{B}(\vec{x},t) &= B_1(\hat{b}_1\cos\phi
77 : * -\hat{e}\sin\phi) + \vec{B_0}
78 : * \f}
79 : *
80 : * Note that the phase speed is not the characteristic Alfvén speed
81 : * \f$c_A\f$, which is the speed in the limiting case where the total magnetic
82 : * field is parallel to the direction of propagation \cite DelZanna2007pk :
83 : *
84 : * \f[c_A^2 = \frac{b^2}{\rho_0 h + b^2}\f]
85 : *
86 : * Where \f$b^2\f$ is the invariant quantity \f$B^2 - E^2\f$, given by:
87 : *
88 : * \f[b^2 = B_0^2 + B_1^2 - B_0^2 v_f^2\f]
89 : */
90 1 : class AlfvenWave : public evolution::initial_data::InitialData,
91 : public AnalyticSolution,
92 : public hydro::TemperatureInitialization<AlfvenWave>,
93 : public MarkAsAnalyticSolution {
94 : public:
95 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
96 :
97 : /// The wave number of the profile.
98 1 : struct WaveNumber {
99 0 : using type = double;
100 0 : static constexpr Options::String help = {"The wave number of the profile."};
101 : };
102 :
103 : /// The constant pressure throughout the fluid.
104 1 : struct Pressure {
105 0 : using type = double;
106 0 : static constexpr Options::String help = {
107 : "The constant pressure throughout the fluid."};
108 0 : static type lower_bound() { return 0.0; }
109 : };
110 :
111 : /// The constant rest mass density throughout the fluid.
112 1 : struct RestMassDensity {
113 0 : using type = double;
114 0 : static constexpr Options::String help = {
115 : "The constant rest mass density throughout the fluid."};
116 0 : static type lower_bound() { return 0.0; }
117 : };
118 :
119 : /// The constant electron fraction throughout the fluid.
120 1 : struct ElectronFraction {
121 0 : using type = double;
122 0 : static constexpr Options::String help = {
123 : "The constant electron fraction throughout the fluid."};
124 0 : static type lower_bound() { return 0.0; }
125 0 : static type upper_bound() { return 1.0; }
126 : };
127 :
128 : /// The adiabatic index for the ideal fluid.
129 1 : struct AdiabaticIndex {
130 0 : using type = double;
131 0 : static constexpr Options::String help = {
132 : "The adiabatic index for the ideal fluid."};
133 0 : static type lower_bound() { return 1.0; }
134 : };
135 :
136 : /// The background static magnetic field vector.
137 1 : struct BackgroundMagneticField {
138 0 : using type = std::array<double, 3>;
139 0 : static std::string name() { return "BkgdMagneticField"; }
140 0 : static constexpr Options::String help = {
141 : "The background magnetic field [B0^x, B0^y, B0^z]."};
142 : };
143 :
144 : /// The sinusoidal magnetic field vector associated with
145 : /// the Alfvén wave, perpendicular to the background
146 : /// magnetic field vector.
147 1 : struct WaveMagneticField {
148 0 : using type = std::array<double, 3>;
149 0 : static constexpr Options::String help = {
150 : "The wave magnetic field [B1^x, B1^y, B1^z]."};
151 : };
152 :
153 0 : using options =
154 : tmpl::list<WaveNumber, Pressure, RestMassDensity, ElectronFraction,
155 : AdiabaticIndex, BackgroundMagneticField, WaveMagneticField>;
156 0 : static constexpr Options::String help = {
157 : "Circularly polarized Alfven wave in Minkowski spacetime."};
158 :
159 0 : AlfvenWave() = default;
160 0 : AlfvenWave(const AlfvenWave& /*rhs*/) = default;
161 0 : AlfvenWave& operator=(const AlfvenWave& /*rhs*/) = default;
162 0 : AlfvenWave(AlfvenWave&& /*rhs*/) = default;
163 0 : AlfvenWave& operator=(AlfvenWave&& /*rhs*/) = default;
164 0 : ~AlfvenWave() override = default;
165 :
166 0 : AlfvenWave(double wavenumber, double pressure, double rest_mass_density,
167 : double electron_fraction, double adiabatic_index,
168 : const std::array<double, 3>& background_magnetic_field,
169 : const std::array<double, 3>& wave_magnetic_field);
170 :
171 0 : auto get_clone() const
172 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
173 :
174 : /// \cond
175 : explicit AlfvenWave(CkMigrateMessage* msg);
176 : using PUP::able::register_constructor;
177 : WRAPPED_PUPable_decl_template(AlfvenWave);
178 : /// \endcond
179 :
180 : /// @{
181 : /// Retrieve hydro variable at `(x, t)`
182 : template <typename DataType>
183 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
184 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
185 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
186 :
187 : template <typename DataType>
188 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
189 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
190 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
191 :
192 : template <typename DataType>
193 1 : auto variables(
194 : const tnsr::I<DataType, 3>& x, double t,
195 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
196 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
197 :
198 : template <typename DataType>
199 1 : auto variables(const tnsr::I<DataType, 3>& x, double /*t*/,
200 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
201 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
202 :
203 : template <typename DataType>
204 1 : auto variables(const tnsr::I<DataType, 3>& x, double /*t*/,
205 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
206 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
207 :
208 : template <typename DataType>
209 1 : auto variables(const tnsr::I<DataType, 3>& x, double /*t*/,
210 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
211 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
212 :
213 : template <typename DataType>
214 1 : auto variables(
215 : const tnsr::I<DataType, 3>& x, double /*t*/,
216 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
217 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
218 :
219 : template <typename DataType>
220 1 : auto variables(const tnsr::I<DataType, 3>& x, double /*t*/,
221 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
222 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
223 :
224 : template <typename DataType>
225 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
226 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
227 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
228 :
229 : template <typename DataType>
230 1 : auto variables(const tnsr::I<DataType, 3>& x, double t,
231 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
232 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
233 : return TemperatureInitialization::variables(
234 : x, t, tmpl::list<hydro::Tags::Temperature<DataType>>{});
235 : }
236 : /// @}
237 :
238 : /// Retrieve a collection of hydro variables at `(x, t)`
239 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
240 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
241 : const tnsr::I<DataType, 3>& x, double t,
242 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
243 : return {tuples::get<Tag1>(variables(x, t, tmpl::list<Tag1>{})),
244 : tuples::get<Tag2>(variables(x, t, tmpl::list<Tag2>{})),
245 : tuples::get<Tags>(variables(x, t, tmpl::list<Tags>{}))...};
246 : }
247 :
248 : /// Retrieve the metric variables
249 : template <typename DataType, typename Tag,
250 : Requires<tmpl::list_contains_v<
251 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
252 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x, double t,
253 : tmpl::list<Tag> /*meta*/) const {
254 : return background_spacetime_.variables(x, t, tmpl::list<Tag>{});
255 : }
256 :
257 : // NOLINTNEXTLINE(google-runtime-references)
258 0 : void pup(PUP::er& /*p*/) override;
259 :
260 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
261 : return equation_of_state_;
262 : }
263 :
264 : protected:
265 0 : friend bool operator==(const AlfvenWave& lhs, const AlfvenWave& rhs);
266 :
267 : // Computes the phase.
268 : template <typename DataType>
269 0 : DataType k_dot_x_minus_vt(const tnsr::I<DataType, 3>& x, double t) const;
270 0 : double wavenumber_ = std::numeric_limits<double>::signaling_NaN();
271 0 : double pressure_ = std::numeric_limits<double>::signaling_NaN();
272 0 : double rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
273 0 : double electron_fraction_ = std::numeric_limits<double>::signaling_NaN();
274 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
275 0 : std::array<double, 3> background_magnetic_field_{
276 : {std::numeric_limits<double>::signaling_NaN(),
277 : std::numeric_limits<double>::signaling_NaN(),
278 : std::numeric_limits<double>::signaling_NaN()}};
279 0 : std::array<double, 3> wave_magnetic_field_{
280 : {std::numeric_limits<double>::signaling_NaN(),
281 : std::numeric_limits<double>::signaling_NaN(),
282 : std::numeric_limits<double>::signaling_NaN()}};
283 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
284 0 : tnsr::I<double, 3> initial_unit_vector_along_background_magnetic_field_{};
285 0 : tnsr::I<double, 3> initial_unit_vector_along_wave_magnetic_field_{};
286 0 : tnsr::I<double, 3> initial_unit_vector_along_wave_electric_field_{};
287 0 : double magnitude_B0_ = std::numeric_limits<double>::signaling_NaN();
288 0 : double magnitude_B1_ = std::numeric_limits<double>::signaling_NaN();
289 0 : double magnitude_E_ = std::numeric_limits<double>::signaling_NaN();
290 0 : double alfven_speed_ = std::numeric_limits<double>::signaling_NaN();
291 0 : double fluid_speed_ = std::numeric_limits<double>::signaling_NaN();
292 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
293 : };
294 :
295 0 : bool operator!=(const AlfvenWave& lhs, const AlfvenWave& rhs);
296 :
297 : } // namespace grmhd::Solutions
|
