Compressible hydrodynamics#
The Euler equations of compressible hydrodynamics take the form:
with \(\rho E = \rho e + \frac{1}{2} \rho U^2\) and \(p = p(\rho, e)\). Note these do not include any dissipation terms, since they are usually negligible in astrophysics.
pyro has several compressible solvers to solve this equation set. The implementations here have flattening at shocks, artificial viscosity, a simple gammalaw equation of state, and (in some cases) a choice of Riemann solvers. Optional constant gravity in the vertical direction is allowed.
Note
All the compressible solvers share the same problems/
directory, which lives in compressible/problems/
. For the
other compressible solvers, we simply use a symboliclink to this
directory in the solver’s directory.
compressible
solver#
pyro.compressible
is based on a directionally unsplit (the corner
transport upwind algorithm) piecewise linear method for the Euler
equations, following [Colella90]. This is overall secondorder
accurate.
The parameters for this solver are:
section:
[compressible]
option
value
description
use_flattening
1
apply flattening at shocks (1)
z0
0.75
flattening z0 parameter
z1
0.85
flattening z1 parameter
delta
0.33
flattening delta parameter
cvisc
0.1
artificial viscosity coefficient
limiter
2
limiter (0 = none, 1 = 2nd order, 2 = 4th order)
grav
0.0
gravitational acceleration (in ydirection)
riemann
HLLC
HLLC or CGF
section:
[driver]
option
value
description
cfl
0.8
section:
[eos]
option
value
description
gamma
1.4
pres = rho ener (gamma  1)
section:
[particles]
option
value
description
do_particles
0
particle_generator
grid
supported problems#
acoustic_pulse
#
The acoustic pulse problem described in McCorquodale & Colella 2011. This uses a uniform background and a small pressure perturbation that drives a low Mach number soundwave. This problem is useful for testing convergence of a compressible solver.
parameters:
name 
default 


1.4 

0.14 
advect
#
A simple advection test. A density perturbation is set with a constant pressure in the domain and a velocity field is set to advect the profile across the domain. This is useful for testing convergence.
bubble
#
A buoyant perturbation (bubble) is placed in an isothermal hydrostatic atmosphere (planeparallel). It will rise and deform (due to shear) parameters:
name 
default 


10.0 

2.0 

2.0 

2.0 

0.25 

5.0 

0.01 
gresho
#
The Gresho vortex problem sets up a toroidal velocity field that is balanced by a radial pressure gradient. This is in equilibrium and the state should remain unchanged in time. This version of the problem is based on Miczek, Roepke, and Edelmann 2014. parameters:
name 
default 


1.0 

0.2 

59.5 

1.0 
hse
#
Initialize an isothermal hydrostatic atmosphere. It should remain static. This is a test of our treatment of the gravitational source term. parameters:
name 
default 


1.0 

1.0 
kh
#
A KelvinHelmholtz shear problem. There are 2 shear layers, with the and an optional vertical bulk velocity. This can be used to test the numerical dissipation in the solver. This setup is based on McNally et al. 2012. parameters:
name 
default 


1.0 

1.0 

2.0 

1.0 

0.0 
logo
#
Generate the pyro logo! The word “pyro” is written in the center of the domain and perturbations are placed in the 4 corners to drive converging shocks inward to scramble the logo.
quad
#
The quadrant problem from ShulzRinne et al. 1993; Lax and Lui 1998. Four different states are initialized in the quadrants of the domain, driving shocks and other hydrodynamic waves at the interfaces. This can be used to test the symmetry of the solver.
parameters:
name 
default 


1.5 

0.0 

0.0 

1.5 

0.532258064516129 

1.206045378311055 

0.0 

0.3 

0.137992831541219 

1.206045378311055 

1.206045378311055 

0.029032258064516 

0.532258064516129 

0.0 

1.206045378311055 

0.3 

0.5 

0.5 
ramp
#
A shock hitting a ramp at an oblique angle. This is based on Woodward & Colella 1984. parameters:
name 
default 


8.0 

7.1447096 

4.125 

116.5 

1.4 

0.0 

0.0 

1.0 
rt
#
A singlemode RayleighTaylor instability. parameters:
name 
default 


1.0 

2.0 

1.0 

0.1 

10.0 
rt2
#
A RT problem with two distinct modes: short wavelength on the left and long wavelength on the right. This allows one to see how the growth rate depends on wavenumber.
parameters:
name 
default 


1.0 

2.0 

1.0 

0.1 

10.0 
sedov
#
The classic Sedov problem. parameters:
name 
default 


0.1 

4 
sod
#
A general shock tube problem for comparing the solver to an exact Riemann solution. parameters:
name 
default 


x 

1.0 

0.125 

0.0 

0.0 

1.0 

0.1 
test
#
A setup intended for unit testing.
compressible_rk
solver#
pyro.compressible_rk
uses a method of lines timeintegration
approach with piecewise linear spatial reconstruction for the Euler
equations. This is overall secondorder accurate.
The parameters for this solver are:
section:
[compressible]
option
value
description
use_flattening
1
apply flattening at shocks (1)
z0
0.75
flattening z0 parameter
z1
0.85
flattening z1 parameter
delta
0.33
flattening delta parameter
cvisc
0.1
artificial viscosity coefficient
limiter
2
limiter (0 = none, 1 = 2nd order, 2 = 4th order)
temporal_method
RK4
integration method (see mesh/integration.py)
grav
0.0
gravitational acceleration (in ydirection)
riemann
HLLC
HLLC or CGF
well_balanced
0
use a wellbalanced scheme to keep the model in hydrostatic equilibrium
section:
[driver]
option
value
description
cfl
0.8
section:
[eos]
option
value
description
gamma
1.4
pres = rho ener (gamma  1)
supported problems#
compressible_rk
uses the problems defined by compressible
.
compressible_fv4
solver#
pyro.compressible_fv4
uses a 4th order accurate method with RK4
time integration, following [McCorquodaleColella11].
The parameter for this solver are:
section:
[compressible]
option
value
description
use_flattening
1
apply flattening at shocks (1)
z0
0.75
flattening z0 parameter
z1
0.85
flattening z1 parameter
delta
0.33
flattening delta parameter
cvisc
0.1
artificial viscosity coefficient
limiter
2
limiter (0 = none, 1 = 2nd order, 2 = 4th order)
temporal_method
RK4
integration method (see mesh/integration.py)
grav
0.0
gravitational acceleration (in ydirection)
riemann
CGF
section:
[driver]
option
value
description
cfl
0.8
section:
[eos]
option
value
description
gamma
1.4
pres = rho ener (gamma  1)
supported problems#
compressible_fv4
uses the problems defined by compressible
.
compressible_sdc
solver#
pyro.compressible_sdc
uses a 4th order accurate method with
spectraldeferred correction (SDC) for the time integration. This
shares much in common with the pyro.compressible_fv4
solver, aside from
how the timeintegration is handled.
The parameters for this solver are:
section:
[compressible]
option
value
description
use_flattening
1
apply flattening at shocks (1)
z0
0.75
flattening z0 parameter
z1
0.85
flattening z1 parameter
delta
0.33
flattening delta parameter
cvisc
0.1
artificial viscosity coefficient
limiter
2
limiter (0 = none, 1 = 2nd order, 2 = 4th order)
temporal_method
RK4
integration method (see mesh/integration.py)
grav
0.0
gravitational acceleration (in ydirection)
riemann
CGF
section:
[driver]
option
value
description
cfl
0.8
section:
[eos]
option
value
description
gamma
1.4
pres = rho ener (gamma  1)
supported problems#
compressible_sdc
uses the problems defined by compressible
.
Example problems#
Note
The 4thorder accurate solver (pyro.compressible_fv4
) requires that
the initialization create cellaverages accurate to 4thorder. To
allow for all the solvers to use the same problem setups, we assume
that the initialization routines initialize cellcenters (which is
fine for 2ndorder accuracy), and the
preevolve()
method will convert
these to cellaverages automatically after initialization.
Sod#
The Sod problem is a standard hydrodynamics problem. It is a onedimensional 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 gammalaw 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 onedimensional, we run it in narrow domains in the x or ydirections. It can be run as:
pyro_sim.py compressible sod inputs.sod.x
pyro_sim.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 xdirection, 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.
Sedov#
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_sim.py compressible sedov inputs.sedov
The video 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 gammalaw 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 2d
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.
quad#
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 appeared in several places (and a detailed investigation is online by Pawel Artymowicz). It is run as:
pyro_sim.py compressible quad inputs.quad
rt#
The RayleighTaylor 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_sim.py compressible rt inputs.rt
bubble#
The bubble problem initializes a hot spot in a stratified domain and watches it buoyantly rise and roll up. This is run as:
pyro_sim.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#
Explorations#
Measure the growth rate of the RayleighTaylor instability for different wavenumbers.
There are multiple Riemann solvers in the compressible algorithm. Run the same problem with the different Riemann solvers and look at the differences. Toro’s text is a good book to help understand what is happening.
Run the problems with and without limiting—do you notice any overshoots?
Extensions#
Limit on the characteristic variables instead of the primitive variables. What changes do you see? (the notes show how to implement this change.)
Add passively advected species to the solver.
Add an external heating term to the equations.
Add 2d axisymmetric coordinates (rz) to the solver. This is discussed in the notes. Run the Sedov problem with the explosion on the symmetric axis—now the solution will behave like the spherical sedov explosion instead of the cylindrical explosion.
Swap the piecewise linear reconstruction for piecewise parabolic (PPM). The notes and the Miller and Colella paper provide a good basis for this. Research the Roe Riemann solver and implement it in pyro.