compressible package¶
The pyro compressible hydrodynamics solver. This implements the second-order (piecewise-linear), unsplit method of Colella 1990.
Subpackages¶
- compressible.problems package
- Submodules
- compressible.problems.acoustic_pulse module
- compressible.problems.advect module
- compressible.problems.bubble module
- compressible.problems.hse module
- compressible.problems.kh module
- compressible.problems.logo module
- compressible.problems.quad module
- compressible.problems.ramp module
- compressible.problems.rt module
- compressible.problems.rt2 module
- compressible.problems.sedov module
- compressible.problems.sod module
- compressible.problems.test module
Submodules¶
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
-
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
- ccdata : CellCenterData2d object
The data object
compressible.derives module¶
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.
-
compressible.eos.
dens
(gamma, pres, eint)[source]¶ Given the pressure and the specific internal energy, return the density
Parameters: - gamma : float
The ratio of specific heats
- pres : float
The pressure
- eint : float
The specific internal energy
Returns: - out : float
The density
compressible.interface module¶
-
compressible.interface.
artificial_viscosity
[source]¶ Compute the artifical 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 ofavisco_x[i,j]
Y
is the location ofavisco_y[i,j]
Parameters: - ng : int
The number of ghost cells
- dx, dy : float
Cell spacings
- cvisc : float
viscosity parameter
- u, v : ndarray
x- and y-velocities
Returns: - out : ndarray, ndarray
Artificial viscosity in the x- and y-directions
-
compressible.interface.
consFlux
[source]¶ Calculate the conservative flux.
Parameters: - idir : int
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
- gamma : float
Adiabatic index
- idens, ixmom, iymom, iener, irhoX : int
The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.
- nspec : int
The number of species
- U_state : ndarray
Conserved state vector.
Returns: - out : ndarray
Conserved flux
-
compressible.interface.
riemann_cgf
[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: - idir : int
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
- ng : int
The number of ghost cells
- nspec : int
The number of species
- idens, ixmom, iymom, iener, irhoX : int
The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.
- lower_solid, upper_solid : int
Are we at lower or upper solid boundaries?
- gamma : float
Adiabatic index
- U_l, U_r : ndarray
Conserved state on the left and right cell edges.
Returns: - out : ndarray
Conserved flux
-
compressible.interface.
riemann_hllc
[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: - idir : int
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
- ng : int
The number of ghost cells
- nspec : int
The number of species
- idens, ixmom, iymom, iener, irhoX : int
The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector.
- lower_solid, upper_solid : int
Are we at lower or upper solid boundaries?
- gamma : float
Adiabatic index
- U_l, U_r : ndarray
Conserved state on the left and right cell edges.
Returns: - out : ndarray
Conserved flux
-
compressible.interface.
riemann_prim
[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: - idir : int
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
- ng : int
The number of ghost cells
- nspec : int
The number of species
- irho, iu, iv, ip, iX : int
The indices of the density, x-velocity, y-velocity, pressure and species fractions in the state vector.
- lower_solid, upper_solid : int
Are we at lower or upper solid boundaries?
- gamma : float
Adiabatic index
- q_l, q_r : ndarray
Primitive state on the left and right cell edges.
Returns: - out : ndarray
Primitive flux
-
compressible.interface.
states
[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 andV_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: - idir : int
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
- ng : int
The number of ghost cells
- dx : float
The cell spacing
- dt : float
The timestep
- irho, iu, iv, ip, ix : int
Indices of the density, x-velocity, y-velocity, pressure and species in the state vector
- nspec : int
The number of species
- gamma : float
Adiabatic index
- qv : ndarray
The primitive state vector
- dqv : ndarray
Spatial derivitive of the state vector
Returns: - out : ndarray, ndarray
State vector predicted to the left and right edges
compressible.simulation module¶
-
class
compressible.simulation.
Simulation
(solver_name, problem_name, rp, timers=None, data_class=<class 'mesh.patch.CellCenterData2d'>)[source]¶ Bases:
simulation_null.NullSimulation
The main simulation class for the corner transport upwind compressible hydrodynamics solver
-
initialize
(extra_vars=None, ng=4)[source]¶ Initialize the grid and variables for compressible flow and set the initial conditions for the chosen problem.
-
-
class
compressible.simulation.
Variables
(myd)[source]¶ Bases:
object
a container class for easy access to the different compressible variable by an integer key
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
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.
-
compressible.unsplit_fluxes.
unsplit_fluxes
(my_data, my_aux, rp, ivars, solid, tc, dt)[source]¶ unsplitFluxes returns the fluxes through the x and y interfaces by doing an unsplit reconstruction of the interface values and then solving the Riemann problem through all the interfaces at once
currently we assume a gamma-law EOS
The runtime parameter grav is assumed to be the gravitational acceleration in the y-direction
Parameters: - my_data : CellCenterData2d object
The data object containing the grid and advective scalar that we are advecting.
- rp : RuntimeParameters object
The runtime parameters for the simulation
- vars : Variables object
The Variables object that tells us which indices refer to which variables
- tc : TimerCollection object
The timers we are using to profile
- dt : float
The timestep we are advancing through.
Returns: - out : ndarray, ndarray
The fluxes on the x- and y-interfaces