pyro.compressible package

pyro.compressible package#

The pyro compressible hydrodynamics solver. This implements the second-order (piecewise-linear), unsplit method of Colella 1990.

Subpackages#

Submodules#

pyro.compressible.BC module#

compressible-specific boundary conditions. Here, in particular, we implement an HSE BC in the vertical direction.

Note: the pyro BC routines operate on a single variable at a time, so some work will necessarily be repeated.

Also note: we may come in here with the aux_data (source terms), so we’ll do a special case for them

pyro.compressible.BC.inflow_post_bc(var, g)[source]#
pyro.compressible.BC.inflow_pre_bc(var, g)[source]#
pyro.compressible.BC.user(bc_name, bc_edge, variable, ccdata)[source]#

A hydrostatic boundary. This integrates the equation of HSE into the ghost cells to get the pressure and density under the assumption that the specific internal energy is constant.

Upon exit, the ghost cells for the input variable will be set

Parameters:
bc_name{‘hse’}

The descriptive name for the boundary condition – this allows for pyro to have multiple types of user-supplied boundary conditions. For this module, it needs to be ‘hse’.

bc_edge{‘ylb’, ‘yrb’}

The boundary to update: ylb = lower y boundary; yrb = upper y boundary.

variable{‘density’, ‘x-momentum’, ‘y-momentum’, ‘energy’}

The variable whose ghost cells we are filling

ccdataCellCenterData2d object

The data object

pyro.compressible.derives module#

pyro.compressible.derives.derive_primitives(myd, varnames)[source]#

derive desired primitive variables from conserved state

pyro.compressible.eos module#

This is a gamma-law equation of state: p = rho e (gamma - 1), where gamma is the constant ratio of specific heats.

pyro.compressible.eos.dens(gamma, p, eint)[source]#

Given the pressure and the specific internal energy, return the density

Parameters:
gammafloat

The ratio of specific heats

pfloat

The pressure

eintfloat

The specific internal energy

Returns:
outfloat

The density

pyro.compressible.eos.pres(gamma, rho, eint)[source]#

Given the density and the specific internal energy, return the pressure

Parameters:
gammafloat

The ratio of specific heats

rhofloat

The density

eintfloat

The specific internal energy

Returns:
outfloat

The pressure

pyro.compressible.eos.rhoe(gamma, p)[source]#

Given the pressure, return (rho * e)

Parameters:
gammafloat

The ratio of specific heats

pfloat

The pressure

Returns:
outfloat

The internal energy density, rho e

pyro.compressible.interface module#

pyro.compressible.interface.artificial_viscosity(ng, dx, dy, Lx, Ly, xmin, ymin, coord_type, cvisc, u, v)[source]#

Compute the artificial viscosity. Here, we compute edge-centered approximations to the divergence of the velocity. This follows directly Colella Woodward (1984) Eq. 4.5

data locations:

j+3/2--+---------+---------+---------+
       |         |         |         |
  j+1  +         |         |         |
       |         |         |         |
j+1/2--+---------+---------+---------+
       |         |         |         |
     j +         X         |         |
       |         |         |         |
j-1/2--+---------+----Y----+---------+
       |         |         |         |
   j-1 +         |         |         |
       |         |         |         |
j-3/2--+---------+---------+---------+
       |    |    |    |    |    |    |
           i-1        i        i+1
     i-3/2     i-1/2     i+1/2     i+3/2

X is the location of avisco_x[i,j] Y is the location of avisco_y[i,j]

Parameters:
ngint

The number of ghost cells

dx, dyfloat

Cell spacings

xmin, yminfloat

Min physical x, y boundary

Lx, Ly: ndarray

Cell size in x, y direction

cviscfloat

viscosity parameter

u, vndarray

x- and y-velocities

Returns:
outndarray, ndarray

Artificial viscosity in the x- and y-directions

pyro.compressible.interface.states(idir, ng, dx, dt, irho, iu, iv, ip, ix, nspec, gamma, qv, dqv)[source]#

predict the cell-centered state to the edges in one-dimension using the reconstructed, limited slopes.

We follow the convection here that V_l[i] is the left state at the i-1/2 interface and V_l[i+1] is the left state at the i+1/2 interface.

We need the left and right eigenvectors and the eigenvalues for the system projected along the x-direction.

Taking our state vector as \(Q = (\rho, u, v, p)^T\), the eigenvalues are \(u - c\), \(u\), \(u + c\).

We look at the equations of hydrodynamics in a split fashion – i.e., we only consider one dimension at a time.

Considering advection in the x-direction, the Jacobian matrix for the primitive variable formulation of the Euler equations projected in the x-direction is:

    / u   r   0   0 \
    | 0   u   0  1/r |
A = | 0   0   u   0  |
    \ 0  rc^2 0   u  /

The right eigenvectors are:

     /  1  \        / 1 \        / 0 \        /  1  \
     |-c/r |        | 0 |        | 0 |        | c/r |
r1 = |  0  |   r2 = | 0 |   r3 = | 1 |   r4 = |  0  |
     \ c^2 /        \ 0 /        \ 0 /        \ c^2 /

In particular, we see from r3 that the transverse velocity (v in this case) is simply advected at a speed u in the x-direction.

The left eigenvectors are:

l1 =     ( 0,  -r/(2c),  0, 1/(2c^2) )
l2 =     ( 1,     0,     0,  -1/c^2  )
l3 =     ( 0,     0,     1,     0    )
l4 =     ( 0,   r/(2c),  0, 1/(2c^2) )

The fluxes are going to be defined on the left edge of the computational zones:

 |             |             |             |
 |             |             |             |
-+------+------+------+------+------+------+--
 |     i-1     |      i      |     i+1     |
              ^ ^           ^
          q_l,i q_r,i  q_l,i+1

q_r,i and q_l,i+1 are computed using the information in zone i,j.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

ngint

The number of ghost cells

dxndarray

The cell spacing

dtfloat

The timestep

irho, iu, iv, ip, ixint

Indices of the density, x-velocity, y-velocity, pressure and species in the state vector

nspecint

The number of species

gammafloat

Adiabatic index

qvndarray

The primitive state vector

dqvndarray

Spatial derivative of the state vector

Returns:
outndarray, ndarray

State vector predicted to the left and right edges

pyro.compressible.riemann module#

pyro.compressible.riemann.consFlux(idir, coord_type, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state)[source]#

Calculate the conservative flux.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

gammafloat

Adiabatic index

idens, ixmom, iymom, iener, irhoXint

The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.

nspecint

The number of species

U_statendarray

Conserved state vector.

Returns:
outndarray

Conserved flux

pyro.compressible.riemann.estimate_wave_speed(rho_l, u_l, p_l, c_l, rho_r, u_r, p_r, c_r, gamma)[source]#
pyro.compressible.riemann.riemann_cgf(idir, ng, idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, gamma, U_l, U_r)[source]#

Solve riemann shock tube problem for a general equation of state using the method of Colella, Glaz, and Ferguson. See Almgren et al. 2010 (the CASTRO paper) for details.

The Riemann problem for the Euler’s equation produces 4 regions, separated by the three characteristics (u - cs, u, u + cs):

u - cs    t    u      u + cs
  \       ^   .       /
   \  *L  |   . *R   /
    \     |  .     /
     \    |  .    /
 L    \   | .   /    R
       \  | .  /
        \ |. /
         \|./
----------+----------------> x

We care about the solution on the axis. The basic idea is to use estimates of the wave speeds to figure out which region we are in, and: use jump conditions to evaluate the state there.

Only density jumps across the u characteristic. All primitive variables jump across the other two. Special attention is needed if a rarefaction spans the axis.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

ngint

The number of ghost cells

nspecint

The number of species

idens, ixmom, iymom, iener, irhoXint

The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gammafloat

Adiabatic index

U_l, U_rndarray

Conserved state on the left and right cell edges.

Returns:
outndarray

Conserved states.

pyro.compressible.riemann.riemann_flux(idir, U_l, U_r, my_data, rp, ivars, lower_solid, upper_solid, tc, return_cons=False)[source]#

This is the general interface that constructs the unsplit fluxes through the idir (1 for x, 2 for y) interfaces using the left and right conserved states by using the riemann solver.

Parameters:
U_l, U_r: ndarray, ndarray

Conserved states in the left and right interface

my_dataCellCenterData2d object

The data object containing the grid and advective scalar that we are advecting.

rpRuntimeParameters object

The runtime parameters for the simulation

ivarsVariables object

The Variables object that tells us which indices refer to which variables

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

tcTimerCollection object

The timers we are using to profile

return_cons: Boolean

If we don’t use HLLC Riemann solver, do we also return conserved states?

Returns:
Fndarray

Fluxes in x or y direction

Optionally:
U: ndarray

Conserved states in x or y direction

pyro.compressible.riemann.riemann_hllc(idir, ng, idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, gamma, U_l, U_r)[source]#

This is the HLLC Riemann solver. The implementation follows directly out of Toro’s book. Note: this does not handle the transonic rarefaction.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

ngint

The number of ghost cells

nspecint

The number of species

idens, ixmom, iymom, iener, irhoXint

The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gammafloat

Adiabatic index

U_l, U_rndarray

Conserved state on the left and right cell edges.

Returns:
outndarray

Conserved flux

pyro.compressible.riemann.riemann_hllc_lowspeed(idir, ng, idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, gamma, U_l, U_r)[source]#

This is the HLLC Riemann solver based on Toro (2009) alternate formulation (Eqs. 10.43 and 10.44) and the low Mach number asymptotic fix of Minoshima & Miyoshi (2021). It is also based on the Quokka implementation.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

ngint

The number of ghost cells

nspecint

The number of species

idens, ixmom, iymom, iener, irhoXint

The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gammafloat

Adiabatic index

U_l, U_rndarray

Conserved state on the left and right cell edges.

Returns:
outndarray

Conserved flux

pyro.compressible.riemann.riemann_prim(idir, ng, irho, iu, iv, ip, iX, nspec, lower_solid, upper_solid, gamma, q_l, q_r)[source]#

this is like riemann_cgf, except that it works on a primitive variable input state and returns the primitive variable interface state

Solve riemann shock tube problem for a general equation of state using the method of Colella, Glaz, and Ferguson. See Almgren et al. 2010 (the CASTRO paper) for details.

The Riemann problem for the Euler’s equation produces 4 regions, separated by the three characteristics \((u - cs, u, u + cs)\):

u - cs    t    u      u + cs
  \       ^   .       /
   \  *L  |   . *R   /
    \     |  .     /
     \    |  .    /
 L    \   | .   /    R
       \  | .  /
        \ |. /
         \|./
----------+----------------> x

We care about the solution on the axis. The basic idea is to use estimates of the wave speeds to figure out which region we are in, and: use jump conditions to evaluate the state there.

Only density jumps across the \(u\) characteristic. All primitive variables jump across the other two. Special attention is needed if a rarefaction spans the axis.

Parameters:
idirint

Are we predicting to the edges in the x-direction (1) or y-direction (2)?

ngint

The number of ghost cells

nspecint

The number of species

irho, iu, iv, ip, iXint

The indices of the density, x-velocity, y-velocity, pressure and species fractions in the state vector.

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gammafloat

Adiabatic index

q_l, q_rndarray

Primitive state on the left and right cell edges.

Returns:
outndarray

Primitive flux

pyro.compressible.simulation module#

class pyro.compressible.simulation.Simulation(solver_name, problem_name, problem_func, rp, *, problem_finalize_func=None, problem_source_func=None, timers=None, data_class=<class 'pyro.mesh.patch.CellCenterData2d'>)[source]#

Bases: NullSimulation

The main simulation class for the corner transport upwind compressible hydrodynamics solver

dovis()[source]#

Do runtime visualization.

evolve()[source]#

Evolve the equations of compressible hydrodynamics through a timestep dt.

initialize(*, extra_vars=None, ng=4)[source]#

Initialize the grid and variables for compressible flow and set the initial conditions for the chosen problem.

method_compute_timestep()[source]#

The timestep function computes the advective timestep (CFL) constraint. The CFL constraint says that information cannot propagate further than one zone per timestep.

We use the driver.cfl parameter to control what fraction of the CFL step we actually take.

write_extras(f)[source]#

Output simulation-specific data to the h5py file f

class pyro.compressible.simulation.Variables(myd)[source]#

Bases: object

a container class for easy access to the different compressible variable by an integer key

pyro.compressible.simulation.cons_to_prim(U, gamma, ivars, myg)[source]#

convert an input vector of conserved variables to primitive variables

pyro.compressible.simulation.get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source=None)[source]#

compute the external sources, including gravity

pyro.compressible.simulation.prim_to_cons(q, gamma, ivars, myg)[source]#

convert an input vector of primitive variables to conserved variables

pyro.compressible.unsplit_fluxes module#

Implementation of the Colella 2nd order unsplit Godunov scheme. This is a 2-dimensional implementation only. We assume that the grid is uniform, but it is relatively straightforward to relax this assumption.

There are several different options for this solver (they are all discussed in the Colella paper).

  • limiter: 0 = no limiting; 1 = 2nd order MC limiter; 2 = 4th order MC limiter

  • riemann: HLLC or CGF (for Colella, Glaz, and Freguson solver)

  • use_flattening: set to 1 to use the multidimensional flattening at shocks

  • delta, z0, z1: flattening parameters (we use Colella 1990 defaults)

The grid indices look like

j+3/2--+---------+---------+---------+
       |         |         |         |
  j+1 _|         |         |         |
       |         |         |         |
       |         |         |         |
j+1/2--+---------XXXXXXXXXXX---------+
       |         X         X         |
    j _|         X         X         |
       |         X         X         |
       |         X         X         |
j-1/2--+---------XXXXXXXXXXX---------+
       |         |         |         |
  j-1 _|         |         |         |
       |         |         |         |
       |         |         |         |
j-3/2--+---------+---------+---------+
       |    |    |    |    |    |    |
           i-1        i        i+1
     i-3/2     i-1/2     i+1/2     i+3/2

We wish to solve

\[U_t + F^x_x + F^y_y = H\]

we want U_{i+1/2}^{n+1/2} – the interface values that are input to the Riemann problem through the faces for each zone.

Taylor expanding yields:

 n+1/2                     dU           dU
U          = U   + 0.5 dx  --  + 0.5 dt --
 i+1/2,j,L    i,j          dx           dt


                           dU             dF^x   dF^y
           = U   + 0.5 dx  --  - 0.5 dt ( ---- + ---- - H )
              i,j          dx              dx     dy


                            dU      dF^x            dF^y
           = U   + 0.5 ( dx -- - dt ---- ) - 0.5 dt ---- + 0.5 dt H
              i,j           dx       dx              dy


                                dt       dU           dF^y
           = U   + 0.5 dx ( 1 - -- A^x ) --  - 0.5 dt ---- + 0.5 dt H
              i,j               dx       dx            dy


                              dt       _            dF^y
           = U   + 0.5  ( 1 - -- A^x ) DU  - 0.5 dt ---- + 0.5 dt H
              i,j             dx                     dy

                   +----------+-----------+  +----+----+   +---+---+
                              |                   |            |

                  this is the monotonized   this is the   source term
                  central difference term   transverse
                                            flux term

There are two components, the central difference in the normal to the interface, and the transverse flux difference. This is done for the left and right sides of all 4 interfaces in a zone, which are then used as input to the Riemann problem, yielding the 1/2 time interface values:

 n+1/2
U
 i+1/2,j

Then, the zone average values are updated in the usual finite-volume way:

 n+1    n     dt    x  n+1/2       x  n+1/2
U    = U    + -- { F (U       ) - F (U       ) }
 i,j    i,j   dx       i-1/2,j        i+1/2,j


              dt    y  n+1/2       y  n+1/2
            + -- { F (U       ) - F (U       ) }
              dy       i,j-1/2        i,j+1/2

Updating U_{i,j}:

  • We want to find the state to the left and right (or top and bottom) of each interface, ex. U_{i+1/2,j,[lr]}^{n+1/2}, and use them to solve a Riemann problem across each of the four interfaces.

  • U_{i+1/2,j,[lr]}^{n+1/2} is comprised of two parts, the computation of the monotonized central differences in the normal direction (eqs. 2.8, 2.10) and the computation of the transverse derivatives, which requires the solution of a Riemann problem in the transverse direction (eqs. 2.9, 2.14).

    • the monotonized central difference part is computed using the primitive variables.

    • We compute the central difference part in both directions before doing the transverse flux differencing, since for the high-order transverse flux implementation, we use these as the input to the transverse Riemann problem.

pyro.compressible.unsplit_fluxes.apply_artificial_viscosity(F_x, F_y, q, my_data, rp, ivars)[source]#

This applies artificial viscosity correction terms to the fluxes.

Parameters:
F_x, F_yndarray, ndarray

Fluxes in x and y interface.

qndarray

Primitive variables

my_dataCellCenterData2d object

The data object containing the grid and advective scalar that we are advecting.

rpRuntimeParameters object

The runtime parameters for the simulation

ivarsVariables object

The Variables object that tells us which indices refer to which variables

dtfloat

The timestep we are advancing through.

Returns:
outndarray, ndarray

Fluxes in x and y interface.

pyro.compressible.unsplit_fluxes.apply_source_terms(U_xl, U_xr, U_yl, U_yr, my_data, my_aux, rp, ivars, tc, dt, *, problem_source=None)[source]#

This function applies source terms including external (gravity), geometric terms, and pressure terms to the left and right interface states (normal conserved states). Both geometric and pressure terms arise purely from geometry.

Parameters:
U_xl, U_xr, U_yl, U_yr: ndarray, ndarray, ndarray, ndarray

Conserved states in the left and right x-interface and left and right y-interface.

my_dataCellCenterData2d object

The data object containing the grid and advective scalar that we are advecting.

my_auxCellCenterData2d object

The data object that carries auxiliary quantities which we need to fill in the ghost cells.

rpRuntimeParameters object

The runtime parameters for the simulation

ivarsVariables object

The Variables object that tells us which indices refer to which variables

tcTimerCollection object

The timers we are using to profile

dtfloat

The timestep we are advancing through.

problem_sourcefunction (optional)

A problem-specific source function to call

Returns:
outndarray, ndarray, ndarray, ndarray

Left and right normal conserved states in x and y interfaces with source terms added.

pyro.compressible.unsplit_fluxes.apply_transverse_flux(U_xl, U_xr, U_yl, U_yr, my_data, rp, ivars, solid, tc, dt)[source]#

This function applies transverse correction terms to the normal conserved states after applying other source terms.

We construct the state perpendicular to the interface by adding the central difference part to the transverse flux difference.

The states that we represent by indices i,j are shown below (1,2,3,4):

j+3/2--+----------+----------+----------+
       |          |          |          |
       |          |          |          |
  j+1 -+          |          |          |
       |          |          |          |
       |          |          |          |    1: U_xl[i,j,:] = U
j+1/2--+----------XXXXXXXXXXXX----------+                      i-1/2,j,L
       |          X          X          |
       |          X          X          |
    j -+        1 X 2        X          |    2: U_xr[i,j,:] = U
       |          X          X          |                      i-1/2,j,R
       |          X    4     X          |
j-1/2--+----------XXXXXXXXXXXX----------+
       |          |    3     |          |    3: U_yl[i,j,:] = U
       |          |          |          |                      i,j-1/2,L
  j-1 -+          |          |          |
       |          |          |          |
       |          |          |          |    4: U_yr[i,j,:] = U
j-3/2--+----------+----------+----------+                      i,j-1/2,R
       |    |     |    |     |    |     |
           i-1         i         i+1
     i-3/2      i-1/2      i+1/2      i+3/2

remember that the fluxes are stored on the left edge, so:

F_x[i,j,:] = F_x
                i-1/2, j

F_y[i,j,:] = F_y
                i, j-1/2
Parameters:
U_xl, U_xr, U_yl, U_yr: ndarray, ndarray, ndarray, ndarray

Conserved states in the left and right x-interface and left and right y-interface.

my_dataCellCenterData2d object

The data object containing the grid and advective scalar that we are advecting.

rpRuntimeParameters object

The runtime parameters for the simulation

ivarsVariables object

The Variables object that tells us which indices refer to which variables

solid: A container class

This is used in Riemann solver to indicate which side has solid boundary

tcTimerCollection object

The timers we are using to profile

dtfloat

The timestep we are advancing through.

Returns:
outndarray, ndarray, ndarray, ndarray

Left and right normal conserved states in x and y interfaces with source terms added.

pyro.compressible.unsplit_fluxes.interface_states(my_data, rp, ivars, tc, dt)[source]#

interface_states returns the normal conserved states in the x and y interfaces. We get the normal fluxes by finding the normal primitive states, Then construct the corresponding conserved states.

Parameters:
my_dataCellCenterData2d object

The data object containing the grid and advective scalar that we are advecting.

rpRuntimeParameters object

The runtime parameters for the simulation

ivarsVariables object

The Variables object that tells us which indices refer to which variables

tcTimerCollection object

The timers we are using to profile

dtfloat

The timestep we are advancing through.

Returns:
outndarray, ndarray, ndarray, ndarray

Left and right normal conserved states in x and y interfaces