template<size_t Dim>
class NewtonianEuler::Solutions::RiemannProblem< Dim >
Analytic solution to the Riemann Problem.
This class implements the exact Riemann solver described in detail in Chapter 4 of [107]. We follow the notation there. The algorithm implemented here allows for 1, 2 and 3D wave propagation along any coordinate axis. Typical initial data for test cases (see [107]) include:
- Sod's Shock Tube (shock on the right, rarefaction on the left):
- \((\rho_L, u_L, p_L) = (1.0, 0.0, 1.0)\)
- \((\rho_R, u_R, p_R) = (0.125, 0.0, 0.1)\)
- Recommended setup for sample run:
- InitialTimeStep: 0.0001
- Final time: 0.2
- DomainCreator along wave propagation (no AMR):
- Interval of length 1
- InitialRefinement: 8
- InitialGridPoints: 2
- "123" problem (two symmetric rarefaction waves):
- \((\rho_L, u_L, p_L) = (1.0, -2.0, 0.4)\)
- \((\rho_R, u_R, p_R) = (1.0, 2.0, 0.4)\)
- Recommended setup for sample run:
- InitialTimeStep: 0.0001
- Final time: 0.15
- DomainCreator along wave propagation (no AMR):
- Interval of length 1
- InitialRefinement: 8
- InitialGridPoints: 2
- Collision of two blast waves (this test is challenging):
- \((\rho_L, u_L, p_L) = (5.99924, 19.5975, 460.894)\)
- \((\rho_R, u_R, p_R) = (5.99242, -6.19633, 46.0950)\)
- Recommended setup for sample run:
- InitialTimeStep: 0.00001
- Final time: 0.012
- DomainCreator along wave propagation (no AMR):
- Interval of length 1
- InitialRefinement: 8
- InitialGridPoints: 2
- Lax problem:
- \((\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)\)
- \((\rho_R, u_R, p_R) = (0.5, 0.0, 0.571)\)
- Recommended setup for sample run:
- InitialTimeStep: 0.00001
- Final time: 0.1
- DomainCreator along wave propagation (no AMR):
- Interval of length 1
- InitialRefinement: 8
- InitialGridPoints: 2
where \(\rho\) is the mass density, \(p\) is the pressure, and \(u\) denotes the normal velocity.
- Note
- Currently the propagation axis must be hard-coded as a
size_t
private member variable propagation_axis_
, which can take one of the three values PropagationAxis::X
, PropagationAxis::Y
, and PropagationAxis::Z
.
Details
The algorithm makes use of the following recipe:
- Given the initial data on both sides of the initial interface of the discontinuity (here called "left" and "right" sides, where a coordinate axis points from left to right), we compute the pressure, \(p_*\), and the normal velocity, \(u_*\), in the so-called star region. This is done in the constructor. Here "normal" refers to the normal direction to the initial interface.
- Given the pressure and the normal velocity in the star region, two
Wave
struct
s are created, which represent the waves propagating at later times on each side of the contact discontinuity. Each Wave
is equipped with two struct
s named Shock
and Rarefaction
which contain functions that compute the primitive variables depending on whether the wave is a shock or a rarefaction.
- If \(p_* > p_K\), the wave is a shock, otherwise the wave is a rarefaction. Here \(K\) stands for \(L\) or \(R\): the left and right initial pressure, respectively. Since this comparison can't be performed at compile time, each
Wave
holds a bool
member is_shock_
which is true
if it is a shock, and false
if it is a rarefaction wave. This variable is used to evaluate the correct functions at run time.
- In order to obtain the primitives at a certain time and spatial location, we evaluate whether the spatial location is on the left of the propagating contact discontinuity \((x < u_* t)\) or on the right \((x > u_* t)\), and we use the corresponding functions for left or right
Wave
s, respectively.
- Note
- The characterization of each propagating wave will only depend on the normal velocity, while the initial jump in the components of the velocity transverse to the wave propagation will be advected at the speed of the contact discontinuity ( \(u_*\)).