swe package

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

Subpackages

• swe.problems package
  • Submodules
  • swe.problems.acoustic_pulse module
  • swe.problems.advect module
  • swe.problems.dam module
  • swe.problems.kh module
  • swe.problems.logo module
  • swe.problems.quad module
  • swe.problems.test module

Submodules

swe.derives module

swe.derives.derive_primitives(myd, varnames)

derive desired primitive variables from conserved state

swe.simulation module

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

Bases: simulation_null.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 swe.simulation.Variables(myd)

Bases: object

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

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

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

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

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

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.

swe.unsplit_fluxes.unsplit_fluxes(my_data, my_aux, 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_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