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/AnalyticData/AnalyticData.hpp"
13 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
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 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::AnalyticData {
29 : /*!
30 : * \brief Initial conditions for relativistic MHD Riemann problems
31 : *
32 : * The standard problems were first collected in \cite Balsara2001. A complete
33 : * Riemann solver for the RMHD equations is presented in \cite Giacomazzo2006
34 : * and can be downloaded from https://www.brunogiacomazzo.org/?page_id=395
35 : *
36 : * The domain is from \f$[-0.5,0.5]^3\f$ with periodic boundary conditions in
37 : * \f$y\f$ and \f$z\f$. The problems are usually run at a resolution of 1600
38 : * finite difference/finite volume grid points and the initial discontinuity is
39 : * located at \f$x=0\f$.
40 : *
41 : * While the problems were originally run in Minkowski space, changing the lapse
42 : * \f$\alpha\f$ to a non-unity constant value and/or the shift \f$\beta^x\f$ to
43 : * a non-zero constant value allows testing some of the metric terms in the
44 : * evolution equations in a fairly simple setup.
45 : *
46 : * Below are the initial conditions for the 5 different Balsara Riemann
47 : * problems. Please note that RP5 has a different final time than the rest, and
48 : * that RP1 has a different adiabatic index than the rest.
49 : *
50 : * RP1:
51 : * - AdiabaticIndex: 2.0
52 : * - LeftDensity: 1.0
53 : * - LeftPressure: 1.0
54 : * - LeftVelocity: [0.0, 0.0, 0.0]
55 : * - LeftMagneticField: [0.5, 1.0, 0.0]
56 : * - RightDensity: 0.125
57 : * - RightPressure: 0.1
58 : * - RightVelocity: [0.0, 0.0, 0.0]
59 : * - RightMagneticField: [0.5, -1.0, 0.0]
60 : * - Lapse: 1.0
61 : * - ShiftX: 0.0
62 : * - Final time: 0.4
63 : *
64 : * RP2:
65 : * - AdiabaticIndex: 1.66666666666666666
66 : * - LeftDensity: 1.0
67 : * - LeftPressure: 30.0
68 : * - LeftVelocity: [0.0, 0.0, 0.0]
69 : * - LeftMagneticField: [5.0, 6.0, 6.0]
70 : * - RightDensity: 1.0
71 : * - RightPressure: 1.0
72 : * - RightVelocity: [0.0, 0.0, 0.0]
73 : * - RightMagneticField: [5.0, 0.7, 0.7]
74 : * - Lapse: 1.0
75 : * - ShiftX: 0.0
76 : * - Final time: 0.4
77 : *
78 : * RP3:
79 : * - AdiabaticIndex: 1.66666666666666666
80 : * - LeftDensity: 1.0
81 : * - LeftPressure: 1000.0
82 : * - LeftVelocity: [0.0, 0.0, 0.0]
83 : * - LeftMagneticField: [10.0, 7.0, 7.0]
84 : * - RightDensity: 1.0
85 : * - RightPressure: 0.1
86 : * - RightVelocity: [0.0, 0.0, 0.0]
87 : * - RightMagneticField: [10.0, 0.7, 0.7]
88 : * - Lapse: 1.0
89 : * - ShiftX: 0.0
90 : * - Final time: 0.4
91 : *
92 : * RP4:
93 : * - AdiabaticIndex: 1.66666666666666666
94 : * - LeftDensity: 1.0
95 : * - LeftPressure: 0.1
96 : * - LeftVelocity: [0.999, 0.0, 0.0]
97 : * - LeftMagneticField: [10.0, 7.0, 7.0]
98 : * - RightDensity: 1.0
99 : * - RightPressure: 0.1
100 : * - RightVelocity: [-0.999, 0.0, 0.0]
101 : * - RightMagneticField: [10.0, -7.0, -7.0]
102 : * - Lapse: 1.0
103 : * - ShiftX: 0.0
104 : * - Final time: 0.4
105 : *
106 : * RP5:
107 : * - AdiabaticIndex: 1.66666666666666666
108 : * - LeftDensity: 1.08
109 : * - LeftPressure: 0.95
110 : * - LeftVelocity: [0.4, 0.3, 0.2]
111 : * - LeftMagneticField: [2.0, 0.3, 0.3]
112 : * - RightDensity: 1.0
113 : * - RightPressure: 1.0
114 : * - RightVelocity: [-0.45, -0.2, 0.2]
115 : * - RightMagneticField: [2.0, -0.7, 0.5]
116 : * - Lapse: 1.0
117 : * - ShiftX: 0.0
118 : * - Final time: 0.55
119 : */
120 1 : class RiemannProblem : public evolution::initial_data::InitialData,
121 : public AnalyticDataBase,
122 : public hydro::TemperatureInitialization<RiemannProblem>,
123 : public MarkAsAnalyticData {
124 : public:
125 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
126 :
127 0 : struct AdiabaticIndex {
128 0 : using type = double;
129 0 : static constexpr Options::String help = {
130 : "The adiabatic index of the ideal fluid"};
131 0 : static type lower_bound() { return 1.0; }
132 : };
133 0 : struct LeftRestMassDensity {
134 0 : using type = double;
135 0 : static std::string name() { return "LeftDensity"; };
136 0 : static constexpr Options::String help = {
137 : "Fluid rest mass density in the left half-domain"};
138 0 : static type lower_bound() { return 0.0; }
139 : };
140 0 : struct RightRestMassDensity {
141 0 : using type = double;
142 0 : static std::string name() { return "RightDensity"; };
143 0 : static constexpr Options::String help = {
144 : "Fluid rest mass density in the right half-domain"};
145 0 : static type lower_bound() { return 0.0; }
146 : };
147 0 : struct LeftPressure {
148 0 : using type = double;
149 0 : static constexpr Options::String help = {
150 : "Fluid pressure in the left half-domain"};
151 0 : static type lower_bound() { return 0.0; }
152 : };
153 0 : struct RightPressure {
154 0 : using type = double;
155 0 : static constexpr Options::String help = {
156 : "Fluid pressure in the right half-domain"};
157 0 : static type lower_bound() { return 0.0; }
158 : };
159 0 : struct LeftSpatialVelocity {
160 0 : using type = std::array<double, 3>;
161 0 : static std::string name() { return "LeftVelocity"; };
162 0 : static constexpr Options::String help = {
163 : "Fluid spatial velocity in the left half-domain"};
164 : };
165 0 : struct RightSpatialVelocity {
166 0 : using type = std::array<double, 3>;
167 0 : static std::string name() { return "RightVelocity"; };
168 0 : static constexpr Options::String help = {
169 : "Fluid spatial velocity in the right half-domain"};
170 : };
171 0 : struct LeftMagneticField {
172 0 : using type = std::array<double, 3>;
173 0 : static constexpr Options::String help = {
174 : "Magnetic field in the left half-domain"};
175 : };
176 0 : struct RightMagneticField {
177 0 : using type = std::array<double, 3>;
178 0 : static constexpr Options::String help = {
179 : "Magnetic field in the right half-domain"};
180 : };
181 0 : struct Lapse {
182 0 : using type = double;
183 0 : static constexpr Options::String help = {
184 : "The value of the lapse. Standard is 1."};
185 0 : static type lower_bound() { return 0.0; }
186 : };
187 0 : struct ShiftX {
188 0 : using type = double;
189 0 : static constexpr Options::String help = {
190 : "The value of the x-component of the shift, beta^x. Standard is 0."};
191 : };
192 :
193 0 : using options =
194 : tmpl::list<AdiabaticIndex, LeftRestMassDensity, RightRestMassDensity,
195 : LeftPressure, RightPressure, LeftSpatialVelocity,
196 : RightSpatialVelocity, LeftMagneticField, RightMagneticField,
197 : Lapse, ShiftX>;
198 :
199 0 : static constexpr Options::String help = {
200 : "Analytic initial data for a GRMHD Riemann problem. The fluid variables "
201 : "are set homogeneously on either half of the domain left and right of "
202 : "x=0."};
203 :
204 0 : RiemannProblem() = default;
205 0 : RiemannProblem(const RiemannProblem& /*rhs*/) = default;
206 0 : RiemannProblem& operator=(const RiemannProblem& /*rhs*/) = default;
207 0 : RiemannProblem(RiemannProblem&& /*rhs*/) = default;
208 0 : RiemannProblem& operator=(RiemannProblem&& /*rhs*/) = default;
209 0 : ~RiemannProblem() override = default;
210 :
211 0 : RiemannProblem(double adiabatic_index, double left_rest_mass_density,
212 : double right_rest_mass_density, double left_pressure,
213 : double right_pressure,
214 : const std::array<double, 3>& left_spatial_velocity,
215 : const std::array<double, 3>& right_spatial_velocity,
216 : const std::array<double, 3>& left_magnetic_field,
217 : const std::array<double, 3>& right_magnetic_field,
218 : double lapse, double shift);
219 :
220 0 : auto get_clone() const
221 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
222 :
223 : /// \cond
224 : explicit RiemannProblem(CkMigrateMessage* msg);
225 : using PUP::able::register_constructor;
226 : WRAPPED_PUPable_decl_template(RiemannProblem);
227 : /// \endcond
228 :
229 : /// @{
230 : /// Retrieve the GRMHD variables at a given position.
231 : template <typename DataType>
232 1 : auto variables(const tnsr::I<DataType, 3>& x,
233 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
234 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
235 :
236 : template <typename DataType>
237 1 : auto variables(const tnsr::I<DataType, 3>& x,
238 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
239 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
240 :
241 : template <typename DataType>
242 1 : auto variables(
243 : const tnsr::I<DataType, 3>& x,
244 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
245 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
246 :
247 : template <typename DataType>
248 1 : auto variables(const tnsr::I<DataType, 3>& x,
249 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
250 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
251 :
252 : template <typename DataType>
253 1 : auto variables(const tnsr::I<DataType, 3>& x,
254 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
255 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
256 :
257 : template <typename DataType>
258 1 : auto variables(const tnsr::I<DataType, 3>& x,
259 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
260 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
261 :
262 : template <typename DataType>
263 1 : auto variables(
264 : const tnsr::I<DataType, 3>& x,
265 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
266 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
267 :
268 : template <typename DataType>
269 1 : auto variables(const tnsr::I<DataType, 3>& x,
270 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
271 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
272 :
273 : template <typename DataType>
274 1 : auto variables(const tnsr::I<DataType, 3>& x,
275 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
276 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
277 :
278 : template <typename DataType>
279 1 : auto variables(const tnsr::I<DataType, 3>& x,
280 : tmpl::list<gr::Tags::Lapse<DataType>> /*meta*/) const
281 : -> tuples::TaggedTuple<gr::Tags::Lapse<DataType>>;
282 :
283 : template <typename DataType>
284 1 : auto variables(const tnsr::I<DataType, 3>& x,
285 : tmpl::list<gr::Tags::Shift<DataType, 3>> /*meta*/) const
286 : -> tuples::TaggedTuple<gr::Tags::Shift<DataType, 3>>;
287 :
288 : template <typename DataType>
289 1 : auto variables(const tnsr::I<DataType, 3>& x,
290 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
291 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
292 : return TemperatureInitialization::variables(
293 : x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
294 : }
295 : /// @}
296 :
297 : /// Retrieve a collection of hydrodynamic variables at position x
298 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
299 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
300 : const tnsr::I<DataType, 3>& x,
301 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
302 : return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
303 : tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
304 : tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
305 : }
306 :
307 : /// Retrieve the metric variables
308 : template <typename DataType, typename Tag>
309 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
310 : tmpl::list<Tag> /*meta*/) const {
311 : constexpr double dummy_time = 0.0;
312 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
313 : }
314 :
315 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
316 : return equation_of_state_;
317 : }
318 :
319 : // NOLINTNEXTLINE(google-runtime-references)
320 0 : void pup(PUP::er& /*p*/) override;
321 :
322 : private:
323 0 : friend bool operator==(const RiemannProblem& lhs, const RiemannProblem& rhs);
324 :
325 0 : friend bool operator!=(const RiemannProblem& lhs, const RiemannProblem& rhs);
326 :
327 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
328 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
329 :
330 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
331 0 : double left_rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
332 0 : double right_rest_mass_density_ =
333 : std::numeric_limits<double>::signaling_NaN();
334 0 : double left_pressure_ = std::numeric_limits<double>::signaling_NaN();
335 0 : double right_pressure_ = std::numeric_limits<double>::signaling_NaN();
336 0 : static constexpr double discontinuity_location_ = 0.0;
337 0 : std::array<double, 3> left_spatial_velocity_{
338 : {std::numeric_limits<double>::signaling_NaN(),
339 : std::numeric_limits<double>::signaling_NaN(),
340 : std::numeric_limits<double>::signaling_NaN()}};
341 0 : std::array<double, 3> right_spatial_velocity_{
342 : {std::numeric_limits<double>::signaling_NaN(),
343 : std::numeric_limits<double>::signaling_NaN(),
344 : std::numeric_limits<double>::signaling_NaN()}};
345 0 : std::array<double, 3> left_magnetic_field_{
346 : {std::numeric_limits<double>::signaling_NaN(),
347 : std::numeric_limits<double>::signaling_NaN(),
348 : std::numeric_limits<double>::signaling_NaN()}};
349 0 : std::array<double, 3> right_magnetic_field_{
350 : {std::numeric_limits<double>::signaling_NaN(),
351 : std::numeric_limits<double>::signaling_NaN(),
352 : std::numeric_limits<double>::signaling_NaN()}};
353 0 : double lapse_ = std::numeric_limits<double>::signaling_NaN();
354 0 : double shift_ = std::numeric_limits<double>::signaling_NaN();
355 : };
356 : } // namespace grmhd::AnalyticData
|
