Demonstrations

This section show simple demonstrations of how to use Flooddrake to solve the shallow water equations. Actual examples of problems that can be solved using Flooddrake can be found in the examples repository.

At the start of each example script, the following must be typed:

from flooddrake import *

Then, the following commands initialize a 2D dam break problem, by setting up a state vector (depth, x momentum, y momentum) and a topography (spatially uniform at 0 for this problem) field. The dam is given as a circle of radius sqrt(0.05) and centred at (0.5, 0.5) in a unit square grid.

# Meshsize
n = 10
mesh = UnitSquareMesh(n, n)

# mixed function space
v_h = FunctionSpace(mesh, "DG", 1)
v_mu = FunctionSpace(mesh, "DG", 1)
v_mv = FunctionSpace(mesh, "DG", 1)
V = v_h * v_mu * v_mv

# setup free surface depth
g = Function(V)
x = SpatialCoordinate(V.mesh())
g.sub(0).interpolate(conditional(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2) < 0.05,
                                 1.0,
                                 0.8))

# setup bed
bed = Function(V)

# setup state
state = State(V, g, bed)

# setup source
source = Function(v_h)

A few things to note from this is that state.w represents the state vector with a modified water depth that allowed g to be the free-surface depth. Basically, state.w.sub(0) represents (g.sub(0) - bed.sub(0)) truncated at 0 (to avoid negative water depth). source represents the source term in the shallow water equations (e.g. rainfall) and is just a Function of the depth FunctionSpace.

Now we can integrate the problem over the time interval [0, 0.75]. We set the maximum time-step to 0.1 and save the final state vector as final_w.

# timestep
solution = Timestepper(V, state.bed, source, MaxTimestep=0.1)

t_visualisation = 0.1
solution.stepper(0, 0.75, state.w, t_visualisation)

t_visualisation is an argument describing at which time points we would like the free-surface depth to be written to the h.pvd file created by Flooddrake. One can define boundary conditions and use them as an argument in the Timestepper class. See help(BoundaryConditions) for more information.

Visualising the second to last frame in h.pvd (at time 0.7) we get: