source: documentation/user_manual/anuga_user_manual.tex @ 2893

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

\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-05-17 04:06:46 +0000 (Wed, 17 May 2006) $
%$LastChangedRevision: 2893 $
%$LastChangedBy: ole $

\input{copyright}


\begin{abstract}

\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 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 problem background, as
well as enough programming experience to adapt the code to
different requirements, as described in this manual,  and to
understand the basic terminology of object-oriented programming.

\section{Structure of This Manual}

FIXME (Ole): Isn't this superseded by the contents section?

This manual is structured as follows:

\begin{itemize}
  \item Background (What \anuga does)
  \item A \emph{Getting Started} section
  \item A detailed description of the public interface
  \item \anuga 's overall architecture, components and file formats
  \item Assumptions
\end{itemize}


\pagebreak
\chapter{Background}

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

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

\section{A Simple Example}

\subsection{Overview}

What follows is a discussion of the structure and operation of a script which we will call \file{runup.py}, with just enough detail to allow
the reader to appreciate what's involved in setting up a scenario
like the one it depicts.

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 many of the basic ideas involved in
setting up a more complex scenario. In the general case 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. 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 be in many typical cases.

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}


%\emph{[More details of the problem background]}

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

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 now.) The above assignment sets up a $10
\times 10$ rectangular mesh, triangulated in a specific way. In
general, the assignment

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

returns:

\begin{itemize}

   \item a list \code{points} of length $N$, where $N = (m + 1)(n + 1)$,
comprising the coordinates \code{(x, y)} of each of the $N$ mesh
points,

   \item a list \code{vertices} of length $2mn$ (each entry specifies the three
vertices of one of the triangles used in the triangulation) , and

   \item a dictionary \code{boundary}, used to tag the triangle edges on
the boundaries. Each key corresponds to a triangle edge on one of
the four boundaries and its value is one of \code{`left'},
\code{`right'}, \code{`top'} and \code{`bottom'}, indicating
which boundary the edge in question belongs to.

\end{itemize}

Since these three variables encapsulate the information needed to
specify the grid, it may be helpful to consider how they look in a
very simple case, not directly related to the example at hand.
Figure \ref{fig:simplemesh} represents a very simple mesh comprising
just 11 points and three triangles. (To avoid confusion, we should
emphasise that this particular mesh is \emph{not} generated by
\code{rectangular}---and is not even rectangular in nature. )


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




\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 uses a Python class \class{Domain}, imported from
\module{shallow\_water}, which is an extension of a more generic
class of the same name in the module \refmodule{pyvolution.domain}
(page \pageref{mod:pyvolution.domain}),
and inherits
some methods from the generic class but has others specific to the
shallow-water scenarios in which it is used. Specific options for
domain are set at this point. One of them is to set the basename for
the output file:

{\small \begin{verbatim}
    domain.set_name('bedslope')
\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
particularly flexible method that allows the user to set quantities
in a variety of ways---using constants, functions, numeric arrays or
expressions involving other quantities, arbitrary data points with
associated values, all of which can be passed as arguments. All
quantities can be initialised using \method{set\_quantity}. For
conserved quantities (\code{stage, xmomentum, ymomentum}) this is
called the \emph{initial condition}, for other quantities that
aren't updated by the equation, the same interface is used to assign
their values. 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.  %[screen shot?]

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 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 just 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
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.)

\subsection{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} 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} 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} Specifies a fixed value at the
      boundary and assigns initial values to the $x$-momentum and $y$-momentum.
    \item \textbf {Time boundary} Like a Dirichlet boundary but with behaviour 
      varying with time.
\end{itemize}

Once the boundary objects have been specified through the
statements above, they can be applied through a statement of the
form

{\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}

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

is the key step that causes the configuration of the domain to
`evolve' in accordance with the model embodied in the code, 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}

%Give details here of the form of the output and explain how it can
%be used with swollen. Include screen shots.//

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


\section{How to Run the Code}

The code can be run in various ways:

\begin{itemize}
  \item{from a Windows command line} as in \code{python runup.py}
  \item{within the Python IDLE environment}
  \item{within emacs}
  \item{from a Linux command line} as in \code{python runup.py}
\end{itemize}


\section{Exploring the model output}

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

The following discussion builds on the concepts introduced through
the \file{runup.py} example and introduces a second
example, \file{run\_sydney\_smf.py}, that follows the same basic
outline, but incorporates more complex features and refers to a
real-life scenario, rather than the artificial illustrative one used
in \file{runup.py}. 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 outputs 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{"runsydneysmf.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}

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

The following remarks may help the reader understand how
\code{pmesh} is used.

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 mesh points are
created inside the polygon through a random process. 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.

While a mesh created inside a polygon offers more flexibility than
one based on a rectangular grid, using \code{pmesh} in the limited
form we have described so far still doesn't provide a way to adapt
to geographic or other features in the landscape, whose presence may
require us to vary the resolution in the neighbourhood of the
features. To cope with this requirement, \code{pmesh} also allows
the user to specify a number of \emph{interior polygons}, which are
triangulated separately, each according to a separately specified
resolution. See Figure \ref{fig:interior meshes}.

\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 inside the bounding polygon, using the specified
resolution, and then creates a separate triangulation inside each of
the interior polygons, using the resolution specified for that
polygon.

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}. The resulting mesh is output to a
\emph{meshfile}\index{meshfile}. 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.)
\code{pmesh} assigns a name to the file by appending the extension
\code{.msh} to the name specified in the input file
\file{project.py}. This name is stored in the variable
\code{meshname}.

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'},
`right0', `right1', `right2', `top', `left1', `left2' and `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 meshfile
\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 bedslopephysical 
above causes Domain to convert a meshfile
\code{meshname} into an instance of the data structure
\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.03)
\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}

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:

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

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

\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{usecache} & Specifies whether caching is to be used\\
\emph{verbose} & If \code{True}, provides detailed terminal output
to the user\\  % \hline
\end{tabular}
%\end{center}

\section{Mesh Generation}
\refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}
\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. 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.
Additionally, the user specifies a list of boundary tags, one for
each edge of the bounding polygon.
\end{funcdesc}


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

% Translate following into layman's language
An instance of the class \class{Mesh} is used to store .  This can
then be used to build the outline of the mesh and then generate the
mesh.
\end{funcdesc}


\begin{funcdesc}  {add_region_from_polygon}{self, polygon, tags=None,
                                max_triangle_area=None, geo_reference=None}
Module: \module{pmesh.mesh.Mesh}

% Translate following into layman's language
This method is used to add a region to a \class{Mesh} instance.  The
region 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{funcdesc}

\begin{funcdesc}  {add_hole_from_polygon}{self, polygon, tags=None,
    geo_reference=None}
Module: \module{pmesh.mesh.Mesh}

% Translate following into layman's language
This method is used to add a `hole'---that is, a region where the
triangular mesh will not be generated---to a \class{Mesh} instance.
The region 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{funcdesc}

\begin{funcdesc}  {generate_mesh}{self,
                      maximum_triangle_area=None,
                      minimum_triangle_angle=28.0,
                      verbose=False}
Module: \module{pmesh.mesh.Mesh}

% Translate following into layman's language
This method is used to generate the triangular mesh.  The  maximal area
of any triangle in the mesh can be specified, along with the minimum
angle of all triangles.
\end{funcdesc}


\begin{funcdesc}  {export_mesh_file}{self,ofile}
Module: \module{pmesh.mesh.Mesh}

% Translate following into layman's language
This method is used to save the mesh to a file. \code{ofile} is the name of the mesh file to be writen,
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{funcdesc}

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

\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{funcdesc} {set\_name}{name}
    Module: \refmodule{pyvolution.domain}, page \pageref{mod:pyvolution.domain}  %\code{pyvolution.domain}

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

\begin{funcdesc} {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{funcdesc}

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

    Specifies the directory used for data, 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{funcdesc}

\begin{funcdesc} {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{funcdesc}

\begin{funcdesc} {set_time}{time=0.0}
    Module: \module{pyvolution.domain}

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

\begin{funcdesc} {set_default_order}{??}
\end{funcdesc}


%%%%%%
\section{Initial Conditions}

\begin{funcdesc}{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)
\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}


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.


\end{funcdesc}







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

\begin{funcdesc}{tsunami_slump}{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}
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.

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 from NetCDF file and returns
a callable object. 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.

To get access to the model time stored within the file funcion use 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, either as a sequence of numbers or
defined at the vertices of a triangular mesh (such
as those stored in \code{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}

\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{funcdesc}{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{funcdesc}

\begin{funcdesc} {get_boundary_tags}{}
Module: \module{pyvolution.mesh}
\end{funcdesc}

%%%
\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 the conserved
quantities.
\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}

The boundary values are obtained from a file and interpolated. The
file is assumed to contain a time series and possibly also spatial
information. The conserved quantities are given as a function of
time.
\end{classdesc}


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





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

  \begin{funcdesc}{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). 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{funcdesc}



\subsection{Diagnostics}

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


  \end{funcdesc}


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


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


  \begin{classdesc}{Geospatial_data}{???}
    Module: \module{geospatial_data.geo_spatial_data}
    Creates a georeferenced geospatial data object from either arrays or 
    a file (pts or xya).

    Objects of this class can be used with \method{set\_quantity}.
    
    FIXME (Ole): Describe methods such as get_attributes() etc
  \end{classdesc}





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

\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)\\
TSH & Triangular meshes (e.g. created from \code{pmesh})\\
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}
\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{Meshfile Formats}

A meshfile is a file that has a specific format suited to specifying
mesh data for \anuga. A meshfile 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 meshfile can be generated from the function
\function{create_mesh_from_regions} (see ) and used to initialise a
domain.

A meshfile describes 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 meshfile may contain `holes' and/or `regions'. A hole
or region is defined by specifying a point and a number of segments
that enclose that point. 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 meshfile 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 lines of this type:

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

        Here the first number specifies the zone (in this case zone 56)  and other numbers specify the
        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 has the following format...

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 and converts it to DEM format (ASC or
  ERS)
  \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}

  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_defualts, co_argcount)} and the function will be recomputed.

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

    {\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). 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.

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}\\
  \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}\\
  \code{-hmin <float>} & Height below which transparency is set to
                                     zero\\
  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
                                     up, default is overhead)\\
  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
  \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}

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


  \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 a single point---or the members of a list of points---
  are inside the specified polygon. If the first argument is a single point,
  returns \code{True} if the point is inside the polygon, or \code{False}
  otherwise. If the first argument is a list of points, 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}{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}



\section{coordinate_transforms}

\section{geo_spatial_data}

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

TBA

\section{pmesh GUI}

\section{alpha_shape}


\section{utilities/numerical_tools} Do now.

\begin{itemize}
  \item \indexedcode{ensure_numeric}
  \item \indexedcode{mean}
  \item
\end{itemize}


\chapter{Modules available in \anuga}


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

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

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


\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 "best" you have available. Onshore DEMs are typically available in 25m, 100m and 250m grids. Note, offshore data is often sparse,
or non-existant.

\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 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 against 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.

\subsection{Boundary Conditions}

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

\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{itemize}
    \item \indexedbold{\anuga} Name of software (joint development between ANU and GA)

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

    \item \indexedbold{Dirichlet boundary} - A Dirichlet boundary condition imposed on a differential equation
 which specifies the values the solution is to take on the boundary of the domain.

    \item \indexedbold{elevation} - refers to bathymetry and topography

    \item \indexedbold{bathymetry} - offshore elevation

    \item \indexedbold{topography} - onshore elevation

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

    \item \indexedbold{forcing term}

    \item \indexedbold{IDLE} - Development environment shipped with Python

    \item \indexedbold{Manning friction coefficient}

    \item \indexedbold{mesh}    - Triangulation of domain

    \item \indexedbold{meshfile}  [generic word for either .tsh or
    .msh file]

    \item \indexedbold{points file}  [generic word for either .pts or
    .xya file]

    \item \indexedbold{grid} - evenly spaced mesh

    \item \indexedbold{NetCDF}

    \item \indexedbold{pmesh} does this really need to be here? it's a class/module?

    \item \indexedbold{pyvolution} does this really need to be here? it's a class/module?

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

    \item \indexedbold{reflective boundary}

    \item \indexedbold{smoothing} is this really needed?

    \item \indexedbold{stage}

%    \item \indexedbold{try this}

    \item \indexedbold{swollen} - visualisation tool

    \item \indexedbold{time boundary} - defined in the manual (flog from there)

    \item \indexedbold{transmissive boundary} - defined in the manual (flog from there)

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

    \item \indexedbold{ymomentum}  - conserved quantity

    \item \indexedbold{resolution} -  The maximal area of a triangular cell in a mesh

    \item \indexedbold{polygon} - A sequence of points in the plane. (Arbitrary polygons can be created
    in this way.)
    \anuga represents a polygon in one of two ways. One way is to represent it as a
    list whose members are either Python tuples
    or Python lists of length 2. The unit square, for example, would be represented by the
    list
    [ [0,0], [1,0], [1,1], [0,1] ]. The alternative is to represent it as an
    $N \times 2$ Numeric array, where $N$ is the number of points.

    NOTE: More can be read in the module utilities/polygon.py ....

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

    \item \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.


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

    \item \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.

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

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

    \item \indexedbold{finite volume} - The method evaluates the terms in the shallow water wave equationas 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.


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

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


\end{itemize}

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



\end{document}