Building a 1D Hydro Code

Introduction

I've created this page to document the important aspects of building a hydro code. It contains general information that I have learned from reading Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction by E.F. Toro. The page contains what I feel is the most important information from my own personal experience as I build my first hydro code.

The Euler Equations

These are the so called Euler Equations. They are also known as the fluid equations in conservative-law form. These are conceptually simple and make intuitive sense.

Name of Law	Formula
Conservation of Mass	frac{\partial}{\partial t}\rho + \frac{\partial}{\partial x} (\rho v) =
Conservation of Momentum	frac{\partial}{\partial t} (\rho v) + \frac{\partial}{\partial x} (\rho v^2 + p) =
Conservation of Energy	frac{\partial}{\partial t} E + \frac{\partial}{\partial x} (v(E+p)) =

Where rh is mass density, v is velocity (only in x-direction for 1D), p is pressure, and E is total energy per unit volume. E is further defined as = \rho (\frac{1}{2} v^2 + e where e is the specific internal energy. The specific internal energy depends on the equation of state. For an ideal gas = \frac{p}{\rho(\gamma - 1) where gamm is the ratio of specific heats. Basically, these laws state that in a given volume, the change in a conserved quantity must be equal to the flux through the boundaries of that volume. In other words, the conserved quantity is in front of the time derivative and its flux is in front of the spatial derivative.

The Euler equations can also be written in integral form for a general control volume [x_1,x_2] \ \mathrm{x} \ [t_1,t_2] . To simplify the notation, let U be a vector containing the conserved quantities and F be the vector of the corresponding fluxes. Now the Euler equations are:

\int^{x_2}_{x_1} U(x,t_2) \ \mathrm{d}x = \int_{x_1}^{x_2} U(x,t_1) \ \mathrm{d}x + \int_{t_1}^{t_2}F(U(x_1,t)) \ \mathrm{d}t - \int_{t_1}^{t_2} F(U(x_2,t)) \ \mathrm{d}t \hspace{1 in} (1)

In this form, it might be easier to see that the change in a conserved quantity is equal to its flux coming into the volume minus its flux going out of the volume.

The Riemann Problem

The Riemann problem is essentially the Euler equations with discrete initial conditions. The initial conditions consist of a left state and a right state separated by a discontinuity. You can turn the differential form of the Euler equations into an eigenvalue problem by writing the flux term as a Jacobian matrix multiplied by the conserved variable vector. The Jacobian matrix involves various derivatives of the conserved quantities. If you solved this eigensystem, you would find three eigenvalues and three corresponding eigenvectors. These eigenvalues correspond to wave speeds of v, v + a, and v - a where v is the velocity and a is the sound speed. The sound speed is defined as:

a=\sqrt{\frac{\gamma p}{\rho}}

An analysis of the eigenvectors would show that there are essentially three types of waves; contact, shock, and rarefaction. The contact wave is always in the middle, and the left or right waves can be either shock or rarefaction waves. The waves divide x-t space into four distinct regions: left data, star left, star right, and right data. Furthermore, velocity and pressure inside the star region is constant while density can be different in the star left and star right regions. Data across a shock wave has a discontinuous jump while data across a rarefaction has a smooth transition. This means that a rarefaction wave actually has some "thickness" to it. The region inside a rarefaction is known as the fan state.

This is essentially the problem set-up. Now we need to actually solve the problem using some type of numerical method.

The Exact Riemann Solver

A description of an exact Riemann solver can be quite tedious and lengthy, so I will leave out many details in this section. In solving the Riemann problem, we will use what are known as primitive variables. The primitive variables are density, velocity, and pressure. It is important to keep in mind that, apart from density, these are different than the conserved variables. Basically, the unknowns are in the star region, so if we can find these quantities then the problem is solved.

1. Find the pressure and velocity in the star region

This pressure is computed through an iteration scheme. You start with a guess* for the star pressure. This can be found using approximate Riemann solvers. Then, you calculate pressure functions to get a better value for the pressure. If the difference between the guess pressure and better pressure is below some tolerance then you have found the star pressure. If it is not below the tolerance, the better pressure becomes the new guess pressure and the process repeats. The star velocity is then very straightforward, as it is computed directly from the star pressure.

2. Sample the wave speed

Now that we have the solution for the star region, we can get the solution everywhere. This is done by sampling the wave speed S as opposed to sampling the position x. There are essentially 10 possible solutions for any given S. Let W be a vector containing the correct values for the primitive variables.

	Left side of contact discontinuity	Right side of contact discontinuity
Data	W_L	W_R
Star (Shock)	W_{\mathrm{*}L}^{SHO}	W_{\mathrm{*}R}^{SHO}
Star (Rarefaction)	W_{\mathrm{*}L}^{RAR}	W_{\mathrm{*}R}^{RAR}
Rarefaction Fan	W_{L}^{FAN}	W_{R}^{FAN}

A sampling procedure involves checking several logical statements to determine which state S is in. With the Riemann problem now solved, we have the correct values for the primitive variables.

The HLLC Solver

The HLLC solver is different in that it computes an approximation for the fluxes directly rather than using a Riemann problem solution to calculate the exact fluxes. The basic idea is to assume a wave configuration of three waves which separate four constant regions. The two middle regions form the star region, so the HLLC solver needs values for the pressure and velocity in the star region. This can be done using the same guess procedure* mentioned in the exact Riemann solver description. These values are then used to compute various wave speeds, and the wave speeds are in turn used to compute the fluxes.

The Roe Solver

The Roe solver is different from all other solvers presented on this page because it makes an approximation to the Riemann problem itself and then solves that problem exactly. The Roe method takes the Riemann problem and linearizes it by assuming a constant coefficient matrix instead of a Jacobian matrix as mentioned in the Riemann problem section. The eigenvalues and eigenvectors of this new linear system are then computed. Then, a parameter known as the wave strength is computed. The wave strengths are found from projecting the difference of the data across a cell boundary onto the eigenvectors. Lastly, the fluxes can be computed using corresponding eigenvalues, wave strengths, and eigenvectors.

The original Roe method has trouble at sonic points inside rarefactions. It will produce waves known as rarefaction shocks which violate an entropy condition. Values from the star region can be used to fix this problem. A simple way to obtain these values is to use the same guess procedure* from the other two Riemann solvers. These star region values define new eigenvalues that are to be used near sonic points. Thus the entropy problem is resolved.

* A Note on the "Guess" Procedure

Since all three Riemann solvers use this procedure, I thought it was important to explain it in more detail. The guess procedure is essentially an adaptive non-iterative Riemann solver. It uses three different approximate Riemann solvers: the Primitive Variable Riemann Solver (PVRS), the Two Rarefaction Riemann Solver (TRRS), and the Two Shock Riemann Solver (TSRS). These solvers are all easy to implement and involve simple approximations. The program chooses which solver to use based on the initial conditions. The PVRS is similar to the Roe solver in that it linearizes the Riemann problem, but it then makes simpler approximations to compute the solution. The TRRS and TSRS are both based on the exact Riemann solver. These basically assume a particular wave configuration and use those simple solutions as opposed to a sampling procedure.

Discretisation

Discretisation is the process of taking a domain and dividing it into cells. The "mesh size" \Delta x is defined as the domain length divided by the number of cells. Each cell contains data corresponding to a specific problem initialization thus creating many Riemann problems. In other words, each neighboring pair of cells is now a "local" Riemann problem. Given two neighboring cells, the left cell is the left data state, the right cell is the right data state, and the intercell boundary is the discontinuity position. Now each local Riemann problem can be solved and that solution can be used to calculate the fluxes required to update the cells to the next time step.

Godunov's First Order Upwind Scheme

First, we want to use the integral form of the Euler equations (1). Next, we define a cell average for cell i at time step n + 1:

U_{i}^{n+1} = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \tilde{U}(x,t^{n+1}) \ \mathrm{d}x$ \hspace{1 in}(2

Where \tilde{U} is the global solution as opposed to the local solution U. However, we already know that the global solution can be written in terms of its fluxes. You just use (1) with a specific control volume:

\int^{x_{i+\frac{1}{2}}}_{x_{i - \frac{1}{2}}} \tilde{U}(x,t^{n+1}) \ \mathrm{d}x = \int^{x_{i+\frac{1}{2}}}_{x_{i - \frac{1}{2}}} \tilde{U}(x,t^{n}) \ \mathrm{d}x + \int_{t^{n}}^{t^{n+1}}F(\tilde{U}(x_{i-\frac{1}{2}},t)) \ \mathrm{d}t - \int_{t^{n}}^{t^{n+1}} F(\tilde{U}(x_{i+\frac{1}{2}},t)) \ \mathrm{d}t \hspace{1 in} (3)

The trick is to write the global solution in terms of the local solution. This can be done through a transformation of coordinates:

\tilde{U}(x,t) = U_{i+\frac{1}{2}}(\bar{x},\bar{t}) \hspace{1 in} (4)

Where

\bar{x}=x-x_{i+\frac{1}{2}} \ \ \mathrm{and} \ \ \bar{t}=t-t^n \hspace{1 in} (5)

Now we can combine (4) and (5) to get:

\tilde{U}(x_{i-\frac{1}{2}},t) = U_{i-\frac{1}{2}}(0) = constant \hspace{1 in} (6)

\tilde{U}(x_{i+\frac{1}{2}},t) = U_{i+\frac{1}{2}}(0) = constant \hspace{1 in} (7)

Substitute (6) and (7) into (3), and then substituting that result into (2) gives:

U_{i}^{n+1} = \frac{1}{\Delta x} \int^{x_{i+\frac{1}{2}}}_{x_{i - \frac{1}{2}}} \tilde{U}(x,t^{n}) \ \mathrm{d}x + \frac{1}{\Delta x} \int_{t^{n}}^{t^{n+1}}F(U_{i-\frac{1}{2}}(0,t)) \ \mathrm{d}t - \frac{1}{\Delta x} \int_{t^{n}}^{t^{n+1}} F(U_{i+\frac{1}{2}}(0,t)) \ \mathrm{d}t \hspace{1 in} (8)

Just looking at the right hand side of (8)…the first term can be rewritten just by using the definition of cell average that we defined in (2), and since the flux integrands are now constants, those integrals can be simplified. The final equation reads:

U_{i}^{n+1} = U_{i}^{n} + \frac{\Delta t}{\Delta x}[F(U_{i-\frac{1}{2}}(0))-F(U_{i+\frac{1}{2}}(0))] \hspace{1 in} (9)

In summary, the necessary quantities to update a cell are: its initial quantities, the mesh size, the time step size, and the fluxes through its boundaries. Initial quantities are known, the mesh size is chosen by the user, and the fluxes are calculated from solving the Riemann problem. However, in order to calculate the fluxes for the first and last cell we need more information.

Boundary Conditions

When you discretize a domain, you have to create "ghost" cells outside your domain so that you can calculate the necessary fluxes. These ghost cells are assigned values depending on the boundary conditions, and this assignment is imposed every time step. The three most common boundary conditions are transparent, reflective, and periodic.

Transparent

Transparent boundary conditions are the simplest of the three. A ghost cell is simply given the values of its neighboring cell.

Reflective

Think of reflective boundary conditions as a rigid wall. A ghost cell is given the same values as its neighboring cell, except its velocity is negated.

Periodic

For periodic boundary conditions, I like to think of the old Atari game Asteroids. If material leaves through the right boundary it comes back through the left boundary. To achieve this, the first ghost cell is given the same values as the last domain cell, and the last ghost cell is given the same values as the first domain cell.

Now that we have boundary conditions, we can get all the fluxes for our domain cells. However, we still need to choose an appropriate time step size to make our updating scheme complete.

The CFL Condition

CFL stands for Courant-Friedrichs-Lewy. The CFL coefficient (often just referred to as the CFL) is used to calculate the time step size. The time step size along with the CFL condition is:

\Delta t = \frac{CFL \ \Delta x}{S_{max}}

Where 0<CFL \le 1. S_{max} is the maximum wave speed within the domain at the current computational time. This condition ensures that the time step is not large enough to allow a wave to pass through an entire cell within that time. A CFL of 1 gives the maximum time step size, which is the most efficient, but can often lead to inaccuracies if S_{max} is not reasonably known. Therefore, it is common to use a CFL of 0.9. A popular estimate for the maximum wave speed is:

S_{max} = max[|v_i| + a_i]

Where v is the particle velocity and a is the sound speed. Like the boundary condition routine, the CFL condition routine must be called every time step since S_{max} is not constant.

Now we have everything we need to write a hydro code.

The MUSCL-Hancock Scheme

The MUSCL-Hancock scheme achieves second order accuracy by adding a few more steps to the process. The essential idea here is that instead of using the initial data directly we are going to alter the data slightly. This process is known as reconstruction. Instead of having constant data across a cell, we give the data a slope. In other words, the data in each cell is now represented by a sloped line rather than a flat line. This methodology can achieve second order accuracy, and other methods can give higher order accuracy. For example, representing the data as a quadratic could achieve third order accuracy and so on.

1. Compute a limited slope

The first thing one usually thinks to use for a slope is the difference of data states across a cell boundary. However, this will produce what are known as spurious oscillations, because a simple difference can overestimate the desired slope. This is why most methods use some type of limited slope. This can be computed directly or via a slope limiter. The slope limiter is another parameter that is just multiplied by the original slope to keep it from becoming too large.

2. Reconstruct the data

Once a suitable slope is found, it can be used to compute what are known as boundary extrapolated values. These are the extreme points of the data located at the intercell boundaries. This data is then evolved in time by half of a time step via a formula similar to (9). Now we have the fully reconstructed, evolved data. This data is then used to solve the Riemann problem, and the process proceeds exactly as the first order Godunov method.

Program Outline

These are the basic steps in my program SEEQUOD (Solver for Euler EQUations in One Dimension). SEEQUOD currently uses either an exact Riemann solver, an HLLC solver, or a Roe solver. SEEQUOD can be easily modified to use either the aforementioned first order Godunov scheme or the second order MUSCL-Hancock scheme. These steps should be general enough to apply to many codes that are used to solve the Euler Equations.

• Read all input data necessary for problem
• Initialize domain cells
• Begin time stepping procedure
  • Apply boundary conditions for "ghost" cells
  • Impose CFL condition to get appropriate time step
  • Begin flux calculation procedure
    • Calculate a limited slope (MUSCL only)
    • Reconstruct the data (MUSCL only)
    • Solve local Riemann problem (exact, HLLC, or Roe)
    • Use solution from Riemann solver to compute fluxes
    • Repeat flux calculation procedure until fluxes are computed for all cells
  • Use fluxes to update conserved variables
  • Use updated conserved variables to update primitive variables
  • Repeat time stepping procedure until final time is reached
• Output data to a file

The MUSCL only steps can easily be omitted to obtain the first order Godunov results.

Test Results

Five different tests were run in accordance with the tests in Ch. 6 of Toro's aforementioned book. All tests were run with a domain length of 1.0 and 100 computing cells. For the first five time steps a CFL of 0.18 was used to give the various waves time to develop. For the time steps following, a more natural CFL of 0.9 was used for better efficiency. Here is a table of all the input data for each test:

Test Number	Discontinuity Position	Output Time	\rho_L	v_L	p_L	\rho_R	v_R	p_R
1	0.3	0.2	1.0	0.75	1.0	0.125	0.0	0.1
2	0.5	0.15	1.0	-2.0	0.4	1.0	2.0	0.4
3	0.5	0.012	1.0	0.0	1000.0	1.0	0.0	0.01
4	0.4	0.035	5.99924	19.5975	460.894	5.99242	-6.19633	46.0950
5	0.8	0.012	1.0	-19.59745	1000.0	1.0	-19.59745	0.01

Using the exact Riemann solver, density, velocity, pressure, and internal energy were plotted for every test. These plots correspond to the first five groups of plots.

Test 1 was also run with the HLLC solver. The density plot was made to compare it to the results of the exact solver. Notice that the HLLC solver actually does a better job of resolving the left rarefaction. The exact solution has a much more noticeable entropy "glitch". This is what we call that gap inside the left rarefaction.

Test 1 was also run with the Roe solver in two different ways. First as a "barebones" Roe solver, and second with an entropy fix. This entropy fix is required for transonic rarefactions, which is the case for the left rarefaction in Test 1. You can see that without this entropy fix, the gap inside that left rarefaction is very large and completely inaccurate.

The tests were also run with the MUSCL-Hancock scheme with the HLLC solver. A comparison with the first order plots for test 1 shows that the MUSCL-Hancock method does indeed give more accurate results. The MUSCL-Hancock scheme produced similar results for the other tests except for test 2. For some reason the MUSCL-Hancock scheme cannot handle the two rarefaction problem no matter what the CFL, slope limiter, etc.

For all the plots, the numerical results (points) are alongside the exact solution (line).

All plots appear to be the same as those presented by Toro. Therefore, the codes SEEQUOD and SEEQUODexact must be correct.

Test 1 (Exact)

Test 2 (Exact)

Test 3 (Exact)

Test 4 (Exact)

Test 5 (Exact)

Test 1 (HLLC vs. Exact)

HLLC	Exact

Test 1 (Roe without entropy fix vs. Roe with entropy fix)

Roe without entropy fix	Roe with entropy fix

Test 1 (Godunov vs. MUSCL-Hancock)

Godunov	MUSCL-Hancock

Test 2 (MUSCL-Hancock)