Question: One of the famous and convenient test cases for shock wave modeling is the 1D Sod's shock tube. This is a Riemann problem for the compressible Euler equations of gas dynamics. The initial set up has been done by dividing domain half with diaphragm. The left and right sides of the domain have initial values chosen to generate a shock wave. The left region has initial pressure and density values higher then those in the right region. Velocities are zero on both sides of the diaphragm.

In my case, I know the initial Mach number of the shock wave and the right hand side initial values. I wish to know how to determine the initial values for the left side when only the Mach number and initial right state are given.

Answer: 
I have [implemented the solution derived below using PyClaw in an IPython notebook](http://nbviewer.ipython.org/gist/ketch/d31fd8d2d7739e59b6c6). If you download that, you can adjust the initial values and see the computed solutions.

General setup
=============

In the solution of a Riemann problem for the 1D Euler equations, there are generally 3 waves. Two of them are genuinely nonlinear; the other is a contact discontinuity, which only carries a jump in the density. My understanding is that you'd like to set up initial left and right states such that the velocity is initially zero everywhere, the rightmost wave is a shock moving at Mach number $M$ relative to the right-state sound speed, and the leftmost wave is a rarefaction (you didn't specify the last requirement, but it would make the problem analogous to the Sod shock tube). This situation is depicted in the figure below, which also establishes the notation I'll use.

![Riemann solution states](https://i.stack.imgur.com/Yeg0Q.png)

Finding the state behind the shock
==================================

Let
$$\mu = \frac{2(M^2-1)}{M(\gamma+1)}.$$

Using the Rankine-Hugoniot conditions, we find that the state just to the right of the contact is given by the following:

\begin{align}
\rho\_r^\*& = \rho\_r \frac{M}{M-\mu} \\
u & = \mu c\_r \\
p & = p\_r \frac{(2M^2 -1)\gamma+1}{\gamma+1}.
\end{align}

If you were willing to have a non-zero velocity for the left state, you could stop here and use the values just derived for the left state. If you want $u\_l=0$, read on.

Finding the left state
======================

Next, we know that the left state and the $[\rho\_l^\*, u, p]$ state must be connected by a left-going rarefaction. The relevant Riemann invariants imply

$$\rho\_l = \frac{4\gamma p\_l}{(\gamma-1)^2 u^2} \left(1- \left(\frac{p}{p\_l}\right)^{\frac{\gamma-1}{2\gamma}}\right)^2$$

and

$$\rho\_l^\* = \left(\frac{p}{p\_l}\right)^{1/\gamma} \rho\_l.$$

So we can prescribe any value we like for $p\_l$ and compute the remaining values $\rho\_l,\rho\_l^\*$ from the equations above. The only catch is that the resulting rarefaction should satisfy the entropy condition, which is true if $p\_l>p$.

Notes
=====

I worked this out using Chapter 14 of [Randall Leveque's text on finite volume methods](http://depts.washington.edu/clawpack/book.html). 

There is a 1-parameter family of possible left states (which seems obvious once you check the number of unknowns and constraints).

If you didn't restrict the left-state velocity to be zero, it would be possible to set up a problem where only the shock wave appears by taking the left state to be $(\rho\_r^\*, u, p)^T$. But since you want zero initial velocities, it's necessary to use the full procedure described above.
If I recall my compressibles, you want to use the [normal shock equations](http://en.wikipedia.org/wiki/Normal_shock_tables). You'll also need to pick a $\gamma$.