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