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. In the near future, it will contain useful 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) = 0
Conservation of Momentum	\frac{\partial}{\partial t} (\rho v) + \frac{\partial}{\partial x} (\rho v^2 + p) = 0
Conservation of Energy	\frac{\partial}{\partial t} E + \frac{\partial}{\partial x} (v(E+p)) = 0

Where \rho 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 E = \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 e = \frac{p}{\rho(\gamma - 1)} where \gamma 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.

Solving the Riemann Problem

The Riemann problem is essentially the Euler equations with discrete initial conditions. Initial data has a left state and a right state, separated by a discontinuity.

Discretisation

Discretisation is the process of taking a domain and dividing it into 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

Program Outline

These are the basic steps in my program SEEQUOD (Solver for Euler EQUations in One Dimension). SEEQUOD currently uses an exact Riemann solver, but will later be updated to use approximate Riemann solvers as well. Also, SEEQUOD currently uses a first order Godunov upwind scheme, but again will later be updated to use higher order update schemes. 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
    • Solve local Riemann problem
    • 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

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 to obtain more accurate results. 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

Density, velocity, pressure, and internal energy were plotted for every test. The following images show the numerical results (points) 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

Test 2

Test 3

Test 4

Test 5