# Compressible gas dynamics

The Euler system of conservation laws for a compressible gas can be written

${\partial }_{t}\left(\begin{array}{c}\hfill \rho \hfill \\ \hfill E\hfill \\ \hfill {w}_{x}\hfill \\ \hfill {w}_{y}\hfill \end{array}\right)+{\nabla }_{x}\cdot \left(\begin{array}{c}\hfill {w}_{x}\hfill \\ \hfill \frac{{w}_{x}}{\rho }\left(E+p\right)\hfill \\ \hfill \frac{{w}_{x}^{2}}{\rho }+p\hfill \\ \hfill \frac{{w}_{y}{w}_{x}}{\rho }\hfill \end{array}\right)+{\nabla }_{y}\cdot \left(\begin{array}{c}\hfill {w}_{y}\hfill \\ \hfill \frac{{w}_{y}}{\rho }\left(E+p\right)\hfill \\ \hfill \frac{{w}_{y}{w}_{x}}{\rho }\hfill \\ \hfill \frac{{w}_{y}^{2}}{\rho }+p\hfill \end{array}\right)=0$

with $\rho$ the gas density, $E$ the total energy, $\mathbf{\text{w}}$ the gas momentum and $p$ the pressure given by the equation of state

$p=\left(\gamma -1\right)\left(E-\rho {\mathbf{\text{u}}}^{2}/2\right)$

with $\gamma$ the polytropic exponent. This system can be solved using the generic solver for systems of conservation laws.

``````#include "conservation.h"
``````

The conserved scalars are the gas density $\rho$ and the total energy $E$. The only conserved vector is the momentum $\mathbf{\text{w}}$. The constant $\gamma$ is represented by gammao here, with a default value of 1.4.

``````scalar ρ[], E[];
vector w[];
scalar * scalars = {ρ, E};
vector * vectors = {w};
double gammao = 1.4 ;``````

The system is entirely defined by the `flux()` function called by the generic solver for conservation laws. The parameter passed to the function is the array `s` which contains the state variables for each conserved field, in the order of their definition above (i.e. scalars then vectors).

``````void flux (const double * s, double * f, double e[2])
{``````

We first recover each value ($\rho$, $E$, ${w}_{x}$ and ${w}_{y}$) and then compute the corresponding fluxes (`f[0]`, `f[1]`, `f[2]` and `f[3]`).

``````  double ρ = s[0], E = s[1], wn = s[2], w2 = 0.;
for (int i = 2; i < 2 + dimension; i++)
w2 += sq(s[i]);
double un = wn/ρ, p = (gammao - 1.)*(E - 0.5*w2/ρ);

f[0] = wn;
f[1] = un*(E + p);
f[2] = un*wn + p;
for (int i = 3; i < 2 + dimension; i++)
f[i] = un*s[i];``````

The minimum and maximum eigenvalues for the Euler system are the characteristic speeds $u±\sqrt{\left(}\gamma p/\rho \right)$.

``````  double c = sqrt(gammao*p/ρ);
e[0] = un - c; // min
e[1] = un + c; // max
}``````