Compressible hydrodynamics

Compressible hydrodynamics#

The Euler equations of compressible hydrodynamics express conservation of mass, momentum, and energy. For the conserved state, \(\Uc = (\rho, \rho \Ub, \rho E)^\intercal\), the conserved system can be written as:

\[\Uc_t + \nabla \cdot {\bf F}(\Uc) = {\bf S}\]

where \({\bf F}\) are the fluxes and \({\bf S}\) are any source terms. In terms of components, they take the form (with a gravity source term):

\[\begin{split}\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \Ub) &= 0 \\ \frac{\partial (\rho \Ub)}{\partial t} + \nabla \cdot (\rho \Ub \Ub) + \nabla p &= \rho {\bf g} \\ \frac{\partial (\rho E)}{\partial t} + \nabla \cdot [(\rho E + p ) \Ub] &= \rho \Ub \cdot {\bf g}\end{split}\]

with \(\rho E = \rho e + \frac{1}{2} \rho |\Ub|^2\). The system is closed with an equation of state of the form:

\[p = p(\rho, e)\]

Note

The Euler equations 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 gamma-law 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 symbolic-link 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 second-order 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 y-direction)

    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

acoustic_pulse.rho0

1.4

acoustic_pulse.drho0

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 (plane-parallel). It will rise and deform (due to shear) parameters:

name

default

bubble.dens_base

10.0

bubble.scale_height

2.0

bubble.x_pert

2.0

bubble.y_pert

2.0

bubble.r_pert

0.25

bubble.pert_amplitude_factor

5.0

bubble.dens_cutoff

0.01

convection#

A heat source in a layer at some height above the bottom will drive convection in an adiabatically stratified atmosphere. parameters:

name

default

convection.dens_base

10.0

convection.scale_height

4.0

convection.y_height

2.0

convection.thickness

0.25

convection.e_rate

0.1

convection.dens_cutoff

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

gresho.rho0

1.0

gresho.r

0.2

gresho.mach

0.1

gresho.t_r

1.0

heating#

A test of the energy sources. This uses a uniform domain and slowly adds heat to the center over time. parameters:

name

default

heating.rho_ambient

1.0

heating.p_ambient

10.0

heating.r_src

0.1

heating.e_rate

0.1

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

hse.dens0

1.0

hse.h

1.0

kh#

A Kelvin-Helmholtz 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

kh.rho_1

1.0

kh.u_1

-1.0

kh.rho_2

2.0

kh.u_2

1.0

kh.bulk_velocity

0.0

plume#

A heat source at a point creates a plume that buoynantly rises in an adiabatically stratified atmosphere. parameters:

name

default

plume.dens_base

10.0

plume.scale_height

4.0

plume.x_pert

2.0

plume.y_pert

2.0

plume.r_pert

0.25

plume.e_rate

0.1

plume.dens_cutoff

0.01

quad#

The quadrant problem from Shulz-Rinne 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

quadrant.rho1

1.5

quadrant.u1

0.0

quadrant.v1

0.0

quadrant.p1

1.5

quadrant.rho2

0.532258064516129

quadrant.u2

1.206045378311055

quadrant.v2

0.0

quadrant.p2

0.3

quadrant.rho3

0.137992831541219

quadrant.u3

1.206045378311055

quadrant.v3

1.206045378311055

quadrant.p3

0.029032258064516

quadrant.rho4

0.532258064516129

quadrant.u4

0.0

quadrant.v4

1.206045378311055

quadrant.p4

0.3

quadrant.cx

0.5

quadrant.cy

0.5

ramp#

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

name

default

ramp.rhol

8.0

ramp.ul

7.1447096

ramp.vl

-4.125

ramp.pl

116.5

ramp.rhor

1.4

ramp.ur

0.0

ramp.vr

0.0

ramp.pr

1.0

rt#

A single-mode Rayleigh-Taylor instability. parameters:

name

default

rt.dens1

1.0

rt.dens2

2.0

rt.amp

1.0

rt.sigma

0.1

rt.p0

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

rt2.dens1

1.0

rt2.dens2

2.0

rt2.amp

1.0

rt2.sigma

0.1

rt2.p0

10.0

rt_multimode#

A multi-mode Rayleigh-Taylor instability. parameters:

name

default

rt_multimode.dens1

1.0

rt_multimode.dens2

2.0

rt_multimode.amp

1.0

rt_multimode.sigma

0.1

rt_multimode.nmodes

10

rt_multimode.p0

10.0

sedov#

The classic Sedov problem. parameters:

name

default

sedov.r_init

0.1

sedov.nsub

4

sod#

A general shock tube problem for comparing the solver to an exact Riemann solution. parameters:

name

default

sod.direction

x

sod.dens_left

1.0

sod.dens_right

0.125

sod.u_left

0.0

sod.u_right

0.0

sod.p_left

1.0

sod.p_right

0.1

test#

A setup intended for unit testing.

compressible_rk solver#

pyro.compressible_rk uses a method of lines time-integration approach with piecewise linear spatial reconstruction for the Euler equations. This is overall second-order 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 y-direction)

    riemann

    HLLC

    HLLC or CGF

    well_balanced

    0

    use a well-balanced 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 y-direction)

    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 spectral-deferred correction (SDC) for the time integration. This shares much in common with the pyro.compressible_fv4 solver, aside from how the time-integration 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 y-direction)

    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.