# Simple Script Example

Here we discuss the structure and operation of a
script called `runup.py`

(which is available in the `examples`

directory of `anuga_core`

.

This example carries out the solution of the shallow-water wave equation in the simple case of a configuration comprising a flat bed, sloping at a fixed angle in one direction and having a constant depth across each line in the perpendicular direction.

The example demonstrates the basic ideas involved in setting up a complex scenario. In general the user specifies the geometry (bathymetry and topography), the initial water level, boundary conditions such as tide, and any forcing terms that may drive the system such as rainfall, abstraction of water, wind stress or atmospheric pressure gradients. Frictional resistance from the different terrains in the model is represented by predefined forcing terms. In this example, the boundary is reflective on three sides and a time dependent wave on one side.

The present example represents a simple scenario and does not include any forcing terms, nor is the data taken from a file as it would typically be.

The conserved quantities involved in the problem are stage (absolute height of water surface), \(x\)-momentum and \(y\)-momentum. Other quantities involved in the computation are the friction and elevation.

Water depth can be obtained through the equation:

```
depth = stage - elevation
```

## Outline of the Program

In outline, `runup.py`

performs the following steps:

Sets up a triangular mesh.

Sets certain parameters governing the mode of operation of the model, specifying, for instance, where to store the model output.

Inputs various quantities describing physical measurements, such as the elevation, to be specified at each mesh point (vertex).

Sets up the boundary conditions.

Carries out the evolution of the model through a series of time steps and outputs the results, providing a results file that can be viewed.

## The Code

For reference we include below the complete code listing for
`runup.py`

. Subsequent paragraphs provide a
‘commentary’ that describes each step of the program and explains it
significance.

```
"""Simple water flow example using ANUGA
Water driven up a linear slope and time varying boundary,
similar to a beach environment
"""
#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
import anuga
from math import sin, pi, exp
#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
domain = anuga.rectangular_cross_domain(10, 10) # Create domain
#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography(x, y):
return -x/2 # linear bed slope
#return x*(-(2.0-x)*.5) # curved bed slope
domain.set_quantity('elevation', topography) # Use function for elevation
domain.set_quantity('friction', 0.1) # Constant friction
domain.set_quantity('stage', -0.4) # Constant negative initial stage
#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
Br = anuga.Reflective_boundary(domain) # Solid reflective wall
Bw = anuga.Time_boundary(domain=domain, # Time dependent boundary
function=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
# Associate boundary tags with boundary objects
domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep=0.1, finaltime=10.0):
domain.print_timestepping_statistics()
```

## Establishing the Domain

The very first thing to do is import the various modules, of which the
`anuga`

module is the most important:

```
import anuga
```

Then we need to set up the triangular mesh to be used for the scenario. This is carried out through the statement:

```
domain = anuga.rectangular_cross_domain(10, 5, len1=10.0, len2=5.0)
```

The above assignment sets up a \(10 \times 5\) rectangular mesh,
triangulated in a regular way with boundary tags
`left`

, `right`

, `top`

or `bottom`

.

It is also possible to set up a domain from “first principles”
using `points`

, `vertices`

and `boundary`

via the assignment:

```
points, vertices, boundary = anuga.rectangular_cross(10, 5, len1=10.0, len2=5.0)
domain = anuga.Domain(points, vertices, boundary)
```

- where:
`points`

is a list giving the coordinates of each mesh point,`vertices`

is a list specifying the three vertices of each triangle, and`boundary`

is a dictionary that stores the edges on the boundary and associates with each a symbolic tag. The edges are represented as pairs (i, j) where i refers to the triangle id and j to the edge id of that triangle. Edge ids are enumerated from 0 to 2 based on the id of the vertex opposite.

This creates an instance of the `Domain`

class, which
represents the domain of the simulation. Specific options are set at
this point, including the basename for the output file and the
directory to be used for data:

```
domain.set_name('runup')
domain.set_datadir('.')
```

In addition, the following statement could be used to state that
quantities `stage`

, `xmomentum`

and `ymomentum`

are
to be stored at every timestep and `elevation`

only once at
the beginning of the simulation:

```
domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
```

However, this is not necessary, as the above is the default behaviour.

## Initial Conditions

The next task is to specify a number of quantities that we wish to
set for each mesh point. The class `{Domain`

has a method
`set_quantity`

, used to specify these quantities. It is a
flexible method that allows the user to set quantities in a variety
of ways – using constants, functions, numeric arrays, expressions
involving other quantities, or arbitrary data points with associated
values, all of which can be passed as arguments. All quantities can
be initialised using `set_quantity`

. For a conserved
quantity (such as `stage, xmomentum, ymomentum`

) this is called
an *initial condition*. However, other quantities that aren’t
updated by the evolution procedure are also assigned values using the same
interface. The code in the present example demonstrates a number of
forms in which we can invoke `set_quantity`

.

### Elevation

The elevation, or height of the bed, is set using a function defined through the statements below, which is specific to this example and specifies a particularly simple initial configuration for demonstration purposes:

```
def topography(x, y):
return -x/2
```

This simply associates an elevation with each point `(x, y)`

of
the plane. It specifies that the bed slopes linearly in the
`x`

direction, with slope \(-\frac{1}{2}\), and is constant in
the `y`

direction.

Once the function `topography`

is specified, the quantity
`elevation`

is assigned through the statement:

```
domain.set_quantity('elevation', topography)
```

Note

If using a function to set a quantity such as `elevation`

it must be vector
compatible. For example, using `math.sqrt`

will not work, but `numpy.sqrt`

will work.

### Friction

The assignment of the friction quantity (a forcing term)
demonstrates another way we can use `set_quantity`

to set
quantities – namely, assign them to a constant numerical value:

```
domain.set_quantity('friction', 0.1)
```

This specifies that the Manning friction coefficient is set to 0.1 at every mesh point.

### Stage

The stage (the height of the water surface) is related to the elevation and the depth at any time by the equation:

```
stage = elevation + depth
```

For this example, we simply assign a constant value to `stage`

,
using the statement:

```
domain.set_quantity('stage', -0.4)
```

which specifies that the surface level is set to a height of \(-0.4\), i.e. 0.4 units (metres) below the zero level.

Although it is not necessary for this example, it may be useful to
digress here and mention a variant to this requirement, which allows
us to illustrate another way to use `set_quantity`

– namely,
incorporating an expression involving other quantities. Suppose,
instead of setting a constant value for the stage, we wished to
specify a constant value for the *depth*. For such a case we
need to specify that `stage`

is everywhere obtained by adding
that value to the value already specified for `elevation`

. We
would do this by means of the statements:

```
h = 0.05 # Constant depth
domain.set_quantity('stage', expression='elevation + %f' % h)
```

That is, the value of `stage`

is set to `h = 0.05`

plus
the value of `elevation`

already defined.

The reader will probably appreciate that this capability to
incorporate expressions into statements using `set_quantity`

greatly expands its power.

## Boundary Conditions

The boundary conditions are specified as follows:

```
Br = anuga.Reflective_boundary(domain)
Bw = anuga.Time_boundary(domain=domain, f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
```

The effect of these statements is to set up a selection of different alternative boundary conditions and store them in variables that can be assigned as needed. Each boundary condition specifies the behaviour at a boundary in terms of the behaviour in neighbouring elements. The boundary conditions introduced here may be briefly described as follows:

Reflective boundary: Returns same

`stage`

as in its neighbour volume but momentum vector reversed 180 degrees (reflected). Specific to the shallow water equation as it works with the momentum quantities assumed to be the second and third conserved quantities. A reflective boundary condition models a solid wall.Time boundary: Set a boundary varying with time.

Before describing how these boundary conditions are assigned,
we recall that a mesh is specified using
three variables `points`

, `vertices`

and `boundary`

.
In the code we are discussing, these three variables are returned by
the function `rectangular_cross`

. The example given in
Section ref{sec:realdataexample} illustrates another way of
assigning the values, by means of the function
`create_domain_from_regions`

.

These variables store the data determining the mesh as follows. (You
may find that the example given in Section ref{sec:meshexample}
helps to clarify the following discussion, even though that example
is a *non-rectangular* mesh.):

`points`

stores a list of 2-tuples giving the coordinates of the mesh points.

`vertices`

stores a list of 3-tuples of numbers, representing vertices of triangles in the mesh. In this list, the triangle whose vertices are`points[i]`

,`points[j]`

,`points[k]`

is represented by the 3-tuple`(i, j, k)`

.The variable

`boundary`

is a Python dictionary that not only stores the edges that make up the boundary but also assigns symbolic tags to these edges to distinguish different parts of the boundary. An edge with endpoints`points[i]`

and`points[j]`

is represented by the 2-tuple`(i, j)`

. The keys for the dictionary are the 2-tuples`(i, j)`

corresponding to boundary edges in the mesh, and the values are the tags are used to label them. In the present example, the value`boundary[(i, j)]`

assigned to`(i, j)]`

is one of the four tags`left`

,`right`

,`top`

or`bottom`

, depending on whether the boundary edge represented by`(i, j)`

occurs at the left, right, top or bottom of the rectangle bounding the mesh. The function`rectangular_cross`

automatically assigns these tags to the boundary edges when it generates the mesh.

The tags provide the means to assign different boundary conditions
to an edge depending on which part of the boundary it belongs to.
(In Section Real Example we describe an example that
uses different boundary tags – in general, the possible tags are entirely
selectable by the user when generating the mesh and not
limited to ‘left’, ‘right’, ‘top’ and ‘bottom’ as in this example.)
All segments in bounding polygon must be tagged. If a tag is not supplied,
the default tag name `exterior`

will be assigned by ANUGA.

Using the boundary objects described above, we assign a boundary condition to each part of the boundary by means of a statement like:

```
domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
```

It is critical that all tags are associated with a boundary condition in this statement. If not the program will halt with a statement like:

```
Traceback (most recent call last):
File "mesh_test.py", line 114, in ?
domain.set_boundary({'west': Bi, 'east': Bo, 'north': Br, 'south': Br})
File "X:\inundation\sandpits\onielsen\anuga_core\source\anuga\
abstract_2d_finite_volumes\domain.py", line 505, in set_boundary
raise msg
ERROR (domain.py): Tag "exterior" has not been bound to a boundary object.
All boundary tags defined in domain must appear in the supplied dictionary.
The tags are: ['ocean', 'east', 'north', 'exterior', 'south']
```

The command `set_boundary`

stipulates that, in the current example, the right
boundary varies with time, as defined by the lambda function, while the other
boundaries are all reflective.

## Evolution

The final statement:

```
for t in domain.evolve(yieldstep=0.1, duration=10.0):
domain.print_timestepping_statistics()
```

causes `domain`

we have just setup to *evolve*, over a series of
steps indicated by the values of `yieldstep`

and
`duration`

, which can be altered as required (an alternative
to `duration`

is `finaltime`

) The value of
`yieldstep`

provides the time interval between successive yields to the evolve loop.
Behind the scenes more *inner* time steps are generally taken.

By default, the current state of the evolution is stored a each yield step.

Time between output can also be controlled by
the argument `outputstep`

which needs to an integer multiple of the `yieldstep`

## Output

The output is a NetCDF file with the extension `.sww`

. It
contains stage and momentum information and can be used with the
ANUGA viewer `anuga_viewer`

to generate a visual
display.

## Exploring the Model Output

The following figures are screenshots from the anuga viewer visualisation
tool `anuga_viewer`

.

The first figure shows the domain with water surface as specified by the initial condition, \(t=0\).

The second figure shows the flow at time \(t=2.3\) and the last figure show the flow at time \(t=4\) where the system has been evolved and the wave is encroaching on the previously dry bed.

Online documentation is available
for the `anuga_viewer`