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"
21 #include "Utilities/TaggedTuple.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 Options::String help = {
145  "The adiabatic index of the fluid."};
146  };
147 
148  /// Initial position of the discontinuity
150  using type = double;
151  static constexpr Options::String help = {
152  "The initial position of the discontinuity."};
153  };
154 
155  /// The mass density on the left of the initial discontinuity
157  using type = double;
158  static constexpr Options::String help = {"The left mass density."};
159  };
160 
161  /// The velocity on the left of the initial discontinuity
162  struct LeftVelocity {
164  static constexpr Options::String help = {"The left velocity."};
165  };
166 
167  /// The pressure on the left of the initial discontinuity
168  struct LeftPressure {
169  using type = double;
170  static constexpr Options::String help = {"The left pressure."};
171  };
172 
173  /// The mass density on the right of the initial discontinuity
175  using type = double;
176  static constexpr Options::String help = {"The right mass density."};
177  };
178 
179  /// The velocity on the right of the initial discontinuity
180  struct RightVelocity {
182  static constexpr Options::String help = {"The right velocity."};
183  };
184 
185  /// The pressure on the right of the initial discontinuity
186  struct RightPressure {
187  using type = double;
188  static constexpr Options::String help = {"The right pressure."};
189  };
190 
191  /// The tolerance for solving for \f$p_*\f$.
193  using type = double;
194  static constexpr Options::String help = {
195  "The tolerance for the numerical solution for p star"};
196  static type suggested_value() noexcept { return 1.e-9; }
197  };
198 
199  // Any of the two states that constitute the initial data, including
200  // some derived quantities that are used repeatedly at each evaluation of the
201  // different waves of the solution.
202  /// Holds initial data on a side of the discontinuity and related quantities
203  struct InitialData {
204  InitialData() = default;
205  InitialData(const InitialData& /*rhs*/) = default;
206  InitialData& operator=(const InitialData& /*rhs*/) = default;
207  InitialData(InitialData&& /*rhs*/) noexcept = default;
208  InitialData& operator=(InitialData&& /*rhs*/) noexcept = default;
209  ~InitialData() = default;
210 
211  InitialData(double mass_density, const std::array<double, Dim>& velocity,
212  double pressure, double adiabatic_index,
213  size_t propagation_axis) noexcept;
214 
215  // clang-tidy: no runtime references
216  void pup(PUP::er& /*p*/) noexcept; // NOLINT
217 
218  double mass_density_ = std::numeric_limits<double>::signaling_NaN();
219  std::array<double, Dim> velocity_ =
221  double pressure_ = std::numeric_limits<double>::signaling_NaN();
222  double sound_speed_ = std::numeric_limits<double>::signaling_NaN();
223  double normal_velocity_ = std::numeric_limits<double>::signaling_NaN();
224 
225  // Data-dependent constants A and B in Eqns. (4.8) of Toro.
226  double constant_a_ = std::numeric_limits<double>::signaling_NaN();
227  double constant_b_ = std::numeric_limits<double>::signaling_NaN();
228 
229  friend bool operator==(const InitialData& lhs,
230  const InitialData& rhs) noexcept {
231  return lhs.mass_density_ == rhs.mass_density_ and
232  lhs.velocity_ == rhs.velocity_ and lhs.pressure_ == rhs.pressure_;
233  }
234  };
235 
236  using options = tmpl::list<AdiabaticIndex, InitialPosition, LeftMassDensity,
239 
240  static constexpr Options::String help = {
241  "Riemann Problem in 1, 2 or 3D along any coordinate axis."};
242 
243  RiemannProblem() = default;
244  RiemannProblem(const RiemannProblem& /*rhs*/) = delete;
245  RiemannProblem& operator=(const RiemannProblem& /*rhs*/) = delete;
246  RiemannProblem(RiemannProblem&& /*rhs*/) noexcept = default;
247  RiemannProblem& operator=(RiemannProblem&& /*rhs*/) noexcept = default;
248  ~RiemannProblem() = default;
249 
251  double adiabatic_index, double initial_position, double left_mass_density,
252  const std::array<double, Dim>& left_velocity, double left_pressure,
253  double right_mass_density, const std::array<double, Dim>& right_velocity,
254  double right_pressure,
255  double pressure_star_tol = PressureStarTol::suggested_value()) noexcept;
256 
257  /// Retrieve a collection of hydrodynamic variables at position `x`
258  /// and time `t`
259  template <typename DataType, typename... Tags>
260  tuples::TaggedTuple<Tags...> variables(
261  const tnsr::I<DataType, Dim, Frame::Inertial>& x, double t,
262  tmpl::list<Tags...> /*meta*/) const noexcept {
263  const Wave left(left_initial_data_, pressure_star_, velocity_star_,
264  adiabatic_index_, Side::Left);
265  const Wave right(right_initial_data_, pressure_star_, velocity_star_,
266  adiabatic_index_, Side::Right);
267 
268  tnsr::I<DataType, Dim, Frame::Inertial> x_shifted(x);
269  x_shifted.get(propagation_axis_) -= initial_position_;
270 
271  return {tuples::get<Tags>(
272  variables(x_shifted, t, tmpl::list<Tags>{}, left, right))...};
273  }
274 
275  const EquationsOfState::IdealFluid<false>& equation_of_state() const
276  noexcept {
277  return equation_of_state_;
278  }
279 
280  // clang-tidy: no runtime references
281  void pup(PUP::er& /*p*/) noexcept; // NOLINT
282 
283  // Retrieve these member variables for testing purposes.
284  constexpr std::array<double, 2> diagnostic_star_region_values() const
285  noexcept {
286  return make_array(pressure_star_, velocity_star_);
287  }
288 
289  private:
290  /// @{
291  /// Retrieve hydro variable at `(x, t)`
292  template <typename DataType>
293  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
294  double t, tmpl::list<Tags::MassDensity<DataType>> /*meta*/,
295  const Wave& left, const Wave& right) const noexcept
297 
298  template <typename DataType>
299  auto variables(
300  const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted, double t,
301  tmpl::list<Tags::Velocity<DataType, Dim>> /*meta*/,
302  const Wave& left, const Wave& right) const noexcept
304 
305  template <typename DataType>
306  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
307  double t, tmpl::list<Tags::Pressure<DataType>> /*meta*/,
308  const Wave& left, const Wave& right) const noexcept
310 
311  template <typename DataType>
312  auto variables(const tnsr::I<DataType, Dim, Frame::Inertial>& x_shifted,
313  double t,
314  tmpl::list<Tags::SpecificInternalEnergy<DataType>> /*meta*/,
315  const Wave& left, const Wave& right) const noexcept
317  /// @}
318 
319  // Any of the two waves propagating on each side of the contact discontinuity.
320  // Depending on whether p_* is larger or smaller
321  // than the initial pressure on the corresponding side, the wave is a
322  // shock (larger) or a rarefaction (smaller) wave. Since p_*
323  // is computed at run time, the characterization of the wave
324  // must also be done at run time. Shock and rarefaction waves will provide
325  // different functions to retrieve the primitive variables at a given (x, t).
326  // Here normal velocity means velocity along the wave propagation.
327  struct Wave {
328  Wave(const InitialData& data, double pressure_star, double velocity_star,
329  double adiabatic_index, const Side& side) noexcept;
330 
331  double mass_density(double x_shifted, double t) const noexcept;
332 
333  double normal_velocity(double x_shifted, double t,
334  double velocity_star) const noexcept;
335 
336  double pressure(double x_shifted, double t, double pressure_star) const
337  noexcept;
338 
339  private:
340  // p_* over initial pressure on the corresponding side
341  double pressure_ratio_ = std::numeric_limits<double>::signaling_NaN();
342  // false if rarefaction wave
343  bool is_shock_ = true;
344 
345  InitialData data_{};
346  Shock shock_{};
347  Rarefaction rarefaction_{};
348  };
349 
350  struct Shock {
351  Shock(const InitialData& data, double pressure_ratio,
352  double adiabatic_index, const Side& side) noexcept;
353 
354  double mass_density(double x_shifted, double t,
355  const InitialData& data) const noexcept;
356  double normal_velocity(double x_shifted, double t, const InitialData& data,
357  double velocity_star) const noexcept;
358 
359  double pressure(double x_shifted, double t, const InitialData& data,
360  double pressure_star) const noexcept;
361 
362  private:
363  double direction_ = std::numeric_limits<double>::signaling_NaN();
364  double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
365  double shock_speed_ = std::numeric_limits<double>::signaling_NaN();
366  };
367 
368  struct Rarefaction {
369  Rarefaction(const InitialData& data, double pressure_ratio,
370  double velocity_star, double adiabatic_index,
371  const Side& side) noexcept;
372 
373  double mass_density(double x_shifted, double t,
374  const InitialData& data) const noexcept;
375  double normal_velocity(double x_shifted, double t, const InitialData& data,
376  double velocity_star) const noexcept;
377 
378  double pressure(double x_shifted, double t, const InitialData& data,
379  double pressure_star) const noexcept;
380 
381  private:
382  double direction_ = std::numeric_limits<double>::signaling_NaN();
383  double gamma_mm_ = std::numeric_limits<double>::signaling_NaN();
384  double gamma_pp_ = std::numeric_limits<double>::signaling_NaN();
385  double mass_density_star_ = std::numeric_limits<double>::signaling_NaN();
386  double sound_speed_star_ = std::numeric_limits<double>::signaling_NaN();
387  double head_speed_ = std::numeric_limits<double>::signaling_NaN();
388  double tail_speed_ = std::numeric_limits<double>::signaling_NaN();
389  };
390 
391  template <size_t SpatialDim>
392  friend bool
393  operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
394  const RiemannProblem<SpatialDim>& lhs,
395  const RiemannProblem<SpatialDim>& rhs) noexcept;
396 
397  double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
398  double initial_position_ = std::numeric_limits<double>::signaling_NaN();
399  size_t propagation_axis_ = PropagationAxis::X;
400  InitialData left_initial_data_{};
401  InitialData right_initial_data_{};
402 
403  double pressure_star_tol_ = std::numeric_limits<double>::signaling_NaN();
404  // the pressure in the star region, p_*
405  double pressure_star_ = std::numeric_limits<double>::signaling_NaN();
406  // the velocity in the star region, u_*
407  double velocity_star_ = std::numeric_limits<double>::signaling_NaN();
408 
409  EquationsOfState::IdealFluid<false> equation_of_state_{};
410 };
411 
412 template <size_t Dim>
413 bool operator!=(const RiemannProblem<Dim>& lhs,
414  const RiemannProblem<Dim>& rhs) noexcept;
415 
416 } // namespace Solutions
417 } // namespace NewtonianEuler
NewtonianEuler::Solutions::RiemannProblem::InitialData
Holds initial data on a side of the discontinuity and related quantities.
Definition: RiemannProblem.hpp:203
NewtonianEuler::Solutions::RiemannProblem::LeftMassDensity
The mass density on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:156
std::rel_ops::operator!=
T operator!=(T... args)
Options.hpp
MakeArray.hpp
NewtonianEuler::Solutions::RiemannProblem::RightVelocity
The velocity on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:180
NewtonianEuler
Items related to evolving the Newtonian Euler system.
Definition: EvolveNewtonianEulerFwd.hpp:8
NewtonianEuler::Solutions::RiemannProblem::RightMassDensity
The mass density on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:174
NewtonianEuler::Solutions::RiemannProblem::AdiabaticIndex
The adiabatic index of the fluid.
Definition: RiemannProblem.hpp:142
make_array
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
cstddef
NewtonianEuler::Solutions::RiemannProblem
Analytic solution to the Riemann Problem.
Definition: EvolveNewtonianEulerFwd.hpp:19
EquationsOfState::IdealFluid< false >
array
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
NewtonianEuler::Solutions::RiemannProblem::variables
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:260
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
Side
Side
Definition: Side.hpp:17
NewtonianEuler::Solutions::RiemannProblem::InitialPosition
Initial position of the discontinuity.
Definition: RiemannProblem.hpp:149
limits
NewtonianEuler::Solutions::RiemannProblem::LeftPressure
The pressure on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:168
NewtonianEuler::Solutions::RiemannProblem::PressureStarTol
The tolerance for solving for .
Definition: RiemannProblem.hpp:192
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
Tensor.hpp
NewtonianEuler::Solutions::RiemannProblem::RightPressure
The pressure on the right of the initial discontinuity.
Definition: RiemannProblem.hpp:186
NewtonianEuler::Solutions::RiemannProblem::LeftVelocity
The velocity on the left of the initial discontinuity.
Definition: RiemannProblem.hpp:162
std::numeric_limits
TMPL.hpp