source: documentation/user_manual/anuga_user_manual.tex @ 3241

% Complete documentation on the extended LaTeX markup used for Python
% documentation is available in ``Documenting Python'', which is part
% of the standard documentation for Python.  It may be found online
% at:
%
%     http://www.python.org/doc/current/doc/doc.html


%labels
%Sections and subsections \label{sec: }
%Chapters \label{ch: }
%Equations \label{eq: }
%Figures \label{fig: }



\documentclass{manual}

\usepackage{graphicx}
\usepackage{datetime}

\input{definitions}

\title{\anuga User Manual}
\author{Geoscience Australia and the Australian National University}

% Please at least include a long-lived email address;
% the rest is at your discretion.
\authoraddress{Geoscience Australia \\
  Email: \email{ole.nielsen@ga.gov.au}
}

%Draft date

% update before release!
% Use an explicit date so that reformatting
% doesn't cause a new date to be used.  Setting
% the date to \today can be used during draft
% stages to make it easier to handle versions.


\longdate       % Make date format long using datetime.sty
%\settimeformat{xxivtime} % 24 hour Format
\settimeformat{oclock} % Verbose
\date{\today, \ \currenttime}
%\hyphenation{set\_datadir}

\ifhtml
\date{\today} % latex2html does not know about datetime
\fi





\release{1.0}   % release version; this is used to define the
                % \version macro

\makeindex          % tell \index to actually write the .idx file
\makemodindex       % If this contains a lot of module sections.

\setcounter{tocdepth}{3}
\setcounter{secnumdepth}{3}


\begin{document}
\maketitle


% This makes the contents more accessible from the front page of the HTML.
\ifhtml
\chapter*{Front Matter\label{front}}
\fi

%Subversion keywords:
%
%$LastChangedDate: 2006-06-19 05:58:37 +0000 (Mon, 19 Jun 2006) $
%$LastChangedRevision: 3179 $
%$LastChangedBy: sexton $

\input{copyright}


\begin{abstract}
\label{def:anuga}

\noindent \anuga\index{\anuga} is a hydrodynamic modelling tool that
allows users to model realistic flow problems in complex geometries.
Examples include dam breaks or the effects of natural hazards such
as riverine flooding, storm surges and tsunami.

The user must specify a study area represented by a mesh of triangular
cells, the topography and bathymetry, frictional resistance, initial
values for water level (called \emph{stage}\index{stage} within \anuga),
boundary
conditions and forces such as windstress or pressure gradients if
applicable.

\anuga tracks the evolution of water depth and horizontal momentum
within each cell over time by solving the shallow water wave equation
governing equation using a finite-volume method.

\anuga cannot model details of breaking waves, flow under ceilings such
as pipes, turbulence and vortices, vertical convection or viscous
flows.

\anuga also incorporates a mesh generator, called \code{pmesh}, that
allows the user to set up the geometry of the problem interactively as
well as tools for interpolation and surface fitting, and a number of
auxiliary tools for visualising and interrogating the model output.

Most \anuga components are written in the object-oriented programming
language Python and most users will interact with \anuga by writing
small Python programs based on the \anuga library
functions. Computationally intensive components are written for
efficiency in C routines working directly with the Numerical Python
structures.


\end{abstract}

\tableofcontents


\chapter{Introduction}


\section{Purpose}

The purpose of this user manual is to introduce the new user to the
inundation software, describe what it can do and give step-by-step
instructions for setting up and running hydrodynamic simulations.

\section{Scope}

This manual covers only what is needed to operate the software after
installation and configuration. It does not includes instructions
for installing the software or detailed API documentation, both of
which will be covered in separate publications and by documentation
in the source code.

\section{Audience}

Readers are assumed to be familiar with the operating environment
and have a general understanding of the subject matter, as well as
enough programming experience to adapt the code to different
requirements and to understand the basic terminology of
object-oriented programming.

\pagebreak
\chapter{Background}


Modelling the effects on the built environment of natural hazards such
as riverine flooding, storm surges and tsunami is critical for
understanding their economic and social impact on our urban
communities.  Geoscience Australia and the Australian National
University are developing a hydrodynamic inundation modelling tool
called \anuga to help simulate the impact of these hazards.

The core of \anuga is the fluid dynamics module, called pyvolution,
which is based on a finite-volume method for solving the shallow water
wave equation.  The study area is represented by a mesh of triangular
cells.  By solving the governing equation within each cell, water
depth and horizontal momentum are tracked over time.

A major capability of pyvolution is that it can model the process of
wetting and drying as water enters and leaves an area.  This means
that it is suitable for simulating water flow onto a beach or dry land
and around structures such as buildings.  Pyvolution is also capable
of modelling hydraulic jumps due to the ability of the finite-volume
method to accommodate discontinuities in the solution.

To set up a particular scenario the user specifies the geometry
(bathymetry and topography), the initial water level (stage),
boundary conditions such as tide, and any forcing terms that may
drive the system such as wind stress or atmospheric pressure
gradients. Gravity and frictional resistance from the different
terrains in the model are represented by predefined forcing terms.

A mesh generator, called pmesh, allows the user to set up the geometry
of the problem interactively and to identify boundary segments and
regions using symbolic tags.  These tags may then be used to set the
actual boundary conditions and attributes for different regions
(e.g. the Manning friction coefficient) for each simulation.

Most \anuga components are written in the object-oriented programming
language Python.  Software written in Python can be produced quickly
and can be readily adapted to changing requirements throughout its
lifetime.  Computationally intensive components are written for
efficiency in C routines working directly with the Numerical Python
structures.  The animation tool developed for \anuga is based on
OpenSceneGraph, an Open Source Software (OSS) component allowing high
level interaction with sophisticated graphics primitives.
See \cite{nielsen2005} for more background on \anuga.

\chapter{Restrictions and limitations on \anuga}

Although a powerful and flexible tool for hydrodynamic modelling, \anuga has a
number of limitations that any potential user need to be aware of. They are

\begin{itemize}
  \item The mathematical model is the 2D shallow water wave equation. 
  As such it cannot resolve vertical convection and consequently not breaking 
  waves or 3D turbulence (e.g.\ vorticity).
  \item The finite volume is a very robust and flexible numerical technique, 
  but it is not the fastest method around. If the geometry is sufficiently 
  simple and if there is no need for wetting or drying, a finite-difference 
  method may be able to solve the problem faster than \anuga.
  %\item Mesh resolutions near coastlines with steep gradients need to be...   
  \item Frictional resistance is implemented using Manning's formula, but 
  \anuga has not yet been validated in regard to bottom roughness
\end{itemize}



\chapter{Getting Started}
\label{ch:getstarted}

This section is designed to assist the reader to get started with
\anuga by working through some examples. Two examples are discussed;
the first is a simple example to illustrate many of the ideas, and
the second is a more realistic example.

\section{A Simple Example}
\label{sec:simpleexample}

\subsection{Overview}

What follows is a discussion of the structure and operation of a
script called \file{runup.py}.

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 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

\begin{tabular}{rcrcl}
  \code{depth} &=& \code{stage} &$-$& \code{elevation}
\end{tabular}


\subsection{Outline of the Program}

In outline, \file{runup.py} performs the following steps:

\begin{enumerate}

   \item Sets up a triangular mesh.

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

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

   \item Sets up the boundary conditions.

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

\end{enumerate}

\subsection{The Code}

%FIXME: we are using the \code function here.
%This should be used wherever possible
For reference we include below the complete code listing for
\file{runup.py}. Subsequent paragraphs provide a
`commentary' that describes each step of the program and explains it
significance.

\verbatiminput{examples/runup.py}
%\verbatiminput{examples/bedslope.py}

\subsection{Establishing the Mesh}\index{mesh, establishing}

The first task is to set up the triangular mesh to be used for the
scenario. This is carried out through the statement:

{\small \begin{verbatim}
    points, vertices, boundary = rectangular(10, 10)
\end{verbatim}}

The function \function{rectangular} is imported from a module
\module{mesh\_factory} defined elsewhere. (\anuga also contains
several other schemes that can be used for setting up meshes, but we
shall not discuss these.) The above assignment sets up a $10 \times
10$ rectangular mesh, triangulated in a regular way. The assignment

{\small \begin{verbatim}
    points, vertices, boundary = rectangular(m, n)
\end{verbatim}}

returns:

\begin{itemize}

   \item a list \code{points} giving the coordinates of each mesh point,

   \item a list \code{vertices} specifying the three vertices of each triangle, and

   \item a dictionary \code{boundary} that stores the edges on
   the boundary and associates each with one of the symbolic tags \code{`left'}, \code{`right'},
   \code{`top'} or \code{`bottom'}.

\end{itemize}

(For more details on symbolic tags, see page
\pageref{ref:tagdescription}.)

An example of a general unstructured mesh and the associated data
structures \code{points}, \code{vertices} and \code{boundary} is
given in Section \ref{sec:meshexample}.




\subsection{Initialising the Domain}

These variables are then used to set up a data structure
\code{domain}, through the assignment:

{\small \begin{verbatim}
    domain = Domain(points, vertices, boundary)
\end{verbatim}}

This creates an instance of the \class{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:

{\small \begin{verbatim}
    domain.set_name('bedslope')
\end{verbatim}}

{\small \begin{verbatim}
    domain.set_datadir('.')
\end{verbatim}}

In addition, the following statement now specifies that the
quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are
to be stored:

{\small \begin{verbatim}
    domain.set_quantities_to_be_stored(['stage', 'xmomentum',
    'ymomentum'])
\end{verbatim}}


\subsection{Initial Conditions}

The next task is to specify a number of quantities that we wish to
set for each mesh point. The class \class{Domain} has a method
\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 \method{set\_quantity}. For a conserved
quantity (such as \code{stage, xmomentum, ymomentum}) this is called
an \emph{initial condition}. However, other quantities that aren't
updated by the equation are also assigned values using the same
interface. The code in the present example demonstrates a number of
forms in which we can invoke \method{set\_quantity}.


\subsubsection{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:

{\small \begin{verbatim}
    def f(x,y):
        return -x/2
\end{verbatim}}

This simply associates an elevation with each point \code{(x, y)} of
the plane.  It specifies that the bed slopes linearly in the
\code{x} direction, with slope $-\frac{1}{2}$,  and is constant in
the \code{y} direction.

Once the function \function{f} is specified, the quantity
\code{elevation} is assigned through the simple statement:

{\small \begin{verbatim}
    domain.set_quantity('elevation', f)
\end{verbatim}}


\subsubsection{Friction}

The assignment of the friction quantity (a forcing term)
demonstrates another way we can use \method{set\_quantity} to set
quantities---namely, assign them to a constant numerical value:

{\small \begin{verbatim}
    domain.set_quantity('friction', 0.1)
\end{verbatim}}

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

\subsubsection{Stage}

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

{\small \begin{verbatim}
    stage = elevation + depth
\end{verbatim}}


For this example, we simply assign a constant value to \code{stage},
using the statement

{\small \begin{verbatim}
    domain.set_quantity('stage', -.4)
\end{verbatim}}

which specifies that the surface level is set to a height of $-0.4$,
i.e. 0.4 units (m) 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 \method{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 \emph{depth}. For such a case we
need to specify that \code{stage} is everywhere obtained by adding
that value to the value already specified for \code{elevation}. We
would do this by means of the statements:

{\small \begin{verbatim}
    h = 0.05 # Constant depth
    domain.set_quantity('stage', expression = 'elevation + %f' %h)
\end{verbatim}}

That is, the value of \code{stage} is set to $\code{h} = 0.05$ plus
the value of \code{elevation} already defined.

The reader will probably appreciate that this capability to
incorporate expressions into statements using \method{set\_quantity}
greatly expands its power.) See Section \ref{sec:Initial Conditions} for more
details.

\subsection{Boundary Conditions}\index{boundary conditions}

The boundary conditions are specified as follows:

{\small \begin{verbatim}
    Br = Reflective_boundary(domain)

    Bt = Transmissive_boundary(domain)

    Bd = Dirichlet_boundary([0.2,0.,0.])

    Bw = Time_boundary(domain=domain,
                f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
\end{verbatim}}

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:

\begin{itemize}
    \item \textbf{Reflective boundary}\label{def:reflective boundary} Returns same \code{stage} as
      as present 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.
    \item \textbf{Transmissive boundary}\label{def:transmissive boundary} Returns same conserved quantities as
      those present in its neighbour volume. This is one way of modelling
      outflow from a domain, but it should be used with caution if flow is
      not steady state as replication of momentum at the boundary
      may cause occasional spurious effects. If this occurs,
      consider using e.g. a Dirichlet boundary condition.
    \item \textbf{Dirichlet boundary}\label{def:dirichlet boundary} Specifies
      constant values for stage, $x$-momentum and $y$-momentum at the boundary.
    \item \textbf{Time boundary}\label{def:time boundary} Like a Dirichlet
      boundary but with behaviour varying with time.
\end{itemize}

\label{ref:tagdescription}Before describing how these boundary
conditions are assigned, we recall that a mesh is specified using
three variables \code{points}, \code{vertices} and \code{boundary}.
In the code we are discussing, these three variables are returned by
the function \code{rectangular}; however, the example given in
Section \ref{sec:realdataexample} illustrates another way of
assigning the values, by means of the function
\code{create\_mesh\_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 \emph{non-rectangular} mesh.)

\begin{itemize}
\item The variable \code{points} stores a list of 2-tuples giving the
coordinates of the mesh points.

\item The variable \code{vertices} stores a list of 3-tuples of
numbers, representing vertices of triangles in the mesh. In this
list, the triangle whose vertices are \code{points[i]},
\code{points[j]}, \code{points[k]} is represented by the 3-tuple
\code{(i, j, k)}.

\item The variable \code{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 \code{points[i]} and
\code{points[j]} is represented by the 2-tuple \code{(i, j)}. The
keys for the dictionary are the 2-tuples \code{(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 \code{boundary[(i,
j)]} assigned to \code{(i, j)]} is one of the four tags
\code{`left'}, \code{`right'}, \code{`top'} or \code{`bottom'},
depending on whether the boundary edge represented by \code{(i, j)}
occurs at the left, right, top or bottom of the rectangle bounding
the mesh. The function \code{rectangular} automatically assigns
these tags to the boundary edges when it generates the mesh.
\end{itemize}

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 \ref{sec:realdataexample} we describe an example that
uses different boundary tags---in general, the possible tags are not
limited to `left', `right', `top' and `bottom', but can be specified
by the user.)

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

{\small \begin{verbatim}
    domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
\end{verbatim}}

This statement 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.

The reader may wish to experiment by varying the choice of boundary
types for one or more of the boundaries. (In the case of \code{Bd}
and \code{Bw}, the three arguments in each case represent the
\code{stage}, $x$-momentum and $y$-momentum, respectively.)

{\small \begin{verbatim}
    Bw = Time_boundary(domain=domain,
                       f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
\end{verbatim}}



\subsection{Evolution}\index{evolution}

The final statement \nopagebreak[3]
{\small \begin{verbatim}
    for t in domain.evolve(yieldstep = 0.1, duration = 4.0):
        print domain.timestepping_statistics()
\end{verbatim}}

causes the configuration of the domain to `evolve', over a series of
steps indicated by the values of \code{yieldstep} and
\code{duration}, which can be altered as required.  The value of
\code{yieldstep} controls the time interval between successive model
outputs.  Behind the scenes more time steps are generally taken.


\subsection{Output}

The output is a NetCDF file with the extension \code{.sww}. It
contains stage and momentum information and can be used with the
\code{swollen} (see Section \ref{sec:swollen}) visualisation package
to generate a visual display. See Section \ref{sec:file formats}
(page \pageref{sec:file formats}) for more on NetCDF and other file
formats.

The following is a listing of the screen output seen by the user
when this example is run:

\verbatiminput{examples/runupoutput.txt}


\section{How to Run the Code}

The code can be run in various ways:

\begin{itemize}
  \item{from a Windows or Unix command line} as in\ \code{python runup.py}
  \item{within the Python IDLE environment}
  \item{within emacs}
  \item{within Windows, by double-clicking the \code{runup.py}
  file.}
\end{itemize}


\section{Exploring the Model Output}

The following figures are screenshots from the \anuga visualisation
tool \code{swollen}. Figure \ref{fig:runupstart} shows the domain
with water surface as specified by the initial condition, $t=0$.
Figure \ref{fig:bedslope2} shows later snapshots for $t=2.3$ and
$t=4$ where the system has been evolved and the wave is encroaching
on the previously dry bed.  All figures are screenshots from an
interactive animation tool called Swollen which is part of \anuga.
Swollen is described in more detail is Section \ref{sec:swollen}.

\begin{figure}[hbt]

  \centerline{\includegraphics[width=75mm, height=75mm]
    {examples/runupstart.eps}}

  \caption{Bedslope example viewed with Swollen}
  \label{fig:runupstart}
\end{figure}


\begin{figure}[hbt]

  \centerline{
    \includegraphics[width=75mm, height=75mm]{examples/runupduring.eps}
    \includegraphics[width=75mm, height=75mm]{examples/runupend.eps}
   }

  \caption{Bedslope example viewed with Swollen}
  \label{fig:bedslope2}
\end{figure}




\clearpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{An Example with Real Data}
\label{sec:realdataexample} The following discussion builds on the
concepts introduced through the \file{runup.py} example and
introduces a second example, \file{run\_sydney.py}.  This refers to
a real-life scenario, in which the domain of interest surrounds the
Sydney region, and predominantly covers Sydney Harbour. A
hypothetical tsunami wave is generated by a submarine mass failure
situated on the edge of the continental shelf.

\subsection{Overview}
As in the case of \file{runup.py}, the actions carried
out by the program can be organised according to this outline:

\begin{enumerate}

   \item Set up a triangular mesh.

   \item Set certain parameters governing the mode of
operation of the model---specifying, for instance, where to store the
model output.

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

   \item Set up the boundary conditions.

   \item Carry out the evolution of the model through a series of time
steps and output the results, providing a results file that can be
visualised.

\end{enumerate}



\subsection{The Code}

Here is the code for \file{run\_sydney\_smf.py}:

\verbatiminput{examples/runsydney.py}

In discussing the details of this example, we follow the outline
given above, discussing each major step of the code in turn.

\subsection{Establishing the Mesh}\index{mesh, establishing}

One obvious way that the present example differs from
\file{runup.py} is in the use of a more complex method to
create the mesh. Instead of imposing a mesh structure on a
rectangular grid, the technique used for this example involves
building mesh structures inside polygons specified by the user,
using a mesh-generator referred to as \code{pmesh}.

In its simplest form, \code{pmesh} creates the mesh within a single
polygon whose vertices are at geographical locations specified by
the user. The user specifies the \emph{resolution}---that is, the
maximal area of a triangle used for triangulation---and a triangular
mesh is created inside the polygon using a mesh generation engine.
On any given platform, the same mesh will be returned. Figure
\ref{fig:pentagon} shows a simple example of this, in which the
triangulation is carried out within a pentagon.


\begin{figure}[hbt]

  \caption{Mesh points are created inside the polygon}
  \label{fig:pentagon}
\end{figure}

Boundary tags are not restricted to \code{`left'}, \code{`right'},
\code{`bottom'} and \code{`top'}, as in the case of
\file{runup.py}. Instead the user specifies a list of
tags appropriate to the configuration being modelled.

In addition, \code{pmesh} provides a way to adapt to geographic or
other features in the landscape, whose presence may require an
increase in resolution. This is done by allowing the user to specify
a number of \emph{interior polygons}, each with a specified
resolution, see Figure \ref{fig:interior meshes}. It is also
possible to specify one or more `holes'---that is, areas bounded by
polygons in which no triangulation is required.

\begin{figure}[hbt]



  \caption{Interior meshes with individual resolution}
  \label{fig:interior meshes}
\end{figure}

In its general form, \code{pmesh} takes for its input a bounding
polygon and (optionally) a list of interior polygons. The user
specifies resolutions, both for the bounding polygon and for each of
the interior polygons. Given this data, \code{pmesh} first creates a
triangular mesh with varying resolution.

The function used to implement this process is
\function{create\_mesh\_from\_regions}. Its arguments include the
bounding polygon and its resolution, a list of boundary tags, and a
list of pairs \code{[polygon, resolution]}, specifying the interior
polygons and their resolutions.

In practice, the details of the polygons used are read from a
separate file \file{project.py}. Here is a complete listing of
\file{project.py}:

\verbatiminput{examples/project.py}

The resulting mesh is output to a \emph{mesh file}\index{mesh
file}\label{def:mesh file}. This term is used to describe a file of
a specific format used to store the data specifying a mesh. (There
are in fact two possible formats for such a file: it can either be a
binary file, with extension \code{.msh}, or an ASCII file, with
extension \code{.tsh}. In the present case, the binary file format
\code{.msh} is used. See Section \ref{sec:file formats} (page
\pageref{sec:file formats}) for more on file formats.)

The statements

{\small \begin{verbatim}
    interior_res = 5000%
    interior_regions = [[project.harbour_polygon_2, interior_res],
                    [project.botanybay_polygon_2, interior_res]]
\end{verbatim}}

are used to read in the specific polygons \code{project.harbour\_polygon\_2} and
\code{botanybay\_polygon\_2} from \file{project.py} and assign a
common resolution of 5000 to each. The statement

{\small \begin{verbatim}
    create_mesh_from_regions(project.diffpolygonall,
                         boundary_tags= {'bottom': [0],
                                         'right1': [1],
                                         'right0': [2],
                                         'right2': [3],
                                         'top': [4],
                                         'left1': [5],
                                         'left2': [6],
                                         'left3': [7]},
                         maximum_triangle_area=100000,
                         filename=meshname,
                         interior_regions=interior_regions)
\end{verbatim}}

is then used to create the mesh, taking the bounding polygon to be
the polygon \code{diffpolygonall} specified in \file{project.py}.
The argument \code{boundary\_tags} assigns a dictionary, whose keys
are the names of the boundary tags used for the bounding
polygon---\code{`bottom'}, \code{`right0'}, \code{`right1'},
\code{`right2'}, \code{`top'}, \code{`left1'}, \code{`left2'} and
\code{`left3'}--- and whose values identify the indices of the
segments associated with each of these tags. (The value associated
with each boundary tag is a one-element list.)


\subsection{Initialising the Domain}

As with \file{runup.py}, once we have created the mesh, the next
step is to create the data structure \code{domain}. We did this for
\file{runup.py} by inputting lists of points and triangles and
specifying the boundary tags directly. However, in the present case,
we use a method that works directly with the mesh file
\code{meshname}, as follows:


{\small \begin{verbatim}
    domain = Domain(meshname, use_cache=True, verbose=True)
\end{verbatim}}

Providing a filename instead of the lists used in \file{runup.py}
above causes \code{Domain} to convert a mesh file \code{meshname}
into an instance of \code{Domain}, allowing us to use methods like
\method{set\_quantity} to set quantities and to apply other
operations.

%(In principle, the
%second argument of \function{pmesh\_to\_domain\_instance} can be any
%subclass of \class{Domain}, but for applications involving the
%shallow-water wave equation, the second argument of
%\function{pmesh\_to\_domain\_instance} can always be set simply to
%\class{Domain}.)

The following statements specify a basename and data directory, and
identify quantities to be stored. For the first two, values are
taken from \file{project.py}.

{\small \begin{verbatim}
    domain.set_name(project.basename)
    domain.set_datadir(project.outputdir)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum',
        'ymomentum'])
\end{verbatim}}


\subsection{Initial Conditions}
Quantities for \file{runsydney.py} are set
using similar methods to those in \file{runup.py}. However,
in this case, many of the values are read from the auxiliary file
\file{project.py} or, in the case of \code{elevation}, from an
ancillary points file.



\subsubsection{Stage}

For the scenario we are modelling in this case, we use a callable
object \code{tsunami\_source}, assigned by means of a function
\function{slump\_tsunami}. This is similar to how we set elevation in
\file{runup.py} using a function---however, in this case the
function is both more complex and more interesting.

The function returns the water displacement for all \code{x} and
\code{y} in the domain. The water displacement is a double Gaussian
function that depends on the characteristics of the slump (length,
thickness, slope, etc), its location (origin) and the depth at that
location.


\subsubsection{Friction}

We assign the friction exactly as we did for \file{runup.py}:

{\small \begin{verbatim}
    domain.set_quantity('friction', 0.0)
\end{verbatim}}


\subsubsection{Elevation}

The elevation is specified by reading data from a file:

{\small \begin{verbatim}
    domain.set_quantity('elevation',
                        filename = project.combineddemname + '.pts',
                        use_cache = True,
                        verbose = True)
\end{verbatim}}

However, before this step can be executed, some preliminary steps
are needed to prepare the file from which the data is taken. Two
source files are used for this data---their names are specified in
the file \file{project.py}, in the variables \code{coarsedemname}
and \code{finedemname}. They contain `coarse' and `fine' data,
respectively---that is, data sampled at widely spaced points over a
large region and data sampled at closely spaced points over a
smaller subregion. The data in these files is combined through the
statement

{\small \begin{verbatim}
combine_rectangular_points_files(project.finedemname + '.pts',
                                 project.coarsedemname + '.pts',
                                 project.combineddemname + '.pts')
\end{verbatim}}

The effect of this is simply to combine the datasets by eliminating
any coarse data associated with points inside the smaller region
common to both datasets. The name to be assigned to the resulting
dataset is also derived from the name stored in the variable
\code{combinedname} in the file \file{project.py}.

\subsection{Boundary Conditions}\index{boundary conditions}

Setting boundaries follows a similar pattern to the one used for
\file{runup.py}, except that in this case we need to associate a
boundary type with each of the
boundary tag names introduced when we established the mesh. In place of the four
boundary types introduced for \file{runup.py}, we use the reflective
boundary for each of the
eight tagged segments defined by \code{create_mesh_from_regions}:

{\small \begin{verbatim}
    Br = Reflective_boundary(domain)
    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
                          'right2': Br, 'top': Br, 'left1': Br,
                          'left2': Br, 'left3': Br} )
\end{verbatim}}

\subsection{Evolution}

With the basics established, the running of the `evolve' step is
very similar to the corresponding step in \file{runup.py}. Here,
the simulation is run for five hours (18000 seconds) with
the output stored every two minutes (120 seconds).

{\small \begin{verbatim}
    import time t0 = time.time()

    for t in domain.evolve(yieldstep = 120, duration = 18000):
        print domain.timestepping_statistics()
        print domain.boundary_statistics(tags = 'bottom')

    print 'That took %.2f seconds' %(time.time()
\end{verbatim}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{\anuga Public Interface}
\label{ch:interface}

This chapter gives an overview of the features of \anuga available
to the user at the public interface. These are grouped under the
following headings, which correspond to the outline of the examples
described in Chapter \ref{ch:getstarted}:

\begin{itemize}
    \item Establishing the Mesh
    \item Initialising the Domain
    \item Specifying the Quantities
    \item Initial Conditions
    \item Boundary Conditions
    \item Forcing Functions
    \item Evolution
\end{itemize}

The listings are intended merely to give the reader an idea of what
each feature is, where to find it and how it can be used---they do
not give full specifications; for these the reader
may consult the code. The code for every function or class contains
a documentation string, or `docstring', that specifies the precise
syntax for its use. This appears immediately after the line
introducing the code, between two sets of triple quotes.

Each listing also describes the location of the module in which
the code for the feature being described can be found. All modules
are in the folder \file{inundation} or one of its subfolders, and the
location of each module is described relative to \file{inundation}. Rather
than using pathnames, whose syntax depends on the operating system,
we use the format adopted for importing the function or class for
use in Python code. For example, suppose we wish to specify that the
function \function{create\_mesh\_from\_regions} is in a module called
\module{mesh\_interface} in a subfolder of \module{inundation} called
\code{pmesh}. In Linux or Unix syntax, the pathname of the file
containing the function, relative to \file{inundation}, would be

\begin{center}
%    \code{pmesh/mesh\_interface.py}
    \code{pmesh}$\slash$\code{mesh\_interface.py}
\end{center}

while in Windows syntax it would be

\begin{center}
    \code{pmesh}$\backslash$\code{mesh\_interface.py}
\end{center}

Rather than using either of these forms, in this chapter we specify
the location simply as \code{pmesh.mesh\_interface}, in keeping with
the usage in the Python statement for importing the function,
namely:
\begin{center}
    \code{from pmesh.mesh\_interface import create\_mesh\_from\_regions}
\end{center}

Each listing details the full set of parameters for the class or
function; however, the description is generally limited to the most
important parameters and the reader is again referred to the code
for more details.

The following parameters are common to many functions and classes
and are omitted from the descriptions given below:

%\begin{center}
\begin{tabular}{ll}  %\hline
%\textbf{Name } & \textbf{Description}\\
%\hline
\emph{use\_cache} & Specifies whether caching is to be used for improved performance. See Section \ref{sec:caching} for details on the underlying caching functionality\\
\emph{verbose} & If \code{True}, provides detailed terminal output
to the user\\  % \hline
\end{tabular}
%\end{center}

\section{Mesh Generation}

Before discussing the part of the interface relating to mesh
generation, we begin with a description of a simple example of a
mesh and use it to describe how mesh data is stored.

\label{sec:meshexample} Figure \ref{fig:simplemesh} represents a
very simple mesh comprising just 11 points and 10 triangles.


\begin{figure}[h]
  \begin{center}
    \includegraphics[width=90mm, height=90mm]{triangularmesh.eps}
  \end{center}

  \caption{A simple mesh}
  \label{fig:simplemesh}
\end{figure}


The variables \code{points}, \code{vertices} and \code{boundary}
represent the data displayed in Figure \ref{fig:simplemesh} as
follows. The list \code{points} stores the coordinates of the
points, and may be displayed schematically as in Table
\ref{tab:points}.


\begin{table}
  \begin{center}
    \begin{tabular}[t]{|c|cc|} \hline
      index & \code{x} & \code{y}\\  \hline
      0 & 1 & 1\\
      1 & 4 & 2\\
      2 & 8 & 1\\
      3 & 1 & 3\\
      4 & 5 & 5\\
      5 & 8 & 6\\
      6 & 11 & 5\\
      7 & 3 & 6\\
      8 & 1 & 8\\
      9 & 4 & 9\\
      10 & 10 & 7\\  \hline
    \end{tabular}
  \end{center}

  \caption{Point coordinates for mesh in
    Figure \protect \ref{fig:simplemesh}}
  \label{tab:points}
\end{table}

The list \code{vertices} specifies the triangles that make up the
mesh. It does this by specifying, for each triangle, the indices
(the numbers shown in the first column above) that correspond to the
three points at its vertices, taken in an anti-clockwise order
around the triangle. Thus, in the example shown in Figure
\ref{fig:simplemesh}, the variable \code{vertices} contains the
entries shown in Table \ref{tab:vertices}. The starting point is
arbitrary so triangle $(0,1,3)$ is considered the same as $(1,3,0)$
and $(3,0,1)$.


\begin{table}
  \begin{center}
    \begin{tabular}{|c|ccc|} \hline
      index & \multicolumn{3}{c|}{\code{vertices}}\\ \hline
      0 & 0 & 1 & 3\\
      1 & 1 & 2 & 4\\
      2 & 2 & 5 & 4\\
      3 & 2 & 6 & 5\\
      4 & 4 & 5 & 9\\
      5 & 4 & 9 & 7\\
      6 & 3 & 4 & 7\\
      7 & 7 & 9 & 8\\
      8 & 1 & 4 & 3\\
      9 & 5 & 10 & 9\\  \hline
    \end{tabular}
  \end{center}

  \caption{Vertices for mesh in Figure \protect \ref{fig:simplemesh}}
  \label{tab:vertices}
\end{table}

Finally, the variable \code{boundary} identifies the boundary
triangles and associates a tag with each.

\refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}\label{sec:meshgeneration}

\begin{funcdesc}  {create\_mesh\_from\_regions}{bounding_polygon,
                             boundary_tags,
                             maximum_triangle_area,
                             filename=None,
                             interior_regions=None,
                             poly_geo_reference=None,
                             mesh_geo_reference=None,
                             minimum_triangle_angle=28.0}
Module: \module{pmesh.mesh\_interface}

This function allows a user to initiate the automatic creation of a
mesh inside a specified polygon (input \code{bounding_polygon}).
Among the parameters that can be set are the \emph{resolution}
(maximal area for any triangle in the mesh) and the minimal angle
allowable in any triangle. The user can specify a number of internal
polygons within each of which a separate mesh is to be created,
generally with a smaller resolution. \code{interior_regions} is a
paired list containing the interior polygon and its resolution.
Additionally, the user specifies a list of boundary tags, one for
each edge of the bounding polygon.

\textbf{WARNING}. Note that the dictionary structure used for the
parameter \code{boundary\_tags} is different from that used for the
variable \code{boundary} that occurs in the specification of a mesh.
In the case of \code{boundary}, the tags are the \emph{values} of
the dictionary, whereas in the case of \code{boundary_tags}, the
tags are the \emph{keys} and the \emph{value} corresponding to a
particular tag is a list of numbers identifying boundary edges
labelled with that tag. Because of this, it is theoretically
possible to assign the same edge to more than one tag. However, an
attempt to do this will cause an error.
\end{funcdesc}




\begin{classdesc}  {Mesh}{userSegments=None,
                 userVertices=None,
                 holes=None,
                 regions=None,
                 geo_reference=None}
Module: \module{pmesh.mesh}

A class used to build a mesh outline and generate a two-dimensional
triangular mesh. The mesh outline is used to describe features on the
mesh, such as the mesh boundary. Many of this classes methods are used
to build a mesh outline, such as \code{add\_vertices} and
\code{add\_region\_from\_polygon}.

\end{classdesc}


\subsection{Key Methods of Class Mesh}


\begin{methoddesc} {add\_hole}{x,y, geo_reference=None}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to build the mesh outline.  It defines a hole,
when the boundary of the hole has already been defined, by selecting a
point within the boundary. 

\end{methoddesc}


\begin{methoddesc}  {add\_hole\_from\_polygon}{self, polygon, tags=None,
    geo_reference=None}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to add a `hole' within a region ---that is, to
define a interior region where the triangular mesh will not be
generated---to a \class{Mesh} instance. The region boundary is described by
the polygon passed in.  Additionally, the user specifies a list of
boundary tags, one for each edge of the bounding polygon.
\end{methoddesc}


\begin{methoddesc} {add\_region}{x,y, geo_reference=None}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to build the mesh outline.  It defines a region,
when the boundary of the region has already been defined, by selecting
a point within the boundary.  A region instance is returned.  This can
be used to set the resolution.

\end{methoddesc}

\begin{methoddesc}  {add\_region\_from\_polygon}{self, polygon, tags=None,
                                max_triangle_area=None, geo_reference=None}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to build the mesh outline.  It adds a region to a
\class{Mesh} instance.  Regions are commonly used to describe an area
with an increased density of triangles, by setting
\code{max_triangle_area}.  The
region boundary is described by the input \code{polygon}.  Additionally, the
user specifies a list of boundary tags, one for each edge of the
bounding polygon.

\end{methoddesc}





\begin{methoddesc} {add\_vertices}{point_data}

Add user vertices. The point_data can be a list of (x,y) values, a numeric
array or a geospatial_data instance.
        
\end{methoddesc}

\begin{methoddesc}  {export\_mesh_file}{self,ofile}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to save the mesh to a file. \code{ofile} is the
name of the mesh file to be written, including the extension.  Use
the extension \code{.msh} for the file to be in NetCDF format and
\code{.tsh} for the file to be ASCII format.
\end{methoddesc}

\begin{methoddesc}  {generate\_mesh}{self,
                      maximum_triangle_area=None,
                      minimum_triangle_angle=28.0,
                      verbose=False}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to generate the triangular mesh.  The  maximal
area of any triangle in the mesh can be specified, which is used to
control the triangle density, along with the
minimum angle in any triangle.
\end{methoddesc}



\begin{methoddesc}  {import_ungenerate_file}{self,ofile, tag=None}
Module: \module{pmesh.mesh},  Class: \class{Mesh}

This method is used to import a polygon file in the ungenerate
format, which is used by arcGIS. The polygons from the file are converted to
vertices and segments. \code{ofile} is the name of the polygon file.
\code{tag} is the tag given to all the polygon's segments.

This function can be used to import building footprints.
\end{methoddesc}

%%%%%%
\section{Initialising the Domain}

%Include description of the class Domain and the module domain.

%FIXME (Ole): This is also defined in a later chapter
%\declaremodule{standard}{pyvolution.domain}

\begin{classdesc} {Domain} {source=None,
                 triangles=None,
                 boundary=None,
                 conserved_quantities=None,
                 other_quantities=None,
                 tagged_elements=None,
                 geo_reference=None,
                 use_inscribed_circle=False,
                 mesh_filename=None,
                 use_cache=False,
                 verbose=False,
                 full_send_dict=None,
                 ghost_recv_dict=None,
                 processor=0,
                 numproc=1}
Module: \refmodule{pyvolution.domain}

This class is used to create an instance of a data structure used to
store and manipulate data associated with a mesh. The mesh is
specified either by assigning the name of a mesh file to
\code{source} or by
\end{classdesc}

\begin{funcdesc}  {pmesh\_to\_domain\_instance}{file_name, DomainClass, use_cache = False, verbose = False}
Module: \module{pyvolution.pmesh2domain}

Once the initial mesh file has been created, this function is
applied to convert it to a domain object---that is, to a member of
the special Python class \class{Domain} (or a subclass of
\class{Domain}), which provides access to properties and methods
that allow quantities to be set and other operations to be carried
out.

\code{file\_name} is the name of the mesh file to be converted,
including the extension. \code{DomainClass} is the class to be
returned, which must be a subclass of \class{Domain} having the same
interface as \class{Domain}---in practice, it can usually be set
simply to \class{Domain}.

This is now superseded by Domain(mesh_filename).
\end{funcdesc}


\subsection{Key Methods of Domain}

\begin{methoddesc} {set\_name}{name}
    Module: \refmodule{pyvolution.domain}, page \pageref{mod:pyvolution.domain}  %\code{pyvolution.domain}

    Assigns the name \code{name} to the domain.
\end{methoddesc}

\begin{methoddesc} {get\_name}{}
    Module: \module{pyvolution.domain}

    Returns the name assigned to the domain by \code{set\_name}. If no name has been
    assigned, returns \code{`domain'}.
\end{methoddesc}

\begin{methoddesc} {set\_datadir}{name}
    Module: \module{pyvolution.domain}

    Specifies the directory used for SWW files, assigning it to the pathname \code{name}. The default value, before
    \code{set\_datadir} has been run, is the value \code{default\_datadir}
    specified in \code{config.py}.

    Since different operating systems use different formats for specifying pathnames,
    it is necessary to specify path separators using the Python code \code{os.sep}, rather than
    the operating-specific ones such as `$\slash$' or `$\backslash$'.
    For this to work you will need to include the statement \code{import os}
    in your code, before the first appearance of \code{set\_datadir}.

    For example, to set the data directory to a subdirectory
    \code{data} of the directory \code{project}, you could use
    the statements:

    {\small \begin{verbatim}
        import os
        domain.set_datadir{'project' + os.sep + 'data'}
    \end{verbatim}}
\end{methoddesc}

\begin{methoddesc} {get\_datadir}{}
    Module: \module{pyvolution.domain}

    Returns the data directory set by \code{set\_datadir} or,
    if \code{set\_datadir} has not
    been run, returns the value \code{default\_datadir} specified in
    \code{config.py}.
\end{methoddesc}

\begin{methoddesc} {set\_time}{time=0.0}
    Module: \module{pyvolution.domain}

    Sets the initial time, in seconds, for the simulation. The
    default is 0.0.
\end{methoddesc}

\begin{methoddesc} {set\_default\_order}{n}
    Sets the default (spatial) order to the value specified by
    \code{n}, which must be either 1 or 2. (Assigning any other value
    to \code{n} will cause an error.)
\end{methoddesc}


%%%%%%
\section{Initial Conditions}
\label{sec:Initial Conditions}
In standard usage of partial differential equations, initial conditions
refers to the values associated to the system variables (the conserved
quantities here) for \code{time = 0}. In setting up a scenario script
as described in Sections \ref{sec:simpleexample} and \ref{sec:realdataexample},
\code{set_quantity} is used to define the initial conditions of variables
other than the conserved quantities, such as friction. Here, we use the terminology
of initial conditions to refer to initial values for variables which need
prescription to solve the shallow water wave equation. Further, it must be noted
that \code{set_quantity} does not necessarily have to be used in the initial
condition setting; it can be used at any time throughout the simulation.

\begin{methoddesc}{set\_quantity}{name,
    numeric = None,
    quantity = None,
    function = None,
    geospatial_data = None,
    filename = None,
    attribute_name = None,
    alpha = None,
    location = 'vertices',
    indices = None,
    verbose = False,
    use_cache = False}
  Module: \module{pyvolution.domain}
  (see also \module{pyvolution.quantity.set\_values})

This function is used to assign values to individual quantities for a
domain. It is very flexible and can be used with many data types: a
statement of the form \code{domain.set\_quantity(name, x)} can be used
to define a quantity having the name \code{name}, where the other
argument \code{x} can be any of the following:

\begin{itemize}
\item a number, in which case all vertices in the mesh gets that for
the quantity in question.
\item a list of numbers or a Numeric array ordered the same way as the mesh vertices.
\item a function (e.g.\ see the samples introduced in Chapter 2)
\item an expression composed of other quantities and numbers, arrays, lists (for
example, a linear combination of quantities, such as
\code{domain.set\_quantity('stage','elevation'+x))}
\item the name of a file from which the data can be read. In this case, the optional argument attribute\_name will select which attribute to use from the file. If left out, set\_quantity will pick one. This is useful in cases where there is only one attribute.
\item a geospatial dataset (See ?????). Optional argument attribute\_name applies here as with files.
\end{itemize}


Exactly one of the arguments
  numeric, quantity, function, points, filename
must be present.


Set quantity will look at the type of the second argument (\code{numeric}) and
determine what action to take.

Values can also be set using the appropriate keyword arguments.
If x is a function, for example, \code{domain.set\_quantity(name, x)}, \code{domain.set\_quantity(name, numeric=x)}, and \code{domain.set\_quantity(name, function=x)}
are all equivalent.


Other optional arguments are
\begin{itemize}
\item \code{indices} which is a list of ids of triangles to which set\_quantity should apply its assignment of values.
\item \code{location} determines which part of the triangles to assign to. Options are 'vertices' (default), 'edges', and 'centroids'.
\end{itemize}


\end{methoddesc}







%%%
\anuga provides a number of predefined initial conditions to be used
with \code{set\_quantity}.

\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
                x0=0.0, y0=0.0, alpha=0.0,
                gravity=9.8, gamma=1.85,
                massco=1, dragco=1, frictionco=0, psi=0,
                dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
                domain=None,
                verbose=False}
Module: \module{pyvolution.smf}

This function returns a callable object representing an initial water
displacement generated by a submarine sediment failure. These failures can take the form of
a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead.

The arguments include as a minimum, the slump or slide length, the water depth to the centre of sediment
mass, and the bathymetric slope. Other slump or slide parameters can be included if they are known.
\end{funcdesc}


%%%
\begin{funcdesc}{file\_function}{filename,
    domain = None,
    quantities = None,
    interpolation_points = None,
    verbose = False,
    use_cache = False}
Module: \module{pyvolution.util}

Reads the time history of spatial data for
specified interpolation points from a NetCDF file (\code{filename})
and returns
a callable object. \code{filename} could be a \code{sww} file.
Returns interpolated values based on the input
file using the underlying \code{interpolation\_function}.

\code{quantities} is either the name of a single quantity to be
interpolated or a list of such quantity names. In the second case, the resulting
function will return a tuple of values---one for each quantity.

\code{interpolation\_points} is a list of absolute UTM coordinates
for points at which values are sought.

The model time stored within the file function can be accessed using
the method \code{f.get\_time()}
\end{funcdesc}

%%%
\begin{classdesc}{Interpolation\_function}{self,
    time,
    quantities,
    quantity_names = None,
    vertex_coordinates = None,
    triangles = None,
    interpolation_points = None,
    verbose = False}
Module: \module{pyvolution.least\_squares}

Given a time series (i.e. a series of values associated with
different times), whose values are either just numbers or a set of
numbers defined at the vertices of a triangular mesh (such as those
stored in SWW files), \code{Interpolation\_function} is used to
create a callable object that interpolates a value for an arbitrary
time \code{t} within the model limits and possibly a point \code{(x,
y)} within a mesh region.

The actual time series at which data is available is specified by
means of an array \code{time} of monotonically increasing times. The
quantities containing the values to be interpolated are specified in
an array---or dictionary of arrays (used in conjunction with the
optional argument \code{quantity\_names}) --- called
\code{quantities}. The optional arguments \code{vertex\_coordinates}
and \code{triangles} represent the spatial mesh associated with the
quantity arrays. If omitted the function created by
\code{Interpolation\_function} will be a function of \code{t} only.

Since, in practice, values need to be computed at specified points,
the syntax allows the user to specify, once and for all, a list
\code{interpolation\_points} of points at which values are required.
In this case, the function may be called using the form \code{f(t,
id)}, where \code{id} is an index for the list
\code{interpolation\_points}.

\end{classdesc}

%%%
%\begin{funcdesc}{set\_region}{functions}
%[Low priority. Will be merged into set\_quantity]

%Module:\module{pyvolution.domain}
%\end{funcdesc}



%%%%%%
\section{Boundary Conditions}\index{boundary conditions}

\anuga provides a large number of predefined boundary conditions,
represented by objects such as \code{Reflective\_boundary(domain)} and
\code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples
in Chapter 2. Alternatively, you may prefer to ``roll your own'',
following the method explained in Section \ref{sec:roll your own}.

These boundary objects may be used with the function \code{set\_boundary} described below
to assign boundary conditions according to the tags used to label boundary segments.

\begin{methoddesc}{set\_boundary}{boundary_map}
Module: \module{pyvolution.domain}

This function allows you to assign a boundary object (corresponding to a
pre-defined or user-specified boundary condition) to every boundary segment that
has been assigned a particular tag.

This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects
and whose keys are the symbolic tags.

\end{methoddesc}

\begin{methoddesc} {get\_boundary\_tags}{}
Module: \module{pyvolution.mesh}

Returns a list of the available boundary tags.
\end{methoddesc}

%%%
\subsection{Predefined boundary conditions}

\begin{classdesc}{Reflective\_boundary}{Boundary}
Module: \module{pyvolution.shallow\_water}

Reflective boundary returns same conserved quantities as those present in
the neighbouring volume but reflected.

This class is specific to the shallow water equation as it works with the
momentum quantities assumed to be the second and third conserved quantities.
\end{classdesc}

%%%
\begin{classdesc}{Transmissive\_boundary}{domain = None}
Module: \module{pyvolution.generic\_boundary\_conditions}

A transmissive boundary returns the same conserved quantities as
those present in the neighbouring volume.

The underlying domain must be specified when the boundary is instantiated.
\end{classdesc}

%%%
\begin{classdesc}{Dirichlet\_boundary}{conserved_quantities=None}
Module: \module{pyvolution.generic\_boundary\_conditions}

A Dirichlet boundary returns constant values for each of conserved
quantities. In the example of \code{Dirichlet\_boundary([0.2, 0.0, 0.0])},
the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and
\code{ymomentum} at the boundary are set to 0.0. The list must contain
a value for each conserved quantity.
\end{classdesc}

%%%
\begin{classdesc}{Time\_boundary}{domain = None, f = None}
Module: \module{pyvolution.generic\_boundary\_conditions}

A time-dependent boundary returns values for the conserved
quantities as a function \code{f(t)} of time. The user must specify
the domain to get access to the model time.
\end{classdesc}

%%%
\begin{classdesc}{File\_boundary}{Boundary}
Module: \module{pyvolution.generic\_boundary\_conditions}

This method may be used if the user wishes to apply a SWW file or
a time series file to a boundary segment or segments.
The boundary values are obtained from a file and interpolated to the
appropriate segments for each conserved quantity.
\end{classdesc}



%%%
\begin{classdesc}{Transmissive\_Momentum\_Set\_Stage\_boundary}{Boundary}
Module: \module{pyvolution.shallow\_water}

This boundary returns same momentum conserved quantities as
those present in its neighbour volume but sets stage as in a Time\_boundary.
The underlying domain must be specified when boundary is instantiated

This type of boundary is useful when stage is known at the boundary as a
function of time, but momenta (or speeds) aren't.

This class is specific to the shallow water equation as it works with the
momentum quantities assumed to be the second and third conserved quantities.
\end{classdesc}



\subsection{User-defined boundary conditions}
\label{sec:roll your own}
[How to roll your own] TBA





\section{Forcing Functions}

\anuga provides a number of predefined forcing functions to be used with .....

%\begin{itemize}


%  \item \indexedcode{}
%  [function, arguments]

%  \item \indexedcode{}

%\end{itemize}



\section{Evolution}\index{evolution}

  \begin{methoddesc}{evolve}{yieldstep = None, finaltime = None, duration = None, skip_initial_step = False}

  Module: \module{pyvolution.domain}

  This function (a method of \class{domain}) is invoked once all the
  preliminaries have been completed, and causes the model to progress
  through successive steps in its evolution, storing results and
  outputting statistics whenever a user-specified period
  \code{yieldstep} is completed (generally during this period the
  model will evolve through several steps internally
  as the method forces the water speed to be calculated
  on successive new cells). The user
  specifies the total time period over which the evolution is to take
  place, by specifying values (in seconds) for either \code{duration}
  or \code{finaltime}, as well as the interval in seconds after which
  results are to be stored and statistics output.

  You can include \method{evolve} in a statement of the type:

  {\small \begin{verbatim}
      for t in domain.evolve(yieldstep, finaltime):
          <Do something with domain and t>
  \end{verbatim}}

  \end{methoddesc}



\subsection{Diagnostics}


  \begin{funcdesc}{statistics}{}
  Module: \module{pyvolution.domain}

  \end{funcdesc}

  \begin{funcdesc}{timestepping\_statistics}{}
  Module: \module{pyvolution.domain}

  Returns a string of the following type for each
  timestep:

  \code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12
  (12)}

  Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
  the number of first-order steps, respectively.
  \end{funcdesc}


  \begin{funcdesc}{boundary\_statistics}{quantities = None, tags = None}
  Module: \module{pyvolution.domain}

  Returns a string of the following type when \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}:

  {\small \begin{verbatim}
 Boundary values at time 0.5000:
    top:
        stage in [ -0.25821218,  -0.02499998]
    bottom:
        stage in [ -0.27098821,  -0.02499974]
  \end{verbatim}}

  \end{funcdesc}


  \begin{funcdesc}{get\_quantity}{name, location='vertices', indices = None}
  Module: \module{pyvolution.domain}
  Allow access to individual quantities and their methods

  \end{funcdesc}


  \begin{funcdesc}{get\_values}{location='vertices', indices = None}
  Module: \module{pyvolution.quantity}

  Extract values for quantity as an array

  \end{funcdesc}


  \begin{funcdesc}{get\_integral}{}
  Module: \module{pyvolution.quantity}

  Return computed integral over entire domain for this quantity

  \end{funcdesc}


\section{Other}

  \begin{funcdesc}{domain.create\_quantity\_from\_expression}{???}

  Handy for creating derived quantities on-the-fly.
  See \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
  \end{funcdesc}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{\anuga System Architecture}

From pyvolution/documentation

\section{File Formats}
\label{sec:file formats}

\anuga makes use of a number of different file formats. The
following table lists all these formats, which are described in more
detail in the paragraphs below.

\bigskip

\begin{center}

\begin{tabular}{|ll|}  \hline

\textbf{Extension} & \textbf{Description} \\
\hline\hline

\code{.sww} & NetCDF format for storing model output
\code{f(t,x,y)}\\

\code{.tms} & NetCDF format for storing time series \code{f(t)}\\

\code{.xya} & ASCII format for storing arbitrary points and
associated attributes\\

\code{.pts} & NetCDF format for storing arbitrary points and
associated attributes\\

\code{.asc} & ASCII format of regular DEMs as output from ArcView\\

\code{.prj} & Associated ArcView file giving more metadata for
\code{.asc} format\\

\code{.ers} & ERMapper header format of regular DEMs for ArcView\\

\code{.dem} & NetCDF representation of regular DEM data\\

\code{.tsh} & ASCII format for storing meshes and associated
boundary and region info\\

\code{.msh} & NetCDF format for storing meshes and associated
boundary and region info\\

\code{.nc} & Native ferret NetCDF format\\

\code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
%\caption{File formats used by \anuga}
\end{tabular}


\end{center}

The above table shows the file extensions used to identify the
formats of files. However, typically, in referring to a format we
capitalise the extension and omit the initial full stop---thus, we
refer, for example, to `SWW files' or `PRJ files'.

\bigskip

A typical dataflow can be described as follows:

\subsection{Manually Created Files}

\begin{tabular}{ll}
ASC, PRJ & Digital elevation models (gridded)\\
NC & Model outputs for use as boundary conditions (e.g. from MOST)
\end{tabular}

\subsection{Automatically Created Files}

\begin{tabular}{ll}
ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert
DEMs to native \code{.pts} file\\

NC $\rightarrow$ SWW & Convert MOST boundary files to
boundary \code{.sww}\\

PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\

TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using
Swollen\\

TSH + Boundary SWW $\rightarrow$ SWW & Simulation using
\code{pyvolution}\\

Polygonal mesh outline $\rightarrow$ & TSH or MSH
\end{tabular}




\bigskip

\subsection{SWW and TMS Formats}

The SWW and TMS formats are both NetCDF formats, and are of key
importance for \anuga.

An SWW file is used for storing \anuga output and therefore pertains
to a set of points and a set of times at which a model is evaluated.
It contains, in addition to dimension information, the following
variables:

\begin{itemize}
    \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays
    \item \code{elevation}, a Numeric array storing bed-elevations
    \item \code{volumes}, a list specifying the points at the vertices of each of the
    triangles
    % Refer here to the example to be provided in describing the simple example
    \item \code{time}, a Numeric array containing times for model
    evaluation
\end{itemize}


The contents of an SWW file may be viewed using the visualisation
tool \code{swollen}, which creates an on-screen geometric
representation. See section \ref{sec:swollen} (page
\pageref{sec:swollen}) in Appendix \ref{ch:supportingtools} for more
on \code{swollen}.

Alternatively, there are tools, such as \code{ncdump}, that allow
you to convert an NetCDF file into a readable format such as the
Class Definition Language (CDL). The following is an excerpt from a
CDL representation of the output file \file{bedslope.sww} generated
from running the simple example \file{runup.py} of
Chapter \ref{ch:getstarted}:

\verbatiminput{examples/bedslopeexcerpt.cdl}

The SWW format is used not only for output but also serves as input
for functions such as \function{file\_boundary} and
\function{file\_function}, described in Chapter \ref{ch:interface}.

A TMS file is used to store time series data that is independent of
position.


\subsection{Mesh File Formats}

A mesh file is a file that has a specific format suited to
triangular meshes and their outlines. A mesh file can have one of
two formats: it can be either a TSH file, which is an ASCII file, or
an MSH file, which is a NetCDF file. A mesh file can be generated
from the function \function{create\_mesh\_from\_regions} (see
Section \ref{sec:meshgeneration}) and used to initialise a domain.

A mesh file can define the outline of the mesh---the vertices and
line segments that enclose the region in which the mesh is
created---and the triangular mesh itself, which is specified by
listing the triangles and their vertices, and the segments, which
are those sides of the triangles that are associated with boundary
conditions.

In addition, a mesh file may contain `holes' and/or `regions'. A
hole represents an area where no mesh is to be created, while a
region is a labelled area used for defining properties of a mesh,
such as friction values.  A hole or region is specified by a point
and bounded by a number of segments that enclose that point.

A mesh file can also contain a georeference, which describes an
offset to be applied to $x$ and $y$ values---eg to the vertices.


\subsection{Formats for Storing Arbitrary Points and Attributes}

An XYA file is used to store data representing arbitrary numerical
attributes associated with a set of points.

The format for an XYA file is:\\
%\begin{verbatim}

            first line:     \code{[attribute names]}\\
            other lines:  \code{x y [attributes]}\\

            for example:\\
            \code{elevation, friction}\\
            \code{0.6, 0.7, 4.9, 0.3}\\
            \code{1.9, 2.8, 5, 0.3}\\
            \code{2.7, 2.4, 5.2, 0.3}

        The first two columns are always implicitly assumed to be $x$, $y$ coordinates.
        Use the same delimiter for the attribute names and the data.

        An XYA file can optionally end with a description of the georeference:

            \code{\#geo reference}\\
            \code{56}\\
            \code{466600.0}\\
            \code{8644444.0}

        Here the first number specifies the UTM zone (in this case zone 56)  and other numbers specify the
        easting and northing
        coordinates (in this case (466600.0, 8644444.0)) of the lower left corner.

A PTS file is a NetCDF representation of the data held in an XYA
file. If the data is associated with a set of $N$ points, then the
data is stored using an $N \times 2$ Numeric array of float
variables for the points and an $N \times 1$ Numeric array for each
attribute.

%\end{verbatim}

\subsection{ArcView Formats}

Files of the three formats ASC, PRJ and ERS are all associated with
data from ArcView.

An ASC file is an ASCII representation of DEM output from ArcView.
It contains a header with the following format:

\begin{tabular}{l l}
\code{ncols}      &   \code{753}\\
\code{nrows}      &   \code{766}\\
\code{xllcorner}  &   \code{314036.58727982}\\
\code{yllcorner}  & \code{6224951.2960092}\\
\code{cellsize}   & \code{100}\\
\code{NODATA_value} & \code{-9999}
\end{tabular}

The remainder of the file contains the elevation data for each grid point
in the grid defined by the above information.

A PRJ file is an ArcView file used in conjunction with an ASC file
to represent metadata for a DEM.


\subsection{DEM Format}

A DEM file is a NetCDF representation of regular DEM data.


\subsection{Other Formats}




\subsection{Basic File Conversions}

  \begin{funcdesc}{sww2dem}{basename_in, basename_out = None,
            quantity = None,
            timestep = None,
            reduction = None,
            cellsize = 10,
            NODATA_value = -9999,
            easting_min = None,
            easting_max = None,
            northing_min = None,
            northing_max = None,
            expand_search = False,
            verbose = False,
            origin = None,
            datum = 'WGS84',
            format = 'ers'}
  Module: \module{pyvolution.data\_manager}

  Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
  ERS) of a desired grid size \code{cellsize} in metres.
  The easting and northing values are used if the user wished to clip the output
  file to a specified rectangular area. The \code{reduction} input refers to a function
  to reduce the quantities over all time step of the SWW file, example, maximum.
  \end{funcdesc}


  \begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
            easting_min=None, easting_max=None,
            northing_min=None, northing_max=None,
            use_cache=False, verbose=False}
  Module: \module{pyvolution.data\_manager}

  Takes DEM data (a NetCDF file representation of data from a regular Digital
  Elevation Model) and converts it to PTS format.
  \end{funcdesc}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Basic \anuga Assumptions}

(From pyvolution/documentation)


Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
If one wished to recreate scenarios prior to that date it must be done
using some relative time (e.g. 0).


All spatial data relates to the WGS84 datum (or GDA94) and has been
projected into UTM with false easting of 500000 and false northing of
1000000 on the southern hemisphere (0 on the northern).

It is assumed that all computations take place within one UTM zone.

DEMs, meshes and boundary conditions can have different origins within
one UTM zone. However, the computation will use that of the mesh for
numerical stability.


%OLD
%The dataflow is: (See data_manager.py and from scenarios)
%
%
%Simulation scenarios
%--------------------%
%%
%
%Sub directories contain scrips and derived files for each simulation.
%The directory ../source_data contains large source files such as
%DEMs provided externally as well as MOST tsunami simulations to be used
%as boundary conditions.
%
%Manual steps are:
%  Creation of DEMs from argcview (.asc + .prj)
%  Creation of mesh from pmesh (.tsh)
%  Creation of tsunami simulations from MOST (.nc)
%%
%
%Typical scripted steps are%
%
%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
%                   native dem and pts formats%
%
%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
%                  as boundary condition%
%
%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
%                   smoothing. The outputs are tsh files with elevation data.%
%
%  run_simulation.py: Use the above together with various parameters to
%                     run inundation simulation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix

\chapter{Supporting Tools}
\label{ch:supportingtools}

This section describes a number of supporting tools, supplied with \anuga, that offer a
variety of types of functionality and enhance the basic capabilities of \anuga.

\section{caching}
\label{sec:caching}

The \code{cache} function is used to provide supervised caching of function 
results. A Python function call of the form

      {\small \begin{verbatim}
      result = func(arg1,...,argn)
      \end{verbatim}}

  can be replaced by

      {\small \begin{verbatim}
      from caching import cache
      result = cache(func,(arg1,...,argn))
      \end{verbatim}}

  which returns the same output but reuses cached
  results if the function has been computed previously in the same context.
  \code{result} and the arguments can be simple types, tuples, list, dictionaries or
  objects, but not unhashable types such as functions or open file objects.
  The function \code{func} may be a member function of an object or a module.

  This type of caching is particularly useful for computationally intensive
  functions with few frequently used combinations of input arguments. Note that
  if the inputs or output are very large caching may not save time because
  disc access may dominate the execution time.

  If the function definition changes after a result has been cached, this will be
  detected by examining the functions \code{bytecode (co\_code, co\_consts,
  func\_defaults, co\_argcount)} and the function will be recomputed.
  However, caching will not detect changes in modules used by \code{func}.
  In this case cache must be cleared manually.

  Options are set by means of the function \code{set\_option(key, value)},
  where \code{key} is a key associated with a
  Python dictionary \code{options}. This dictionary stores settings such as the name of
  the directory used, the maximum
  number of cached files allowed, and so on.

  The \code{cache} function allows the user also to specify a list of dependent files. If any of these
  have been changed, the function is recomputed and the results stored again.

  %Other features include support for compression and a capability to \ldots


   \textbf{USAGE:} \nopagebreak

    {\small \begin{verbatim}
    result = cache(func, args, kwargs, dependencies, cachedir, verbose,
                   compression, evaluate, test, return_filename)
    \end{verbatim}}


\section{swollen}
\label{sec:swollen}
 The output generated by \anuga may be viewed by
means of the visualisation tool \code{swollen}, which takes the
\code{SWW} file output by \anuga and creates a visual representation
of the data. Examples may be seen in Figures \ref{fig:runupstart}
and \ref{fig:bedslope2}. To view an \code{SWW} file with
\code{swollen} in the Windows environment, you can simply drag the
icon representing the file over an icon on the desktop for the
\code{swollen} executable file (or a shortcut to it), or set up a
file association to make files with the extension \code{.sww} open
with \code{swollen}. Alternatively, you can operate \code{swollen}
from the command line, in both Windows and Linux environments.

On successful operation, you will see an interactive moving-picture
display. You can use keys and the mouse to slow down, speed up or
stop the display, change the viewing position or carry out a number
of other simple operations. Help is also displayed when you press
the \code{h} key.

The main keys operating the interactive screen are:\\

\begin{center}
\begin{tabular}{|ll|}   \hline

\code{w} & toggle wireframe \\

space bar & start/stop\\

up/down arrows & increase/decrease speed\\

left/right arrows & direction in time \emph{(when running)}\\
& step through simulation \emph{(when stopped)}\\

left mouse button & rotate\\

middle mouse button & pan\\

right mouse button & zoom\\  \hline

\end{tabular}
\end{center}

\vfill

The following table describes how to operate swollen from the command line:

Usage: \code{swollen [options] swwfile \ldots}\\  \nopagebreak
Options:\\  \nopagebreak
\begin{tabular}{ll}
  \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\
                                    & \code{HEAD\_MOUNTED\_DISPLAY}\\
  \code{--rgba} & Request a RGBA colour buffer visual\\
  \code{--stencil} & Request a stencil buffer visual\\
  \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
                                    & overridden by environmental variable\\
  \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD\_BUFFER | HORIZONTAL\_SPLIT |}\\
                                    & \code{VERTICAL\_SPLIT | LEFT\_EYE | RIGHT\_EYE |}\\
                                     & \code{ON | OFF} \\
  \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
  \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
\end{tabular}

\begin{tabular}{ll}
  \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
  \code{-help} & Display this information\\
  \code{-hmax <float>} & Height above which transparency is set to
                                     \code{alphamax}\\
\end{tabular}

\begin{tabular}{ll}

  \code{-hmin <float>} & Height below which transparency is set to
                                     zero\\
\end{tabular}

\begin{tabular}{ll}
  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
                                     up, default is overhead)\\
\end{tabular}

\begin{tabular}{ll}
  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\

\end{tabular}

\begin{tabular}{ll}
  \code{-movie <dirname>} & Save numbered images to named directory and
                                     quit\\

  \code{-nosky} & Omit background sky\\


  \code{-scale <float>} & Vertical scale factor\\
  \code{-texture <file>} & Image to use for bedslope topography\\
  \code{-tps <rate>} & Timesteps per second\\
  \code{-version} & Revision number and creation (not compile)
                                     date\\
\end{tabular}

\section{utilities/polygons}

  \declaremodule{standard}{utilities.polygon}
  \refmodindex{utilities.polygon}

  \begin{classdesc}{Polygon\_function}{regions, default = 0.0, geo_reference = None}
  Module: \code{utilities.polygon}

  Creates a callable object that returns one of a specified list of values when
  evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the
  point belongs to. The parameter \code{regions} is a list of pairs
  \code{(P, v)}, where each \code{P} is a polygon and each \code{v}
  is either a constant value or a function of coordinates \code{x}
  and \code{y}, specifying the return value for a point inside \code{P}. The
  optional parameter \code{default} may be used to specify a value
  for a point not lying inside any of the specified polygons. When a
  point lies in more than one polygon, the return value is taken to
  be the value for whichever of these polygon appears later in the
  list.
  %FIXME (Howard): CAN x, y BE VECTORS?

  \end{classdesc}

  \begin{funcdesc}{read\_polygon}{filename}
  Module: \code{utilities.polygon}

  Reads the specified file and returns a polygon. Each
  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
  as coordinates of one vertex of the polygon.
  \end{funcdesc}

  \begin{funcdesc}{populate\_polygon}{polygon, number_of_points, seed = None, exclude = None}
  Module: \code{utilities.polygon}

  Populates the interior of the specified polygon with the specified number of points,
  selected by means of a uniform distribution function.
  \end{funcdesc}

  \begin{funcdesc}{point\_in\_polygon}{polygon, delta=1e-8}
  Module: \code{utilities.polygon}

  Returns a point inside the specified polygon and close to the edge. The distance between
  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
  second argument \code{delta}, which is taken as $10^{-8}$ by default.
  \end{funcdesc}

  \begin{funcdesc}{inside\_polygon}{points, polygon, closed = True, verbose = False}
  Module: \code{utilities.polygon}

  Used to test whether the members of a list of points
  are inside the specified polygon. Returns a Numeric
  array comprising the indices of the points in the list that lie inside the polygon.
  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
  Points on the edges of the polygon are regarded as inside if
  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
  \end{funcdesc}

  \begin{funcdesc}{outside\_polygon}{points, polygon, closed = True, verbose = False}
  Module: \code{utilities.polygon}

  Exactly like \code{inside\_polygon}, but with the words `inside' and `outside' interchanged.
  \end{funcdesc}

  \begin{funcdesc}{is\_inside\_polygon}{point, polygon, closed=True, verbose=False}
  Module: \code{utilities.polygon}

  Returns \code{True} if \code{point} is inside \code{polygon} or
  \code{False} otherwise. Points on the edges of the polygon are regarded as inside if
  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
  \end{funcdesc}

  \begin{funcdesc}{is\_outside\_polygon}{point, polygon, closed=True, verbose=False}
  Module: \code{utilities.polygon}

  Exactly like \code{is\_outside\_polygon}, but with the words `inside' and `outside' interchanged.
  \end{funcdesc}

  \begin{funcdesc}{point\_on\_line}{x, y, x0, y0, x1, y1}
  Module: \code{utilities.polygon}

  Returns \code{True} or \code{False}, depending on whether the point with coordinates
  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
  and \code{x1, y1} (extended if necessary at either end).
  \end{funcdesc}

  \begin{funcdesc}{separate\_points\_by\_polygon}{points, polygon, closed = True, verbose = False}
    \indexedcode{separate\_points\_by\_polygon}
  Module: \code{utilities.polygon}

  \end{funcdesc}

  \begin{funcdesc}{polygon\_area}{polygon}
  Module: \code{utilities.polygon}

  Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html)
  \end{funcdesc}

  \begin{funcdesc}{plot\_polygons}{polygons, figname, verbose = False}
  Module: \code{utilities.polygon}

  Plots each polygon contained in input polygon list, e.g.
 \code{polygons = [poly1, poly2, poly3]} where \code{poly1 = [[x11,y11],[x12,y12],[x13,y13]]}
 etc.  Each polygon is closed for plotting purposes and subsequent plot saved to \code{figname}.
  Returns list containing the minimum and maximum of \code{x} and \code{y},
  i.e. \code{[x_{min}, x_{max}, y_{min}, y_{max}]}.
  \end{funcdesc}

\section{coordinate\_transforms}

\section{geospatial\_data}

This describes a class that represents arbitrary point data in UTM
coordinates along with named attribute values.

%FIXME (Ole): This gives a LaTeX error
%\declaremodule{standard}{geospatial_data}
%\refmodindex{geospatial_data}



\begin{classdesc}{Geospatial\_data}
  {data_points = None,
    attributes = None,
    geo_reference = None,
    default_attribute_name = None,
    file_name = None}
Module: \code{geospatial\_data}

This class is used to store a set of data points and associated
attributes, allowing these to be manipulated by methods defined for
the class.

The data points are specified either by reading them from a NetCDF
or XYA file, identified through the parameter \code{file\_name}, or
by providing their \code{x}- and \code{y}-coordinates in metres,
either as a sequence of 2-tuples of floats or as an $M \times 2$
Numeric array of floats, where $M$ is the number of points.
Coordinates are interpreted relative to the origin specified by the
object \code{geo\_reference}, which contains data indicating the UTM
zone, easting and northing. If \code{geo\_reference} is not
specified, a default is used.

Attributes are specified through the parameter \code{attributes},
set either to a list or array of length $M$ or to a dictionary whose
keys are the attribute names and whose values are lists or arrays of
length $M$. One of the attributes may be specified as the default
attribute, by assigning its name to \code{default\_attribute\_name}.
If no value is specified, the default attribute is taken to be the
first one.
\end{classdesc}


\begin{methoddesc}{import\_points\_file}{delimiter = None, verbose = False}

\end{methoddesc}


\begin{methoddesc}{export\_points\_file}{ofile, absolute=False}

\end{methoddesc}


\begin{methoddesc}{get\_data\_points}{absolute = True}

\end{methoddesc}


\begin{methoddesc}{set\_attributes}{attributes}

\end{methoddesc}


\begin{methoddesc}{get\_attributes}{attribute_name = None}

\end{methoddesc}


\begin{methoddesc}{get\_all\_attributes}{}

\end{methoddesc}


\begin{methoddesc}{set\_default\_attribute\_name}{default_attribute_name}

\end{methoddesc}


\begin{methoddesc}{set\_geo\_reference}{geo_reference}

\end{methoddesc}


\begin{methoddesc}{add}{}

\end{methoddesc}


\section{pmesh GUI}

\section{alpha\_shape}
\emph{Alpha shapes} are used to generate close-fitting boundaries
around sets of points. The alpha shape algorithm produces a shape
that approximates to the `shape formed by the points'---or the shape
that would be seen by viewing the points from a coarse enough
resolution. For the simplest types of point sets, the alpha shape
reduces to the more precise notion of the convex hull. However, for
many sets of points the convex hull does not provide a close fit and
the alpha shape usually fits more closely to the original point set,
offering a better approximation to the shape being sought.

In \anuga, an alpha shape is used to generate a polygonal boundary
around a set of points before mesh generation. The algorithm uses a
parameter $\alpha$ that can be adjusted to make the resultant shape
resemble the shape suggested by intuition more closely. An alpha
shape can serve as an initial boundary approximation that the user
can adjust as needed.

The following paragraphs describe the class used to model an alpha
shape and some of the important methods and attributes associated
with instances of this class.

\begin{classdesc}{Alpha\_Shape}{points, alpha = None}
Module: \code{alpha\_shape}

To instantiate this class the user supplies the points from which
the alpha shape is to be created (in the form of a list of 2-tuples
\code{[[x1, y1],[x2, y2]}\ldots\code{]}, assigned to the parameter
\code{points}) and, optionally, a value for the parameter
\code{alpha}. The alpha shape is then computed and the user can then
retrieve details of the boundary through the attributes defined for
the class.
\end{classdesc}


\begin{funcdesc}{alpha\_shape\_via\_files}{point_file, boundary_file, alpha= None}
Module: \code{alpha\_shape}

This function reads points from the specified point file
\code{point\_file}, computes the associated alpha shape (either
using the specified value for \code{alpha} or, if no value is
specified, automatically setting it to an optimal value) and outputs
the boundary to a file named \code{boundary\_file}. This output file
lists the coordinates \code{x, y} of each point in the boundary,
using one line per point.
\end{funcdesc}


\begin{methoddesc}{set\_boundary\_type}{self,raw_boundary=True,
                          remove_holes=False,
                          smooth_indents=False,
                          expand_pinch=False,
                          boundary_points_fraction=0.2}
Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}

This function sets flags that govern the operation of the algorithm
that computes the boundary, as follows:

\code{raw\_boundary = True} returns raw boundary, i.e. the regular edges of the
                alpha shape.\\
\code{remove\_holes = True} removes small holes (`small' is defined by
\code{boundary\_points\_fraction})\\
\code{smooth\_indents = True} removes sharp triangular indents in
boundary\\
\code{expand\_pinch = True} tests for pinch-off and
corrects---preventing a boundary vertex from having more than two edges.
\end{methoddesc}


\begin{methoddesc}{get\_boundary}{}
Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}

Returns a list of tuples representing the boundary of the alpha
shape. Each tuple represents a segment in the boundary by providing
the indices of its two endpoints.
\end{methoddesc}


\begin{methoddesc}{write\_boundary}{file_name}
Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}

Writes the list of 2-tuples returned by \code{get\_boundary} to the
file \code{file\_name}, using one line per tuple.
\end{methoddesc}

\section{Numerical Tools}

The following table describes some useful numerical functions that
may be found in the module \module{utilities.numerical\_tools}:

\begin{tabular}{|p{8cm} p{8cm}|}  \hline
\code{angle(v1, v2=None)} & Angle between two-dimensional vectors
\code{v1} and \code{v2}, or between \code{v1} and the $x$-axis if
\code{v2} is \code{None}. Value is in range $0$ to $2\pi$. \\

\code{normal\_vector(v)} & Normal vector to \code{v}.\\

\code{mean(x)} & Mean value of a vector \code{x}.\\

\code{cov(x, y=None)} & Covariance of vectors \code{x} and \code{y}.
If \code{y} is \code{None}, returns \code{cov(x, x)}.\\

\code{err(x, y=0, n=2, relative=True)} & Relative error of
$\parallel$\code{x}$-$\code{y}$\parallel$ to
$\parallel$\code{y}$\parallel$ (2-norm if \code{n} = 2 or Max norm
if \code{n} = \code{None}). If denominator evaluates to zero or if
\code{y}
is omitted or if \code{relative = False}, absolute error is returned.\\

\code{norm(x)} & 2-norm of \code{x}.\\

\code{corr(x, y=None)} & Correlation of \code{x} and \code{y}. If
\code{y} is \code{None} returns autocorrelation of \code{x}.\\

\code{ensure\_numeric(A, typecode = None)} & Returns a Numeric array
for any sequence \code{A}. If \code{A} is already a Numeric array it
will be returned unaltered. Otherwise, an attempt is made to convert
it to a Numeric array. (Needed because \code{array(A)} can
cause memory overflow.)\\

\code{histogram(a, bins, relative=False)} & Standard histogram. If
\code{relative} is \code{True}, values will be normalised against
the total and thus represent frequencies rather than counts.\\

\code{create\_bins(data, number\_of\_bins = None)} & Safely create
bins for use with histogram. If \code{data} contains only one point
or is constant, one bin will be created. If \code{number\_of\_bins}
is omitted, 10 bins will be created.\\  \hline

\end{tabular}


\chapter{Modules available in \anuga}


\section{\module{pyvolution.general\_mesh} }
\declaremodule[pyvolution.generalmesh]{}{pyvolution.general\_mesh}
\label{mod:pyvolution.generalmesh}

\section{\module{pyvolution.neighbour\_mesh} }
\declaremodule[pyvolution.neighbourmesh]{}{pyvolution.neighbour\_mesh}
\label{mod:pyvolution.neighbourmesh}

\section{\module{pyvolution.domain} --- Generic module for 2D triangular domains for finite-volume computations of conservation laws}
\declaremodule{}{pyvolution.domain}
\label{mod:pyvolution.domain}


\section{\module{pyvolution.quantity}}
\declaremodule{}{pyvolution.quantity}
\label{mod:pyvolution.quantity}


Class Quantity - Implements values at each triangular element

To create:

   Quantity(domain, vertex_values)

   domain: Associated domain structure. Required.

   vertex_values: N x 3 array of values at each vertex for each element.
                  Default None

   If vertex_values are None Create array of zeros compatible with domain.
   Otherwise check that it is compatible with dimenions of domain.
   Otherwise raise an exception



\section{\module{pyvolution.shallow\_water} --- 2D triangular domains for finite-volume
computations of the shallow water wave equation. This module contains a specialisation
of class Domain from module domain.py consisting of methods specific to the Shallow Water
Wave Equation
}
\declaremodule[pyvolution.shallowwater]{}{pyvolution.shallow\_water}
\label{mod:pyvolution.shallowwater}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Frequently Asked Questions}


\section{General Questions}

\subsubsection{What is \anuga?}

\subsubsection{Why is it called \anuga?}

\subsubsection{How do I obtain a copy of \anuga?}

\subsubsection{What developments are expected for \anuga in the future?}

\subsubsection{Are there any published articles about \anuga that I can reference?}

\section{Modelling Questions}

\subsubsection{Which type of problems are \anuga good for?}

\subsubsection{Which type of problems are beyond the scope of \anuga?}

\subsubsection{Can I start the simulation at an arbitrary time?}
Yes, using \code{domain.set\_time()} you can specify an arbitrary
starting time. This is for example useful in conjunction with a
file\_boundary, which may start hours before anything hits the model
boundary. By assigning a later time for the model to start,
computational resources aren't wasted.

\subsubsection{Can I change values for any quantity during the simulation?}
Yes, using \code{domain.set\_quantity()} inside the domain.evolve
loop you can change values of any quantity. This is for example
useful if you wish to let the system settle for a while before
assigning an initial condition. Another example would be changing
the values for elevation to model e.g. erosion.

\subsubsection{Can I change boundary conditions during the simulation?}
Not sure, but it would be nice :-)

\subsubsection{Why does a file\_function return a list of numbers when evaluated?}
Currently, file\_function works by returning values for the conserved
quantities \code{stage}, \code{xmomentum} and \code{ymomentum} at a given point in time
and space as a triplet. To access e.g.\ \code{stage} one must specify element 0 of the
triplet returned by file\_function.

\subsubsection{Which diagnostics are available to troubleshoot a simulation?}

\subsubsection{How do I use a DEM in my simulation?}
You use \code{dem2pts} to convert your DEM to the required .pts format. This .pts file is then called
when setting the elevation data to the mesh in \code{domain.set_quantity}

\subsubsection{What sort of DEM resolution should I use?}
Try and work with the \emph{best} you have available. Onshore DEMs
are typically available in 25m, 100m and 250m grids. Note, offshore
data is often sparse, or non-existent.

\subsubsection{What sort of mesh resolution should I use?}
The mesh resolution should be commensurate with your DEM - it does not make sense to put in place a mesh which is finer than your DEM. As an example,
if your DEM is on a 25m grid, then the cell resolution should be of the order of 315$m^2$ (this represents half the area of the square grid). Ideally,
you need a fine mesh over regions where the DEM changes rapidly, and other areas of significant interest, such as the coast.


\subsubsection{How do I tag interior polygons?}
At the moment create_mesh_from_regions does not allow interior
polygons with symbolic tags. If tags are needed, the interior
polygons must be created subsequently. For example, given a filename
of polygons representing solid walls (in Arc Ungenerate format) can
be tagged as such using the code snippet:
\begin{verbatim}
  # Create mesh outline with tags
  mesh = create_mesh_from_regions(bounding_polygon,
                                  boundary_tags=boundary_tags)
  # Add buildings outlines with tags set to 'wall'. This would typically
  # bind to a Reflective boundary
  mesh.import_ungenerate_file(buildings_filename, tag='wall')

  # Generate and write mesh to file
  mesh.generate_mesh(maximum_triangle_area=max_area)
  mesh.export_mesh_file(mesh_filename)
\end{verbatim}

Note that a mesh object is returned from \code{create_mesh_from_regions}
when file name is omitted.

\subsubsection{How often should I store the output?}
This will depend on what you are trying to answer with your model and how much memory you have available on your machine. If you need
to look in detail at the evolution, then you will need to balance your storage requirements and the duration of the simulation.
If the SWW file exceeds 1Gb, another SWW file will be created until the end of the simulation. As an example, to store all the conserved
quantities on a mesh with approximately 300000 triangles on a 2 min interval for 5 hours will result in approximately 350Mb SWW file
(as for the \file{run\_sydney\_smf.py} example).

\subsection{Boundary Conditions}

\subsubsection{How do I create a Dirichlet boundary condition?}

A Dirichlet boundary condition sets a constant value for the
conserved quantities at the boundaries. A list containing
the constant values for stage, xmomentum and ymomentum is constructed
and used in the function call, e.g. \code{Dirichlet_boundary([0.2,0.,0.])}

\subsubsection{How do I know which boundary tags are available?}
The method \code{domain.get\_boundary\_tags()} will return a list of
available tags for use with
\code{domain.set\_boundary\_condition()}.





\chapter{Glossary}

\begin{tabular}{|lp{10cm}|c|}  \hline
%\begin{tabular}{|llll|}  \hline
    \emph{Term} & \emph{Definition} & \emph{Page}\\  \hline

    \indexedbold{\anuga} & Name of software (joint development between ANU and
    GA) & \pageref{def:anuga}\\

    \indexedbold{bathymetry} & offshore elevation &\\

    \indexedbold{conserved quantity} & conserved (stage, x and y
    momentum) & \\

%    \indexedbold{domain} & The domain of a function is the set of all input values to the
%    function.&\\

    \indexedbold{Digital Elevation Model (DEM)} & DEMs are digital files consisting of points of elevations,
sampled systematically at equally spaced intervals.& \\

    \indexedbold{Dirichlet boundary} & A boundary condition imposed on a differential equation
 that specifies the values the solution is to take on the boundary of the
 domain. & \pageref{def:dirichlet boundary}\\

    \indexedbold{edge} & A triangular cell within the computational mesh can be depicted
    as a set of vertices joined by lines (the edges). & \\

    \indexedbold{elevation} & refers to bathymetry and topography &\\

    \indexedbold{evolution} & integration of the shallow water wave equations
    over time &\\

    \indexedbold{finite volume method} & The method evaluates the terms in the shallow water
    wave equation as fluxes at the surfaces of each finite volume. Because the
    flux entering a given volume is identical to that leaving the adjacent volume,
    these methods are conservative. Another advantage of the finite volume method is
    that it is easily formulated to allow for unstructured meshes. The method is used
    in many computational fluid dynamics packages. & \\

    \indexedbold{forcing term} & &\\

    \indexedbold{flux} & the amount of flow through the volume per unit
    time & \\

    \indexedbold{grid} & Evenly spaced mesh & \\

    \indexedbold{latitude} & The angular distance on a mericlear north and south of the
    equator, expressed in degrees and minutes. & \\

    \indexedbold{longitude} & The angular distance east or west, between the meridian
    of a particular place on Earth and that of the Prime Meridian (located in Greenwich,
    England) expressed in degrees or time.& \\

    \indexedbold{Manning friction coefficient} & &\\

    \indexedbold{mesh} & Triangulation of domain &\\

    \indexedbold{mesh file} & A TSH or MSH file & \pageref{def:mesh file}\\

    \indexedbold{NetCDF} & &\\

    \indexedbold{northing} & A rectangular (x,y) coordinate measurement of distance
    north from a north-south reference line, usually a meridian used as the axis of
    origin within a map zone or projection. Northing is a UTM (Universal Transverse
    Mercator) coordinate. & \\

    \indexedbold{points file} & A PTS or XYA file & \\  \hline

    \end{tabular}

    \begin{tabular}{|lp{10cm}|c|}  \hline

    \indexedbold{polygon} & A sequence of points in the plane. \anuga represents a polygon
    either as a list consisting of Python tuples or lists of length 2 or as an $N \times 2$
    Numeric array, where $N$ is the number of points.

    The unit square, for example, would be represented either as
    \code{[ [0,0], [1,0], [1,1], [0,1] ]} or as \code{array( [0,0], [1,0], [1,1],
    [0,1] )}.

    NOTE: For details refer to the module \module{utilities/polygon.py}. &
    \\     \indexedbold{resolution} &  The maximal area of a triangular cell in a
    mesh & \\


    \indexedbold{reflective boundary} & Models a solid wall. Returns same conserved
    quantities as those present in the neighbouring volume but reflected. Specific to the
    shallow water equation as it works with the momentum quantities assumed to be the
    second and third conserved quantities. & \pageref{def:reflective boundary}\\

    \indexedbold{stage} & &\\

%    \indexedbold{try this}

    \indexedbold{swollen} & visualisation tool used with \anuga & \pageref{sec:swollen}\\

    \indexedbold{time boundary} & Returns values for the conserved
quantities as a function of time. The user must specify
the domain to get access to the model time. & \pageref{def:time boundary}\\

    \indexedbold{topography} & onshore elevation &\\

    \indexedbold{transmissive boundary} & & \pageref{def:transmissive boundary}\\

    \indexedbold{vertex} & A point at which edges meet. & \\

    \indexedbold{xmomentum} & conserved quantity (note, two-dimensional SWW equations say
    only \code{x} and \code{y} and NOT \code{z}) &\\

    \indexedbold{ymomentum}  & conserved quantity & \\  \hline

    \end{tabular}


%The \code{\e appendix} markup need not be repeated for additional
%appendices.


%
%  The ugly "%begin{latexonly}" pseudo-environments are really just to
%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
%  not really valuable.
%
%  If you don't want the Module Index, you can remove all of this up
%  until the second \input line.
%

%begin{latexonly}
%\renewcommand{\indexname}{Module Index}
%end{latexonly}
\input{mod\jobname.ind}        % Module Index
%
%begin{latexonly}
%\renewcommand{\indexname}{Index}
%end{latexonly}
\input{\jobname.ind}            % Index



\begin{thebibliography}{99}
\bibitem[nielsen2005]{nielsen2005} 
{\it Hydrodynamic modelling of coastal inundation}.
Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman.
In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
Modelling and Simulation. Modelling and Simulation Society of Australia and
New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.\\
http://www.mssanz.org.au/modsim05/papers/nielsen.pdf




\end{thebibliography}{99}

\end{document}