RiemannProblem.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <limits>
9 
10 #include "DataStructures/DataVector.hpp"
13 #include "Evolution/Systems/NewtonianEuler/Sources/NoSource.hpp"
14 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
15 #include "Options/Options.hpp"
16 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
17 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
18 #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
19 #include "Utilities/MakeArray.hpp"
20 #include "Utilities/TMPL.hpp"
22 
23 /// \cond
24 namespace PUP {
25 class er; // IWYU pragma: keep
26 } // namespace PUP
27 /// \endcond
28 
29 namespace NewtonianEuler {
30 namespace Solutions {
31 
32 /*!
33  * \brief Analytic solution to the Riemann Problem
34  *
35  * This class implements the exact Riemann solver described in detail in
36  * Chapter 4 of \cite Toro2009. We follow the notation there.
37  * The algorithm implemented here allows for 1, 2 and 3D wave propagation
38  * along any coordinate axis. Typical initial data for test cases
39  * (see \cite Toro2009) include:
40  *
41  * - Sod's Shock Tube (shock on the right, rarefaction on the left):
42  * - \f$(\rho_L, u_L, p_L) = (1.0, 0.0, 1.0)\f$
43  * - \f$(\rho_R, u_R, p_R) = (0.125, 0.0, 0.1)\f$
44  * - Recommended setup for sample run:
45  * - InitialTimeStep: 0.0001
46  * - Final time: 0.2
47  * - DomainCreator along wave propagation (no AMR):
48  * - Interval of length 1
49  * - InitialRefinement: 8
50  * - InitialGridPoints: 2
51  *
52  * - "123" problem (two symmetric rarefaction waves):
53  * - \f$(\rho_L, u_L, p_L) = (1.0, -2.0, 0.4)\f$
54  * - \f$(\rho_R, u_R, p_R) = (1.0, 2.0, 0.4)\f$
55  * - Recommended setup for sample run:
56  * - InitialTimeStep: 0.0001
57  * - Final time: 0.15
58  * - DomainCreator along wave propagation (no AMR):
59  * - Interval of length 1
60  * - InitialRefinement: 8
61  * - InitialGridPoints: 2
62  *
63  * - Collision of two blast waves (this test is challenging):
64  * - \f$(\rho_L, u_L, p_L) = (5.99924, 19.5975, 460.894)\f$
65  * - \f$(\rho_R, u_R, p_R) = (5.99242, -6.19633, 46.0950)\f$
66  * - Recommended setup for sample run:
67  * - InitialTimeStep: 0.00001
68  * - Final time: 0.012
69  * - DomainCreator along wave propagation (no AMR):
70  * - Interval of length 1
71  * - InitialRefinement: 8
72  * - InitialGridPoints: 2
73  *
74  * - Lax problem:
75  * - \f$(\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)\f$
76  * - \f$(\rho_R, u_R, p_R) = (0.5, 0.0, 0.571)\f$
77  * - Recommended setup for sample run:
78  * - InitialTimeStep: 0.00001
79  * - Final time: 0.1
80  * - DomainCreator along wave propagation (no AMR):
81  * - Interval of length 1
82  * - InitialRefinement: 8
83  * - InitialGridPoints: 2
84  *
85  * where \f$\rho\f$ is the mass density, \f$p\f$ is the pressure, and
86  * \f$u\f$ denotes the normal velocity.
87  *
88  * \note Currently the propagation axis must be hard-coded as a `size_t`
89  * private member variable `propagation_axis_`, which can take one of
90  * the three values `PropagationAxis::X`, `PropagationAxis::Y`, and
91  * `PropagationAxis::Z`.
92  *
93  * \details The algorithm makes use of the following recipe:
94  *
95  * - Given the initial data on both sides of the initial interface of the
96  * discontinuity (here called "left" and "right" sides, where a coordinate
97  * axis points from left to right), we compute the pressure,
98  * \f$p_*\f$, and the normal velocity,
99  * \f$u_*\f$, in the so-called star region. This is done in the constructor.
100  * Here "normal" refers to the normal direction to the initial interface.
101  *
102  * - Given the pressure and the normal velocity in the star region, two
103  * `Wave` `struct`s are created, which represent the waves propagating
104  * at later times on each side of the contact discontinuity. Each `Wave`
105  * is equipped with two `struct`s named `Shock` and `Rarefaction` which
106  * contain functions that compute the primitive variables depending on whether
107  * the wave is a shock or a rarefaction.
108  *
109  * - If \f$p_* > p_K\f$, the wave is a shock, otherwise the wave is a
110  * rarefaction. Here \f$K\f$ stands for \f$L\f$ or \f$R\f$: the left
111  * and right initial pressure, respectively. Since this comparison can't be
112  * performed at compile time, each `Wave` holds a `bool` member `is_shock_`
113  * which is `true` if it is a shock, and `false` if it is a
114  * rarefaction wave. This variable is used to evaluate the correct functions
115  * at run time.
116  *
117  * - In order to obtain the primitives at a certain time and spatial location,
118  * we evaluate whether the spatial location is on the left of the propagating
119  * contact discontinuity \f$(x < u_* t)\f$ or on the right \f$(x > u_* t)\f$,
120  * and we use the corresponding functions for left or right `Wave`s,
121  * respectively.
122  *
123  * \note The characterization of each propagating wave will only
124  * depend on the normal velocity, while the initial jump in the components of
125  * the velocity transverse to the wave propagation will be advected at the
126  * speed of the contact discontinuity (\f$u_*\f$).
127  */
128 template <size_t Dim>
129 class RiemannProblem : public MarkAsAnalyticSolution {
130  enum class Side { Left, Right };
131  enum PropagationAxis { X = 0, Y = 1, Z = 2 };
132 
133  struct Wave;
134  struct Shock;
135  struct Rarefaction;
136 
137  public:
138  using equation_of_state_type = EquationsOfState::IdealFluid<false>;
139  using source_term_type = Sources::NoSource;
140 
141  /// The adiabatic index of the fluid.
142  struct AdiabaticIndex {
143  using type = double;
144  static constexpr OptionString help = {"The adiabatic index of the fluid."};
145  };
146 
147  /// Initial position of the discontinuity
149  using type = double;
150  static constexpr OptionString help = {
151  "The initial position of the discontinuity."};
152  };
153 
154  /// The mass density on the left of the initial discontinuity
156  using type = double;
157  static constexpr OptionString help = {"The left mass density."};
158  };
159 
160  /// The velocity on the left of the initial discontinuity
161  struct LeftVelocity {
163  static constexpr OptionString help = {"The left velocity."};
164  };
165 
166  /// The pressure on the left of the initial discontinuity
167  struct LeftPressure {
168  using type = double;
169  static constexpr OptionString help = {"The left pressure."};
170  };
171 
172  /// The mass density on the right of the initial discontinuity
174  using type = double;
175  static constexpr OptionString help = {"The right mass density."};
176  };
177 
178  /// The velocity on the right of the initial discontinuity
179  struct RightVelocity {
181  static constexpr OptionString help = {"The right velocity."};
182  };
183 
184  /// The pressure on the right of the initial discontinuity
185  struct RightPressure {
186  using type = double;
187  static constexpr OptionString help = {"The right pressure."};
188  };
189 
190  /// The tolerance for solving for \f$p_*\f$.
192  using type = double;
193  static constexpr OptionString help = {
194  "The tolerance for the numerical solution for p star"};
195  static type default_value() noexcept { return 1.e-9; }
196  };
197 
198  // Any of the two states that constitute the initial data, including
199  // some derived quantities that are used repeatedly at each evaluation of the
200  // different waves of the solution.
201  /// Holds initial data on a side of the discontinuity and related quantities
202  struct InitialData {
203  InitialData() = default;
204  InitialData(const InitialData& /*rhs*/) = default;
205  InitialData& operator=(const InitialData& /*rhs*/) = default;
206  InitialData(InitialData&& /*rhs*/) noexcept = default;
207  InitialData& operator=(InitialData&& /*rhs*/) noexcept = default;
208  ~InitialData() = default;
209 
210  InitialData(double mass_density, std::array<double, Dim> velocity,
211  double pressure, double adiabatic_index,
212  size_t propagation_axis) noexcept;
213 
214  // clang-tidy: no runtime references
215  void pup(PUP::er& /*p*/) noexcept; // NOLINT
216 
217  double mass_density_ = std::numeric_limits<double>::signaling_NaN();
218  std::array<double, Dim> velocity_ =
220  double pressure_ = std::numeric_limits<double>::signaling_NaN();
221  double sound_speed_ = std::numeric_limits<double>::signaling_NaN();
222  double normal_velocity_ = std::numeric_limits<double>::signaling_NaN();
223 
224  // Data-dependent constants A and B in Eqns. (4.8) of Toro.
225  double constant_a_ = std::numeric_limits<double>::signaling_NaN();
226  double constant_b_ = std::numeric_limits<double>::signaling_NaN();
227 
228  friend bool operator==(const InitialData& lhs,
229  const InitialData& rhs) noexcept {
230  return lhs.mass_density_ == rhs.mass_density_ and
231  lhs.velocity_ == rhs.velocity_ and lhs.pressure_ == rhs.pressure_;
232  }
233  };
234 
235  using options = tmpl::list<AdiabaticIndex, InitialPosition, LeftMassDensity,
238 
239  static constexpr OptionString help = {
240  "Riemann Problem in 1, 2 or 3D along any coordinate axis."};
241 
242  RiemannProblem() = default;
243  RiemannProblem(const RiemannProblem& /*rhs*/) = delete;
244  RiemannProblem& operator=(const RiemannProblem& /*rhs*/) = delete;
245  RiemannProblem(RiemannProblem&& /*rhs*/) noexcept = default;
246  RiemannProblem& operator=(RiemannProblem&& /*rhs*/) noexcept = default;
247  ~RiemannProblem() = default;
248 
250  double adiabatic_index, double initial_position, double left_mass_density,
251  std::array<double, Dim> left_velocity, double left_pressure,
252  double right_mass_density, std::array<double, Dim> right_velocity,
253  double right_pressure,
254  double pressure_star_tol = PressureStarTol::default_value()) noexcept;
255 
256  /// Retrieve a collection of hydrodynamic variables at position `x`
257  /// and time `t`
258  template <typename DataType, typename... Tags>
260  const tnsr::I<DataType, Dim, Frame::Inertial>& x, double t,
261  tmpl::list<Tags...> /*meta*/) const noexcept {
262  Wave left(left_initial_data_, pressure_star_, velocity_star_,
263  adiabatic_index_, Side::Left);
264  Wave right(right_initial_data_, pressure_star_, velocity_star_,
265  adiabatic_index_, Side::Right);
266 
267  tnsr::I<DataType, Dim, Frame::Inertial> x_shifted(x);
268  x_shifted.get(propagation_axis_) -= initial_position_;
269 
270  return {tuples::get<Tags>(
271  variables(x_shifted, t, tmpl::list<Tags>{}, left, right))...};
272  }
273 
274  const EquationsOfState::IdealFluid<false>& equation_of_state() const
275  noexcept {
276  return equation_of_state_;
277  }
278 
279  // clang-tidy: no runtime references
280  void pup(PUP::er& /*p*/) noexcept; // NOLINT
281 
282  // Retrieve these member variables for testing purposes.
283  constexpr std::array<double, 2> diagnostic_star_region_values() const
284  noexcept {
285  return make_array(pressure_star_, velocity_star_);
286  }
287 
288  private:
289  // @{
290  /// Retrieve hydro variable at `(x, t)`
291  template <typename DataType>
292  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
293  double t, tmpl::list<Tags::MassDensity<DataType>> /*meta*/,
294  const Wave& left, const Wave& right) const noexcept
296 
297  template <typename DataType>
298  auto variables(
299  const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, double t,
301  const Wave& left, const Wave& right) const noexcept
303 
304  template <typename DataType>
305  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
306  double t, tmpl::list<Tags::Pressure<DataType>> /*meta*/,
307  const Wave& left, const Wave& right) const noexcept
309 
310  template <typename DataType>
311  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
312  double t,
313  tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/,
314  const Wave& left, const Wave& right) const noexcept
316  // @}
317 
318  // Any of the two waves propagating on each side of the contact discontinuity.
319  // Depending on whether p_* is larger or smaller
320  // than the initial pressure on the corresponding side, the wave is a
321  // shock (larger) or a rarefaction (smaller) wave. Since p_*
322  // is computed at run time, the characterization of the wave
323  // must also be done at run time. Shock and rarefaction waves will provide
324  // different functions to retrieve the primitive variables at a given (x, t).
325  // Here normal velocity means velocity along the wave propagation.
326  struct Wave {
327  Wave(const InitialData& data, double pressure_star, double velocity_star,
328  double adiabatic_index, const Side& side) noexcept;
329 
330  double mass_density(double x_shifted, double t) const noexcept;
331 
332  double normal_velocity(double x_shifted, double t,
333  double velocity_star) const noexcept;
334 
335  double pressure(double x_shifted, double t, double pressure_star) const
336  noexcept;
337 
338  private:
339  // p_* over initial pressure on the corresponding side
340  double pressure_ratio_ = std::numeric_limits<double>::signaling_NaN();
341  // false if rarefaction wave
342  bool is_shock_ = true;
343 
344  InitialData data_{};
345  Shock shock_{};
346  Rarefaction rarefaction_{};
347  };
348 
349  struct Shock {
350  Shock(const InitialData& data, double pressure_ratio,
351  double adiabatic_index, const Side& side) noexcept;
352 
353  double mass_density(double x_shifted, double t,
354  const InitialData& data) const noexcept;
355  double normal_velocity(double x_shifted, double t, const InitialData& data,
356  double velocity_star) const noexcept;
357 
358  double pressure(double x_shifted, double t, const InitialData& data,
359  double pressure_star) const noexcept;
360 
361  private:
362  double direction_ = std::numeric_limits<double>::signaling_NaN();
363  double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
364  double shock_speed_ = std::numeric_limits<double>::signaling_NaN();
365  };
366 
367  struct Rarefaction {
368  Rarefaction(const InitialData& data, double pressure_ratio,
369  double velocity_star, double adiabatic_index,
370  const Side& side) noexcept;
371 
372  double mass_density(double x_shifted, double t,
373  const InitialData& data) const noexcept;
374  double normal_velocity(double x_shifted, double t, const InitialData& data,
375  double velocity_star) const noexcept;
376 
377  double pressure(double x_shifted, double t, const InitialData& data,
378  double pressure_star) const noexcept;
379 
380  private:
381  double direction_ = std::numeric_limits<double>::signaling_NaN();
382  double gamma_mm_ = std::numeric_limits<double>::signaling_NaN();
383  double gamma_pp_ = std::numeric_limits<double>::signaling_NaN();
384  double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
385  double sound_speed_star_ = std::numeric_limits<double>::signaling_NaN();
386  double head_speed_ = std::numeric_limits<double>::signaling_NaN();
387  double tail_speed_ = std::numeric_limits<double>::signaling_NaN();
388  };
389 
390  template <size_t SpatialDim>
391  friend bool
392  operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
393  const RiemannProblem<SpatialDim>& lhs,
394  const RiemannProblem<SpatialDim>& rhs) noexcept;
395 
396  double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
397  double initial_position_ = std::numeric_limits<double>::signaling_NaN();
398  size_t propagation_axis_ = PropagationAxis::X;
399  InitialData left_initial_data_{};
400  InitialData right_initial_data_{};
401 
402  double pressure_star_tol_ = std::numeric_limits<double>::signaling_NaN();
403  // the pressure in the star region, p_*
404  double pressure_star_ = std::numeric_limits<double>::signaling_NaN();
405  // the velocity in the star region, u_*
406  double velocity_star_ = std::numeric_limits<double>::signaling_NaN();
407 
408  EquationsOfState::IdealFluid<false> equation_of_state_{};
409 };
410 
411 template <size_t Dim>
412 bool operator!=(const RiemannProblem<Dim>& lhs,
413  const RiemannProblem<Dim>& rhs) noexcept;
414 
415 } // namespace Solutions
416 } // namespace NewtonianEuler
The mass density of the fluid.
Definition: Tags.hpp:29
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
T signaling_NaN(T... args)
Analytic solution to the Riemann Problem.
Definition: EvolveNewtonianEulerFwd.hpp:11
The mass density on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:173
Defines function make_array.
The adiabatic index of the fluid.
Definition: RiemannProblem.hpp:142
Holds initial data on a side of the discontinuity and related quantities.
Definition: RiemannProblem.hpp:202
The macroscopic or flow velocity of the fluid.
Definition: Tags.hpp:59
Defines classes and functions for making classes creatable from input files.
Side
A label for the side of a manifold.
Definition: Side.hpp:17
The velocity on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:179
The tolerance for solving for .
Definition: RiemannProblem.hpp:191
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:29
The specific internal energy of the fluid.
Definition: Tags.hpp:68
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:273
constexpr std::array< T, Size > make_array(Args &&... args) noexcept(noexcept(MakeArray_detail::MakeArray< Size==0 >::template apply< T >(std::make_index_sequence<(Size==0 ? Size :Size - 1)>{}, std::forward< Args >(args)...)))
Create a std::array<T, Size>{{T(args...), T(args...), ...}}
Definition: MakeArray.hpp:68
Definition: DataBoxTag.hpp:29
The fluid pressure.
Definition: Tags.hpp:75
Defines classes for Tensor.
Defines a list of useful type aliases for tensors.
Wraps the template metaprogramming library used (brigand)
Initial position of the discontinuity.
Definition: RiemannProblem.hpp:148
The mass density on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:155
The velocity on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:161
tuples::TaggedTuple< Tags... > variables(const tnsr::I< DataType, Dim, Frame::Inertial > &x, double t, tmpl::list< Tags... >) const noexcept
Retrieve a collection of hydrodynamic variables at position x and time t
Definition: RiemannProblem.hpp:259
The pressure on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:167
The pressure on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:185
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8