hydro by example

A simple python-based tutorial on computational methods for hydrodynamics









low Mach





∗ Compressible hydrodynamics solver

pyro has two compressible solvers:

The basic theory is outlined in Chapters 6 and 7 of the notes:
notes on computational hydro

The implementation here has flattening at shocks, artificial viscosity, a simple gamma-law equation of state, and two different Riemann solvers. Optional constant gravity in the vertical direction is allowed.

The main parameters that affect this solver are:

cfl the advective CFL number (what fraction of a zone can we cross in a single timestep)
use_flattening do we flatten the profiles at shocks? (0=no, 1=yes)
z0, z1, delta the parameters that affect the flattening algorithm
cvisc the coefficient for the artifical viscosity
limiter what type of limiting to use in reconstructing the slopes. 0 means use an unlimited second-order centered difference. 1 is the MC limiter, and 2 is the 4th-order MC limiter
riemann which Riemann solver do we use? "HLLE" for the HLLE solver, or "CGF" for the Colella, Glaz, and Ferguson solver
grav the gravitational acceleration (vertical/y direction)
temporal_method the MOL integration method to use (RK2, TVD2, TVD3, RK4) (compressible_rk only)
gamma the constant adiabatic index for the gas

∗ Examples


The Sod problem is a standard hydrodynamics problem. It is a one-dimensional shock tube (two states separated by an interface), that exhibits all three hydrodynamic waves: a shock, contact, and rarefaction. Furthermore, there are exact solutions for a gamma-law equation of state, so we can check our solution against these exact solutions. See Toro's book for details on this problem and the exact Riemann solver.

Because it is one-dimensional, we run it in narrow domains in the x- or y-directions. It can be run as:

./pyro.py compressible sod inputs.sod.x
./pyro.py compressible sod inputs.sod.y

A simple script, sod_compare.py in analysis/ will read a pyro output file and plot the solution over the exact Sod solution. Below we see the result for a Sod run with 128 points in the x-direction, gamma = 1.4, and run until t = 0.2 s.

We see excellent agreement for all quantities. The shock wave is very steep, as expected. The contact wave is smeared out over ~5 zones—this is discussed in the notes above, and can be improved in the PPM method with contact steepening.


The Sedov blast wave problem is another standard test with an analytic solution (Sedov 1959). A lot of energy is point into a point in a uniform medium and a blast wave propagates outward. The Sedov problem is run as:

./pyro.py compressible sedov inputs.sedov

The image below shows the output from a 128 x 128 grid with the energy put in a radius of 0.0125 surrounding the center of the domain. A gamma-law EOS with gamma = 1.4 is used, and we run until 0.1

We see some grid effects because it is hard to initialize a small circular explosion on a rectangular grid. To compare to the analytic solution, we need to radially bin the data. Since this is a 2-d explosion, the physical geometry it represents is a cylindrical blast wave, so we compare to Sedov's cylindrical solution. The radial binning is done with the sedov_compare.py script in analysis/

This shows good agreement with the analytic solution.


The quad problem sets up different states in four regions of the domain and watches the complex interfaces that develop as shocks interact. This problem has appear in several places (and a detailed investigation is online by Pawel Artymowicz). It is run as:

./pyro.py compressible quad inputs.quad


The Rayleigh-Taylor problem puts a dense fluid over a lighter one and perturbs the interface with a sinusoidal velocity. Hydrostatic boundary conditions are used to ensure any initial pressure waves can escape the domain. It is run as:

./pyro.py compressible er inputs.rt


The bubble problem initializes a hot spot in a stratified domain and watches it buoyantly rise and roll up. This is run as:

./pyro.py compressible bubble inputs.bubble

The shock at the top of the domain is because we cut off the stratified atmosphere at some low density and the resulting material above that rains down on our atmosphere. Also note the acoustic signal propagating outward from the bubble (visible in the U and e panels).

∗ Exercises



∗ Going further

The compressible algorithm presented here is essentially the single-grid hydrodynamics algorithm used in the Castro code—an adaptive mesh radiation hydrodynamics code developed at CCSE/LBNL. Castro is freely available for download.

A simple, pure Fortran, 1-d compressible hydrodynamics code that does piecewise constant, linear, or parabolic (PPM) reconstruction is also available. See the hydro1d page.