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