pyro.swe package

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

Subpackages

• pyro.swe.problems package
  • Submodules
  • pyro.swe.problems.acoustic_pulse module
    • finalize()
    • init_data()
  • pyro.swe.problems.advect module
    • finalize()
    • init_data()
  • pyro.swe.problems.dam module
    • finalize()
    • init_data()
  • pyro.swe.problems.kh module
    • finalize()
    • init_data()
  • pyro.swe.problems.logo module
    • finalize()
    • init_data()
  • pyro.swe.problems.quad module
    • finalize()
    • init_data()
  • pyro.swe.problems.test module
    • finalize()
    • init_data()

Submodules

pyro.swe.derives module

pyro.swe.derives.derive_primitives(myd, varnames)

derive desired primitive variables from conserved state

pyro.swe.interface module

pyro.swe.interface.consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_state)

Calculate the conserved flux for the shallow water equations. In the x-direction, this is given by:

    /      hu       \
F = | hu^2 + gh^2/2 |
    \      huv      /

Parameters:

idirint

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

gfloat

Graviational acceleration

ih, ixmom, iymom, ihXint

The indices of the height, x-momentum, y-momentum, height*species fraction in the conserved state vector.

nspecint

The number of species

U_statendarray

Conserved state vector.

Returns:

outndarray

Conserved flux

pyro.swe.interface.riemann_hllc(idir, ng, ih, ixmom, iymom, ihX, nspec, lower_solid, upper_solid, g, U_l, U_r)

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

ih, ixmom, iymom, ihXint

The indices of the height, x-momentum, y-momentum and height*species fractions in the conserved state vector.

nspecint

The number of species

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gfloat

Gravitational acceleration

U_l, U_rndarray

Conserved state on the left and right cell edges.

Returns:

outndarray

Conserved flux

pyro.swe.interface.riemann_roe(idir, ng, ih, ixmom, iymom, ihX, nspec, lower_solid, upper_solid, g, U_l, U_r)

This is the Roe Riemann solver with entropy fix. The implementation follows Toro’s SWE book and the clawpack 2d SWE Roe solver.

Parameters:

idirint

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

ngint

The number of ghost cells

ih, ixmom, iymom, ihXint

The indices of the height, x-momentum, y-momentum and height*species fractions in the conserved state vector.

nspecint

The number of species

lower_solid, upper_solidint

Are we at lower or upper solid boundaries?

gfloat

Gravitational acceleration

U_l, U_rndarray

Conserved state on the left and right cell edges.

Returns:

outndarray

Conserved flux

pyro.swe.interface.states(idir, ng, dx, dt, ih, iu, iv, ix, nspec, g, qv, dqv)

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   0   0 \
    | g   u   0 |
A = \ 0   0   u /

The right eigenvectors are:

     /  h  \       /  0  \      /  h  \
r1 = | -c  |  r2 = |  0  | r3 = |  c  |
     \  0  /       \  1  /      \  0  /

The left eigenvectors are:

l1 =     ( 1/(2h),  -h/(2hc),  0 )
l2 =     ( 0,          0,  1 )
l3 =     ( -1/(2h), -h/(2hc),  0 )

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

dxfloat

The cell spacing

dtfloat

The timestep

ih, iu, iv, ixint

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

nspecint

The number of species

gfloat

Gravitational acceleration

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.swe.simulation module

class pyro.swe.simulation.Simulation(solver_name, problem_name, rp, timers=None, data_class=<class 'pyro.mesh.patch.CellCenterData2d'>)

Bases: NullSimulation

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

dovis()

Do runtime visualization.

evolve()

Evolve the equations of swe hydrodynamics through a timestep dt.

initialize(extra_vars=None, ng=4)

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

method_compute_timestep()

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.

class pyro.swe.simulation.Variables(myd)

Bases: object

a container class for easy access to the different swe variables by an integer key

pyro.swe.simulation.cons_to_prim(U, ivars, myg)

Convert an input vector of conserved variables \(U = (h, hu, hv, {hX})\) to primitive variables \(q = (h, u, v, {X})\).

pyro.swe.simulation.prim_to_cons(q, ivars, myg)

Convert an input vector of primitive variables \(q = (h, u, v, {X})\) to conserved variables \(U = (h, hu, hv, {hX})\)

pyro.swe.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 Roe-fix

• 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.swe.unsplit_fluxes.unsplit_fluxes(my_data, rp, ivars, solid, tc, dt)

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

The runtime parameter g is assumed to be the gravitational acceleration in the y-direction

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

varsVariables 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

The fluxes on the x- and y-interfaces