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