source: anuga_core/documentation/user_manual/anuga_user_manual.tex @ 7342

Last change on this file since 7342 was 7342, checked in by ole, 15 years ago

Refactored stofrage of quantities to use new concept of static and dynamic quantities.
This is in preparation for flexibly allowing quantities such as elevation or friction
to be time dependent.

All tests and the okushiri validation pass.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 207.6 KB
Line 
1% Complete documentation on the extended LaTeX markup used for Python
2% documentation is available in ''Documenting Python'', which is part
3% of the standard documentation for Python.  It may be found online
4% at:
5%
6%     http://www.python.org/doc/current/doc/doc.html
7
8%labels
9%Sections and subsections \label{sec: }
10%Chapters \label{ch: }
11%Equations \label{eq: }
12%Figures \label{fig: }
13
14% Is latex failing with;
15% 'modanuga_user_manual.ind' not found?
16% try this command-line
17%   makeindex modanuga_user_manual.idx
18% To produce the modanuga_user_manual.ind file.
19
20%%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%%
21%
22% ensure_geospatial
23% ensure_absolute
24% set_geo_reference
25
26\documentclass{manual}
27
28\usepackage{graphicx}
29\usepackage[english]{babel}
30\usepackage{datetime}
31\usepackage[hang,small,bf]{caption}
32
33\input{definitions}
34
35%%%%%
36% Set penalties for widows, etc, very high
37%%%%%
38
39\widowpenalty=10000
40\clubpenalty=10000
41\raggedbottom
42
43\title{\anuga User Manual}
44\author{Geoscience Australia and the Australian National University}
45
46% Please at least include a long-lived email address;
47% the rest is at your discretion.
48\authoraddress{Geoscience Australia \\
49  Email: \email{ole.nielsen@ga.gov.au}
50}
51
52%Draft date
53
54% update before release!
55% Use an explicit date so that reformatting
56% doesn't cause a new date to be used.  Setting
57% the date to \today can be used during draft
58% stages to make it easier to handle versions.
59
60\longdate       % Make date format long using datetime.sty
61%\settimeformat{xxivtime} % 24 hour Format
62\settimeformat{oclock} % Verbose
63\date{\today, \ \currenttime}
64%\hyphenation{set\_datadir}
65
66\ifhtml
67  \date{\today} % latex2html does not know about datetime
68\fi
69
70\input{version} % Get version info - this file may be modified by
71                % update_anuga_user_manual.py - if not a dummy
72                % will be used.
73
74%\release{1.0}   % release version; this is used to define the
75%                % \version macro
76
77\makeindex          % tell \index to actually write the .idx file
78\makemodindex       % If this contains a lot of module sections.
79
80\setcounter{tocdepth}{3}
81\setcounter{secnumdepth}{3}
82
83%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84
85\begin{document}
86\maketitle
87
88% This makes the contents more accessible from the front page of the HTML.
89\ifhtml
90  \chapter*{Front Matter\label{front}}
91\fi
92
93%Subversion keywords:
94%
95%$LastChangedDate: 2009-08-08 15:20:42 +0000 (Sat, 08 Aug 2009) $
96%$LastChangedRevision: 7342 $
97%$LastChangedBy: ole $
98
99\input{copyright}
100
101\begin{abstract}
102\label{def:anuga}
103
104\noindent \anuga\index{\anuga} is a hydrodynamic modelling tool that
105allows users to model realistic flow problems in complex 2D geometries.
106Examples include dam breaks or the effects of natural hazards such
107as riverine flooding, storm surges and tsunami.
108
109The user must specify a study area represented by a mesh of triangular
110cells, the topography and bathymetry, frictional resistance, initial
111values for water level (called \emph{stage}\index{stage} within \anuga),
112boundary conditions and forces such as rainfall, stream flows, windstress or pressure gradients if applicable.
113
114\anuga tracks the evolution of water depth and horizontal momentum
115within each cell over time by solving the shallow water wave equation
116governing equation using a finite-volume method.
117
118\anuga also incorporates a mesh generator that
119allows the user to set up the geometry of the problem interactively as
120well as tools for interpolation and surface fitting, and a number of
121auxiliary tools for visualising and interrogating the model output.
122
123Most \anuga components are written in the object-oriented programming
124language Python and most users will interact with \anuga by writing
125small Python programs based on the \anuga library
126functions. Computationally intensive components are written for
127efficiency in C routines working directly with Python numpy structures.
128
129\end{abstract}
130
131\tableofcontents
132
133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134
135\chapter{Introduction}
136
137
138\section{Purpose}
139
140The purpose of this user manual is to introduce the new user to the
141inundation software system, describe what it can do and give step-by-step
142instructions for setting up and running hydrodynamic simulations.
143The stable release of \anuga and this manual are available on sourceforge ati
144\url{http://sourceforge.net/projects/anuga}. A snapshot of work in progress is
145available through the \anuga software repository at
146\url{https://datamining.anu.edu.au/svn/ga/anuga_core}
147where the more adventurous reader might like to go.
148
149This manual describes \anuga version 1.0. To check for later versions of this manual
150go to \url{https://datamining.anu.edu.au/anuga}.
151
152\section{Scope}
153
154This manual covers only what is needed to operate the software after
155installation and configuration. It does not includes instructions
156for installing the software or detailed API documentation, both of
157which will be covered in separate publications and by documentation
158in the source code.
159
160The latest installation instructions may be found at:
161\url{http://datamining.anu.edu.au/\~{}ole/anuga/user_manual/anuga_installation_guide.pdf}.
162
163\section{Audience}
164
165Readers are assumed to be familiar with the Python Programming language and
166its object oriented approach.
167Python tutorials include
168\url{http://docs.python.org/tut},
169\url{http://www.sthurlow.com/python}, and
170%\url{http://datamining.anu.edu.au/\%7e ole/work/teaching/ctac2006/exercise1.pdf}.
171\url{http://datamining.anu.edu.au/\~{}ole/work/teaching/ctac2006/exercise1.pdf}.
172
173Readers also need to have a general understanding of scientific modelling,
174as well as enough programming experience to adapt the code to different
175requirements.
176
177\pagebreak
178
179%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180
181\chapter{Background}
182
183Modelling the effects on the built environment of natural hazards such
184as riverine flooding, storm surges and tsunami is critical for
185understanding their economic and social impact on our urban
186communities.  Geoscience Australia and the Australian National
187University are developing a hydrodynamic inundation modelling tool
188called \anuga to help simulate the impact of these hazards.
189
190The core of \anuga is the fluid dynamics module, called \code{shallow\_water},
191which is based on a finite-volume method for solving the Shallow Water
192Wave Equation.  The study area is represented by a mesh of triangular
193cells.  By solving the governing equation within each cell, water
194depth and horizontal momentum are tracked over time.
195
196A major capability of \anuga is that it can model the process of
197wetting and drying as water enters and leaves an area.  This means
198that it is suitable for simulating water flow onto a beach or dry land
199and around structures such as buildings.  \anuga is also capable
200of modelling hydraulic jumps due to the ability of the finite-volume
201method to accommodate discontinuities in the solution\footnote{
202While \anuga works with discontinuities in the conserved quantities stage,
203xmomentum and ymomentum, it does not allow discontinuities in the bed elevation.}.
204
205To set up a particular scenario the user specifies the geometry
206(bathymetry and topography), the initial water level (stage),
207boundary conditions such as tide, and any forcing terms that may
208drive the system such as rainfall, abstraction of water, wind stress or atmospheric pressure
209gradients. Gravity and frictional resistance from the different
210terrains in the model are represented by predefined forcing terms.
211See section \ref{sec:forcing terms} for details on forcing terms available in \anuga.
212
213The built-in mesh generator, called \code{graphical\_mesh\_generator},
214allows the user to set up the geometry
215of the problem interactively and to identify boundary segments and
216regions using symbolic tags.  These tags may then be used to set the
217actual boundary conditions and attributes for different regions
218(e.g.\ the Manning friction coefficient) for each simulation.
219
220Most \anuga components are written in the object-oriented programming
221language Python.  Software written in Python can be produced quickly
222and can be readily adapted to changing requirements throughout its
223lifetime.  Computationally intensive components are written for
224efficiency in C routines working directly with Python numeric
225structures.  The animation tool developed for \anuga is based on
226OpenSceneGraph, an Open Source Software (OSS) component allowing high
227level interaction with sophisticated graphics primitives.
228See \cite{nielsen2005} for more background on \anuga.
229
230\chapter{Restrictions and limitations on \anuga}
231\label{ch:limitations}
232
233Although a powerful and flexible tool for hydrodynamic modelling, \anuga has a
234number of limitations that any potential user needs to be aware of. They are:
235
236\begin{itemize}
237  \item The mathematical model is the 2D shallow water wave equation.
238  As such it cannot resolve vertical convection and consequently not breaking
239  waves or 3D turbulence (e.g.\ vorticity).
240  %\item The surface is assumed to be open, e.g.\ \anuga cannot model
241  %flow under ceilings or in pipes
242  \item All spatial coordinates are assumed to be UTM (meters). As such,
243  \anuga is unsuitable for modelling flows in areas larger than one UTM zone
244  (6 degrees wide).
245  \item Fluid is assumed to be inviscid -- i.e.\ no kinematic viscosity included.
246  \item The finite volume is a very robust and flexible numerical technique,
247  but it is not the fastest method around. If the geometry is sufficiently
248  simple and if there is no need for wetting or drying, a finite-difference
249  method may be able to solve the problem faster than \anuga.
250%\item Mesh resolutions near coastlines with steep gradients need to be...
251  \item Frictional resistance is implemented using Manning's formula, but
252  \anuga has not yet been fully validated in regard to bottom roughness.
253%\item \anuga contains no tsunami-genic functionality relating to earthquakes.
254\end{itemize}
255
256%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
257
258\chapter{Getting Started}
259\label{ch:getstarted}
260
261This section is designed to assist the reader to get started with
262\anuga by working through some examples. Two examples are discussed;
263the first is a simple example to illustrate many of the concepts, and
264the second is a more realistic example.
265
266
267\section{A Simple Example}
268\label{sec:simpleexample}
269
270\subsection{Overview}
271
272What follows is a discussion of the structure and operation of a
273script called \file{runup.py}.
274
275This example carries out the solution of the shallow-water wave
276equation in the simple case of a configuration comprising a flat
277bed, sloping at a fixed angle in one direction and having a
278constant depth across each line in the perpendicular direction.
279
280The example demonstrates the basic ideas involved in setting up a
281complex scenario. In general the user specifies the geometry
282(bathymetry and topography), the initial water level, boundary
283conditions such as tide, and any forcing terms that may drive the
284system such as rainfall, abstraction of water, wind stress or atmospheric pressure gradients.
285Frictional resistance from the different terrains in the model is
286represented by predefined forcing terms. In this example, the
287boundary is reflective on three sides and a time dependent wave on
288one side.
289
290The present example represents a simple scenario and does not
291include any forcing terms, nor is the data taken from a file as it
292would typically be.
293
294The conserved quantities involved in the
295problem are stage (absolute height of water surface),
296$x$-momentum and $y$-momentum. Other quantities
297involved in the computation are the friction and elevation.
298
299Water depth can be obtained through the equation:
300
301\begin{verbatim}
302depth = stage - elevation
303\end{verbatim}
304
305\subsection{Outline of the Program}
306
307In outline, \file{runup.py} performs the following steps:
308\begin{enumerate}
309   \item Sets up a triangular mesh.
310   \item Sets certain parameters governing the mode of
311         operation of the model, specifying, for instance,
312         where to store the model output.
313   \item Inputs various quantities describing physical measurements, such
314         as the elevation, to be specified at each mesh point (vertex).
315   \item Sets up the boundary conditions.
316   \item Carries out the evolution of the model through a series of time
317         steps and outputs the results, providing a results file that can
318         be viewed.
319\end{enumerate}
320
321\subsection{The Code}
322
323For reference we include below the complete code listing for
324\file{runup.py}. Subsequent paragraphs provide a
325'commentary' that describes each step of the program and explains it
326significance.
327
328\label{ref:runup_py_code}
329\verbatiminput{demos/runup.py}
330
331\subsection{Establishing the Mesh}\index{mesh, establishing}
332
333The first task is to set up the triangular mesh to be used for the
334scenario. This is carried out through the statement:
335
336\begin{verbatim}
337points, vertices, boundary = rectangular_cross(10, 10)
338\end{verbatim}
339
340The function \function{rectangular_cross} is imported from a module
341\module{mesh\_factory} defined elsewhere. (\anuga also contains
342several other schemes that can be used for setting up meshes, but we
343shall not discuss these.) The above assignment sets up a $10 \times
34410$ rectangular mesh, triangulated in a regular way. The assignment:
345
346\begin{verbatim}
347points, vertices, boundary = rectangular_cross(m, n)
348\end{verbatim}
349
350returns:
351
352\begin{itemize}
353   \item a list \code{points} giving the coordinates of each mesh point,
354   \item a list \code{vertices} specifying the three vertices of each triangle, and
355   \item a dictionary \code{boundary} that stores the edges on
356         the boundary and associates each with one of the symbolic tags \code{'left'}, \code{'right'},
357         \code{'top'} or \code{'bottom'}.
358\end{itemize}
359
360(For more details on symbolic tags, see page
361\pageref{ref:tagdescription}.)
362
363An example of a general unstructured mesh and the associated data
364structures \code{points}, \code{vertices} and \code{boundary} is
365given in Section \ref{sec:meshexample}.
366
367\subsection{Initialising the Domain}
368
369These variables are then used to set up a data structure
370\code{domain}, through the assignment:
371
372\begin{verbatim}
373domain = Domain(points, vertices, boundary)
374\end{verbatim}
375
376This creates an instance of the \class{Domain} class, which
377represents the domain of the simulation. Specific options are set at
378this point, including the basename for the output file and the
379directory to be used for data:
380
381\begin{verbatim}
382domain.set_name('runup')
383domain.set_datadir('.')
384\end{verbatim}
385
386In addition, the following statement could be used to state that
387quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are
388to be stored at every timestep and \code{elevation} only once at
389the beginning of the simulation:
390
391\begin{verbatim}
392domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
393\end{verbatim}
394
395However, this is not necessary, as the above is the default behaviour.
396
397\subsection{Initial Conditions}
398
399The next task is to specify a number of quantities that we wish to
400set for each mesh point. The class \class{Domain} has a method
401\method{set\_quantity}, used to specify these quantities. It is a
402flexible method that allows the user to set quantities in a variety
403of ways -- using constants, functions, numeric arrays, expressions
404involving other quantities, or arbitrary data points with associated
405values, all of which can be passed as arguments. All quantities can
406be initialised using \method{set\_quantity}. For a conserved
407quantity (such as \code{stage, xmomentum, ymomentum}) this is called
408an \emph{initial condition}. However, other quantities that aren't
409updated by the equation are also assigned values using the same
410interface. The code in the present example demonstrates a number of
411forms in which we can invoke \method{set\_quantity}.
412
413\subsubsection{Elevation}
414
415The elevation, or height of the bed, is set using a function
416defined through the statements below, which is specific to this
417example and specifies a particularly simple initial configuration
418for demonstration purposes:
419
420\begin{verbatim}
421def topography(x, y):
422    return -x/2
423\end{verbatim}
424
425This simply associates an elevation with each point \code{(x, y)} of
426the plane.  It specifies that the bed slopes linearly in the
427\code{x} direction, with slope $-\frac{1}{2}$,  and is constant in
428the \code{y} direction.
429
430Once the function \function{topography} is specified, the quantity
431\code{elevation} is assigned through the simple statement:
432
433\begin{verbatim}
434domain.set_quantity('elevation', topography)
435\end{verbatim}
436
437NOTE: If using function to set \code{elevation} it must be vector
438compatible. For example, using square root will not work.
439
440\subsubsection{Friction}
441
442The assignment of the friction quantity (a forcing term)
443demonstrates another way we can use \method{set\_quantity} to set
444quantities -- namely, assign them to a constant numerical value:
445
446\begin{verbatim}
447domain.set_quantity('friction', 0.1)
448\end{verbatim}
449
450This specifies that the Manning friction coefficient is set to 0.1
451at every mesh point.
452
453\subsubsection{Stage}
454
455The stage (the height of the water surface) is related to the
456elevation and the depth at any time by the equation:
457
458\begin{verbatim}
459stage = elevation + depth
460\end{verbatim}
461
462For this example, we simply assign a constant value to \code{stage},
463using the statement:
464
465\begin{verbatim}
466domain.set_quantity('stage', -0.4)
467\end{verbatim}
468
469which specifies that the surface level is set to a height of $-0.4$,
470i.e.\ 0.4 units (metres) below the zero level.
471
472Although it is not necessary for this example, it may be useful to
473digress here and mention a variant to this requirement, which allows
474us to illustrate another way to use \method{set\_quantity} -- namely,
475incorporating an expression involving other quantities. Suppose,
476instead of setting a constant value for the stage, we wished to
477specify a constant value for the \emph{depth}. For such a case we
478need to specify that \code{stage} is everywhere obtained by adding
479that value to the value already specified for \code{elevation}. We
480would do this by means of the statements:
481
482\begin{verbatim}
483h = 0.05    # Constant depth
484domain.set_quantity('stage', expression='elevation + %f' % h)
485\end{verbatim}
486
487That is, the value of \code{stage} is set to $\code{h} = 0.05$ plus
488the value of \code{elevation} already defined.
489
490The reader will probably appreciate that this capability to
491incorporate expressions into statements using \method{set\_quantity}
492greatly expands its power. See Section \ref{sec:initial conditions} for more
493details.
494
495\subsection{Boundary Conditions}\index{boundary conditions}
496
497The boundary conditions are specified as follows:
498
499\begin{verbatim}
500Br = Reflective_boundary(domain)
501Bt = Transmissive_boundary(domain)
502Bd = Dirichlet_boundary([0.2, 0.0, 0.0])
503Bw = Time_boundary(domain=domain,
504                   f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
505\end{verbatim}
506
507The effect of these statements is to set up a selection of different
508alternative boundary conditions and store them in variables that can be
509assigned as needed. Each boundary condition specifies the
510behaviour at a boundary in terms of the behaviour in neighbouring
511elements. The boundary conditions introduced here may be briefly described as
512follows:
513\begin{itemize}
514    \item \textbf{Reflective boundary}\label{def:reflective boundary}
515          Returns same \code{stage} as in its neighbour volume but momentum
516          vector reversed 180 degrees (reflected).
517          Specific to the shallow water equation as it works with the
518          momentum quantities assumed to be the second and third conserved
519          quantities. A reflective boundary condition models a solid wall.
520    \item \textbf{Transmissive boundary}\label{def:transmissive boundary} 
521          Returns same conserved quantities as
522          those present in its neighbour volume. This is one way of modelling
523          outflow from a domain, but it should be used with caution if flow is
524          not steady state as replication of momentum at the boundary
525          may cause numerical instabilities propagating into the domain and
526          eventually causing \anuga to crash. If this occurs,
527          consider using e.g.\ a Dirichlet boundary condition with a stage value
528          less than the elevation at the boundary.
529    \item \textbf{Dirichlet boundary}\label{def:dirichlet boundary} Specifies
530          constant values for stage, $x$-momentum and $y$-momentum at the boundary.
531    \item \textbf{Time boundary}\label{def:time boundary} Like a Dirichlet
532          boundary but with behaviour varying with time.
533\end{itemize}
534
535\label{ref:tagdescription}Before describing how these boundary
536conditions are assigned, we recall that a mesh is specified using
537three variables \code{points}, \code{vertices} and \code{boundary}.
538In the code we are discussing, these three variables are returned by
539the function \code{rectangular}. The example given in
540Section \ref{sec:realdataexample} illustrates another way of
541assigning the values, by means of the function
542\code{create_mesh_from_regions}.
543
544These variables store the data determining the mesh as follows. (You
545may find that the example given in Section \ref{sec:meshexample}
546helps to clarify the following discussion, even though that example
547is a \emph{non-rectangular} mesh.)
548\begin{itemize}
549    \item The variable \code{points} stores a list of 2-tuples giving the
550          coordinates of the mesh points.
551    \item The variable \code{vertices} stores a list of 3-tuples of
552          numbers, representing vertices of triangles in the mesh. In this
553          list, the triangle whose vertices are \code{points[i]},
554          \code{points[j]}, \code{points[k]} is represented by the 3-tuple
555          \code{(i, j, k)}.
556    \item The variable \code{boundary} is a Python dictionary that
557          not only stores the edges that make up the boundary but also assigns
558          symbolic tags to these edges to distinguish different parts of the
559          boundary. An edge with endpoints \code{points[i]} and
560          \code{points[j]} is represented by the 2-tuple \code{(i, j)}. The
561          keys for the dictionary are the 2-tuples \code{(i, j)} corresponding
562          to boundary edges in the mesh, and the values are the tags are used
563          to label them. In the present example, the value \code{boundary[(i, j)]}
564          assigned to \code{(i, j)]} is one of the four tags
565          \code{'left'}, \code{'right'}, \code{'top'} or \code{'bottom'},
566          depending on whether the boundary edge represented by \code{(i, j)}
567          occurs at the left, right, top or bottom of the rectangle bounding
568          the mesh. The function \code{rectangular} automatically assigns
569          these tags to the boundary edges when it generates the mesh.
570\end{itemize}
571
572The tags provide the means to assign different boundary conditions
573to an edge depending on which part of the boundary it belongs to.
574(In Section \ref{sec:realdataexample} we describe an example that
575uses different boundary tags -- in general, the possible tags are entirely selectable by the user when generating the mesh and not
576limited to 'left', 'right', 'top' and 'bottom' as in this example.)
577All segments in bounding polygon must be tagged. If a tag is not supplied, the default tag name 'exterior' will be assigned by \anuga.
578
579Using the boundary objects described above, we assign a boundary
580condition to each part of the boundary by means of a statement like:
581
582\begin{verbatim}
583domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
584\end{verbatim}
585
586It is critical that all tags are associated with a boundary condition in this statement.
587If not the program will halt with a statement like:
588
589\begin{verbatim}
590Traceback (most recent call last):
591  File "mesh_test.py", line 114, in ?
592    domain.set_boundary({'west': Bi, 'east': Bo, 'north': Br, 'south': Br})
593  File "X:\inundation\sandpits\onielsen\anuga_core\source\anuga\abstract_2d_finite_volumes\domain.py", line 505, in set_boundary
594    raise msg
595ERROR (domain.py): Tag "exterior" has not been bound to a boundary object.
596All boundary tags defined in domain must appear in the supplied dictionary.
597The tags are: ['ocean', 'east', 'north', 'exterior', 'south']
598\end{verbatim}
599
600The command \code{set_boundary} stipulates that, in the current example, the right
601boundary varies with time, as defined by the lambda function, while the other
602boundaries are all reflective.
603
604The reader may wish to experiment by varying the choice of boundary
605types for one or more of the boundaries. (In the case of \code{Bd}
606and \code{Bw}, the three arguments in each case represent the
607\code{stage}, $x$-momentum and $y$-momentum, respectively.)
608
609\begin{verbatim}
610Bw = Time_boundary(domain=domain, f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
611\end{verbatim}
612
613\subsection{Evolution}\index{evolution}
614
615The final statement:
616
617\begin{verbatim}
618for t in domain.evolve(yieldstep=0.1, duration=10.0):
619    print domain.timestepping_statistics()
620\end{verbatim}
621
622causes the configuration of the domain to 'evolve', over a series of
623steps indicated by the values of \code{yieldstep} and
624\code{duration}, which can be altered as required.  The value of
625\code{yieldstep} controls the time interval between successive model
626outputs.  Behind the scenes more time steps are generally taken.
627
628\subsection{Output}
629
630The output is a NetCDF file with the extension \code{.sww}. It
631contains stage and momentum information and can be used with the
632\anuga viewer \code{anuga\_viewer} to generate a visual
633display (see Section \ref{sec:animate}). See Section \ref{sec:file formats}
634(page \pageref{sec:file formats}) for more on NetCDF and other file
635formats.
636
637The following is a listing of the screen output seen by the user
638when this example is run:
639
640\verbatiminput{examples/runupoutput.txt}
641
642
643\section{How to Run the Code}
644
645The code can be run in various ways:
646\begin{itemize}
647  \item{from a Windows or Unix command line} as in\ \code{python runup.py}
648  \item{within the Python IDLE environment}
649  \item{within emacs}
650  \item{within Windows, by double-clicking the \code{runup.py}
651  file.}
652\end{itemize}
653
654
655\section{Exploring the Model Output}
656
657The following figures are screenshots from the \anuga visualisation
658tool \code{animate}. Figure \ref{fig:runupstart} shows the domain
659with water surface as specified by the initial condition, $t=0$.
660Figure \ref{fig:runup2} shows later snapshots for $t=2.3$ and
661$t=4$ where the system has been evolved and the wave is encroaching
662on the previously dry bed.
663
664\code{animate} is described in more detail in Section \ref{sec:animate}.
665
666\begin{figure}[htp]
667  \centerline{\includegraphics[width=75mm, height=75mm]
668    {graphics/bedslopestart.jpg}}
669  \caption{Runup example viewed with the \anuga viewer}
670  \label{fig:runupstart}
671\end{figure}
672
673\begin{figure}[htp]
674  \centerline{
675    \includegraphics[width=75mm, height=75mm]{graphics/bedslopeduring.jpg}
676    \includegraphics[width=75mm, height=75mm]{graphics/bedslopeend.jpg}
677   }
678  \caption{Runup example viewed with ANGUA viewer}
679  \label{fig:runup2}
680\end{figure}
681
682\clearpage
683
684
685\section{A slightly more complex example}
686\label{sec:channelexample}
687
688\subsection{Overview}
689
690The next example is about waterflow in a channel with varying boundary conditions and
691more complex topographies. These examples build on the
692concepts introduced through the \file{runup.py} in Section \ref{sec:simpleexample}.
693The example will be built up through three progressively more complex scripts.
694
695\subsection{Overview}
696
697As in the case of \file{runup.py}, the actions carried
698out by the program can be organised according to this outline:
699\begin{enumerate}
700   \item Set up a triangular mesh.
701   \item Set certain parameters governing the mode of
702         operation of the model -- specifying, for instance, where to store the
703         model output.
704   \item Set up initial conditions for various quantities such as the elevation, to be specified at each mesh point (vertex).
705   \item Set up the boundary conditions.
706   \item Carry out the evolution of the model through a series of time
707         steps and output the results, providing a results file that can be
708         viewed.
709\end{enumerate}
710
711\subsection{The Code}
712
713Here is the code for the first version of the channel flow \file{channel1.py}:
714
715\verbatiminput{demos/channel1.py}
716
717In discussing the details of this example, we follow the outline
718given above, discussing each major step of the code in turn.
719
720\subsection{Establishing the Mesh}\index{mesh, establishing}
721
722In this example we use a similar simple structured triangular mesh as in \file{runup.py}
723for simplicity, but this time we will use a symmetric one and also
724change the physical extent of the domain. The assignment:
725
726\begin{verbatim}
727points, vertices, boundary = rectangular_cross(m, n, len1=length, len2=width)
728\end{verbatim}
729
730returns an \code{mxn} mesh similar to the one used in the previous example, except that now the
731extent in the x and y directions are given by the value of \code{length} and \code{width}
732respectively.
733
734Defining \code{m} and \code{n} in terms of the extent as in this example provides a convenient way of
735controlling the resolution: By defining \code{dx} and \code{dy} to be the desired size of each
736hypothenuse in the mesh we can write the mesh generation as follows:
737
738\begin{verbatim}
739length = 10.0
740width = 5.0
741dx = dy = 1           # Resolution: Length of subdivisions on both axes
742
743points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
744                                               len1=length, len2=width)
745\end{verbatim}
746
747which yields a mesh of length=10m, width=5m with 1m spacings. To increase the resolution,
748as we will later in this example, one merely decreases the values of \code{dx} and \code{dy}.
749
750The rest of this script is similar to the previous example on page \pageref{ref:runup_py_code}.
751% except for an application of the 'expression' form of \code{set\_quantity} where we use
752% the value of \code{elevation} to define the (dry) initial condition for \code{stage}:
753%\begin{verbatim}
754%  domain.set_quantity('stage', expression='elevation')
755%\end{verbatim}
756
757
758\section{Model Output}
759
760The following figure is a screenshot from the \anuga visualisation
761tool \code{animate} of output from this example.
762
763\begin{figure}[htp]
764  \centerline{\includegraphics[height=75mm]
765    {graphics/channel1.png}}%
766  \caption{Simple channel example viewed with the \anuga viewer.}
767  \label{fig:channel1}
768\end{figure}
769
770\subsection{Changing boundary conditions on the fly}
771\label{sec:change boundary}
772
773Here is the code for the second version of the channel flow \file{channel2.py}:
774
775\verbatiminput{demos/channel2.py}
776
777This example differs from the first version in that a constant outflow boundary condition has
778been defined:
779
780\begin{verbatim}
781Bo = Dirichlet_boundary([-5, 0, 0])    # Outflow
782\end{verbatim}
783
784and that it is applied to the right hand side boundary when the water level there exceeds 0m.
785
786\begin{verbatim}
787for t in domain.evolve(yieldstep=0.2, finaltime=40.0):
788    domain.write_time()
789
790    if domain.get_quantity('stage').get_values(interpolation_points=[[10, 2.5]]) > 0:
791        print 'Stage > 0: Changing to outflow boundary'
792        domain.set_boundary({'right': Bo})
793\end{verbatim}
794
795\label{sec:change boundary code}
796The \code{if} statement in the timestepping loop (\code{evolve}) gets the quantity
797\code{stage} and obtains the interpolated value at the point (10m,
7982.5m) which is on the right boundary. If the stage exceeds 0m a
799message is printed and the old boundary condition at tag 'right' is
800replaced by the outflow boundary using the method:
801
802\begin{verbatim}
803domain.set_boundary({'right': Bo})
804\end{verbatim}
805
806This type of dynamically varying boundary could for example be
807used to model the breakdown of a sluice door when water exceeds a certain level.
808
809\subsection{Output}
810
811The text output from this example looks like this:
812
813\begin{verbatim}
814...
815Time = 15.4000, delta t in [0.03789902, 0.03789916], steps=6 (6)
816Time = 15.6000, delta t in [0.03789896, 0.03789908], steps=6 (6)
817Time = 15.8000, delta t in [0.03789891, 0.03789903], steps=6 (6)
818Stage > 0: Changing to outflow boundary
819Time = 16.0000, delta t in [0.02709050, 0.03789898], steps=6 (6)
820Time = 16.2000, delta t in [0.03789892, 0.03789904], steps=6 (6)
821...
822\end{verbatim}
823
824\subsection{Flow through more complex topograhies}
825
826Here is the code for the third version of the channel flow \file{channel3.py}:
827
828\verbatiminput{demos/channel3.py}
829
830This example differs from the first two versions in that the topography
831contains obstacles.
832
833This is accomplished here by defining the function \code{topography} as follows:
834
835\begin{verbatim}
836def topography(x,y):
837    """Complex topography defined by a function of vectors x and y."""
838
839    z = -x/10
840
841    N = len(x)
842    for i in range(N):
843        # Step
844        if 10 < x[i] < 12:
845            z[i] += 0.4 - 0.05*y[i]
846
847        # Constriction
848        if 27 < x[i] < 29 and y[i] > 3:
849            z[i] += 2
850
851        # Pole
852        if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
853            z[i] += 2
854
855    return z
856\end{verbatim}
857
858In addition, changing the resolution to \code{dx = dy = 0.1} creates a finer mesh resolving the new features better.
859
860A screenshot of this model at time 15s is:
861\begin{figure}[htp]
862  \centerline{\includegraphics[height=75mm]
863    {graphics/channel3.png}}
864  \caption{More complex flow in a channel}
865  \label{fig:channel3}
866\end{figure}
867
868
869\section{An Example with Real Data}
870
871\label{sec:realdataexample} The following discussion builds on the
872concepts introduced through the \file{runup.py} example and
873introduces a second example, \file{runcairns.py}.  This refers to
874a {\bf hypothetical} scenario using real-life data,
875in which the domain of interest surrounds the
876Cairns region. Two scenarios are given; firstly, a
877hypothetical tsunami wave is generated by a submarine mass failure
878situated on the edge of the continental shelf, and secondly, a fixed wave
879of given amplitude and period is introduced through the boundary.
880
881{\bf
882Each scenario has been designed to generate a tsunami which will
883inundate the Cairns region. To achieve this, suitably large
884parameters were chosen and were not based on any known tsunami sources
885or realistic amplitudes.
886}
887
888\subsection{Overview}
889As in the case of \file{runup.py}, the actions carried
890out by the program can be organised according to this outline:
891\begin{enumerate}
892   \item Set up a triangular mesh.
893
894   \item Set certain parameters governing the mode of
895         operation of the model -- specifying, for instance, where to store the
896         model output.
897
898   \item Input various quantities describing physical measurements, such
899         as the elevation, to be specified at each mesh point (vertex).
900
901   \item Set up the boundary conditions.
902
903   \item Carry out the evolution of the model through a series of time
904         steps and output the results, providing a results file that can be
905         visualised.
906\end{enumerate}
907
908\subsection{The Code}
909
910Here is the code for \file{runcairns.py}:
911
912\verbatiminput{demos/cairns/runcairns.py}
913
914In discussing the details of this example, we follow the outline
915given above, discussing each major step of the code in turn.
916
917\subsection{Establishing the Mesh}\index{mesh, establishing}
918
919One obvious way that the present example differs from
920\file{runup.py} is in the use of a more complex method to
921create the mesh. Instead of imposing a mesh structure on a
922rectangular grid, the technique used for this example involves
923building mesh structures inside polygons specified by the user,
924using a mesh-generator.
925
926The mesh-generator creates the mesh within a single
927polygon whose vertices are at geographical locations specified by
928the user. The user specifies the \emph{resolution} -- that is, the
929maximal area of a triangle used for triangulation -- and a triangular
930mesh is created inside the polygon using a mesh generation engine.
931On any given platform, the same mesh will be returned each time the
932script is run.
933
934Boundary tags are not restricted to \code{'left'}, \code{'bottom'},
935\code{'right'} and \code{'top'}, as in the case of
936\file{runup.py}. Instead the user specifies a list of
937tags appropriate to the configuration being modelled.
938
939In addition, the mesh-generator provides a way to adapt to geographic or
940other features in the landscape, whose presence may require an
941increase in resolution. This is done by allowing the user to specify
942a number of \emph{interior polygons}, each with a specified
943resolution. It is also
944possible to specify one or more 'holes' -- that is, areas bounded by
945polygons in which no triangulation is required.
946
947In its general form, the mesh-generator takes for its input a bounding
948polygon and (optionally) a list of interior polygons. The user
949specifies resolutions, both for the bounding polygon and for each of
950the interior polygons. Given this data, the mesh-generator first creates a
951triangular mesh with varying resolution.
952
953The function used to implement this process is
954\function{create\_domain\_from\_regions} which creates a Domain object as
955well as a mesh file.  Its arguments include the
956bounding polygon and its resolution, a list of boundary tags, and a
957list of pairs \code{[polygon, resolution]} specifying the interior
958polygons and their resolutions.
959
960The resulting mesh is output to a \emph{mesh file}\index{mesh
961file}\label{def:mesh file}. This term is used to describe a file of
962a specific format used to store the data specifying a mesh. (There
963are in fact two possible formats for such a file: it can either be a
964binary file, with extension \code{.msh}, or an ASCII file, with
965extension \code{.tsh}. In the present case, the binary file format
966\code{.msh} is used. See Section \ref{sec:file formats} (page
967\pageref{sec:file formats}) for more on file formats.
968
969In practice, the details of the polygons used are read from a
970separate file \file{project.py}. Here is a complete listing of
971\file{project.py}:
972
973\verbatiminput{demos/cairns/project.py}
974
975Figure \ref{fig:cairns3d} illustrates the landscape of the region
976for the Cairns example. Understanding the landscape is important in
977determining the location and resolution of interior polygons. The
978supporting data is found in the ASCII grid, \code{cairns.asc}, which
979has been sourced from the publically available Australian Bathymetry
980and Topography Grid 2005, \cite{grid250}. The required resolution
981for inundation modelling will depend on the underlying topography and
982bathymetry; as the terrain becomes more complex, the desired resolution
983would decrease to the order of tens of metres.
984
985\clearpage
986
987\begin{figure}[htp]
988  \centerline{\includegraphics[scale=0.5]{graphics/cairns3.jpg}}
989  \caption{Landscape of the Cairns scenario.}
990  \label{fig:cairns3d}
991\end{figure}
992
993The following statements are used to read in the specific polygons
994from \code{project.cairns} and assign a defined resolution to
995each polygon.
996
997\begin{verbatim}
998islands_res = 100000
999cairns_res = 100000
1000shallow_res = 500000
1001interior_regions = [[project.poly_cairns,  cairns_res],
1002                    [project.poly_island0, islands_res],
1003                    [project.poly_island1, islands_res],
1004                    [project.poly_island2, islands_res],
1005                    [project.poly_island3, islands_res],
1006                    [project.poly_shallow, shallow_res]]
1007\end{verbatim}
1008
1009Figure \ref{fig:cairnspolys}
1010illustrates the polygons used for the Cairns scenario.
1011
1012\clearpage
1013
1014\begin{figure}[htp]
1015  \centerline{\includegraphics[scale=0.5]
1016      {graphics/cairnsmodel.jpg}}
1017  \caption{Interior and bounding polygons for the Cairns example.}
1018  \label{fig:cairnspolys}
1019\end{figure}
1020
1021The statement:
1022
1023\begin{verbatim}
1024remainder_res = 10000000
1025domain = create_domain_from_regions(project.bounding_polygon,
1026                                    boundary_tags={'top': [0],
1027                                                   'ocean_east': [1],
1028                                                   'bottom': [2],
1029                                                   'onshore': [3]},
1030                                    maximum_triangle_area=project.default_res,
1031                                    mesh_filename=project.meshname,
1032                                    interior_regions=project.interior_regions,
1033                                    use_cache=True,
1034                                    verbose=True)
1035\end{verbatim}
1036
1037is then used to create the mesh, taking the bounding polygon to be
1038the polygon \code{bounding\_polygon} specified in \file{project.py}.
1039The argument \code{boundary\_tags} assigns a dictionary, whose keys
1040are the names of the boundary tags used for the bounding
1041polygon -- \code{'top'}, \code{'ocean\_east'}, \code{'bottom'}, and
1042\code{'onshore'} -- and whose values identify the indices of the
1043segments associated with each of these tags.
1044The polygon may be arranged either clock-wise or counter clock-wise and the
1045indices refer to edges in the order they appear: Edge 0 connects vertex 0 and vertex 1, edge 1 connects vertex 1 and 2; and so forth.
1046(Here, the values associated with each boundary tag are one-element lists, but they can have as many indices as there are edges)
1047If polygons intersect, or edges coincide (or are even very close) the resolution may be undefined in some regions.
1048Use the underlying mesh interface for such cases
1049(see Chapter \ref{sec:mesh interface}).
1050If a segment is omitted in the tags definition an Exception is raised.
1051
1052Note that every point on each polygon defining the mesh will be used as vertices in triangles.
1053Consequently, polygons with points very close together will cause triangles with very small
1054areas to be generated irrespective of the requested resolution.
1055Make sure points on polygons are spaced to be no closer than the smallest resolution requested.
1056
1057\subsection{Initialising the Domain}
1058
1059Since we used \code{create_domain_from_regions} to create the mesh file, we do not need to
1060create the domain explicitly, as the above function does both mesh and domain creation.
1061
1062The following statements specify a basename and data directory, and
1063sets a minimum storable height, which helps with visualisation and post-processing
1064if one wants to remove water less than 1cm deep (for instance).
1065
1066\begin{verbatim}
1067domain.set_name('cairns_' + project.scenario) # Name of SWW file
1068domain.set_datadir('.')                       # Store SWW output here
1069domain.set_minimum_storable_height(0.01)      # Store only depth > 1cm
1070\end{verbatim}
1071
1072\subsection{Initial Conditions}
1073
1074Quantities for \file{runcairns.py} are set
1075using similar methods to those in \file{runup.py}. However,
1076in this case, many of the values are read from the auxiliary file
1077\file{project.py} or, in the case of \code{elevation}, from an
1078auxiliary points file.
1079
1080\subsubsection{Stage}
1081
1082The stage is initially set to 0.0 (i.e.\ Mean Sea Level) by the following statements:
1083
1084\begin{verbatim}
1085tide = 0.0
1086domain.set_quantity('stage', tide)
1087\end{verbatim}
1088
1089It could also take the value of the highest astronomical tide.
1090
1091%For the scenario we are modelling in this case, we use a callable
1092%object \code{tsunami_source}, assigned by means of a function
1093%\function{slide\_tsunami}. This is similar to how we set elevation in
1094%\file{runup.py} using a function -- however, in this case the
1095%function is both more complex and more interesting.
1096
1097%The function returns the water displacement for all \code{x} and
1098%\code{y} in the domain. The water displacement is a double Gaussian
1099%function that depends on the characteristics of the slide (length,
1100%width, thickness, slope, etc), its location (origin) and the depth at that
1101%location. For this example, we choose to apply the slide function
1102%at a specified time into the simulation. {\bf Note, the parameters used
1103%in this example have been deliberately chosen to generate a suitably
1104%large amplitude tsunami which would inundate the Cairns region.}
1105
1106\subsubsection{Friction}
1107
1108We assign the friction exactly as we did for \file{runup.py}:
1109
1110\begin{verbatim}
1111domain.set_quantity('friction', 0.0)
1112\end{verbatim}
1113
1114\subsubsection{Elevation}
1115
1116The elevation is specified by reading data from a file with a name derived from
1117\code{project.demname} with the \code{.pts} extension:
1118
1119\begin{verbatim}
1120domain.set_quantity('elevation',
1121                    filename=project.demname + '.pts',
1122                    use_cache=True,
1123                    verbose=True,
1124                    alpha=0.1)
1125\end{verbatim}
1126
1127The \code{alpha} parameter controls how smooth the elevation surface
1128should be.  See section \ref{class:alpha_shape}, page \pageref{class:alpha_shape}.
1129
1130Setting \code{cache=True} allows \anuga to save the result in order
1131to make subsequent runs faster.
1132
1133Using \code{verbose=True} tells the function to write diagnostics to
1134the screen.
1135
1136\subsection{Boundary Conditions}\index{boundary conditions}
1137
1138Setting boundaries follows a similar pattern to the one used for
1139\file{runup.py}, except that in this case we need to associate a
1140boundary type with each of the
1141boundary tag names introduced when we established the mesh. In place of the four
1142boundary types introduced for \file{runup.py}, we use the reflective
1143boundary for each of the tagged segments defined by \code{create_domain_from_regions}:
1144
1145\begin{verbatim}
1146Bd = Dirichlet_boundary([tide,0,0]) # Mean water level
1147Bs = Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary
1148
1149if project.scenario == 'fixed_wave':
1150    # Huge 50m wave starting after 60 seconds and lasting 1 hour.
1151    Bw = Time_boundary(domain=domain,
1152                       function=lambda t: [(60<t<3660)*50, 0, 0])
1153    domain.set_boundary({'ocean_east': Bw,
1154                         'bottom': Bs,
1155                         'onshore': Bd,
1156                         'top': Bs})
1157
1158if project.scenario == 'slide':
1159    # Boundary conditions for slide scenario
1160    domain.set_boundary({'ocean_east': Bd,
1161                         'bottom': Bd,
1162                         'onshore': Bd,
1163                         'top': Bd})
1164\end{verbatim}
1165
1166Note that we use different boundary conditions depending on the \code{scenario}
1167defined in \file{project.py}.
1168
1169It is not a requirement in \anuga to have this code structure, just an example of
1170how the script can take different actions depending on a variable.
1171
1172\subsection{Evolution}
1173
1174With the basics established, the running of the 'evolve' step is
1175very similar to the corresponding step in \file{runup.py}, except we have different \code{evolve}
1176loops for the two scenarios.
1177
1178For the slide scenario, the simulation is run for an intial 60 seconds, at which time
1179the slide occurs.  We use the function \function{tsunami_source} to adjust \code{stage}
1180values.  We then run the simulation until 5000 seconds with the output stored
1181every ten seconds:
1182
1183\begin{verbatim}
1184if project.scenario == 'slide':
1185    # Initial run without any event
1186    for t in domain.evolve(yieldstep=10, finaltime=60):
1187        print domain.timestepping_statistics()
1188        print domain.boundary_statistics(tags='ocean_east')
1189
1190    # Add slide to water surface
1191    if allclose(t, 60):
1192        domain.add_quantity('stage', tsunami_source)
1193
1194    # Continue propagating wave
1195    for t in domain.evolve(yieldstep=10, finaltime=5000,
1196                           skip_initial_step=True):
1197        print domain.timestepping_statistics()
1198        print domain.boundary_statistics(tags='ocean_east')
1199
1200if project.scenario == 'fixed_wave':
1201    # Save every two mins leading up to wave approaching land
1202    for t in domain.evolve(yieldstep=120, finaltime=5000):
1203        print domain.timestepping_statistics()
1204        print domain.boundary_statistics(tags='ocean_east')
1205
1206    # Save every 30 secs as wave starts inundating ashore
1207    for t in domain.evolve(yieldstep=10, finaltime=10000,
1208                           skip_initial_step=True):
1209        print domain.timestepping_statistics()
1210        print domain.boundary_statistics(tags='ocean_east')
1211\end{verbatim}
1212
1213For the fixed wave scenario, the simulation is run to 10000 seconds,
1214with the first half of the simulation stored at two minute intervals,
1215and the second half of the simulation stored at ten second intervals.
1216This functionality is especially convenient as it allows the detailed
1217parts of the simulation to be viewed at higher time resolution.
1218
1219This also demonstrates the ability of \anuga to dynamically override values.  The
1220\code{method add_quantity()} works like \code{set_quantity()} except that it adds the new
1221surface to what exists already.  In this case it adds the initial shape of the water
1222displacement to the water level.
1223
1224\section{Exploring the Model Output}
1225
1226Now that the scenario has been run, the user can view the output in a number of ways.
1227As described earlier, the user may run \code{animate} to view a three-dimensional representation
1228of the simulation.
1229
1230The user may also be interested in a maximum inundation map. This simply shows the
1231maximum water depth over the domain and is achieved with the function \code{sww2dem}
1232described in Section \ref{sec:basicfileconversions}).
1233\file{ExportResults.py} demonstrates how this function can be used:
1234
1235\verbatiminput{demos/cairns/ExportResults.py}
1236
1237The script generates a maximum water depth ASCII grid at a defined
1238resolution (here 100 m$^2$) which can then be viewed in a GIS environment, for
1239example. The parameters used in the function are defined in \file{project.py}.
1240Figures \ref{fig:maxdepthcairnsslide} and \ref{fig:maxdepthcairnsfixedwave} show
1241the maximum water depth within the defined region for the slide and fixed wave scenario
1242respectively. {\bf Note, these inundation maps have been based on purely hypothetical
1243scenarios and were designed explicitly for demonstration purposes only.}
1244The user could develop a maximum absolute momentum or other expressions which can be
1245derived from the quantities.
1246It must be noted here that depth is more meaningful when the elevation is positive
1247(\code{depth} = \code{stage} $-$ \code{elevation}) as it describes the water height
1248above the available elevation. When the elevation is negative, depth is meauring the
1249water height from the sea floor. With this in mind, maximum inundation maps are
1250typically "clipped" to the coastline. However, the data input here did not contain a
1251coastline.
1252
1253\clearpage
1254
1255\begin{figure}[htp]
1256  \centerline{\includegraphics[scale=0.5]{graphics/slidedepth.jpg}}
1257  \caption{Maximum inundation map for the Cairns slide scenario. \bf Note, this
1258           inundation map has been based on a purely hypothetical scenario which was
1259           designed explictiy for demonstration purposes only.}
1260  \label{fig:maxdepthcairnsslide}
1261\end{figure}
1262
1263\clearpage
1264
1265\begin{figure}[htp]
1266  \centerline{\includegraphics[scale=0.5]{graphics/fixedwavedepth.jpg}}
1267  \caption{Maximum inundation map for the Cairns fixed wave scenario.
1268           \bf Note, this inundation map has been based on a purely hypothetical scenario which was
1269           designed explictiy for demonstration purposes only.}
1270  \label{fig:maxdepthcairnsfixedwave}
1271\end{figure}
1272
1273\clearpage
1274
1275The user may also be interested in interrogating the solution at a particular spatial
1276location to understand the behaviour of the system through time. To do this, the user
1277must first define the locations of interest. A number of locations have been
1278identified for the Cairns scenario, as shown in Figure \ref{fig:cairnsgauges}.
1279
1280\begin{figure}[htp]
1281  \centerline{\includegraphics[scale=0.5]{graphics/cairnsgauges.jpg}}
1282  \caption{Point locations to show time series information for the Cairns scenario.}
1283  \label{fig:cairnsgauges}
1284\end{figure}
1285
1286These locations
1287must be stored in either a .csv or .txt file. The corresponding .csv file for
1288the gauges shown in Figure \ref{fig:cairnsgauges} is \file{gauges.csv}:
1289
1290\verbatiminput{demos/cairns/gauges.csv}
1291
1292Header information has been included to identify the location in terms of eastings and
1293northings, and each gauge is given a name. The elevation column can be zero here.
1294This information is then passed to the function \code{sww2csv_gauges} (shown in
1295\file{GetTimeseries.py} which generates the csv files for each point location. The CSV files
1296can then be used in \code{csv2timeseries_graphs} to create the timeseries plot for each desired
1297quantity. \code{csv2timeseries_graphs} relies on \code{pylab} to be installed which is not part
1298of the standard \code{anuga} release, however it can be downloaded and installed from \code{http://matplotlib.sourceforge.net/}
1299
1300\verbatiminput{demos/cairns/GetTimeseries.py}
1301
1302Here, the time series for the quantities stage, depth and speed will be generated for
1303each gauge defined in the gauge file. As described earlier, depth is more meaningful
1304for onshore gauges, and stage is more appropriate for offshore gauges.
1305
1306As an example output,
1307Figure \ref{fig:reef} shows the time series for the quantity stage for the
1308Elford Reef location for each scenario (the elevation at this location is negative,
1309therefore stage is the more appropriate quantity to plot). Note the large negative stage value when the slide was
1310introduced. This is due to the double gaussian form of the initial surface
1311displacement of the slide. By contrast, the time series for depth is shown for the onshore location of the Cairns
1312Airport in Figure \ref{fig:airportboth}.
1313
1314\begin{figure}[htp]
1315  \centerline{\includegraphics[scale=0.5]{graphics/gaugeElfordReefstage.png}}
1316  \caption{Time series information of the quantity stage for the Elford Reef location for the
1317           fixed wave and slide scenario.}
1318  \label{fig:reef}
1319\end{figure}
1320
1321\begin{figure}[htp]
1322  \centerline{\includegraphics[scale=0.5]{graphics/gaugeCairnsAirportdepth.png}}
1323  \caption{Time series information of the quantity depth for the Cairns Airport
1324           location for the slide and fixed wave scenario.}
1325  \label{fig:airportboth}
1326\end{figure}
1327
1328%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1329
1330\chapter{\anuga Public Interface}
1331\label{ch:interface}
1332
1333This chapter gives an overview of the features of \anuga available
1334to the user at the public interface. These are grouped under the
1335following headings, which correspond to the outline of the examples
1336described in Chapter \ref{ch:getstarted}:
1337\begin{itemize}
1338    \item Establishing the Mesh: Section \ref{sec:establishing the mesh}
1339    \item Initialising the Domain: Section \ref{sec:initialising the domain}
1340%    \item Specifying the Quantities: Section \ref{sec:quantities}
1341    \item Initial Conditions: Section \ref{sec:initial conditions}
1342    \item Boundary Conditions: Section \ref{sec:boundary conditions}
1343    \item Forcing Terms: Section \ref{sec:forcing terms}
1344    \item Evolution: Section \ref{sec:evolution}
1345\end{itemize}
1346
1347The listings are intended merely to give the reader an idea of what
1348each feature is, where to find it and how it can be used -- they do
1349not give full specifications; for these the reader
1350may consult the code. The code for every function or class contains
1351a documentation string, or 'docstring', that specifies the precise
1352syntax for its use. This appears immediately after the line
1353introducing the code, between two sets of triple quotes.
1354
1355Each listing also describes the location of the module in which
1356the code for the feature being described can be found. All modules
1357are in the folder \file{inundation} or one of its subfolders, and the
1358location of each module is described relative to \file{inundation}. Rather
1359than using pathnames, whose syntax depends on the operating system,
1360we use the format adopted for importing the function or class for
1361use in Python code. For example, suppose we wish to specify that the
1362function \function{create\_mesh\_from\_regions} is in a module called
1363\module{mesh\_interface} in a subfolder of \module{inundation} called
1364\code{pmesh}. In Linux or Unix syntax, the pathname of the file
1365containing the function, relative to \file{inundation}, would be:
1366
1367\begin{verbatim}
1368pmesh/mesh_interface.py
1369\end{verbatim}
1370
1371\label{sec:mesh interface}
1372while in Windows syntax it would be:
1373
1374\begin{verbatim}
1375pmesh\mesh_interface.py
1376\end{verbatim}
1377
1378Rather than using either of these forms, in this chapter we specify
1379the location simply as \code{pmesh.mesh_interface}, in keeping with
1380the usage in the Python statement for importing the function,
1381namely:
1382
1383\begin{verbatim}
1384from pmesh.mesh_interface import create_mesh_from_regions
1385\end{verbatim}
1386
1387Each listing details the full set of parameters for the class or
1388function; however, the description is generally limited to the most
1389important parameters and the reader is again referred to the code
1390for more details.
1391
1392The following parameters are common to many functions and classes
1393and are omitted from the descriptions given below:
1394
1395%\begin{tabular}{ll}
1396\begin{tabular}{p{2.0cm} p{14.0cm}}
1397  \emph{use\_cache} & Specifies whether caching is to be used for improved performance.
1398                      See Section \ref{sec:caching} for details on the underlying caching functionality\\
1399  \emph{verbose} & If \code{True}, provides detailed terminal output to the user\\
1400\end{tabular}
1401
1402
1403\section{Mesh Generation}\index{Mesh!generation}
1404\label{sec:establishing the mesh}
1405Before discussing the part of the interface relating to mesh
1406generation, we begin with a description of a simple example of a
1407mesh and use it to describe how mesh data is stored.
1408
1409\label{sec:meshexample} Figure \ref{fig:simplemesh} represents a
1410very simple mesh comprising just 11 points and 10 triangles.
1411
1412\begin{figure}[htp]
1413  \begin{center}
1414    \includegraphics[width=90mm, height=90mm]{triangularmesh.jpg}
1415  \end{center}
1416  \caption{A simple mesh}
1417  \label{fig:simplemesh}
1418\end{figure}
1419
1420\clearpage
1421
1422The variables \code{points}, \code{triangles} and \code{boundary}
1423represent the data displayed in Figure \ref{fig:simplemesh} as
1424follows. The list \code{points} stores the coordinates of the
1425points, and may be displayed schematically as in Table \ref{tab:points}.
1426
1427\begin{table}[htp]
1428  \begin{center}
1429    \begin{tabular}[t]{|c|cc|} \hline
1430      index & \code{x} & \code{y}\\  \hline
1431      0 & 1 & 1\\
1432      1 & 4 & 2\\
1433      2 & 8 & 1\\
1434      3 & 1 & 3\\
1435      4 & 5 & 5\\
1436      5 & 8 & 6\\
1437      6 & 11 & 5\\
1438      7 & 3 & 6\\
1439      8 & 1 & 8\\
1440      9 & 4 & 9\\
1441      10 & 10 & 7\\  \hline
1442    \end{tabular}
1443  \end{center}
1444  \caption{Point coordinates for mesh in Figure \protect \ref{fig:simplemesh}}
1445  \label{tab:points}
1446\end{table}
1447
1448The list \code{triangles} specifies the triangles that make up the
1449mesh. It does this by specifying, for each triangle, the indices
1450(the numbers shown in the first column above) that correspond to the
1451three points at the triangles vertices, taken in an anti-clockwise order
1452around the triangle. Thus, in the example shown in Figure
1453\ref{fig:simplemesh}, the variable \code{triangles} contains the
1454entries shown in Table \ref{tab:triangles}. The starting point is
1455arbitrary so triangle $(0,1,3)$ is considered the same as $(1,3,0)$
1456and $(3,0,1)$.
1457
1458\begin{table}[htp]
1459  \begin{center}
1460    \begin{tabular}{|c|ccc|}
1461      \hline
1462      index & \multicolumn{3}{c|}{\code{points}}\\
1463      \hline
1464      0 & 0 & 1 & 3\\
1465      1 & 1 & 2 & 4\\
1466      2 & 2 & 5 & 4\\
1467      3 & 2 & 6 & 5\\
1468      4 & 4 & 5 & 9\\
1469      5 & 4 & 9 & 7\\
1470      6 & 3 & 4 & 7\\
1471      7 & 7 & 9 & 8\\
1472      8 & 1 & 4 & 3\\
1473      9 & 5 & 10 & 9\\
1474      \hline
1475    \end{tabular}
1476  \end{center}
1477
1478  \caption{Triangles for mesh in Figure \protect \ref{fig:simplemesh}}
1479  \label{tab:triangles}
1480\end{table}
1481
1482Finally, the variable \code{boundary} identifies the boundary
1483triangles and associates a tag with each.
1484
1485% \refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}
1486\label{sec:meshgeneration}
1487
1488\begin{funcdesc}{create_mesh_from_regions}{bounding_polygon,
1489                                           boundary_tags,
1490                                           maximum_triangle_area=None,
1491                                           filename=None,
1492                                           interior_regions=None,
1493                                           interior_holes=None,
1494                                           poly_geo_reference=None,
1495                                           mesh_geo_reference=None,
1496                                           minimum_triangle_angle=28.0,
1497                                           fail_if_polygons_outside=True,
1498                                           use_cache=False,
1499                                           verbose=True}
1500Module: \module{pmesh.mesh\_interface}
1501
1502This function allows a user to initiate the automatic creation of a
1503mesh inside a specified polygon (input \code{bounding_polygon}).
1504Among the parameters that can be set are the \emph{resolution}
1505(maximal area for any triangle in the mesh) and the minimal angle
1506allowable in any triangle. The user can specify a number of internal
1507polygons within each of which the resolution of the mesh can be
1508specified. \code{interior_regions} is a paired list containing the
1509interior polygon and its resolution.  Additionally, the user specifies
1510a list of boundary tags, one for each edge of the bounding polygon.
1511
1512\textbf{WARNING}. Note that the dictionary structure used for the
1513parameter \code{boundary\_tags} is different from that used for the
1514variable \code{boundary} that occurs in the specification of a mesh.
1515In the case of \code{boundary}, the tags are the \emph{values} of
1516the dictionary, whereas in the case of \code{boundary_tags}, the
1517tags are the \emph{keys} and the \emph{value} corresponding to a
1518particular tag is a list of numbers identifying boundary edges
1519labelled with that tag. Because of this, it is theoretically
1520possible to assign the same edge to more than one tag. However, an
1521attempt to do this will cause an error.
1522
1523\textbf{WARNING}. Do not have polygon lines cross or be on-top of each
1524    other. This can result in regions of unspecified resolutions. Do
1525    not have polygon close to each other. This can result in the area
1526    between the polygons having small triangles.  For more control
1527    over the mesh outline use the methods described below.
1528
1529\end{funcdesc}
1530
1531\begin{funcdesc}{create_domain_from_regions}{bounding_polygon,
1532                                             boundary_tags,
1533                                             maximum_triangle_area=None,
1534                                             mesh_filename=None,
1535                                             interior_regions=None,
1536                                             interior_holes=None,
1537                                             poly_geo_reference=None,
1538                                             mesh_geo_reference=None,
1539                                             minimum_triangle_angle=28.0,
1540                                             fail_if_polygons_outside=True,
1541                                             use_cache=False,
1542                                             verbose=True}
1543
1544Module: \module{interface.py}
1545
1546This higher-level function allows a user to create a domain (and associated mesh)
1547inside a specified polygon.
1548
1549\code{bounding_polygon} is a list of points in Eastings and Northings,
1550relative to the zone stated in \code{poly_geo_reference} if specified.
1551Otherwise points are just x, y coordinates with no particular
1552association to any location.
1553
1554\code{boundary_tags} is a dictionary of symbolic tags. For every tag there
1555is a list of indices referring to segments associated with that tag.
1556If a segment is omitted it will be assigned the default tag ''.
1557
1558\code{maximum_triangle_area} is the maximal area per triangle
1559for the bounding polygon, excluding the interior regions.
1560
1561\code{mesh_filename} is the name of the file to contain the generated
1562mesh data.
1563
1564\code{interior_regions} is a list of tuples consisting of (polygon,
1565resolution) for each region to be separately refined. Do not have
1566polygon lines cross or be on-top of each other.  Also do not have
1567polygons close to each other.
1568
1569\code{poly_geo_reference} is the geo_reference of the bounding polygon and
1570the interior polygons.
1571If none, assume absolute.  Please pass one though, since absolute
1572references have a zone.
1573
1574\code{mesh_geo_reference} is the geo_reference of the mesh to be created.
1575If none is given one will be automatically generated.  It will use
1576the lower left hand corner of  bounding_polygon (absolute)
1577as the x and y values for the geo_ref.
1578
1579\code{minimum_triangle_angle} is the minimum angle allowed for each generated triangle.
1580This controls the \emph{slimness} allowed for a triangle.
1581
1582\code{fail_if_polygons_outside} -- if True (the default) an Exception in thrown
1583if interior polygons fall outside the bounding polygon. If False, these
1584will be ignored and execution continues.
1585
1586\textbf{WARNING}. Note that the dictionary structure used for the
1587parameter \code{boundary_tags} is different from that used for the
1588variable \code{boundary} that occurs in the specification of a mesh.
1589In the case of \code{boundary}, the tags are the \emph{values} of
1590the dictionary, whereas in the case of \code{boundary_tags}, the
1591tags are the \emph{keys} and the \emph{value} corresponding to a
1592particular tag is a list of numbers identifying boundary edges
1593labelled with that tag. Because of this, it is theoretically
1594possible to assign the same edge to more than one tag. However, an
1595attempt to do this will cause an error.
1596
1597\textbf{WARNING}. Do not have polygon lines cross or be on-top of each
1598    other. This can result in regions of unspecified resolutions. Do
1599    not have polygon close to each other. This can result in the area
1600    between the polygons having small triangles.  For more control
1601    over the mesh outline use the methods described below.
1602
1603\end{funcdesc}
1604
1605\subsection{Advanced mesh generation}
1606
1607For more control over the creation of the mesh outline, use the
1608methods of the class \class{Mesh}.
1609
1610\begin{classdesc}{Mesh}{userSegments=None,
1611                        userVertices=None,
1612                        holes=None,
1613                        regions=None,
1614                        geo_reference=None}
1615Module: \module{pmesh.mesh}
1616
1617A class used to build a mesh outline and generate a two-dimensional
1618triangular mesh. The mesh outline is used to describe features on the
1619mesh, such as the mesh boundary. Many of this class's methods are used
1620to build a mesh outline, such as \code{add_vertices()} and
1621\code{add_region_from_polygon()}.
1622
1623\code{userSegments} and \code{userVertices} define the outline enclosing the mesh.
1624
1625\code{holes} describes any regions inside the mesh that are not to be included in the mesh.
1626
1627\code{geo_reference} defines the geo_reference to which all point information is relative.
1628If \code{geo_reference} is \code{None} then the default geo_reference is used.
1629\end{classdesc}
1630
1631\subsubsection{Key Methods of Class Mesh}
1632
1633\begin{methoddesc}{\emph{<mesh>}.add_hole}{x, y, geo_reference=None}
1634Module: \module{pmesh.mesh}
1635
1636This method adds a hole to the mesh outline.
1637
1638\code{x} and \code{y} define a point on the already defined hole boundary.
1639
1640If \code{geo_reference} is not supplied the points are assumed to be absolute.
1641\end{methoddesc}
1642
1643\begin{methoddesc}{\emph{<mesh>}.add_hole_from_polygon}{polygon,
1644                                                        segment_tags=None,
1645                                                        geo_reference=None}
1646Module: \module{pmesh.mesh}
1647
1648This method is used to add a 'hole' within a region -- that is, to
1649define a interior region where the triangular mesh will not be
1650generated -- to a \class{Mesh} instance. The region boundary is described by
1651the polygon passed in.  Additionally, the user specifies a list of
1652boundary tags, one for each edge of the bounding polygon.
1653
1654\code{polygon} is the polygon that defines the hole to be added to the \code{<mesh>}.
1655
1656\code{segment_tags} -- ??
1657
1658If \code{geo_reference} is \code{None} then the default \code{geo_reference} is used.
1659\end{methoddesc}
1660
1661\begin{methoddesc}{\emph{<mesh>}.add_points_and_segments}{points,
1662                                                          segments=None,
1663                                                          segment_tags=None}
1664Module: \module{pmesh.mesh}
1665
1666This adds points and segments connecting the points to a mesh.
1667
1668\code{points} is a list of points.
1669
1670\code{segments} is a list of segments.  Each segment is defined by the start and end
1671of the line by its point index, e.g.\ use \code{segments = [[0,1],[1,2]]} to make a
1672polyline between points 0, 1 and 2.
1673
1674\code{segment_tags} may be used to optionally define a tag for each segment.
1675\end{methoddesc}
1676
1677\begin{methoddesc}{\emph{<mesh>}.add_region}{x,y, geo_reference=None, tag=None}
1678Module: \module{pmesh.mesh}
1679
1680This method adds a region to a mesh outline.
1681
1682\code{x} and \code{y} define a point on the already-defined region that is to
1683be added to the mesh.
1684
1685If \code{geo_reference} is not supplied the points data is assumed to be absolute.
1686
1687\code{tag} -- ??
1688
1689A region instance is returned.  This can be used to set the resolution of the added region.
1690\end{methoddesc}
1691
1692\begin{methoddesc}{\emph{<mesh>}.add_region_from_polygon}{polygon,
1693                                                          segment_tags=None,
1694                                                          max_triangle_area=None,
1695                                                          geo_reference=None,
1696                                                          region_tag=None}
1697Module: \module{pmesh.mesh}
1698
1699This method adds a region to a
1700\class{Mesh} instance.  Regions are commonly used to describe an area
1701with an increased density of triangles by setting \code{max_triangle_area}.
1702
1703\code{polygon} describes the region boundary to add to the \code{<mesh>}.
1704
1705\code{segment_tags} specifies a list of segment tags, one for each edge of the
1706bounding polygon.
1707
1708If \code{geo_reference} is not supplied the points data is assumed to be absolute.
1709
1710\code{region_tag} sets the region tag.
1711\end{methoddesc}
1712
1713\begin{methoddesc}{\emph{<mesh>}.add_vertices}{point_data}
1714Module: \module{pmesh.mesh}
1715
1716Add user vertices to a mesh.
1717
1718\code{point_data} is the list of point data, and can be a list of (x,y) values,
1719a numeric array or a geospatial_data instance.
1720\end{methoddesc}
1721
1722\begin{methoddesc}{\emph{<mesh>}.auto_segment}{alpha=None,
1723                                               raw_boundary=True,
1724                                               remove_holes=False,
1725                                               smooth_indents=False,
1726                                               expand_pinch=False}
1727Module: \module{pmesh.mesh}
1728
1729Add segments between some of the user vertices to give the vertices an
1730outline.  The outline is an alpha shape. This method is
1731useful since a set of user vertices need to be outlined by segments
1732before generate_mesh is called.
1733
1734\code{alpha} determines the $smoothness$ of the alpha shape.
1735
1736\code{raw_boundary}, if \code{True} instructs the function to return the raw
1737boundary, i.e.\ the regular edges of the alpha shape.
1738
1739\code{remove_holes}, if \code{True} enables a filter to remove small holes
1740(small is defined by  boundary_points_fraction).
1741
1742\code{smooth_indents}, if \code{True} removes sharp triangular indents
1743in the boundary.
1744
1745\code{expand_pinch}, if \code{True} tests for pinch-off and corrects
1746(i.e.\ a boundary vertex with more than two edges).
1747\end{methoddesc}
1748
1749\begin{methoddesc}{\emph{<mesh>}.export_mesh_file}{ofile}
1750Module: \module{pmesh.mesh}
1751
1752This method is used to save a mesh to a file.
1753
1754\code{ofile} is the name of the mesh file to be written, including the extension.
1755Use the extension \code{.msh} for the file to be in NetCDF format and
1756\code{.tsh} for the file to be ASCII format.
1757\end{methoddesc}
1758
1759\begin{methoddesc}{\emph{<mesh>}.generate_mesh}{maximum_triangle_area="",
1760                                                minimum_triangle_angle=28.0,
1761                                                verbose=True}
1762Module: \module{pmesh.mesh}
1763
1764This method is used to generate the triangular mesh.
1765
1766\code{maximum_triangle_area} sets the maximum area of any triangle in the mesh.
1767
1768\code{minimum_triangle_angle} sets the minimum area of any triangle in the mesh.
1769
1770These two parameters can be used to control the triangle density.
1771\end{methoddesc}
1772
1773\begin{methoddesc}{\emph{<mesh>}.import_ungenerate_file}{ofile,
1774                                                         tag=None,
1775                                                         region_tag=None}
1776Module: \module{pmesh.mesh}
1777
1778This method is used to import a polygon file in the ungenerate format,
1779which is used by arcGIS. The polygons from the file are converted to
1780vertices and segments.
1781
1782\code{ofile} is the name of the polygon file.
1783
1784\code{tag} is the tag given to all the polygon's segments.
1785If \code{tag} is not supplied then the segment will not effect the water
1786flow, it will only effect the mesh generation.
1787
1788\code{region_tag} is the tag given to all the polygon's segments.  If
1789it is a string the tag will be assigned to all regions.  If it
1790is a list the first value in the list will be applied to the first
1791polygon etc.
1792
1793This function can be used to import building footprints.
1794\end{methoddesc}
1795
1796
1797\section{Initialising the Domain}\index{Initialising the Domain}
1798\label{sec:initialising the domain}
1799
1800\begin{classdesc}{Domain}{source=None,
1801                          triangles=None,
1802                          boundary=None,
1803                          conserved_quantities=None,
1804                          other_quantities=None,
1805                          tagged_elements=None,
1806                          geo_reference=None,
1807                          use_inscribed_circle=False,
1808                          mesh_filename=None,
1809                          use_cache=False,
1810                          verbose=False,
1811                          full_send_dict=None,
1812                          ghost_recv_dict=None,
1813                          processor=0,
1814                          numproc=1,
1815                          number_of_full_nodes=None,
1816                          number_of_full_triangles=None}
1817Module: \refmodule{abstract_2d_finite_volumes.domain}
1818
1819This class is used to create an instance of a structure used to
1820store and manipulate data associated with a mesh. The mesh is
1821specified either by assigning the name of a mesh file to
1822\code{source} or by specifying the points, triangle and boundary of the
1823mesh.
1824\end{classdesc}
1825
1826\subsection{Key Methods of Domain}
1827
1828\begin{methoddesc}{\emph{<domain>}.set_name}{name}
1829Module: \refmodule{abstract_2d_finite_volumes.domain},
1830page \pageref{mod:domain}
1831
1832\code{name} is used to name the domain.  The \code{name} is also used to identify the output SWW file.
1833If no name is assigned to a domain, the assumed name is \code{'domain'}.
1834\end{methoddesc}
1835
1836\begin{methoddesc}{\emph{<domain>}.get_name}{}
1837Module: \module{abstract_2d_finite_volumes.domain}
1838
1839Returns the name assigned to the domain by \code{set_name()}. If no name has been
1840assigned, returns \code{'domain'}.
1841\end{methoddesc}
1842
1843\begin{methoddesc}{\emph{<domain>}.set_datadir}{path}
1844Module: \module{abstract_2d_finite_volumes.domain}
1845
1846\code{path} specifies the path to the directory used to store SWW files.
1847
1848Before this method is used to set the SWW directory path, the assumed directory
1849path is \code{default_datadir} specified in \code{config.py}.
1850
1851Since different operating systems use different formats for specifying pathnames
1852it is necessary to specify path separators using the Python code \code{os.sep} rather than
1853the operating-specific ones such as '$\slash$' or '$\backslash$'.
1854For this to work you will need to include the statement \code{import os}
1855in your code, before the first use of \code{set_datadir()}.
1856
1857For example, to set the data directory to a subdirectory
1858\code{data} of the directory \code{project}, you could use
1859the statements:
1860
1861\begin{verbatim}
1862import os
1863domain.set_datadir{'project' + os.sep + 'data'}
1864\end{verbatim}
1865\end{methoddesc}
1866
1867\begin{methoddesc}{\emph{<domain>}.get_datadir}{}
1868Module: \module{abstract_2d_finite_volumes.domain}
1869
1870Returns the path to the directory where SWW files will be stored.
1871
1872If the path has not previously been set with \code{set_datadir()} this method
1873will return the value \code{default_datadir} specified in \code{config.py}.
1874\end{methoddesc}
1875
1876\begin{methoddesc}{\emph{<domain>}.set_minimum_allowed_height}{minimum_allowed_height}
1877Module: \module{shallow_water.shallow_water_domain}
1878
1879Set the minimum depth (in metres) that will be recognised in
1880the numerical scheme (including limiters and flux computations)
1881
1882\code{minimum_allowed_height} is the new minimum allowed height value.
1883
1884Default value is $10^{-3}$ metre, but by setting this to a greater value,
1885e.g.\ for large scale simulations, the computation time can be
1886significantly reduced.
1887\end{methoddesc}
1888
1889\begin{methoddesc}{\emph{<domain>}.set_minimum_storable_height}{minimum_storable_height}
1890Module: \module{shallow_water.shallow_water_domain}
1891
1892Sets the minimum depth that will be recognised when writing
1893to an SWW file. This is useful for removing thin water layers
1894that seems to be caused by friction creep.
1895
1896\code{minimum_storable_height} is the new minimum storable height value.
1897\end{methoddesc}
1898
1899\begin{methoddesc}{\emph{<domain>}.set_maximum_allowed_speed}{maximum_allowed_speed}
1900Module: \module{shallow_water.shallow_water_domain}
1901
1902Set the maximum particle speed that is allowed in water
1903shallower than \code{minimum_allowed_height}. This is useful for
1904controlling speeds in very thin layers of water and at the same time
1905allow some movement avoiding pooling of water.
1906
1907\code{maximum_allowed_speed} sets the maximum allowed speed value.
1908\end{methoddesc}
1909
1910\begin{methoddesc}{\emph{<domain>}.set_time}{time=0.0}
1911Module: \module{abstract_2d_finite_volumes.domain}
1912
1913\code{time} sets the initial time, in seconds, for the simulation. The
1914default is 0.0.
1915\end{methoddesc}
1916
1917\begin{methoddesc}{\emph{<domain>}.set_default_order}{n}
1918Module: \module{abstract_2d_finite_volumes.domain}
1919
1920Sets the default (spatial) order to the value specified by
1921\code{n}, which must be either 1 or 2. (Assigning any other value
1922to \code{n} will cause an error.)
1923\end{methoddesc}
1924
1925\begin{methoddesc}{\emph{<domain>}.set_store_vertices_uniquely}{flag, reduction=None}
1926Module: \module{shallow_water.shallow_water_domain}
1927
1928Decide whether vertex values should be stored uniquely as
1929computed in the model or whether they should be reduced to one
1930value per vertex using averaging.
1931
1932\code{flag} may be \code{True} (meaning allow surface to be discontinuous) or \code{False} (meaning smooth vertex values).
1933
1934\code{reduction} defines the smoothing operation if \code{flag} is \code{False}.  If not
1935supplied, \code{reduction} is assumed to be \code{mean}.
1936
1937Triangles stored in the SWW file can be discontinuous reflecting
1938the internal representation of the finite-volume scheme
1939(this is a feature allowing for arbitrary steepness of the water surface gradient
1940as well as the momentum gradients).
1941However, for visual purposes and also for use with \code{Field_boundary}
1942(and \code{File_boundary}), it is often desirable to store triangles
1943with values at each vertex point as the average of the potentially
1944discontinuous numbers found at vertices of different triangles sharing the
1945same vertex location.
1946
1947Storing one way or the other is controlled in \anuga through the method
1948\code{<domain>.store_vertices_uniquely()}. Options are:
1949\begin{itemize}
1950  \item \code{<domain>.store_vertices_uniquely(True)}: Allow discontinuities in the SWW file
1951  \item \code{<domain>.store_vertices_uniquely(False)}: (Default).
1952  Average values
1953  to ensure continuity in SWW file. The latter also makes for smaller
1954  SWW files.
1955\end{itemize}
1956
1957Note that when model data in the SWW file are averaged (i.e.\ not stored uniquely),
1958then there will most likely be a small discrepancy between values extracted from the SWW
1959file and the same data stored in the model domain. This must be borne in mind when comparing
1960data from the SWW files with that of the model internally.
1961\end{methoddesc}
1962
1963
1964\begin{methoddesc}{\emph{<domain>}.set_quantities_to_be_stored}{quantity_dictionary}
1965Module: \module{shallow_water.shallow_water_domain}
1966
1967Selects quantities that is to be stored in the sww files.
1968The argument can be None, in which case nothing is stored.
1969
1970Otherwise, the argument must be a dictionary where the keys are names of quantities
1971already defined within ANUGA and the values are either 1 or 2. If the value is 1, the quantity
1972will be stored once at the beginning of the simulation, if the value is 2 it will be stored
1973at each timestep. The ANUGA default is equivalent to the call
1974\begin{verbatim} 
1975domain.set_quantities_to_be_stored({'elevation': 1,
1976                                    'stage': 2,
1977                                    'xmomentum': 2,
1978                                    'ymomentum': 2})
1979\end{verbatim} 
1980\end{methoddesc}
1981
1982
1983% Structural methods
1984\begin{methoddesc}{\emph{<domain>}.get_nodes}{absolute=False}
1985Module: \module{abstract_2d_finite_volumes.domain}
1986
1987Return x,y coordinates of all nodes in the domain mesh.  The nodes are ordered
1988in an \code{Nx2} array where N is the number of nodes.  This is the same format
1989they were provided in the constructor i.e.\ without any duplication.
1990
1991\code{absolute} is a boolean which determines whether coordinates
1992are to be made absolute by taking georeference into account.
1993Default is \code{False} as many parts of \anuga expect relative coordinates.
1994\end{methoddesc}
1995
1996\begin{methoddesc}{\emph{<domain>}.get_vertex_coordinates}{absolute=False}
1997Module: \module{abstract_2d_finite_volumes.domain}
1998
1999\label{pg:get vertex coordinates}
2000Return vertex coordinates for all triangles as a \code{3*Mx2} array
2001where the jth vertex of the ith triangle is located in row 3*i+j and
2002M is the number of triangles in the mesh.
2003
2004\code{absolute} is a boolean which determines whether coordinates
2005are to be made absolute by taking georeference into account.
2006Default is \code{False} as many parts of \anuga expect relative coordinates.
2007\end{methoddesc}
2008
2009\begin{methoddesc}{\emph{<domain>}.get_centroid_coordinates}{absolute=False}
2010Module: \module{abstract_2d_finite_volumes.domain}
2011
2012Return centroid coordinates for all triangles as an \code{Mx2} array.
2013   
2014\code{absolute} is a boolean which determines whether coordinates
2015are to be made absolute by taking georeference into account.
2016Default is \code{False} as many parts of \anuga expect relative coordinates.
2017\end{methoddesc}
2018
2019\begin{methoddesc}{\emph{<domain>}.get_triangles}{indices=None}
2020Module: \module{abstract_2d_finite_volumes.domain}
2021
2022Return an \code{Mx3} integer array where M is the number of triangles.
2023Each row corresponds to one triangle and the three entries are
2024indices into the mesh nodes which can be obtained using the method
2025\code{get_nodes()}.
2026
2027\code{indices}, if specified, is the set of triangle \code{id}s of interest.
2028\end{methoddesc}
2029
2030\begin{methoddesc}{\emph{<domain>}.get_disconnected_triangles}{}
2031Module: \module{abstract_2d_finite_volumes.domain}
2032
2033Get the domain mesh based on nodes obtained from \code{get_vertex_coordinates()}.
2034
2035Returns an \code{Mx3} array of integers where each row corresponds to
2036a triangle. A triangle is a triplet of indices into
2037point coordinates obtained from \code{get_vertex_coordinates()} and each
2038index appears only once.
2039
2040This provides a mesh where no triangles share nodes
2041(hence the name disconnected triangles) and different
2042nodes may have the same coordinates.
2043
2044This version of the mesh is useful for storing meshes with
2045discontinuities at each node and is e.g.\ used for storing
2046data in SWW files.
2047
2048The triangles created will have the format:
2049
2050\begin{verbatim}
2051[[0,1,2],
2052 [3,4,5],
2053 [6,7,8],
2054 ...
2055 [3*M-3 3*M-2 3*M-1]]
2056\end{verbatim}
2057\end{methoddesc}
2058
2059
2060\section{Initial Conditions}\index{Initial Conditions}
2061\label{sec:initial conditions}
2062In standard usage of partial differential equations, initial conditions
2063refers to the values associated to the system variables (the conserved
2064quantities here) for \code{time = 0}. In setting up a scenario script
2065as described in Sections \ref{sec:simpleexample} and \ref{sec:realdataexample},
2066\code{set_quantity} is used to define the initial conditions of variables
2067other than the conserved quantities, such as friction. Here, we use the terminology
2068of initial conditions to refer to initial values for variables which need
2069prescription to solve the shallow water wave equation. Further, it must be noted
2070that \code{set_quantity()} does not necessarily have to be used in the initial
2071condition setting; it can be used at any time throughout the simulation.
2072
2073\begin{methoddesc}{\emph{<domain>}.set_quantity}{numeric=None,
2074                                                 quantity=None,
2075                                                 function=None,
2076                                                 geospatial_data=None,
2077                                                 filename=None,
2078                                                 attribute_name=None,
2079                                                 alpha=None,
2080                                                 location='vertices',
2081                                                 polygon=None,
2082                                                 indices=None,
2083                                                 smooth=False,
2084                                                 verbose=False,
2085                                                 use_cache=False}
2086Module: \module{abstract_2d_finite_volumes.domain} \\
2087(This method passes off to \module{abstract_2d_finite_volumes.quantity.set_values()})
2088
2089This function is used to assign values to individual quantities for a
2090domain. It is very flexible and can be used with many data types: a
2091statement of the form \code{\emph{<domain>}.set_quantity(name, x)} can be used
2092to define a quantity having the name \code{name}, where the other
2093argument \code{x} can be any of the following:
2094
2095\begin{itemize}
2096  \item a number, in which case all vertices in the mesh gets that for
2097        the quantity in question
2098  \item a list of numbers or a numeric array ordered the same way as the mesh vertices
2099  \item a function (e.g.\ see the samples introduced in Chapter 2)
2100  \item an expression composed of other quantities and numbers, arrays, lists (for
2101        example, a linear combination of quantities, such as
2102        \code{\emph{<domain>}.set_quantity('stage','elevation'+x))}
2103  \item the name of a file from which the data can be read. In this case, the optional
2104        argument \code{attribute_name} will select which attribute to use from the file. If left out,
2105        \code{set_quantity()} will pick one. This is useful in cases where there is only one attribute
2106  \item a geospatial dataset (See Section \ref{sec:geospatial}).
2107        Optional argument \code{attribute_name} applies here as with files
2108\end{itemize}
2109
2110Exactly one of the arguments \code{numeric}, \code{quantity}, \code{function},
2111\code{geospatial_data} and \code{filename} must be present.
2112
2113\code{set_quantity()} will look at the type of the \code{numeric} and
2114determine what action to take.
2115
2116Values can also be set using the appropriate keyword arguments.
2117If \code{x} is a function, for example, \code{domain.set_quantity(name, x)}, \code{domain.set_quantity(name, numeric=x)},
2118and \code{domain.set_quantity(name, function=x)} are all equivalent.
2119
2120Other optional arguments are:
2121\begin{itemize}
2122  \item \code{indices} which is a list of ids of triangles to which \code{set_quantity()}
2123        should apply its assignment of values.
2124  \item \code{location} determines which part of the triangles to assign to.
2125        Options are 'vertices' (the default), 'edges', 'unique vertices', and 'centroids'.
2126        If 'vertices' is used, edge and centroid values are automatically computed as the
2127        appropriate averages. This option ensures continuity of the surface.
2128        If, on the other hand, 'centroids' is used, vertex and edge values will be set to the
2129        same value effectively creating a piecewise constant surface with possible
2130        discontinuities at the edges.
2131\end{itemize}
2132
2133\anuga provides a number of predefined initial conditions to be used
2134with \code{set_quantity()}. See for example callable object \code{slump_tsunami} below.
2135\end{methoddesc}
2136
2137\begin{methoddesc}{\emph{<domain>}.add_quantity}{numeric=None,
2138                                                 quantity=None,
2139                                                 function=None,
2140                                                 geospatial_data=None,
2141                                                 filename=None,
2142                                                 attribute_name=None,
2143                                                 alpha=None,
2144                                                 location='vertices',
2145                                                 polygon=None,
2146                                                 indices=None,
2147                                                 smooth=False,
2148                                                 verbose=False,
2149                                                 use_cache=False}
2150Module: \module{abstract_2d_finite_volumes.domain} \\
2151(passes off to \module{abstract_2d_finite_volumes.domain.set_quantity()})
2152
2153\label{add quantity}
2154This function is used to \emph{add} values to individual quantities for a
2155domain. It has the same syntax as \code{\emph{<domain>}.set_quantity(name, x)}.
2156
2157A typical use of this function is to add structures to an existing elevation model:
2158
2159\begin{verbatim} 
2160# Create digital elevation model from points file
2161domain.set_quantity('elevation', filename='elevation_file.pts, verbose=True)
2162
2163# Add buildings from file
2164building_polygons, building_heights = csv2building_polygons(building_file)
2165
2166B = []
2167for key in building_polygons:
2168    poly = building_polygons[key]
2169    elev = building_heights[key]
2170    B.append((poly, elev))
2171
2172domain.add_quantity('elevation', Polygon_function(B, default=0.0))
2173\end{verbatim}
2174\end{methoddesc}
2175
2176\begin{methoddesc}{\emph{<domain>}.set_region}{tag, quantity, X, location='vertices'}
2177Module: \module{abstract_2d_finite_volumes.domain} \\
2178(see also \module{abstract_2d_finite_volumes.quantity.set_values})
2179
2180This function is used to assign values to individual quantities given
2181a regional tag.   It is similar to \code{set_quantity()}.
2182
2183For example, if in the mesh-generator a regional tag of 'ditch' was
2184used, \code{set_region()} can be used to set elevation of this region to
2185-10m. \code{X} is the constant or function to be applied to the \code{quantity},
2186over the tagged region.  \code{location} describes how the values will be
2187applied.  Options are 'vertices' (the default), 'edges', 'unique
2188vertices', and 'centroids'.
2189
2190This method can also be called with a list of region objects.  This is
2191useful for adding quantities in regions, and having one quantity
2192value based on another quantity. See  \module{abstract_2d_finite_volumes.region} for
2193more details.
2194\end{methoddesc}
2195
2196\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
2197                                radius=None, dphi=0.48, x0=0.0, y0=0.0, alpha=0.0,
2198                                gravity=9.8, gamma=1.85,
2199                                massco=1, dragco=1, frictionco=0,
2200                                dx=None, kappa=3.0, kappad=1.0, zsmall=0.01, scale=None,
2201                                domain=None,
2202                                verbose=False}
2203Module: \module{shallow\_water.smf}
2204
2205This function returns a callable object representing an initial water
2206displacement generated by a submarine sediment failure. These failures can take the form of
2207a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead.
2208
2209\code{length} is the length of the slide or slump.
2210
2211\code{depth} is the water depth to the centre of the sediment mass.
2212
2213\code{slope} is the bathymetric slope.
2214
2215Other slump or slide parameters can be included if they are known.
2216\end{funcdesc}
2217
2218\begin{funcdesc}{\emph{<callable_object>} = file_function}{filename,
2219                                domain=None,
2220                                quantities=None,
2221                                interpolation_points=None,
2222                                time_thinning=1,
2223                                time_limit=None,
2224                                verbose=False,
2225                                use_cache=False,
2226                                boundary_polygon=None}
2227Module: \module{abstract_2d_finite_volumes.util}
2228
2229Reads the time history of spatial data for specified interpolation points from
2230a NetCDF file and returns a callable object.  Values returned from the \code{\emph{<callable_object>}}
2231are interpolated values based on the input file using the underlying \code{interpolation_function}.
2232
2233\code{filename} is the name of the input file. 
2234This would be either an SWW, TMS or STS file.
2235
2236\code{quantities} is either the name of a single quantity to be
2237interpolated or a list of such quantity names. In the second case, the resulting
2238function will return a tuple of values -- one for each quantity.
2239If the NetCDF file uses names other than 'stage', 'xmomentum', and 'ymomentum'
2240the name(s) appearing in the file must be explicitly stated using the
2241quantities keyword. This is for example be the case if a 'tms' file is used
2242to specify time dependent precipitation. In this case the keyword might be called 'rainfall' both in the call to file\_function and in the 'tms' file.
2243
2244\code{interpolation_points} is a list of absolute coordinates or a
2245geospatial object for points at which values are sought.
2246
2247\code{boundary_polygon} is a list of coordinates specifying the vertices of the boundary.
2248This must be the same polygon as used when calling \code{create_mesh_from_regions()}.
2249This argument can only be used when reading boundary data from an STS format file.
2250
2251The model time stored within the file function can be accessed using
2252the method \code{\emph{<callable_object>}.get_time()}
2253
2254The underlying algorithm used is as follows:\\
2255Given a time series (i.e.\ a series of values associated with
2256different times), whose values are either just numbers, a set of
2257numbers defined at the vertices of a triangular mesh (such as those
2258stored in SWW files) or a set of
2259numbers defined at a number of points on the boundary (such as those
2260stored in STS files), \code{Interpolation_function()} is used to
2261create a callable object that interpolates a value for an arbitrary
2262time \code{t} within the model limits and possibly a point \code{(x, y)}
2263within a mesh region.
2264
2265The actual time series at which data is available is specified by
2266means of an array \code{time} of monotonically increasing times. The
2267quantities containing the values to be interpolated are specified in
2268an array -- or dictionary of arrays (used in conjunction with the
2269optional argument \code{quantity_names}) -- called
2270\code{quantities}. The optional arguments \code{vertex_coordinates}
2271and \code{triangles} represent the spatial mesh associated with the
2272quantity arrays. If omitted the function must be created using an STS file
2273or a TMS file.
2274
2275Since, in practice, values need to be computed at specified points,
2276the syntax allows the user to specify, once and for all, a list
2277\code{interpolation_points} of points at which values are required.
2278In this case, the function may be called using the form \code{\emph{<callable_object>}(t, id)},
2279where \code{id} is an index for the list \code{interpolation_points}.
2280\end{funcdesc}
2281
2282\begin{classdesc}{\emph{<callable_object>} = Interpolation_function}{time,
2283                                          quantities,
2284                                          quantity_names=None,
2285                                          vertex_coordinates=None,
2286                                          triangles=None,
2287                                          interpolation_points=None,
2288                                          time_thinning=1,
2289                                          verbose=False,
2290                                          gauge_neighbour_id=None}
2291Module: \module{fit_interpolate.interpolate}
2292
2293Given a time series (i.e.\ a series of values associated with
2294different times) whose values are either just numbers or a set of
2295numbers defined at the vertices of a triangular mesh (such as those
2296stored in SWW files), \code{Interpolation_function} is used to
2297create a callable object that interpolates a value for an arbitrary
2298time \code{t} within the model limits and possibly a point \code{(x, y)}
2299within a mesh region.
2300
2301The actual time series at which data is available is specified by
2302means of an array \code{time} of monotonically increasing times. The
2303quantities containing the values to be interpolated are specified in
2304an array -- or dictionary of arrays (used in conjunction with the
2305optional argument \code{quantity\_names}) -- called
2306\code{quantities}. The optional arguments \code{vertex_coordinates}
2307and \code{triangles} represent the spatial mesh associated with the
2308quantity arrays. If omitted the function created by
2309\code{Interpolation_function} will be a function of \code{t} only.
2310
2311Since, in practice, values need to be computed at specified points,
2312the syntax allows the user to specify, once and for all, a list
2313\code{interpolation_points} of points at which values are required.
2314In this case, the function may be called using the form \code{f(t, id)},
2315where \code{id} is an index for the list \code{interpolation_points}.
2316\end{classdesc}
2317
2318
2319\section{Boundary Conditions}\index{boundary conditions}
2320\label{sec:boundary conditions}
2321
2322\anuga provides a large number of predefined boundary conditions,
2323represented by objects such as \code{Reflective_boundary(domain)} and
2324\code{Dirichlet_boundary([0.2, 0.0, 0.0])}, described in the examples
2325in Chapter 2. Alternatively, you may prefer to ''roll your own'',
2326following the method explained in Section \ref{sec:roll your own}.
2327
2328These boundary objects may be used with the function \code{set_boundary} described below
2329to assign boundary conditions according to the tags used to label boundary segments.
2330
2331\begin{methoddesc}{\emph{<domain>}.set_boundary}{boundary_map}
2332Module: \module{abstract_2d_finite_volumes.domain}
2333
2334This function allows you to assign a boundary object (corresponding to a
2335pre-defined or user-specified boundary condition) to every boundary segment that
2336has been assigned a particular tag.
2337
2338\code{boundary_map} is a dictionary of boundary objects keyed by symbolic tags.
2339\end{methoddesc}
2340
2341\begin{methoddesc} {\emph{<domain>}.get_boundary_tags}{}
2342Module: \module{abstract\_2d\_finite\_volumes.domain}
2343
2344Returns a list of the available boundary tags.
2345\end{methoddesc}
2346
2347\subsection{Predefined boundary conditions}
2348
2349\begin{classdesc}{Reflective_boundary}{domain=None}
2350Module: \module{shallow_water}
2351
2352Reflective boundary returns same conserved quantities as those present in
2353the neighbouring volume but reflected.
2354
2355This class is specific to the shallow water equation as it works with the
2356momentum quantities assumed to be the second and third conserved quantities.
2357\end{classdesc}
2358
2359\begin{classdesc}{Transmissive_boundary}{domain=None}
2360  \label{pg: transmissive boundary}
2361Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
2362
2363A transmissive boundary returns the same conserved quantities as
2364those present in the neighbouring volume.
2365
2366The underlying domain must be specified when the boundary is instantiated.
2367\end{classdesc}
2368
2369\begin{classdesc}{Dirichlet_boundary}{conserved_quantities=None}
2370Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
2371
2372A Dirichlet boundary returns constant values for each of conserved
2373quantities. In the example of \code{Dirichlet_boundary([0.2, 0.0, 0.0])},
2374the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and
2375\code{ymomentum} at the boundary are set to 0.0. The list must contain
2376a value for each conserved quantity.
2377\end{classdesc}
2378
2379\begin{classdesc}{Time_boundary}{domain=None,
2380                                 function=None,
2381                                 default_boundary=None,
2382                                 verbose=False}
2383Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
2384
2385A time-dependent boundary returns values for the conserved
2386quantities as a function of time \code{function(t)}. The user must specify
2387the domain to get access to the model time.
2388
2389Optional argument \code{default_boundary} can be used to specify another boundary object
2390to be used in case model time exceeds the time available in the file used by \code{File_boundary}.
2391The \code{default_boundary} could be a simple Dirichlet condition or
2392even another \code{Time_boundary} typically using data pertaining to another time interval.
2393\end{classdesc}
2394
2395\begin{classdesc}{File_boundary}{filename,
2396                                 domain,
2397                                 time_thinning=1,
2398                                 time_limit=None,
2399                                 boundary_polygon=None,
2400                                 default_boundary=None,
2401                                 use_cache=False,
2402                                 verbose=False}
2403Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
2404
2405This method may be used if the user wishes to apply a SWW file, STS file or
2406a time series file (TMS) to a boundary segment or segments.
2407The boundary values are obtained from a file and interpolated to the
2408appropriate segments for each conserved quantity.
2409
2410Optional argument \code{default_boundary} can be used to specify another boundary object
2411to be used in case model time exceeds the time available in the file used by \code{File_boundary}.
2412The \code{default_boundary} could be a simple Dirichlet condition or
2413even another \code{File_boundary} typically using data pertaining to another time interval.
2414\end{classdesc}
2415
2416\begin{classdesc}{Field_boundary}{filename,
2417                                  domain,
2418                                  mean_stage=0.0,
2419                                  time_thinning=1,
2420                                  time_limit=None,
2421                                  boundary_polygon=None,
2422                                  default_boundary=None,
2423                                  use_cache=False,
2424                                  verbose=False}
2425Module: \module{shallow_water.shallow_water_domain}
2426
2427This method works in the same way as \code{File_boundary} except that it
2428allows for the value of stage to be offset by a constant specified in the
2429keyword argument \code{mean_stage}.
2430
2431This functionality allows for models to be run for a range of tides using
2432the same boundary information (from STS, SWW or TMS files). The tidal value
2433for each run would then be specified in the keyword argument \code{mean_stage}.
2434If \code{mean_stage} = 0.0, \code{Field_boundary} and \code{File_boundary}
2435behave identically.
2436
2437Note that if the optional argument \code{default_boundary} is specified
2438its stage value will be adjusted by \code{mean_stage} just like the values
2439obtained from the file.
2440
2441See \code{File_boundary} for further details.
2442\end{classdesc}
2443
2444\begin{classdesc}{Transmissive_momentum_set_stage_boundary}{domain=None,
2445                                                            function=None}
2446Module: \module{shallow_water.shallow_water_domain}
2447\label{pg: transmissive momentum set stage boundary}
2448
2449This boundary returns the same momentum conserved quantities as
2450those present in its neighbour volume but sets stage as in a \code{Time_boundary}.
2451The underlying domain must be specified when boundary is instantiated.
2452
2453This type of boundary is useful when stage is known at the boundary as a
2454function of time, but momenta (or speeds) aren't.
2455
2456This class is specific to the shallow water equation as it works with the
2457momentum quantities assumed to be the second and third conserved quantities.
2458
2459In some circumstances, this boundary condition may cause numerical instabilities for similar
2460reasons as what has been observed with the simple fully transmissive boundary condition
2461(see Page \pageref{pg: transmissive boundary}).
2462This could for example be the case if a planar wave is reflected out through this boundary.
2463\end{classdesc}
2464
2465\begin{classdesc}{Transmissive_stage_zero_momentum_boundary}{domain=None}
2466Module: \module{shallow_water}
2467\label{pg: transmissive stage zero momentum boundary}
2468
2469This boundary returns same stage conserved quantities as
2470those present in its neighbour volume but sets momentum to zero.
2471The underlying domain must be specified when boundary is instantiated
2472
2473This type of boundary is useful when stage is known at the boundary as a
2474function of time, but momentum should be set to zero. This is for example
2475the case where a boundary is needed in the ocean on the two sides perpendicular
2476to the coast to maintain the wave height of the incoming wave.
2477
2478This class is specific to the shallow water equation as it works with the
2479momentum quantities assumed to be the second and third conserved quantities.
2480
2481This boundary condition should not cause the numerical instabilities seen with the transmissive momentum
2482boundary conditions (see Page \pageref{pg: transmissive boundary} and
2483Page \pageref{pg: transmissive momentum set stage boundary}).
2484\end{classdesc}
2485
2486\begin{classdesc}{Dirichlet_discharge_boundary}{domain=None, stage0=None, wh0=None}
2487Module: \module{shallow_water.shallow_water_domain}
2488
2489\code{stage0} sets the stage.
2490
2491\code{wh0} sets momentum in the inward normal direction.
2492\end{classdesc}
2493
2494\subsection{User-defined boundary conditions}
2495\label{sec:roll your own}
2496
2497All boundary classes must inherit from the generic boundary class
2498\code{Boundary} and have a method called \code{evaluate()} which must
2499take as inputs \code{self}, \code{vol_id} and \code{edge_id} where \code{self} refers to the
2500object itself and \code{vol_id} and \code{edge_id} are integers referring to
2501particular edges. The method must return a list of three floating point
2502numbers representing values for \code{stage},
2503\code{xmomentum} and \code{ymomentum}, respectively.
2504
2505The constructor of a particular boundary class may be used to specify
2506particular values or flags to be used by the \code{evaluate()} method.
2507Please refer to the source code for the existing boundary conditions
2508for examples of how to implement boundary conditions.
2509
2510
2511\section{Forcing Terms}\index{Forcing terms}
2512\label{sec:forcing terms}
2513
2514\anuga provides a number of predefined forcing functions to be used with simulations.
2515Gravity and friction are always calculated using the elevation and friction quantities,
2516but the user may additionally add forcing terms to the list
2517\code{domain.forcing_terms} and have them affect the model.
2518
2519Currently, predefined forcing terms are: \\
2520\begin{classdesc}{General_forcing}{domain,
2521                                  quantity_name,
2522                                  rate=0.0,
2523                                  center=None,
2524                                  radius=None,
2525                                  polygon=None,
2526                                  default_rate=None,
2527                                  verbose=False}
2528Module: \module{shallow_water.shallow_water_domain}
2529
2530This is a general class to modify any quantity according to a given rate of change.
2531Other specific forcing terms are based on this class but it can be used by itself as well (e.g.\ to modify momentum).
2532
2533\code{domain} is the domain being evolved.
2534
2535\code{quantity_name} is the name of the quantity that will be affected by this forcing term.
2536
2537\code{rate} is the rate at which the quantity should change. This can be either a constant or a
2538function of time. Positive values indicate increases, negative values indicate decreases.
2539The parameter \code{rate} can be \code{None} at initialisation but must be specified
2540before a forcing term is applied (i.e.\ simulation has started).
2541The default value is 0.0 -- i.e.\ no forcing.
2542
2543\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
2544
2545\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
2546
2547Note: specifying \code{center}, \code{radius} and \code{polygon} will cause an exception to be thrown.
2548Moreover, if the specified polygon or circle does not lie fully within the mesh boundary an Exception will be thrown.
2549
2550Example:
2551
2552\begin{verbatim}
2553P = [[x0, y0], [x1, y0], [x1, y1], [x0, y1]]    # Square polygon
2554
2555xmom = General_forcing(domain, 'xmomentum', polygon=P)
2556ymom = General_forcing(domain, 'ymomentum', polygon=P)
2557
2558xmom.rate = f
2559ymom.rate = g
2560
2561domain.forcing_terms.append(xmom)
2562domain.forcing_terms.append(ymom)
2563\end{verbatim}
2564
2565Here, \code{f} and \code{g} are assumed to be defined as functions of time providing
2566a time dependent rate of change for xmomentum and ymomentum respectively.
2567\code{P} is assumed to be the polygon, specified as a list of points.
2568\end{classdesc}
2569
2570\begin{classdesc}{Inflow}{domain,
2571                         rate=0.0,
2572                         center=None, radius=None,
2573                         polygon=None,
2574                         default_rate=None,
2575                         verbose=False}
2576Module: \module{shallow_water.shallow_water_domain}
2577
2578This is a general class for inflow and abstraction of water according to a given rate of change.
2579This class will always modify the \code{stage} quantity.
2580
2581Inflow is based on the \code{General_forcing} class so the functionality is similar.
2582
2583\code{domain} is the domain being evolved.
2584
2585\code{rate} is the flow rate ($m^3/s$) at which the quantity should change. This can be either a constant or a
2586function of time. Positive values indicate inflow, negative values indicate outflow.
2587Note: The specified flow will be divided by the area of the inflow region and then applied to update the
2588stage quantity.
2589
2590\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
2591
2592\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
2593
2594Example:
2595
2596\begin{verbatim}
2597hydrograph = Inflow(center=(320, 300), radius=10,
2598                    rate=file_function('QPMF_Rot_Sub13.tms'))
2599
2600domain.forcing_terms.append(hydrograph)
2601\end{verbatim}
2602
2603Here, \code{'QPMF_Rot_Sub13.tms'} is assumed to be a NetCDF file in the TMS format defining a timeseries for a hydrograph.
2604\end{classdesc}
2605
2606\begin{classdesc}{Rainfall}{domain,
2607                           rate=0.0,
2608                           center=None,
2609                           radius=None,
2610                           polygon=None,
2611                           default_rate=None,
2612                           verbose=False}
2613Module: \module{shallow_water.shallow_water_domain}
2614
2615This is a general class for implementing rainfall over the domain, possibly restricted to a given circle or polygon.
2616This class will always modify the \code{stage} quantity.
2617
2618Rainfall is based on the \code{General_forcing} class so the functionality is similar.
2619
2620\code{domain} is the domain being evolved.
2621
2622\code{rate} is the total rain rate over the specified domain.
2623Note: Raingauge Data needs to reflect the time step.
2624For example, if rain gauge is \code{mm} read every \code{dt} seconds, then the input
2625here is as \code{mm/dt} so 10 mm in 5 minutes becomes
262610/(5x60) = 0.0333mm/s.  This parameter can be either a constant or a
2627function of time. Positive values indicate rain being added (or be used for general infiltration),
2628negative values indicate outflow at the specified rate (presumably this could model evaporation or abstraction).
2629
2630\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
2631
2632\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
2633
2634Example:
2635
2636\begin{verbatim}
2637catchmentrainfall = Rainfall(rate=file_function('Q100_2hr_Rain.tms'))
2638domain.forcing_terms.append(catchmentrainfall)
2639\end{verbatim}
2640
2641Here, \code{'Q100_2hr_Rain.tms'} is assumed to be a NetCDF file in the TMS format defining a timeseries for the rainfall.
2642\end{classdesc}
2643
2644\begin{classdesc}{Culvert_flow}{domain,
2645                 culvert_description_filename=None,
2646                 culvert_routine=None,
2647                 end_point0=None,
2648                 end_point1=None,
2649                 enquiry_point0=None,
2650                 enquiry_point1=None,
2651                 type='box',
2652                 width=None,
2653                 height=None,
2654                 length=None,
2655                 number_of_barrels=1,
2656                 trigger_depth=0.01,
2657                 manning=None,
2658                 sum_loss=None,
2659                 use_velocity_head=False,
2660                 use_momentum_jet=False,
2661                 label=None,
2662                 description=None,
2663                 update_interval=None,
2664                 log_file=False,
2665                 discharge_hydrograph=False,
2666                 verbose=False}
2667Module: \module{culvert_flows.culvert_class}
2668
2669This is a general class for implementing flow through a culvert.
2670This class modifies the quantities \code{stage}, \code{xmomentum} and \code{ymomentum} in areas at both ends of the culvert.
2671
2672The \code{Culvert_flow} forcing term uses \code{Inflow} and \code{General_forcing} to update the quantities.
2673The flow direction is determined on-the-fly so openings are referenced simple as opening0 and opening1
2674with either being able to take the role as Inflow or Outflow.
2675
2676The \code{Culvert_flow} class takes as input:
2677\begin{itemize}
2678  \item \code{domain}: a reference to the domain being evolved
2679  \item \code{culvert_description_filename}:
2680  \item \code{culvert_routine}:
2681  \item \code{end_point0}: Coordinates of one opening
2682  \item \code{end_point1}: Coordinates of other opening
2683  \item \code{enquiry_point0}:
2684  \item \code{enquiry_point1}:
2685  \item \code{type}: (default is 'box')
2686  \item \code{width}:
2687  \item \code{height}:
2688  \item \code{length}:
2689  \item \code{number_of_barrels}: Number of identical pipes in the culvert (default is 1)
2690  \item \code{trigger_depth}: (default is 0.01)
2691  \item \code{manning}: Mannings Roughness for Culvert
2692  \item \code{sum_loss}:
2693  \item \code{use_velocity_head}:
2694  \item \code{use_momentum_jet}:
2695  \item \code{label}: Short text naming the culvert
2696  \item \code{description}: Text describing the culvert
2697  \item \code{update_interval}:
2698  \item \code{log_file}:
2699  \item \code{discharge_hydrograph}:
2700\end{itemize}
2701
2702The user can specify different culvert routines. Hower \anuga currently provides only one, namely the
2703\code{boyd_generalised_culvert_model} as used in the example below:
2704
2705\begin{verbatim}
2706from anuga.culvert_flows.culvert_class import Culvert_flow
2707from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
2708
2709culvert1 = Culvert_flow(domain,
2710                        label='Culvert No. 1',
2711                        description='This culvert is a test unit 1.2m Wide by 0.75m High',
2712                        end_point0=[9.0, 2.5],
2713                        end_point1=[13.0, 2.5],
2714                        width=1.20,
2715                        height=0.75,
2716                        culvert_routine=boyd_generalised_culvert_model,
2717                        number_of_barrels=1,
2718                        verbose=True)
2719
2720culvert2 = Culvert_flow(domain,
2721                        label='Culvert No. 2',
2722                        description='This culvert is a circular test with d=1.2m',
2723                        end_point0=[9.0, 1.5],
2724                        end_point1=[30.0, 3.5],
2725                        diameter=1.20,
2726                        invert_level0=7,
2727                        culvert_routine=boyd_generalised_culvert_model,
2728                        number_of_barrels=1,
2729                        verbose=True)
2730
2731domain.forcing_terms.append(culvert1)
2732domain.forcing_terms.append(culvert2)
2733\end{verbatim}
2734\end{classdesc}
2735
2736
2737\section{Evolution}\index{evolution}
2738\label{sec:evolution}
2739
2740\begin{methoddesc}{\emph{<domain>}.evolve}{yieldstep=None,
2741                                           finaltime=None,
2742                                           duration=None,
2743                                           skip_initial_step=False}
2744Module: \module{abstract_2d_finite_volumes.domain}
2745
2746This method is invoked once all the
2747preliminaries have been completed, and causes the model to progress
2748through successive steps in its evolution, storing results and
2749outputting statistics whenever a user-specified period
2750\code{yieldstep} is completed.  Generally during this period the
2751model will evolve through several steps internally
2752as the method forces the water speed to be calculated
2753on successive new cells.
2754
2755\code{yieldstep} is the interval in seconds between yields where results are
2756stored, statistics written and the domain is inspected or possibly modified.
2757If omitted an internal predefined \code{yieldstep} is used.  Internally, smaller
2758timesteps may be taken.
2759
2760\code{duration} is the duration of the simulation in seconds.
2761
2762\code{finaltime} is the time in seconds where simulation should end. This is currently
2763relative time, so it's the same as \code{duration}.  If both \code{duration} and
2764\code{finaltime} are given an exception is thrown.
2765
2766\code{skip_initial_step} is a boolean flag that decides whether the first
2767yield step is skipped or not. This is useful for example to avoid
2768duplicate steps when multiple evolve processes are dove tailed.
2769
2770The code specified by the user in the block following the evolve statement is
2771only executed once every \code{yieldstep} even though
2772\anuga typically will take many more internal steps behind the scenes.
2773
2774You can include \method{evolve} in a statement of the type:
2775
2776\begin{verbatim}
2777for t in domain.evolve(yieldstep, finaltime):
2778    <Do something with domain and t>
2779\end{verbatim}
2780\end{methoddesc}
2781
2782\subsection{Diagnostics}
2783\label{sec:diagnostics}
2784
2785\begin{methoddesc}{\emph{<domain>}.statistics}{}
2786Module: \module{abstract\_2d\_finite\_volumes.domain}
2787
2788Outputs statistics about the mesh within the \code{Domain}.
2789\end{methoddesc}
2790
2791\begin{methoddesc}{\emph{<domain>}.timestepping_statistics}{track_speeds=False, triangle_id=None}
2792Module: \module{abstract_2d_finite_volumes.domain}
2793
2794Returns a string of the following type for each timestep:\\
2795\code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12 (12)}
2796
2797Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
2798the number of first-order steps, respectively.
2799
2800The optional keyword argument \code{track_speeds} will
2801generate a histogram of speeds generated by each triangle if set to \code{True}. The
2802speeds relate to the size of the timesteps used by \anuga and
2803this diagnostics may help pinpoint problem areas where excessive speeds
2804are generated.
2805
2806The optional keyword argument \code{triangle_id} can be used to specify a particular
2807triangle rather than the one with the largest speed.
2808\end{methoddesc}
2809
2810\begin{methoddesc}{\emph{<domain>}.boundary_statistics}{quantities=None,
2811                                                      tags=None}
2812Module: \module{abstract_2d_finite_volumes.domain}
2813
2814Generates output about boundary forcing at each timestep.
2815
2816\code{quantities} names the quantities to be reported -- may be \code{None},
2817a string or a list of strings.
2818
2819\code{tags} names the tags to be reported -- may be either None, a string or a list of strings.
2820
2821When \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}
2822will return a string like:
2823
2824\begin{verbatim}
2825Boundary values at time 0.5000:
2826    top:
2827        stage in [ -0.25821218,  -0.02499998]
2828    bottom:
2829        stage in [ -0.27098821,  -0.02499974]
2830\end{verbatim}
2831\end{methoddesc}
2832
2833\begin{methoddesc}{Q = \emph{<domain>}.get_quantity}{name,
2834                                               location='vertices',
2835                                               indices=None}
2836Module: \module{abstract_2d_finite_volumes.domain}
2837
2838This function returns a Quantity object Q.
2839Access to its values should be done through \code{Q.get_values()} documented on Page \pageref{pg:get values}.
2840
2841\code{name} is the name of the quantity to retrieve.
2842
2843\code{location} is
2844
2845\code{indices} is
2846\end{methoddesc}
2847
2848\begin{methoddesc}{\emph{<domain>}.set_quantities_to_be_monitored}{quantity,
2849                                             polygon=None,
2850                                             time_interval=None}
2851Module: \module{abstract\_2d\_finite\_volumes.domain}
2852
2853Selects quantities and derived quantities for which extrema attained at internal timesteps
2854will be collected.
2855
2856\code{quantity} specifies the quantity or quantities to be monitored and must be either:
2857\begin{itemize}
2858  \item the name of a quantity or derived quantity such as 'stage-elevation',
2859  \item a list of quantity names, or
2860  \item \code{None}.
2861\end{itemize}
2862
2863\code{polygon} can be used to monitor only triangles inside the polygon. Otherwise
2864all triangles will be included.
2865
2866\code{time_interval} will restrict monitoring to time steps in the interval. Otherwise
2867all timesteps will be included.
2868
2869Information can be tracked in the evolve loop by printing \code{quantity_statistics} and
2870collected data will be stored in the SWW file.
2871\end{methoddesc}
2872
2873\begin{methoddesc}{\emph{<domain>}.quantity_statistics}{precision='\%.4f'}
2874Module: \module{abstract_2d_finite_volumes.domain}
2875
2876Reports on extrema attained by selected quantities.
2877
2878Returns a string of the following type for each timestep:
2879
2880\begin{verbatim}
2881Monitored quantities at time 1.0000:
2882  stage-elevation:
2883    values since time = 0.00 in [0.00000000, 0.30000000]
2884    minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333)
2885    maximum attained at time = 0.00000000, location = (0.83333333, 0.16666667)
2886  ymomentum:
2887    values since time = 0.00 in [0.00000000, 0.06241221]
2888    minimum attained at time = 0.00000000, location = (0.33333333, 0.16666667)
2889    maximum attained at time = 0.22472667, location = (0.83333333, 0.66666667)
2890  xmomentum:
2891    values since time = 0.00 in [-0.06062178, 0.47886313]
2892    minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333)
2893    maximum attained at time = 0.35103646, location = (0.83333333, 0.16666667)
2894\end{verbatim}
2895
2896The quantities (and derived quantities) listed here must be selected at model
2897initialisation time using the method \code{domain.set_quantities_to_be_monitored()}.
2898
2899The optional keyword argument \code{precision='\%.4f'} will
2900determine the precision used for floating point values in the output.
2901This diagnostics helps track extrema attained by the selected quantities
2902at every internal timestep.
2903
2904These values are also stored in the SWW file for post-processing.
2905\end{methoddesc}
2906
2907\begin{methoddesc}{\emph{<quantity>}.get_values}{interpolation_points=None,
2908                   location='vertices',
2909                   indices=None,
2910                   use_cache=False,
2911                   verbose=False}
2912\label{pg:get values}
2913Module: \module{abstract_2d_finite_volumes.quantity}
2914
2915Extract values for quantity as a numeric array.
2916
2917\code{interpolation_points} is a list of (x, y) coordinates where the value is
2918sought (using interpolation). If \code{interpolation_points} is given, values
2919for \code{location} and \code{indices} are ignored.
2920Assume either an absolute UTM coordinates or geospatial data object.
2921   
2922\code{location} specifies where values are to be stored.
2923Permissible options are 'vertices', 'edges', 'centroids' or 'unique vertices'.
2924
2925The returned values will have the leading dimension equal to length of the \code{indices} list or
2926N (all values) if \code{indices} is \code{None}.
2927
2928If \code{location} is 'centroids' the dimension of returned
2929values will be a list or a numeric array of length N, N being
2930the number of elements.
2931     
2932If \code{location} is 'vertices' or 'edges' the dimension of
2933returned values will be of dimension \code{Nx3}.
2934
2935If \code{location} is 'unique vertices' the average value at
2936each vertex will be returned and the dimension of returned values
2937will be a 1d array of length "number of vertices"
2938     
2939\code{indices} is the set of element ids that the operation applies to.
2940
2941The values will be stored in elements following their internal ordering.
2942\end{methoddesc}
2943
2944\begin{methoddesc}{\emph{<quantity>}.set_values}{numeric=None,
2945                               quantity=None,
2946                               function=None,
2947                               geospatial_data=None,
2948                               filename=None,
2949                               attribute_name=None,
2950                               alpha=None,
2951                               location='vertices',
2952                               polygon=None,
2953                               indices=None,
2954                               smooth=False,
2955                               verbose=False,
2956                               use_cache=False}
2957Module: \module{abstract_2d_finite_volumes.quantity}
2958
2959Assign values to a quantity object.
2960
2961This method works the same way as \code{set_quantity()} except that it doesn't take
2962a quantity name as the first argument since it is applied directly to the quantity.
2963Use \code{set_values} is used to assign
2964values to a new quantity that has been created but which is
2965not part of the domain's predefined quantities.
2966
2967\code{location} is ??
2968
2969\code{indices} is ??
2970
2971The method \code{set_values()} is always called by \code{set_quantity()}
2972behind the scenes.
2973\end{methoddesc}
2974
2975\begin{methoddesc}{\emph{<quantity>}.get_integral}{}
2976Module: \module{abstract_2d_finite_volumes.quantity}
2977
2978Return the computed integral over the entire domain for the quantity.
2979\end{methoddesc}
2980
2981\begin{methoddesc}{\emph{<quantity>}.get_maximum_value}{indices=None}
2982Module: \module{abstract_2d_finite_volumes.quantity}
2983
2984Return the maximum value of a quantity (on centroids).
2985
2986\code{indices} is the optional set of element \code{id}s that
2987the operation applies to.
2988
2989We do not seek the maximum at vertices as each vertex can
2990have multiple values -- one for each triangle sharing it.
2991\end{methoddesc}
2992
2993\begin{methoddesc}{\emph{<quantity>}.get_maximum_location}{indices=None}
2994Module: \module{abstract_2d_finite_volumes.quantity}
2995
2996Return the location of the maximum value of a quantity (on centroids).
2997
2998\code{indices} is the optional set of element \code{id}s that
2999the operation applies to.
3000
3001We do not seek the maximum at vertices as each vertex can
3002have multiple values -- one for each triangle sharing it.
3003
3004If there are multiple cells with the same maximum value, the
3005first cell encountered in the triangle array is returned.
3006\end{methoddesc}
3007
3008\begin{methoddesc}{\emph{<domain>}.get_wet_elements}{indices=None}
3009Module: \module{shallow_water.shallow_water_domain}
3010
3011Returns the indices for elements where h $>$ minimum_allowed_height
3012
3013\code{indices} is the optional set of element \code{id}s that
3014the operation applies to.
3015\end{methoddesc}
3016
3017\begin{methoddesc}{\emph{<domain>}.get_maximum_inundation_elevation}{indices=None}
3018Module: \module{shallow_water.shallow_water_domain}
3019
3020Return highest elevation where h $>$ 0.
3021
3022\code{indices} is the optional set of element \code{id}s that
3023the operation applies to.
3024
3025Example to find maximum runup elevation:
3026\begin{verbatim}
3027z = domain.get_maximum_inundation_elevation()
3028\end{verbatim}
3029\end{methoddesc}
3030
3031\begin{methoddesc}{\emph{<domain>}.get_maximum_inundation_location}{indices=None}
3032Module: \module{shallow_water.shallow_water_domain}
3033
3034Return location (x,y) of highest elevation where h $>$ 0.
3035
3036\code{indices} is the optional set of element \code{id}s that
3037the operation applies to.
3038
3039Example to find maximum runup location:
3040\begin{verbatim}
3041x, y = domain.get_maximum_inundation_location()
3042\end{verbatim}
3043\end{methoddesc}
3044
3045
3046\section{Queries of SWW model output files}
3047After a model has been run, it is often useful to extract various information from the SWW
3048output file (see Section \ref{sec:sww format}). This is typically more convenient than using the
3049diagnostics described in Section \ref{sec:diagnostics} which rely on the model running -- something
3050that can be very time consuming. The SWW files are easy and quick to read and offer information
3051about the model results such as runup heights, time histories of selected quantities,
3052flow through cross sections and much more.
3053
3054\begin{funcdesc}{elevation = get_maximum_inundation_elevation}{filename,
3055                                     polygon=None,
3056                                     time_interval=None,
3057                                     verbose=False}
3058Module: \module{shallow_water.data_manager}
3059
3060Return the highest elevation where depth is positive ($h > 0$).
3061
3062\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
3063
3064\code{polygon} restricts the query to the points that lie within the polygon.
3065
3066\code {time_interval} restricts the query to within the time interval.
3067
3068Example to find maximum runup elevation:
3069
3070\begin{verbatim}
3071max_runup = get_maximum_inundation_elevation(filename)
3072\end{verbatim}
3073
3074If no inundation is found (within the \code{polygon} and \code{time_interval}, if specified)
3075the return value is \code{None}. This indicates "No Runup" or "Everything is dry".
3076\end{funcdesc}
3077
3078\begin{funcdesc}{(x, y) = get_maximum_inundation_location}{filename,
3079                                    polygon=None,
3080                                    time_interval=None,
3081                                    verbose=False}
3082Module: \module{shallow_water.data_manager}
3083
3084Return location (x,y) of the highest elevation where depth is positive ($h > 0$).
3085
3086\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
3087
3088\code{polygon} restricts the query to the points that lie within the polygon.
3089
3090\code {time_interval} restricts the query to within the time interval.
3091
3092Example to find maximum runup location:
3093
3094\begin{verbatim}
3095max_runup_location = get_maximum_inundation_location(filename)
3096\end{verbatim}
3097
3098If no inundation is found (within the \code{polygon} and \code{time_interval}, if specified)
3099the return value is \code{None}. This indicates "No Runup" or "Everything is dry".
3100is \code{None}. This indicates "No Runup" or "Everything is dry".
3101\end{funcdesc}
3102
3103\begin{funcdesc}{sww2timeseries}{swwfiles,
3104                                 gauge_filename,
3105                                 production_dirs,
3106                                 report=None,
3107                                 reportname=None,
3108                                 plot_quantity=None,
3109                                 generate_fig=False,
3110                                 surface=None,
3111                                 time_min=None,
3112                                 time_max=None,
3113                                 time_thinning=1,
3114                                 time_unit=None,
3115                                 title_on=None,
3116                                 use_cache=False,
3117                                 verbose=False}
3118Module: \module{abstract_2d_finite_volumes.util}
3119
3120Read a set of SWW files and plot the time series for the prescribed quantities
3121at defined gauge locations and prescribed time range.
3122
3123\code{swwfiles} is a dictionary of SWW files with keys of \code{label_id}.
3124
3125\code{gauge_filename} is the path to a file containing gauge data.
3126
3127THIS NEEDS MORE WORK.  WORK ON FUNCTION __DOC__ STRING, IF NOTHING ELSE!
3128\end{funcdesc}
3129
3130\begin{funcdesc}{(time, Q) = get_flow_through_cross_section}{filename, polyline, verbose=False}
3131Module: \module{shallow_water.data_manager}
3132
3133Obtain flow ($m^3/s$) perpendicular to specified cross section.
3134
3135\code{filename} is the path to the SWW file.
3136
3137\code{polyline} is the representation of the desired cross section -- it may contain
3138multiple sections allowing for complex shapes. Assumes absolute UTM coordinates.
3139
3140Returns a tuple \code{time, Q} where:
3141
3142\code{time} is all the stored times in the SWW file.
3143
3144\code{Q} is a hydrograph of total flow across the given segments for all stored times.
3145
3146The normal flow is computed for each triangle intersected by the \code{polyline} and
3147added up.  If multiple segments at different angles are specified the normal flows
3148may partially cancel each other.
3149
3150Example to find flow through cross section:
3151
3152\begin{verbatim}
3153cross_section = [[x, 0], [x, width]]
3154time, Q = get_flow_through_cross_section(filename, cross_section)
3155\end{verbatim}
3156\end{funcdesc}
3157
3158
3159\section{Other}
3160\begin{methoddesc}{quantity = \emph{<domain>}.create_quantity_from_expression}{string}
3161Module: \module{abstract_2d_finite_volumes.domain}
3162
3163Create a new quantity from other quantities in the domain using an arbitrary expression.
3164
3165\code{string} is a string containing an arbitrary quantity expression.
3166
3167Returns the new \code{Quantity} object.
3168
3169Handy for creating derived quantities on-the-fly:
3170
3171\begin{verbatim}
3172Depth = domain.create_quantity_from_expression('stage-elevation')
3173
3174exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
3175Absolute_momentum = domain.create_quantity_from_expression(exp)
3176\end{verbatim}
3177
3178%See also \file{Analytical_solution_circular_hydraulic_jump.py} for an example.
3179\end{methoddesc}
3180
3181%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3182
3183\chapter{\anuga System Architecture}
3184
3185\section{File Formats}
3186\label{sec:file formats}
3187
3188\anuga makes use of a number of different file formats. The
3189following table lists all these formats, which are described in more
3190detail in the paragraphs below.
3191
3192\begin{center}
3193\begin{tabular}{|ll|}  \hline
3194  \textbf{Extension} & \textbf{Description} \\
3195  \hline\hline
3196  \code{.sww} & NetCDF format for storing model output with mesh information \code{f(t,x,y)}\\
3197  \code{.sts} & NetCDF format for storing model ouput \code{f(t,x,y)} without mesh information\\
3198  \code{.tms} & NetCDF format for storing time series \code{f(t)}\\
3199  \code{.csv/.txt} & ASCII format for storing arbitrary points and associated attributes\\
3200  \code{.pts} & NetCDF format for storing arbitrary points and associated attributes\\
3201  \code{.asc} & ASCII format of regular DEMs as output from ArcView\\
3202  \code{.prj} & Associated ArcView file giving more metadata for \code{.asc} format\\
3203  \code{.ers} & ERMapper header format of regular DEMs for ArcView\\
3204  \code{.dem} & NetCDF representation of regular DEM data\\
3205  \code{.tsh} & ASCII format for storing meshes and associated boundary and region info\\
3206  \code{.msh} & NetCDF format for storing meshes and associated boundary and region info\\
3207  \code{.nc} & Native ferret NetCDF format\\
3208  \code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
3209\end{tabular}
3210\end{center}
3211
3212The above table shows the file extensions used to identify the
3213formats of files. However, typically, in referring to a format we
3214capitalise the extension and omit the initial full stop -- thus, we
3215refer to 'SWW files' or 'PRJ files', not 'sww files' or '.prj files'.
3216
3217\bigskip
3218
3219A typical dataflow can be described as follows:
3220
3221SOMETHING MISSING HERE!?
3222
3223\subsection{Manually Created Files}
3224
3225\begin{tabular}{ll}
3226ASC, PRJ & Digital elevation models (gridded)\\
3227NC & Model outputs for use as boundary conditions (e.g.\ from MOST)
3228\end{tabular}
3229
3230\subsection{Automatically Created Files}
3231
3232\begin{tabular}{ll}
3233  ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert DEMs to native \code{.pts} file\\
3234  NC $\rightarrow$ SWW & Convert MOST boundary files to boundary \code{.sww}\\
3235  PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\
3236  TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using \code{animate}\\
3237  TSH + Boundary SWW $\rightarrow$ SWW & Simulation using \code{\anuga}\\
3238  Polygonal mesh outline $\rightarrow$ & TSH or MSH
3239\end{tabular}
3240
3241\bigskip
3242
3243\subsection{SWW, STS and TMS Formats}
3244\label{sec:sww format}
3245The SWW, STS and TMS formats are all NetCDF formats and are of key importance for \anuga.
3246
3247An SWW file is used for storing \anuga output and therefore pertains
3248to a set of points and a set of times at which a model is evaluated.
3249It contains, in addition to dimension information, the following
3250variables:
3251
3252\begin{itemize}
3253  \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
3254  \item \code{elevation}: a numeric array storing bed-elevations
3255  \item \code{volumes}: a list specifying the points at the vertices of each of the triangles
3256    % Refer here to the example to be provided in describing the simple example
3257  \item \code{time}: a numeric array containing times for model evaluation
3258\end{itemize}
3259
3260The contents of an SWW file may be viewed using the anuga viewer \code{animate},
3261which creates an on-screen visialisation.  See section \ref{sec:animate}
3262(page \pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more on \code{animate}.
3263
3264Alternatively, there are tools, such as \code{ncdump}, that allow
3265you to convert a NetCDF file into a readable format such as the
3266Class Definition Language (CDL). The following is an excerpt from a
3267CDL representation of the output file \file{runup.sww} generated
3268from running the simple example \file{runup.py} of Chapter \ref{ch:getstarted}:
3269
3270%FIXME (Ole): Should put in example with nonzero xllcorner, yllcorner
3271\verbatiminput{examples/bedslopeexcerpt.cdl}
3272
3273The SWW format is used not only for output but also serves as input
3274for functions such as \function{file\_boundary} and
3275\function{file\_function}, described in Chapter \ref{ch:interface}.
3276
3277An STS file is used for storing a set of points and associated times.
3278It contains, in addition to dimension information, the following
3279variables:
3280\begin{itemize}
3281  \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
3282  \item \code{permutation}: Original indices of the points as specified by the optional \code{ordering_file} 
3283                            (see the function \code{urs2sts()} in Section \ref{sec:basicfileconversions})
3284  \item \code{elevation}: a numeric array storing bed-elevations
3285    % Refer here to the example to be provided in describing the simple example
3286  \item \code{time}: a numeric array containing times for model evaluation
3287\end{itemize}
3288
3289The only difference between the STS format and the SWW format is the former does
3290not contain a list specifying the points at the vertices of each of the triangles
3291(\code{volumes}). Consequently information (arrays) stored within an STS file such
3292as \code{elevation} can be accessed in exactly the same way as it would be extracted
3293from an SWW file.
3294
3295A TMS file is used to store time series data that is independent of position.
3296
3297\subsection{Mesh File Formats}
3298
3299A mesh file is a file that has a specific format suited to
3300triangular meshes and their outlines. A mesh file can have one of
3301two formats: it can be either a TSH file, which is an ASCII file, or
3302an MSH file, which is a NetCDF file. A mesh file can be generated
3303from the function \function{create_mesh_from_regions()} (see
3304Section \ref{sec:meshgeneration}) and be used to initialise a domain.
3305
3306A mesh file can define the outline of the mesh -- the vertices and
3307line segments that enclose the region in which the mesh is
3308created -- and the triangular mesh itself, which is specified by
3309listing the triangles and their vertices, and the segments, which
3310are those sides of the triangles that are associated with boundary
3311conditions.
3312
3313In addition, a mesh file may contain 'holes' and/or 'regions'. A
3314hole represents an area where no mesh is to be created, while a
3315region is a labelled area used for defining properties of a mesh,
3316such as friction values.  A hole or region is specified by a point
3317and bounded by a number of segments that enclose that point.
3318
3319A mesh file can also contain a georeference, which describes an
3320offset to be applied to $x$ and $y$ values -- e.g.\ to the vertices.
3321
3322\subsection{Formats for Storing Arbitrary Points and Attributes}
3323
3324A CSV/TXT file is used to store data representing
3325arbitrary numerical attributes associated with a set of points.
3326
3327The format for an CSV/TXT file is:\\
3328 \\
3329first line:   \code{[column names]}\\
3330other lines:  \code{[x value], [y value], [attributes]}\\
3331
3332for example:
3333
3334\begin{verbatim}
3335x, y, elevation, friction
33360.6, 0.7, 4.9, 0.3
33371.9, 2.8, 5.0, 0.3
33382.7, 2.4, 5.2, 0.3
3339\end{verbatim}
3340
3341The delimiter is a comma. The first two columns are assumed to
3342be $x$ and $y$ coordinates.
3343
3344A PTS file is a NetCDF representation of the data held in an points CSV
3345file. If the data is associated with a set of $N$ points, then the
3346data is stored using an $N \times 2$ numeric array of float
3347variables for the points and an $N \times 1$ numeric array for each
3348attribute.
3349
3350\subsection{ArcView Formats}
3351
3352Files of the three formats ASC, PRJ and ERS are all associated with
3353data from ArcView.
3354
3355An ASC file is an ASCII representation of DEM output from ArcView.
3356It contains a header with the following format:
3357
3358\begin{tabular}{l l}
3359\code{ncols}      &   \code{753}\\
3360\code{nrows}      &   \code{766}\\
3361\code{xllcorner}  &   \code{314036.58727982}\\
3362\code{yllcorner}  & \code{6224951.2960092}\\
3363\code{cellsize}   & \code{100}\\
3364\code{NODATA_value} & \code{-9999}
3365\end{tabular}
3366
3367The remainder of the file contains the elevation data for each grid point
3368in the grid defined by the above information.
3369
3370A PRJ file is an ArcView file used in conjunction with an ASC file
3371to represent metadata for a DEM.
3372
3373\subsection{DEM Format}
3374
3375A DEM file in \anuga is a NetCDF representation of regular DEM data.
3376
3377\subsection{Other Formats}
3378
3379\subsection{Basic File Conversions}
3380\label{sec:basicfileconversions}
3381
3382\begin{funcdesc}{sww2dem}{(basename_in,
3383            basename_out=None,
3384            quantity=None,
3385            timestep=None,
3386            reduction=None,
3387            cellsize=10,
3388            number_of_decimal_places=None,
3389            NODATA_value=-9999,
3390            easting_min=None,
3391            easting_max=None,
3392            northing_min=None,
3393            northing_max=None,
3394            verbose=False,
3395            origin=None,
3396            datum='WGS84',
3397            format='ers',
3398            block_size=None}
3399Module: \module{shallow_water.data_manager}
3400
3401Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
3402ERS) of a desired grid size \code{cellsize} in metres. The user can select how
3403many decimal places the output data is represented with by using \code{number_of_decimal_places},
3404with the default being 3.
3405
3406The $easting$ and $northing$ values are used if the user wishes to determine the output
3407within a specified rectangular area. The \code{reduction} input refers to a function
3408to reduce the quantities over all time step of the SWW file, e.g.\ maximum.
3409\end{funcdesc}
3410
3411\begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
3412            easting_min=None, easting_max=None,
3413            northing_min=None, northing_max=None,
3414            use_cache=False, verbose=False}
3415  Module: \module{shallow\_water.data\_manager}
3416
3417  Takes DEM data (a NetCDF file representation of data from a regular Digital
3418  Elevation Model) and converts it to PTS format.
3419\end{funcdesc}
3420
3421\begin{funcdesc}{urs2sts}{basename_in, basename_out=None,
3422            weights=None, verbose=False,
3423            origin=None,mean_stage=0.0,
3424            zscale=1.0, ordering_filename=None}
3425  Module: \module{shallow\_water.data\_manager}
3426
3427  Takes URS data (timeseries data in mux2 format) and converts it to STS format.
3428  The optional filename \code{ordering\_filename} specifies the permutation of indices
3429  of points to select along with their longitudes and latitudes. This permutation will also be
3430  stored in the STS file. If absent, all points are taken and the permutation will be trivial,
3431  i.e.\ $0, 1, \ldots, N-1$, where $N$ is the total number of points stored. 
3432\end{funcdesc}
3433
3434\begin{funcdesc}{csv2building\_polygons}{file\_name, floor\_height=3}
3435  Module: \module{shallow\_water.data\_manager}
3436
3437  Convert CSV files of the form:
3438
3439  \begin{verbatim} 
3440easting,northing,id,floors
3441422664.22,870785.46,2,0
3442422672.48,870780.14,2,0
3443422668.17,870772.62,2,0
3444422660.35,870777.17,2,0
3445422664.22,870785.46,2,0
3446422661.30,871215.06,3,1
3447422667.50,871215.70,3,1
3448422668.30,871204.86,3,1
3449422662.21,871204.33,3,1
3450422661.30,871215.06,3,1
3451  \end{verbatim}
3452
3453  to a dictionary of polygons with \code{id} as key.
3454  The associated number of \code{floors} are converted to m above MSL and
3455  returned as a separate dictionary also keyed by \code{id}.
3456   
3457  Optional parameter \code{floor_height} is the height of each building story.
3458
3459  These can e.g.\ be converted to a \code{Polygon_function} for use with \code{add_quantity}
3460  as shown on page \pageref{add quantity}.
3461\end{funcdesc}
3462
3463%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3464
3465\chapter{\anuga mathematical background}
3466\label{cd:mathematical background}
3467
3468
3469\section{Introduction}
3470
3471This chapter outlines the mathematics underpinning \anuga.
3472
3473
3474\section{Model}
3475\label{sec:model}
3476
3477The shallow water wave equations are a system of differential
3478conservation equations which describe the flow of a thin layer of
3479fluid over terrain. The form of the equations are:
3480\[
3481\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
3482x}+\frac{\partial \GG}{\partial y}=\SSS
3483\]
3484where $\UU=\left[ {{\begin{array}{*{20}c}
3485 h & {uh} & {vh} \\
3486\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
3487$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
3488entering the system are bed elevation $z$ and stage (absolute water
3489level) $w$, where the relation $w = z + h$ holds true at all times.
3490The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
3491by
3492\[
3493\EE=\left[ {{\begin{array}{*{20}c}
3494 {uh} \hfill \\
3495 {u^2h+gh^2/2} \hfill \\
3496 {uvh} \hfill \\
3497\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
3498 {vh} \hfill \\
3499 {vuh} \hfill \\
3500 {v^2h+gh^2/2} \hfill \\
3501\end{array} }} \right]
3502\]
3503and the source term (which includes gravity and friction) is given
3504by
3505\[
3506\SSS=\left[ {{\begin{array}{*{20}c}
3507 0 \hfill \\
3508 -{gh(z_{x} + S_{fx} )} \hfill \\
3509 -{gh(z_{y} + S_{fy} )} \hfill \\
3510\end{array} }} \right]
3511\]
3512where $S_f$ is the bed friction. The friction term is modelled using
3513Manning's resistance law
3514\[
3515S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
3516=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
3517\]
3518in which $\eta$ is the Manning resistance coefficient.
3519The model does not currently include consideration of kinematic viscosity or dispersion.
3520
3521As demonstrated in our papers, \cite{ZR1999,nielsen2005} these
3522equations and their implementation in \anuga provide a reliable
3523model of general flows associated with inundation such as dam breaks
3524and tsunamis.
3525
3526
3527\section{Finite Volume Method}
3528\label{sec:fvm}
3529
3530We use a finite-volume method for solving the shallow water wave
3531equations \cite{ZR1999}. The study area is represented by a mesh of
3532triangular cells as in Figure~\ref{fig:mesh} in which the conserved
3533quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
3534in each volume are to be determined. The size of the triangles may
3535be varied within the mesh to allow greater resolution in regions of
3536particular interest.
3537
3538\begin{figure}[htp] \begin{center}
3539  \includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-five}
3540  \caption{Triangular mesh used in our finite volume method. Conserved
3541           quantities $h$, $uh$ and $vh$ are associated with the centroid of
3542           each triangular cell.}
3543  \label{fig:mesh}
3544\end{center} \end{figure}
3545
3546The equations constituting the finite-volume method are obtained by
3547integrating the differential conservation equations over each
3548triangular cell of the mesh. Introducing some notation we use $i$ to
3549refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
3550set of indices referring to the cells neighbouring the $i$th cell.
3551Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
3552the length of the edge between the $i$th and $j$th cells.
3553
3554By applying the divergence theorem we obtain for each volume an
3555equation which describes the rate of change of the average of the
3556conserved quantities within each cell, in terms of the fluxes across
3557the edges of the cells and the effect of the source terms. In
3558particular, rate equations associated with each cell have the form
3559$$
3560 \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
3561$$
3562where
3563\begin{itemize}
3564  \item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
3565  \item $\SSS_i$ is the source term associated with the $i$th cell, and
3566  \item $\HH_{ij}$ is the outward normal flux of material across the \textit{ij}th edge.
3567\end{itemize}
3568
3569%\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
3570%cells
3571%\item $m_{ij}$ is the midpoint of
3572%the \textit{ij}th edge,
3573%\item
3574%$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
3575%normal along the \textit{ij}th edge, and The
3576
3577The flux $\HH_{ij}$ is evaluated using a numerical flux function
3578$\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
3579water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
3580$$
3581H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
3582$$
3583
3584Then
3585$$
3586\HH_{ij}  = \HH(\UU_i(m_{ij}),
3587\UU_j(m_{ij}); \mathbf{n}_{ij})
3588$$
3589where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
3590$\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
3591\textit{ij}th edge. The function $\UU_i(x)$ for $x \in
3592T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
3593neighbouring  cells.
3594
3595We use a second order reconstruction to produce a piece-wise linear
3596function construction of the conserved quantities for  all $x \in
3597T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}). This
3598function is allowed to be discontinuous across the edges of the
3599cells, but the slope of this function is limited to avoid
3600artificially introduced oscillations.
3601
3602Godunov's method (see \cite{Toro1992}) involves calculating the
3603numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
3604solving the corresponding one dimensional Riemann problem normal to
3605the edge. We use the central-upwind scheme of \cite{KurNP2001} to
3606calculate an approximation of the flux across each edge.
3607
3608\begin{figure}[htp] \begin{center}
3609  \includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-reconstruct}
3610  \caption{From the values of the conserved quantities at the centroid
3611           of the cell and its neighbouring cells, a discontinuous piecewise
3612           linear reconstruction of the conserved quantities is obtained.}
3613  \label{fig:mesh:reconstruct}
3614\end{center} \end{figure}
3615
3616In the computations presented in this paper we use an explicit Euler
3617time stepping method with variable timestepping adapted to the
3618observed CFL condition:
3619
3620\begin{equation}
3621  \Delta t = \min_{k,i=[0,1,2]}  \min \left( \frac{r_k}{v_{k,i}}, \frac{r_{n_{k,i}}}{v_{k,i}} \right )
3622  \label{eq:CFL condition}
3623\end{equation}
3624where $r_k$ is the radius of the $k$'th triangle and $v_{k,i}$ is the maximal velocity across
3625edge joining triangle $k$ and it's $i$'th neighbour, triangle $n_{k,i}$, as calculated by the
3626numerical flux function
3627using the central upwind scheme of \cite{KurNP2001}. The symbol $r_{n_{k,i}}$  denotes the radius
3628of the $i$'th neighbour of triangle $k$. The radii are calculated as radii of the inscribed circles
3629of each triangle.
3630
3631
3632\section{Flux limiting}
3633
3634The shallow water equations are solved numerically using a
3635finite volume method on an unstructured triangular grid.
3636The upwind central scheme due to Kurganov and Petrova is used as an
3637approximate Riemann solver for the computation of inviscid flux functions.
3638This makes it possible to handle discontinuous solutions.
3639
3640To alleviate the problems associated with numerical instabilities due to
3641small water depths near a wet/dry boundary we employ a new flux limiter that
3642ensures that unphysical fluxes are never encounted.
3643
3644Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction,
3645$w$ the absolute water level (stage) and
3646$z$ the bed elevation. The latter are assumed to be relative to the
3647same height datum.
3648The conserved quantities tracked by \anuga are momentum in the
3649$x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$)
3650and depth ($h = w-z$).
3651
3652The flux calculation requires access to the velocity vector $(u, v)$
3653where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively.
3654In the presence of very small water depths, these calculations become
3655numerically unreliable and will typically cause unphysical speeds.
3656
3657We have employed a flux limiter which replaces the calculations above with
3658the limited approximations.
3659\begin{equation}
3660  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h},
3661\end{equation}
3662where $h_0$ is a regularisation parameter that controls the minimal
3663magnitude of the denominator. Taking the limits we have for $\hat{u}$
3664\[
3665  \lim_{h \rightarrow 0} \hat{u} =
3666  \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0
3667\]
3668and
3669\[
3670  \lim_{h \rightarrow \infty} \hat{u} =
3671  \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u
3672\]
3673with similar results for $\hat{v}$.
3674
3675The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator)
3676\[
3677  1 - h_0/h^2 = 0
3678\]
3679or
3680\[
3681  h_0 = h^2
3682\]
3683
3684\anuga has a global parameter $H_0$ that controls the minimal depth which
3685is considered in the various equations. This parameter is typically set to
3686$10^{-3}$. Setting
3687\[
3688  h_0 = H_0^2
3689\]
3690provides a reasonable balance between accurracy and stability. In fact,
3691setting $h=H_0$ will scale the predicted speed by a factor of $0.5$:
3692\[
3693  \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = 
3694  \left[ \frac{\mu}{H_0 + H_0^2/H_0} \right] = 
3695  \frac{\mu}{2 H_0} = \frac{\mu}{2 h} = \frac{u}{2}
3696\]
3697In general, for multiples of the minimal depth $N H_0$ one obtains
3698\begin{equation}
3699  \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} =
3700  \frac{\mu}{N H_0 + H_0/N} =   
3701  \frac{\mu}{h (1 + 1/N^2)} 
3702  \label{eq:flux limit multiple} 
3703\end{equation} 
3704which converges quadratically to the true value with the multiple N.
3705
3706Although this equation can be used for any depth, we have restricted its use to depths less than $10 * H_0$ (or 1 cm) to computational resources. 
3707According to Equation \ref{eq:flux limit multiple} this cutoff   
3708affects the calculated velocity by less than 1 \%.
3709
3710%The developed numerical model has been applied to several test cases
3711%as well as to real flows. numeric tests prove the robustness and accuracy of the model.
3712
3713
3714\section{Slope limiting}
3715A multidimensional slope-limiting technique is employed to achieve second-order spatial
3716accuracy and to prevent spurious oscillations. This is using the MinMod limiter and is
3717documented elsewhere.
3718
3719However close to the bed, the limiter must ensure that no negative depths occur.
3720 On the other hand, in deep water, the bed topography should be ignored for the
3721purpose of the limiter.
3722
3723Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
3724let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.
3725Define the minimal depth across all vertices as $\hmin$ as
3726\[
3727  \hmin = \min_i h_i
3728\]
3729
3730Let $\tilde{w_i}$ be the stage obtained from a gradient limiter
3731limiting on stage only. The corresponding depth is then defined as
3732\[
3733  \tilde{h_i} = \tilde{w_i} - z_i
3734\]
3735We would use this limiter in deep water which we will define (somewhat boldly)
3736as
3737\[
3738  \hmin \ge \epsilon
3739\]
3740
3741Similarly, let $\bar{w_i}$ be the stage obtained from a gradient
3742limiter limiting on depth respecting the bed slope.
3743The corresponding depth is defined as
3744\[
3745  \bar{h_i} = \bar{w_i} - z_i
3746\]
3747
3748We introduce the concept of a balanced stage $w_i$ which is obtained as
3749the linear combination
3750
3751\[
3752  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
3753\]
3754or
3755\[
3756  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
3757\]
3758where $\alpha \in [0, 1]$.
3759
3760Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope
3761is ignored we have immediately that
3762\[
3763  \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0
3764\]
3765%where the maximal bed elevation range $dz$ is defined as
3766%\[
3767%  dz = \max_i |z_i - z|
3768%\]
3769
3770If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that
3771no negative depths occur. Formally, we will require that
3772\[
3773  \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i
3774\]
3775or
3776\begin{equation}
3777  \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i
3778  \label{eq:limiter bound}
3779\end{equation}
3780
3781There are two cases:
3782\begin{enumerate}
3783  \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage)
3784  vertex is at least as far away from the bed than the shallow water
3785  (limited using depth). In this case we won't need any contribution from
3786  $\bar{h_i}$ and can accept any $\alpha$.
3787
3788  E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to
3789  \[
3790    \tilde{h_i} > \epsilon
3791  \]
3792  whereas $\alpha=0$ yields
3793  \[
3794    \bar{h_i} > \epsilon
3795  \]
3796  all well and good.
3797  \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is
3798  closer to the bed than the shallow water vertex or even below the bed.
3799  In this case we need to find an $\alpha$ that will ensure a positive depth.
3800  Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one
3801  obtains the bound
3802  \[
3803    \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i
3804  \]
3805\end{enumerate}
3806
3807Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one
3808arrives at the definition
3809\[
3810  \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}}
3811\]
3812which will guarantee that no vertex 'cuts' through the bed. Finally, should
3813$\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting
3814$\alpha=0$ and similarly capping $\alpha$ at 1 just in case.
3815
3816%Furthermore,
3817%dropping the $\epsilon$ ensures that alpha is always positive and also
3818%provides a numerical safety {??)
3819
3820%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3821
3822\chapter{Basic \anuga Assumptions}
3823
3824
3825\section{Time}
3826
3827Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
3828If one wished to recreate scenarios prior to that date it must be done
3829using some relative time (e.g.\ 0).
3830
3831The \anuga domain has an attribute \code{starttime} which is used in cases where the
3832simulation should be started later than the beginning of some input data such as those
3833obtained from boundaries or forcing functions (hydrographs, file_boundary etc).
3834
3835The \code{domain.startime} may be adjusted in \code{file_boundary} in the case the
3836input data does not itself start until a later time.
3837
3838 
3839\section{Spatial data}
3840
3841\subsection{Projection}
3842All spatial data relates to the WGS84 datum (or GDA94) and assumes a projection
3843into UTM with false easting of 500000 and false northing of
38441000000 on the southern hemisphere (0 on the northern hemisphere).
3845All locations must consequently be specified in Cartesian coordinates
3846(eastings, northings) or (x,y) where the unit is metres.
3847Alternative projections can be used, but \anuga does have the concept of a UTM zone
3848that must be the same for all coordinates in the model.
3849
3850\subsection{Internal coordinates}
3851It is important to realise that for numerical precision \anuga uses coordinates that are relative
3852to the lower left node of the rectangle containing the mesh ($x_{\mbox{min}}$, $y_{\mbox{min}}$).
3853This origin is referred to internally as \code{xllcorner}, \code{yllcorner} following the ESRI ASCII grid notation. 
3854The SWW file format also includes \code{xllcorner}, \code{yllcorner} and any coordinate in the file should be adjusted
3855by adding this origin. See Section \ref{sec:sww format}.
3856
3857Throughout the \anuga interface, functions have optional boolean arguments \code{absolute} which controls
3858whether spatial data received is using the internal representation (\code{absolute=False}) or the
3859user coordinate set (\code{absolute=True}). See e.g.\ \code{get_vertex_coordinates()} on \pageref{pg:get vertex coordinates}.
3860
3861DEMs, meshes and boundary conditions can have different origins. However, the internal representation in \anuga 
3862will use the origin of the mesh.
3863
3864\subsection{Polygons}
3865When generating a mesh it is assumed that polygons do not cross.
3866Having polygons that cross can cause bad meshes to be produced or the mesh generation itself may fail.
3867
3868%OLD
3869%The dataflow is: (See data_manager.py and from scenarios)
3870%
3871%
3872%Simulation scenarios
3873%--------------------%
3874%%
3875%
3876%Sub directories contain scrips and derived files for each simulation.
3877%The directory ../source_data contains large source files such as
3878%DEMs provided externally as well as MOST tsunami simulations to be used
3879%as boundary conditions.
3880%
3881%Manual steps are:
3882%  Creation of DEMs from argcview (.asc + .prj)
3883%  Creation of mesh from pmesh (.tsh)
3884%  Creation of tsunami simulations from MOST (.nc)
3885%%
3886%
3887%Typical scripted steps are%
3888%
3889%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
3890%                   native dem and pts formats%
3891%
3892%  prepare_pts.py: Convert netcdf output from MOST to an SWW file suitable
3893%                  as boundary condition%
3894%
3895%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
3896%                   smoothing. The outputs are tsh files with elevation data.%
3897%
3898%  run_simulation.py: Use the above together with various parameters to
3899%                     run inundation simulation.
3900
3901%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3902
3903\appendix
3904
3905
3906\chapter{Supporting Tools}
3907\label{ch:supportingtools}
3908
3909This section describes a number of supporting tools, supplied with \anuga, that offer a
3910variety of types of functionality and can enhance the basic capabilities of \anuga.
3911
3912
3913\section{caching}
3914\label{sec:caching}
3915
3916The \code{cache} function is used to provide supervised caching of function
3917results. A Python function call of the form:
3918
3919\begin{verbatim}
3920result = func(arg1, ..., argn)
3921\end{verbatim}
3922
3923can be replaced by:
3924
3925\begin{verbatim}
3926from caching import cache
3927result = cache(func,(arg1, ..., argn))
3928\end{verbatim}
3929
3930which returns the same output but reuses cached
3931results if the function has been computed previously in the same context.
3932\code{result} and the arguments can be simple types, tuples, list, dictionaries or
3933objects, but not unhashable types such as functions or open file objects.
3934The function \code{func} may be a member function of an object or a module.
3935
3936This type of caching is particularly useful for computationally intensive
3937functions with few frequently used combinations of input arguments. Note that
3938if the inputs or output are very large caching may not save time because
3939disc access may dominate the execution time.
3940
3941If the function definition changes after a result has been cached, this will be
3942detected by examining the functions bytecode and the function will be recomputed.
3943However, caching will not detect changes in modules used by \code{func}.
3944In this case the cache must be cleared manually.
3945
3946Options are set by means of the function \code{set_option(key, value)},
3947where \code{key} is a key associated with a
3948Python dictionary \code{options}. This dictionary stores settings such as the name of
3949the directory used, the maximum number of cached files allowed, and so on.
3950
3951The \code{cache} function allows the user also to specify a list of dependent files. If any of these
3952have been changed, the function is recomputed and the results stored again.
3953
3954%Other features include support for compression and a capability to \ldots
3955
3956USAGE: \nopagebreak
3957
3958\begin{verbatim}
3959result = cache(func, args, kwargs, dependencies, cachedir, verbose,
3960               compression, evaluate, test, return_filename)
3961\end{verbatim}
3962
3963
3964
3965\pagebreak
3966\section{anuga\_viewer}
3967\label{sec:animate}
3968
3969The output generated by \anuga may be viewed by
3970means of the visualisation tool \code{anuga\_viewer}, which takes an
3971SWW file generated by \anuga and creates a visual representation
3972of the data. Examples may be seen in Figures \ref{fig:runupstart}
3973and \ref{fig:runup2}. To view an SWW file with
3974\code{anuga\_viewer} in the Windows environment you have to run it in a command line as in
3975\begin{verbatim}
3976anuga_viewer <swwfile>
3977\end{verbatim} 
3978
3979or if a georeferenced tif file (or jpg cut to the same shape as the domain) is needed
3980
3981\begin{verbatim} 
3982anuga_viewer <swwfile> <tif file>
3983\end{verbatim} 
3984
3985%, you can simply drag the
3986%icon representing the file over an icon on the desktop for the
3987%\code{anuga\_viewer} executable file (or a shortcut to it), or set up a
3988%file association to make files with the extension \code{.sww} open
3989%with \code{animate}. Alternatively, you can operate \code{animate}
3990%from the command line.
3991
3992
3993
3994
3995%% \pagebreak
3996%% \section{\anuga viewer -- animate}
3997%% \label{sec:animate}
3998
3999%% The output generated by \anuga may be viewed by
4000%% means of the visualisation tool \code{animate}, which takes an
4001%% SWW file generated by \anuga and creates a visual representation
4002%% of the data. Examples may be seen in Figures \ref{fig:runupstart}
4003%% and \ref{fig:runup2}. To view an SWW file with
4004%% \code{animate} in the Windows environment, you can simply drag the
4005%% icon representing the file over an icon on the desktop for the
4006%% \code{animate} executable file (or a shortcut to it), or set up a
4007%% file association to make files with the extension \code{.sww} open
4008%% with \code{animate}. Alternatively, you can operate \code{animate}
4009%% from the command line.
4010
4011%% Upon starting the viewer, you will see an interactive moving-picture
4012%% display. You can use the keyboard and mouse to slow down, speed up or
4013%% stop the display, change the viewing position or carry out a number
4014%% of other simple operations. Help is also displayed when you press
4015%% the \code{h} key.
4016
4017%% The main keys operating the interactive screen are:
4018%% \begin{center}
4019%% \begin{tabular}{|ll|}   \hline
4020%%   \code{h} & toggle on-screen help display \\
4021%%   \code{w} & toggle wireframe \\
4022%%   space bar & start/stop\\
4023%%   up/down arrows & increase/decrease speed\\
4024%%   left/right arrows & direction in time \emph{(when running)}\\
4025%%                     & step through simulation \emph{(when stopped)}\\
4026%%   left mouse button & rotate\\
4027%%   middle mouse button & pan\\
4028%%   right mouse button & zoom\\  \hline
4029%% \end{tabular}
4030%% \end{center}
4031
4032%% % \vfill
4033
4034%% The following describes how to operate \code{animate} from the command line:
4035
4036%% Usage: \code{animate [options] swwfile \ldots}\\ \nopagebreak
4037%% where: \\ \nopagebreak
4038%% \begin{tabular}{ll}
4039%%   \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\
4040%%                                     & \code{HEAD\_MOUNTED\_DISPLAY}\\
4041%%   \code{--rgba} & Request a RGBA colour buffer visual\\
4042%%   \code{--stencil} & Request a stencil buffer visual\\
4043%%   \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
4044%%                                     & overridden by environmental variable\\
4045%%   \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD\_BUFFER | HORIZONTAL\_SPLIT |}\\
4046%%                                     & \code{VERTICAL\_SPLIT | LEFT\_EYE | RIGHT\_EYE |}\\
4047%%                                      & \code{ON | OFF} \\
4048%%   \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
4049%%   \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
4050%%   \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
4051%%   \code{-help} & Display this information\\
4052%%   \code{-hmax <float>} & Height above which transparency is set to
4053%%                                      \code{alphamax}\\
4054%%   \code{-hmin <float>} & Height below which transparency is set to
4055%%                                      zero\\
4056%%   \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
4057%%                                      up, default is overhead)\\
4058%%   \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
4059%%   \code{-movie <dirname>} & Save numbered images to named directory and
4060%%                                      quit\\
4061%%   \code{-nosky} & Omit background sky\\
4062%%   \code{-scale <float>} & Vertical scale factor\\
4063%%   \code{-texture <file>} & Image to use for bedslope topography\\
4064%%   \code{-tps <rate>} & Timesteps per second\\
4065%%   \code{-version} & Revision number and creation (not compile)
4066%%                                      date\\
4067%% \end{tabular}
4068
4069
4070\pagebreak
4071\section{utilities/polygons}
4072
4073\declaremodule{standard}{utilities.polygon}
4074\refmodindex{utilities.polygon}
4075
4076\begin{classdesc}{<callable> = Polygon_function}{regions,
4077                                    default=0.0,
4078                                    geo_reference=None}
4079Module: \code{utilities.polygon}
4080
4081Creates a callable object that returns one of a specified list of values when
4082evaluated at a point \code{(x, y)}, depending on which polygon, from a specified list of polygons, the
4083point belongs to.
4084
4085\code{regions} is a list of pairs
4086\code{(P, v)}, where each \code{P} is a polygon and each \code{v}
4087is either a constant value or a function of coordinates \code{x}
4088and \code{y}, specifying the return value for a point inside \code{P}.
4089
4090\code{default} may be used to specify a value (or a function)
4091for a point not lying inside any of the specified polygons.
4092
4093When a point lies in more than one polygon, the return value is taken to
4094be the value for whichever of these polygon appears later in the list.
4095
4096%FIXME (Howard): CAN x, y BE VECTORS?
4097
4098\code{geo_reference} refers to the status of points
4099that are passed into the function. Typically they will be relative to
4100some origin.
4101
4102Typical usage may take the form:
4103
4104\begin{verbatim}
4105set_quantity('elevation',
4106             Polygon_function([(P1, v1), (P2, v2)],
4107                              default=v3,
4108                              geo_reference=domain.geo_reference))
4109\end{verbatim}
4110\end{classdesc}
4111
4112\begin{funcdesc}{<polygon> = read_polygon}{filename, split=','}
4113Module: \code{utilities.polygon}
4114
4115Reads the specified file and returns a polygon.
4116Each line of the file must contain exactly two numbers, separated by a delimiter, which are interpreted
4117as coordinates of one vertex of the polygon.
4118
4119\code{filename} is the path to the file to read.
4120
4121\code{split} sets the delimiter character between the numbers on one line of
4122the file.  If not specified, the delimiter is the ',' character.
4123\end{funcdesc}
4124
4125\label{ref:function_populate_polygon}
4126\begin{funcdesc}{populate_polygon}{polygon, number_of_points, seed=None, exclude=None}
4127Module: \code{utilities.polygon}
4128
4129Populates the interior of the specified polygon with the specified number of points,
4130selected by means of a uniform distribution function.
4131
4132\code{polygon} is the polygon to populate.
4133
4134\code{number_of_points} is the (optional) number of points.
4135
4136\code{seed}is the optional seed for random number generator.
4137
4138\code{exclude} is a list of polygons (inside main polygon) from where points should be excluded.
4139\end{funcdesc}
4140
4141\label{ref:function_point_in_polygon}
4142\begin{funcdesc}{<point> = point_in_polygon}{polygon, delta=1e-8}
4143Module: \code{utilities.polygon}
4144
4145Returns a point inside the specified polygon and close to the edge. The distance between
4146the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
4147second argument \code{delta}, which is taken as $10^{-8}$ by default.
4148\end{funcdesc}
4149
4150\label{ref:function_inside_polygon}
4151\begin{funcdesc}{<array> = inside_polygon}{points, polygon, closed=True, verbose=False}
4152Module: \code{utilities.polygon}
4153
4154Get a set of points that lie inside a given polygon.
4155
4156\code{points} is the list of points to test.
4157
4158\code{polygon} is the polygon to test the points against.
4159
4160\code{closed} specifies if points on the polygon edge are considered to be inside
4161or outside the polygon -- \code{True} means they are inside.
4162
4163Returns a numeric array comprising the indices of the points in the list
4164that lie inside the polygon.  If none of the points are inside, returns
4165\code{zeros((0,), 'l') (ie, an empty numeric array)}.
4166
4167Compare with \code{outside_polygon()}, page \pageref{ref:function_outside_polygon}.
4168\end{funcdesc}
4169
4170\label{ref:function_outside_polygon}
4171\begin{funcdesc}{<array> = outside_polygon}{points, polygon, closed=True, verbose=False}
4172Module: \code{utilities.polygon}
4173
4174Get a set of points that lie outside a given polygon.
4175
4176\code{points} is the list of points to test.
4177
4178\code{polygon} is the polygon to test the points against.
4179
4180\code{closed} specifies if points on the polygon edge are considered to be outside
4181or inside the polygon -- \code{True} means they are outside.
4182
4183Returns a numeric array comprising the indices of the points in the list
4184that lie outside the polygon.  If none of the points are outside, returns
4185\code{zeros((0,), 'l')} (ie, an empty numeric array).
4186
4187Compare with \code{inside_polygon()}, page \pageref{ref:function_inside_polygon}.
4188\end{funcdesc}
4189
4190\label{ref:function_is_inside_polygon}
4191\begin{funcdesc}{<boolean> = is_inside_polygon}{point, polygon, closed=True, verbose=False}
4192Module: \code{utilities.polygon}
4193
4194Determines if a single point is inside a polygon.
4195
4196\code{point} is the point to test.
4197
4198\code{polygon} is the polygon to test \code{point} against.
4199
4200\code{closed} is a flag that forces the function to consider a point on the polygon
4201edge to be inside or outside -- if \code{True} a point on the edge is considered inside the
4202polygon.
4203
4204Returns \code{True} if \code{point} is inside \code{polygon}.
4205
4206Compare with \code{inside_polygon()}, page \pageref{ref:function_inside_polygon}.
4207\end{funcdesc}
4208
4209\label{ref:function_is_outside_polygon}
4210\begin{funcdesc}{<boolean> = is_outside_polygon}{point, polygon, closed=True, verbose=False,
4211%                                                 points_geo_ref=None, polygon_geo_ref=None
4212                                                }
4213Module: \code{utilities.polygon}
4214
4215Determines if a single point is outside a polygon.
4216
4217\code{point} is the point to test.
4218
4219\code{polygon} is the polygon to test \code{point} against.
4220
4221\code{closed}  is a flag that forces the function to consider a point on the polygon
4222edge to be outside or inside -- if \code{True} a point on the edge is considered inside the
4223polygon.
4224
4225%\code{points_geo_ref} is ??
4226%
4227%\code{polygon_geo_ref} is ??
4228
4229Compare with \code{outside_polygon()}, page \pageref{ref:function_outside_polygon}.
4230\end{funcdesc}
4231
4232\label{ref:function_point_on_line}
4233\begin{funcdesc}{<boolean> = point_on_line}{point, line, rtol=1.0e-5, atol=1.0e-8}
4234Module: \code{utilities.polygon}
4235
4236Determine if a point is on a line to some tolerance.  The line is considered to
4237extend past its end-points.
4238
4239\code{point} is the point to test ([$x$, $y$]).
4240
4241\code{line} is the line to test \code{point} against ([[$x1$,$y1$], [$x2$,$y2$]]).
4242
4243\code{rtol} is the relative tolerance to use when testing for coincidence.
4244
4245\code{atol} is the absolute tolerance to use when testing for coincidence.
4246
4247Returns \code{True} if the point is on the line, else \code{False}.
4248\end{funcdesc}
4249
4250\label{ref:function_separate_points_by_polygon}
4251\begin{funcdesc}{(indices, count) = separate_points_by_polygon}{points, polygon,
4252                                             closed=True,
4253                                             check_input=True,
4254                                             verbose=False}
4255\indexedcode{separate_points_by_polygon}
4256Module: \code{utilities.polygon}
4257
4258Separate a set of points into points that are inside and outside a polygon.
4259
4260\code{points} is a list of points to separate.
4261
4262\code{polygon} is the polygon used to separate the points.
4263
4264\code{closed} determines whether points on the polygon edge should be
4265regarded as inside or outside the polygon.  \code{True} means they are inside.
4266
4267\code{check_input} specifies whether the input parameters are checked -- \code{True}
4268means check the input parameters.
4269
4270The function returns a tuple \code{(indices, count)} where \code{indices} is a list of
4271point $indices$ from the input \code{points} list, with the indices of points inside the
4272polygon at the left and indices of points outside the polygon listed at the right.  The
4273\code{count} value is the count (from the left) of the indices of the points $inside$ the
4274polygon.
4275\end{funcdesc}
4276
4277\begin{funcdesc}{<area> = polygon_area}{polygon}
4278Module: \code{utilities.polygon}
4279
4280Returns area of an arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html).
4281\end{funcdesc}
4282
4283\begin{funcdesc}{[$x_{min}$, $x_{max}$, $y_{min}$, $y_{max}$] = plot_polygons}
4284                               {polygons_points, style=None,
4285                                figname=None, label=None, verbose=False}
4286Module: \code{utilities.polygon}
4287
4288Plot a list of polygons to a file.
4289
4290\code{polygons_points} is a list of polygons to plot.
4291
4292\code{style} is a list of style strings to be applied to the corresponding polygon
4293in \code{polygons_points}. A polygon can be closed for plotting purposes by assigning
4294the style string 'line' to it in the appropriate place in the \code{style} list.
4295The default style is 'line'.
4296
4297\code{figname} is the path to the file to save the plot in.  If not specified, use
4298\file{test_image.png}.
4299
4300The function returns a list containing the minimum and maximum of the points in all the
4301input polygons, i.e.\ \code{[$x_{min}$, $x_{max}$, $y_{min}$, $y_{max}$]}.
4302\end{funcdesc}
4303
4304
4305\pagebreak
4306\section{coordinate_transforms}
4307
4308
4309\pagebreak
4310\section{geospatial_data}
4311\label{sec:geospatial}
4312
4313This describes a class that represents arbitrary point data in UTM
4314coordinates along with named attribute values.
4315
4316%FIXME (Ole): This gives a LaTeX error
4317%\declaremodule{standard}{geospatial_data}
4318%\refmodindex{geospatial_data}
4319
4320\begin{classdesc}{Geospatial_data}
4321                {data_points=None,
4322                 attributes=None,
4323                 geo_reference=None,
4324                 default_attribute_name=None,
4325                 file_name=None,
4326                 latitudes=None,
4327                 longitudes=None,
4328                 points_are_lats_longs=False,
4329                 max_read_lines=None,
4330                 load_file_now=True,
4331                 verbose=False}
4332Module: \code{geospatial_data.geospatial_data}
4333
4334This class is used to store a set of data points and associated
4335attributes, allowing these to be manipulated by methods defined for
4336the class.
4337
4338The data points are specified either by reading them from a NetCDF
4339or CSV file, identified through the parameter \code{file_name}, or
4340by providing their \code{x}- and \code{y}-coordinates in metres,
4341either as a sequence of 2-tuples of floats or as an $M \times 2$
4342numeric array of floats, where $M$ is the number of points.
4343
4344Coordinates are interpreted relative to the origin specified by the
4345object \code{geo_reference}, which contains data indicating the UTM
4346zone, easting and northing. If \code{geo_reference} is not
4347specified, a default is used.
4348
4349Attributes are specified through the parameter \code{attributes},
4350set either to a list or array of length $M$ or to a dictionary whose
4351keys are the attribute names and whose values are lists or arrays of
4352length $M$. One of the attributes may be specified as the default
4353attribute, by assigning its name to \code{default_attribute_name}.
4354If no value is specified, the default attribute is taken to be the
4355first one.
4356
4357Note that the \code{Geospatial_data} object currently reads entire datasets
4358into memory i.e.\ no memomry blocking takes place.
4359For this we refer to the \code{set_quantity()} method which will read PTS and CSV
4360files into \anuga using memory blocking allowing large files to be used.
4361\end{classdesc}
4362
4363\begin{methoddesc}{\emph{<Geospatial_data>}.import_points_file}
4364        {file_name, delimiter=None, verbose=False}
4365Module: \code{geospatial_data.geospatial_data}
4366
4367Import a TXT, CSV or PTS points data file into a code{Geospatial_data} object.
4368
4369\code{file_name} is the path to a TXT, CSV or PTS points data file.
4370
4371\code{delimiter} is currently unused.
4372\end{methoddesc}
4373
4374\begin{methoddesc}{\emph{<Geospatial_data>}.export_points_file}{file_name, absolute=True,
4375                                       as_lat_long=False, isSouthHemisphere=True}
4376Module: \code{geospatial_data.geospatial_data}
4377
4378Export a CSV or PTS points data file from a \code{Geospatial_data} object.
4379
4380\code{file_name} is the path to the CSV or PTS points file to write.
4381
4382\code{absolute} determines if the exported data is absolute or relative to the
4383\code{Geospatial_data} object geo_reference.  If \code{True} the exported
4384data is absolute.
4385
4386\code{as_lat_long} exports the points data as latitudes and longitudes if \code{True}.
4387
4388\code{isSouthHemisphere} has effect only if \code{as_lat_long} is \code{True} and causes
4389latitude/longitude values to be for the southern (\code{True}) or northern hemispheres
4390(\code{False}).
4391\end{methoddesc}
4392
4393\begin{methoddesc}{points = \emph{<Geospatial_data>}.get_data_points}
4394        {absolute=True, geo_reference=None,
4395         as_lat_long=False, isSouthHemisphere=True}
4396Module: \code{geospatial_data.geospatial_data}
4397
4398Get the coordinates for all the data points as an $N \times 2$ array.
4399
4400\code{absolute} determines if the exported data is absolute or relative to the
4401\code{Geospatial_data} object geo_reference.  If \code{True} the exported
4402data is absolute.
4403
4404\code{geo_reference} is the geo_reference the points are relative to, if supplied.
4405
4406\code{as_lat_long} exports the points data as latitudes and longitudes if \code{True}.
4407
4408\code{isSouthHemisphere} has effect only if \code{as_lat_long} is \code{True} and causes
4409latitude/longitude values to be for the southern (\code{True}) or northern hemispheres
4410(\code{False}).
4411\end{methoddesc}
4412
4413\begin{methoddesc}{\emph{<Geospatial_data>}.set_attributes}{attributes}
4414Module: \code{geospatial_data.geospatial_data}
4415
4416Set the attributes for a \code{Geospatial_data} object.
4417
4418\code{attributes} is the new value for the object's attributes.  May be a dictionary or \code{None}.
4419\end{methoddesc}
4420
4421\begin{methoddesc}{atributes = \emph{<Geospatial_data>}.get_attributes}{attribute_name=None}
4422Module: \code{geospatial_data.geospatial_data}
4423
4424Get a named attribute from a \code{Geospatial_data} object.
4425
4426\code{attribute_name} is the name of the desired attribute.  If \code{None}, return
4427the default attribute.
4428\end{methoddesc}
4429
4430\begin{methoddesc}{\emph{<Geospatial_data>}.get_all_attributes}{}
4431Module: \code{geospatial_data.geospatial_data}
4432
4433Get all attributes of a \code{Geospatial_data} object.
4434
4435Returns \code{None} or the attributes dictionary (which may be empty).
4436\end{methoddesc}
4437
4438\begin{methoddesc}{\emph{<Geospatial_data>}.set_default_attribute_name}{default_attribute_name}
4439Module: \code{geospatial_data.geospatial_data}
4440
4441Set the default attribute name of a \code{Geospatial_data} object.
4442
4443\code{default_attribute_name} is the new default attribute name.
4444\end{methoddesc}
4445
4446\begin{methoddesc}{\emph{<Geospatial_data>}.set_geo_reference}{geo_reference}
4447Module: \code{geospatial_data.geospatial_data}
4448
4449Set the internal geo_reference of a \code{Geospatial_data} object.
4450
4451\code{geo_reference} is the new internal geo_reference for the object.
4452If \code{None} will use the default geo_reference.
4453
4454If the \code{Geospatial_data} object already has an internal geo_reference
4455then the points data will be changed to use the new geo_reference.
4456\end{methoddesc}
4457
4458\begin{methoddesc}{\emph{<Geospatial_data>}.__add__}{other}
4459Module: \code{geospatial_data.geospatial_data}
4460
4461The \code{__add__()} method is defined so it is possible to add two
4462\code{Geospatial_data} objects.
4463\end{methoddesc}
4464
4465\label{ref:function_clip}
4466\begin{methoddesc}{geospatial = \emph{<Geospatial_data>}.clip}{polygon, closed=True, verbose=False}
4467Module: \code{geospatial_data.geospatial_data}
4468
4469Clip a \code{Geospatial_data} object with a polygon.
4470
4471\code{polygon} is the polygon to clip the \code{Geospatial_data} object with.
4472This may be a list of points, an $N \times 2$ array or a \code{Geospatial_data}
4473object.
4474
4475\code{closed} determines whether points on the \code{polygon} edge are inside (\code{True})
4476or outside (\code{False}) the polygon.
4477
4478Returns a new \code{Geospatial_data} object representing points inside the
4479
4480Compare with \code{clip_outside()}, page \pageref{ref:function_clip_outside}.
4481specified polygon.
4482\end{methoddesc}
4483
4484\label{ref:function_clip_outside}
4485\begin{methoddesc}{geospatial = \emph{<Geospatial_data>}.clip_outside}
4486        {polygon, closed=True, verbose=False}
4487Module: \code{geospatial_data.geospatial_data}
4488
4489Clip a \code{Geospatial_data} object with a polygon.
4490
4491\code{polygon} is the polygon to clip the \code{Geospatial_data} object with.
4492This may be a list of points, an $N \times 2$ array or a \code{Geospatial_data}
4493object.
4494
4495\code{closed} determines whether points on the \code{polygon} edge are inside (\code{True})
4496or outside (\code{False}) the polygon.
4497
4498Returns a new \code{Geospatial_data} object representing points outside the
4499specified polygon.
4500
4501Compare with \code{clip()}, page \pageref{ref:function_clip}.
4502\end{methoddesc}
4503
4504\begin{methoddesc}{(g1, g2) = \emph{<Geospatial_data>}.split}
4505        {factor=0.5, seed_num=None, verbose=False}
4506Module: \code{geospatial_data.geospatial_data}
4507
4508Split a \code{Geospatial_data} object into two objects of predetermined ratios.
4509
4510\code{factor} is the ratio of the size of the first returned object to the
4511original object.  If '0.5' is supplied, the two resulting objects will be
4512of equal size.
4513
4514\code{seed_num}, if supplied, will be the random number generator seed used for
4515the split.
4516
4517Points of the two new geospatial_data object are selected RANDOMLY.
4518
4519Returns two geospatial_data objects that are disjoint sets of the original.
4520\end{methoddesc}
4521
4522\subsection{Miscellaneous Functions}
4523
4524The functions here are not \code{Geospatial_data} object methods, but are used with them.
4525
4526\begin{methoddesc}{X = find_optimal_smoothing_parameter}
4527                  {data_file,
4528                   alpha_list=None,
4529                   mesh_file=None,
4530                   boundary_poly=None,
4531                   mesh_resolution=100000,
4532                   north_boundary=None,
4533                   south_boundary=None,
4534                   east_boundary=None,
4535                   west_boundary=None,
4536                   plot_name='all_alphas',
4537                   split_factor=0.1,
4538                   seed_num=None,
4539                   cache=False,
4540                   verbose=False}
4541Module: \code{geospatial_data.geospatial_data}
4542
4543Calculate the minimum covariance from a set of points in a file.  It does this
4544by removing a small random sample of points from \code{data_file} and creating
4545models with different alpha values from \code{alpha_list} and cross validates
4546the predicted value to the previously removed point data. Returns the
4547alpha value which has the smallest covariance.
4548
4549\code{data_file} is the input data file and must not contain points outside
4550the boundaries defined and is either a PTS, TXT or CSV file.
4551
4552\code{alpha_list} is the list of alpha values to use.
4553
4554\code{mesh_file} is the path to a mesh file to create (if supplied).
4555If \code{None} a mesh file will be created (named \file{temp.msh}).
4556NOTE: if there is a \code{mesh_resolution} defined or any boundaries are defined,
4557any input \code{mesh_file} value is ignored.
4558
4559\code{mesh_resolution} is the maximum area size for a triangle.
4560
4561\code{north_boundary}\\
4562\code{south_boundary}\\
4563\code{east_boundary}\\
4564\code{west_boundary} are the boundary values to use.
4565
4566\code{plot_name} is the path name of the plot file to write.
4567
4568\code{seed_num} is the random number generator seed to use.
4569
4570The function returns a tuple \code{(min_covar, alpha)} where \code{min_covar} is
4571the minumum normalised covariance and \code{alpha} is the alpha value that
4572created it.  A plot file is also written.
4573
4574This is an example of function usage: \nopagebreak
4575
4576\begin{verbatim}
4577convariance_value, alpha = $\backslash$
4578        find_optimal_smoothing_parameter(data_file=fileName,
4579                                         alpha_list=[0.0001, 0.01, 1],
4580                                         mesh_file=None,
4581                                         mesh_resolution=3,
4582                                         north_boundary=5,
4583                                         south_boundary=-5,
4584                                         east_boundary=5,
4585                                         west_boundary=-5,
4586                                         plot_name='all_alphas',
4587                                         seed_num=100000,
4588                                         verbose=False)
4589\end{verbatim}
4590
4591NOTE: The function will not work if the \code{data_file} extent is greater than the
4592\code{boundary_poly} polygon or any of the boundaries, e.g.\ \code{north_boundary}, etc.
4593\end{methoddesc}
4594
4595
4596\pagebreak
4597\section{Graphical Mesh Generator GUI}
4598The program \code{graphical_mesh_generator.py} in the \code{pmesh} module
4599allows the user to set up the mesh of the problem interactively.
4600It can be used to build the outline of a mesh or to visualise a mesh
4601automatically generated.
4602
4603Graphical Mesh Generator will let the user select various modes. The
4604current allowable modes are $vertex$, $segment$, $hole$ or $region$.  The mode
4605describes what sort of object is added or selected in response to
4606mouse clicks.  When changing modes any prior selected objects become
4607deselected.
4608
4609In general the left mouse button will add an object and the right
4610mouse button will select an object.  A selected object can de deleted
4611by pressing the the middle mouse button (scroll bar).
4612
4613
4614\pagebreak
4615\section{class Alpha_Shape}
4616\emph{Alpha shapes} are used to generate close-fitting boundaries
4617around sets of points. The alpha shape algorithm produces a shape
4618that approximates to the 'shape formed by the points' -- or the shape
4619that would be seen by viewing the points from a coarse enough
4620resolution. For the simplest types of point sets, the alpha shape
4621reduces to the more precise notion of the convex hull. However, for
4622many sets of points the convex hull does not provide a close fit and
4623the alpha shape usually fits more closely to the original point set,
4624offering a better approximation to the shape being sought.
4625
4626In \anuga, an alpha shape is used to generate a polygonal boundary
4627around a set of points before mesh generation. The algorithm uses a
4628parameter $\alpha$ that can be adjusted to make the resultant shape
4629resemble the shape suggested by intuition more closely. An alpha
4630shape can serve as an initial boundary approximation that the user
4631can adjust as needed.
4632
4633The following paragraphs describe the class used to model an alpha
4634shape and some of the important methods and attributes associated
4635with instances of this class.
4636
4637\label{class:alpha_shape}
4638\begin{classdesc}{Alpha_Shape}{points, alpha=None}
4639Module: \code{alpha_shape}
4640
4641Instantiate an instance of the \code{Alpha_Shape} class.
4642
4643\code{points} is an $N \times 2$ list of points (\code{[[x1, y1],[x2, y2]}\ldots\code{]}).
4644
4645\code{alpha} is the 'fitting' parameter.
4646\end{classdesc}
4647
4648\begin{funcdesc}{alpha_shape_via_files}{point_file, boundary_file, alpha= None}
4649Module: \code{alpha_shape}
4650
4651This function reads points from the specified point file
4652\code{point_file}, computes the associated alpha shape (either
4653using the specified value for \code{alpha} or, if no value is
4654specified, automatically setting it to an optimal value) and outputs
4655the boundary to a file named \code{boundary_file}. This output file
4656lists the coordinates \code{(x, y)} of each point in the boundary,
4657using one line per point.
4658\end{funcdesc}
4659
4660\label{ref:method_set_boundary_type}
4661\begin{methoddesc}{\emph{<Alpha_shape>}.set_boundary_type}
4662        {raw_boundary=True,
4663         remove_holes=False,
4664         smooth_indents=False,
4665         expand_pinch=False,
4666         boundary_points_fraction=0.2}
4667Module: \code{alpha_shape}
4668
4669This function sets internal state that controls how the \code{Alpha_shape}
4670boundary is presented or exported.
4671
4672\code{raw_boundary} sets the type to $raw$ if \code{True},
4673i.e.\ the regular edges of the alpha shape.
4674
4675\code{remove_holes}, if \code{True} removes small holes ('small' is defined by
4676\code{boundary_points_fraction}).
4677
4678\code{smooth_indents}, if \code{True} removes sharp triangular indents in
4679the boundary.
4680
4681\code{expand_pinch}, if \code{True} tests for pinch-off and
4682corrects -- preventing a boundary vertex from having more than two edges.
4683\end{methoddesc}
4684
4685\label{ref:method_get_boundary}
4686\begin{methoddesc}{boundary = \emph{<Alpha_shape>}.get_boundary}{}
4687Module: \code{alpha_shape}
4688
4689Returns a list of tuples representing the boundary of the alpha
4690shape. Each tuple represents a segment in the boundary by providing
4691the indices of its two endpoints.
4692
4693See \code{set_boundary_type()}, page \pageref{ref:method_set_boundary_type}.
4694\end{methoddesc}
4695
4696\label{ref:method_write_boundary}
4697\begin{methoddesc}{\emph{<Alpha_shape>}.write_boundary}{file_name}
4698Module: \code{alpha_shape}
4699
4700Writes the list of 2-tuples returned by \code{get_boundary()} to the
4701file \code{file_name}, using one line per tuple.
4702
4703See \code{set_boundary_type()}, page \pageref{ref:method_set_boundary_type}. \\
4704See \code{get_boundary()}, page \pageref{ref:method_get_boundary}.
4705\end{methoddesc}
4706
4707
4708\pagebreak
4709\section{Numerical Tools}
4710
4711The following table describes some useful numerical functions that
4712may be found in the module \module{utilities.numerical\_tools}:
4713
4714\begin{tabular}{|p{7.4cm} p{8.6cm}|}
4715  \hline
4716  \code{angle(v1, v2=None)} & Angle between two-dimensional vectors
4717                              \code{v1} and \code{v2}, or between \code{v1} and the $x$-axis if
4718                              \code{v2} is \code{None}. Value is in range $0$ to $2\pi$. \\
4719  & \\
4720  \code{normal\_vector(v)} & Normal vector to \code{v}.\\
4721  & \\
4722  \code{mean(x)} & Mean value of a vector \code{x}.\\
4723  & \\
4724  \code{cov(x, y=None)} & Covariance of vectors \code{x} and \code{y}.
4725                          If \code{y} is \code{None}, returns \code{cov(x, x)}.\\
4726  & \\
4727  \code{err(x, y=0, n=2, relative=True)} & Relative error of $\parallel$\code{x}$-$\code{y}$\parallel$
4728                                           to $\parallel$\code{y}$\parallel$ (2-norm if \code{n} = 2 or
4729                                           Max norm if \code{n} = \code{None}). If denominator evaluates
4730                                           to zero or if \code{y} is omitted or if \code{relative=False},
4731                                           absolute error is returned.\\
4732  & \\
4733  \code{norm(x)} & 2-norm of \code{x}.\\
4734  & \\
4735  \code{corr(x, y=None)} & Correlation of \code{x} and \code{y}. If
4736                           \code{y} is \code{None} returns autocorrelation of \code{x}.\\
4737  & \\
4738  \code{ensure\_numeric(A, typecode=None)} & Returns a numeric array for any sequence \code{A}. If
4739                                             \code{A} is already a numeric array it will be returned
4740                                             unaltered. Otherwise, an attempt is made to convert
4741                                             it to a numeric array. (Needed because \code{array(A)} can
4742                                             cause memory overflow.)\\
4743  & \\
4744  \code{histogram(a, bins, relative=False)} & Standard histogram. If \code{relative} is \code{True},
4745                                              values will be normalised against the total and thus
4746                                              represent frequencies rather than counts.\\
4747  & \\
4748  \code{create\_bins(data, number\_of\_bins=None)} & Safely create bins for use with histogram.
4749                                                     If \code{data} contains only one point
4750                                                     or is constant, one bin will be created.
4751                                                     If \code{number\_of\_bins} is omitted,
4752                                                     10 bins will be created.\\
4753  \hline
4754\end{tabular}
4755
4756
4757\section{Finding the Optimal Alpha Value}
4758
4759The function ????
4760more to come very soon
4761
4762%?% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4763%?%
4764%?% \chapter{Modules available in \anuga}
4765%?%
4766%?%
4767%?% %abstract_2d_finite_volumes
4768%?%
4769%?% \section{\module{abstract_2d_finite_volumes.domain}}
4770%?% Generic module for 2D triangular domains for finite-volume computations of conservation laws
4771%?% \declaremodule[domain]{}{domain}
4772%?% \label{mod:domain}
4773%?%
4774%?%
4775%?% \section{\module{abstract_2d_finite_volumes.ermapper_grids}}
4776%?% \declaremodule[ermappergrids]{}{ermapper_grids}
4777%?% \label{mod:ermapper_grids}
4778%?%
4779%?%
4780%?% \section{\module{abstract_2d_finite_volumes.general_mesh} }
4781%?% \declaremodule[generalmesh]{}{general_mesh}
4782%?% \label{mod:general_mesh}
4783%?%
4784%?%
4785%?% \section{\module{abstract_2d_finite_volumes.generic_boundary_conditions} }
4786%?% \declaremodule[genericboundaryconditions]{}{generic_boundary_conditions}
4787%?% \label{mod:generic_boundary_conditions}
4788%?%
4789%?%
4790%?% \section{\module{abstract_2d_finite_volumes.mesh_factory} }
4791%?% \declaremodule[meshfactory]{}{mesh_factory}
4792%?% \label{mod:mesh_factory}
4793%?%
4794%?%
4795%?% \section{\module{abstract_2d_finite_volumes.mesh_factory} }
4796%?% \declaremodule[meshfactory]{}{mesh_factory}
4797%?% \label{mod:mesh_factory}
4798%?%
4799%?%
4800%?% \section{\module{abstract_2d_finite_volumes.neighbour_mesh} }
4801%?% \declaremodule[neighbourmesh]{}{neighbour_mesh}
4802%?% \label{mod:neighbour_mesh}
4803%?%
4804%?%
4805%?% \section{\module{abstract_2d_finite_volumes.pmesh2domain} }
4806%?% \declaremodule[pmesh2domain]{}{pmesh2domain}
4807%?% \label{mod:pmesh2domain}
4808%?%
4809%?%
4810%?% \section{\module{abstract_2d_finite_volumes.quantity}}
4811%?% \declaremodule{}{quantity}
4812%?% \label{mod:quantity}
4813%?%
4814%?% \begin{verbatim}
4815%?% Class Quantity - Implements values at each triangular element
4816%?%
4817%?% To create:
4818%?%
4819%?%    Quantity(domain, vertex_values)
4820%?%
4821%?%    domain: Associated domain structure. Required.
4822%?%
4823%?%    vertex_values: Nx3 array of values at each vertex for each element.
4824%?%                   Default None
4825%?%
4826%?%    If vertex_values are None Create array of zeros compatible with domain.
4827%?%    Otherwise check that it is compatible with dimenions of domain.
4828%?%    Otherwise raise an exception
4829%?% \end{verbatim}
4830%?%
4831%?%
4832%?% \section{\module{abstract_2d_finite_volumes.region} }
4833%?% \declaremodule[region]{}{region}
4834%?% \label{mod:region}
4835%?%
4836%?%
4837%?% \section{\module{abstract_2d_finite_volumes.util} }
4838%?% \declaremodule[util]{}{util}
4839%?% \label{mod:util}
4840%?%
4841%?%
4842%?% advection
4843%?% alpha_shape
4844%?% caching
4845%?% coordinate_transforms
4846%?% culvert_flows
4847%?% damage_modelling
4848%?% euler
4849%?% fit_interpolate
4850%?% geospatial_data
4851%?% lib
4852%?% load_mesh
4853%?% mesh_engine
4854%?% pmesh
4855%?% SConstruct
4856%?% shallow_water
4857%?% utilities
4858%?%
4859%?%
4860%?% \section{\module{shallow\_water}}
4861%?%
4862%?% 2D triangular domains for finite-volume
4863%?% computations of the shallow water wave equation.
4864%?% This module contains a specialisation of class Domain from module domain.py consisting of methods specific to the Shallow Water
4865%?% Wave Equation
4866%?%
4867%?% \declaremodule[shallowwater]{}{shallow\_water}
4868%?% \label{mod:shallowwater}
4869
4870%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4871
4872\chapter{\anuga Full-scale Validations}
4873
4874
4875\section{Overview}
4876
4877There are some long-running validations that are not included in the small-scale
4878validations that run when you execute the \module{validate_all.py}
4879script in the \module{anuga_validation/automated_validation_test} directory.
4880These validations are not run automatically since they can take a large amount
4881of time and require an internet connection and user input.
4882
4883
4884\section{Patong Beach}
4885
4886The Patong Beach validation is run from the \module{automated_validation_tests/patong_beach_validation}
4887directory.  Just execute the \module{validate_patong.py} script in that directory.
4888This will run a Patong Beach simulation and compare the generated SWW file with a
4889known good copy of that file.
4890
4891The script attempts to refresh the validation data files from master copies held
4892on the Sourceforge project site.  If you don't have an internet connection you
4893may still run the validation, as long as you have the required files.
4894
4895You may download the validation data files by hand and then run the validation.
4896Just go to the \anuga Sourceforge project download page at
4897\module{http://sourceforge.net/project/showfiles.php?group_id=172848} and select
4898the \module{validation_data} package, \module{patong-1.0} release.  You need the
4899\module{data.tgz} file and one or more of the \module{patong.sww.\{BASIC|TRIAL|FINAL\}.tgz}
4900files.
4901
4902The BASIC validation is the quickest and the FINAL validation is the slowest.
4903The \module{validate.py} script will use whatever files you have, BASIC first and
4904FINAL last.
4905
4906%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4907
4908\chapter{Frequently Asked Questions}
4909
4910The Frequently Asked Questions have been move to the online FAQ at: \\
4911\url{https://datamining.anu.edu.au/anuga/wiki/FrequentlyAskedQuestions}
4912
4913%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4914
4915\chapter{Glossary}
4916
4917\begin{tabular}{|lp{10cm}|c|}  \hline
4918  \emph{Term} & \emph{Definition} & \emph{Page}\\
4919  \hline
4920  \indexedbold{\anuga} & Name of software (joint development between ANU and GA) & \pageref{def:anuga}\\
4921  \indexedbold{bathymetry} & offshore elevation & \\
4922  \indexedbold{conserved quantity} & conserved (stage, x and y momentum) & \\
4923%  \indexedbold{domain} & The domain of a function is the set of all input values to the function.& \\
4924  \indexedbold{Digital Elevation Model (DEM)} & DEMs are digital files consisting of points of elevations,
4925                                                sampled systematically at equally spaced intervals.& \\
4926  \indexedbold{Dirichlet boundary} & A boundary condition imposed on a differential equation that specifies
4927                                     the values the solution is to take on the boundary of the domain.
4928                                   & \pageref{def:dirichlet boundary}\\
4929  \indexedbold{edge} & A triangular cell within the computational mesh can be depicted
4930                       as a set of vertices joined by lines (the edges). & \\
4931  \indexedbold{elevation} & refers to bathymetry and topography & \\
4932  \indexedbold{evolution} & integration of the shallow water wave equations over time & \\
4933  \indexedbold{finite volume method} & The method evaluates the terms in the shallow water wave equation as
4934                                       fluxes at the surfaces of each finite volume. Because the flux entering
4935                                       a given volume is identical to that leaving the adjacent volume, these
4936                                       methods are conservative. Another advantage of the finite volume method is
4937                                       that it is easily formulated to allow for unstructured meshes. The method
4938                                       is used in many computational fluid dynamics packages. & \\
4939  \indexedbold{forcing term} & &\\
4940  \indexedbold{flux} & the amount of flow through the volume per unit time & \\
4941  \indexedbold{grid} & Evenly spaced mesh & \\
4942  \indexedbold{latitude} & The angular distance on a mericlear north and south of the equator, expressed in degrees and minutes. & \\
4943  \indexedbold{longitude} & The angular distance east or west, between the meridian of a particular place on Earth
4944                            and that of the Prime Meridian (located in Greenwich, England) expressed in degrees or time.& \\
4945  \indexedbold{Manning friction coefficient} & &\\
4946  \indexedbold{mesh} & Triangulation of domain &\\
4947  \indexedbold{mesh file} & A TSH or MSH file & \pageref{def:mesh file}\\
4948  \indexedbold{NetCDF} & &\\
4949  \indexedbold{node} & A point at which edges meet & \\
4950  \indexedbold{northing} & A rectangular (x,y) coordinate measurement of distance north from a north-south
4951                           reference line, usually a meridian used as the axis of origin within a map zone
4952                           or projection. Northing is a UTM (Universal Transverse Mercator) coordinate. & \\
4953  \indexedbold{points file} & A PTS or CSV file & \\
4954  \hline
4955\end{tabular}
4956
4957\begin{tabular}{|lp{10cm}|c|}
4958  \hline
4959  \indexedbold{polygon} & A sequence of points in the plane. \anuga represents a polygon either as a list
4960                          consisting of Python tuples or lists of length 2 or as an $N \times 2$ numeric array,
4961                          where $N$ is the number of points.
4962
4963                          The unit square, for example, would be represented either as \code{[ [0,0], [1,0], [1,1], [0,1] ]}
4964                          or as \code{array( [0,0], [1,0], [1,1], [0,1] )}.
4965
4966                          NOTE: For details refer to the module \module{utilities/polygon.py}. & \\
4967  \indexedbold{resolution} &  The maximal area of a triangular cell in a mesh & \\
4968  \indexedbold{reflective boundary} & Models a solid wall. Returns same conserved quantities as those present in the
4969                                      neighbouring volume but reflected. Specific to the shallow water equation as
4970                                      it works with the momentum quantities assumed to be the second and third
4971                                      conserved quantities. & \pageref{def:reflective boundary}\\
4972  \indexedbold{stage} & &\\
4973%  \indexedbold{try this}
4974  \indexedbold{animate} & visualisation tool used with \anuga & \pageref{sec:animate}\\
4975  \indexedbold{time boundary} & Returns values for the conserved quantities as a function of time.
4976                                The user must specify the domain to get access to the model time.
4977                              & \pageref{def:time boundary}\\
4978  \indexedbold{topography} & onshore elevation &\\
4979  \indexedbold{transmissive boundary} & & \pageref{def:transmissive boundary}\\
4980  \indexedbold{vertex} & A point at which edges meet. & \\
4981  \indexedbold{xmomentum} & conserved quantity (note, two-dimensional SWW equations say only
4982                            \code{x} and \code{y} and NOT \code{z}) &\\
4983  \indexedbold{ymomentum}  & conserved quantity & \\
4984  \hline
4985\end{tabular}
4986
4987%The \code{\e appendix} markup need not be repeated for additional
4988%appendices.
4989
4990%
4991%  The ugly "%begin{latexonly}" pseudo-environments are really just to
4992%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
4993%  not really valuable.
4994%
4995%  If you don't want the Module Index, you can remove all of this up
4996%  until the second \input line.
4997%
4998
4999%?% %begin{latexonly}
5000%?% %\renewcommand{\indexname}{Module Index}
5001%?% %end{latexonly}
5002%?% \input{mod\jobname.ind}        % Module Index
5003
5004%begin{latexonly}
5005%\renewcommand{\indexname}{Index}
5006%end{latexonly}
5007\input{\jobname.ind}            % Index
5008
5009
5010\begin{thebibliography}{99}
5011%\begin{thebibliography}
5012\bibitem[nielsen2005]{nielsen2005}
5013{\it Hydrodynamic modelling of coastal inundation}.
5014Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman.
5015In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
5016Modelling and Simulation. Modelling and Simulation Society of Australia and
5017New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.\\
5018http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
5019
5020\bibitem[grid250]{grid250}
5021Australian Bathymetry and Topography Grid, June 2005.
5022Webster, M.A. and Petkovic, P.
5023Geoscience Australia Record 2005/12. ISBN: 1-920871-46-2.\\
5024http://www.ga.gov.au/meta/ANZCW0703008022.html
5025
5026\bibitem[ZR1999]{ZR1999}
5027\newblock {Catastrophic Collapse of Water Supply Reservoirs in Urban Areas}.
5028\newblock C.~Zoppou and S.~Roberts.
5029\newblock {\em ASCE J. Hydraulic Engineering}, 125(7):686--695, 1999.
5030
5031\bibitem[Toro1999]{Toro1992}
5032\newblock Riemann problems and the waf method for solving the two-dimensional
5033  shallow water equations.
5034\newblock E.~F. Toro.
5035\newblock {\em Philosophical Transactions of the Royal Society, Series A},
5036  338:43--68, 1992.
5037
5038\bibitem{KurNP2001}
5039\newblock Semidiscrete central-upwind schemes for hyperbolic conservation laws
5040  and hamilton-jacobi equations.
5041\newblock A.~Kurganov, S.~Noelle, and G.~Petrova.
5042\newblock {\em SIAM Journal of Scientific Computing}, 23(3):707--740, 2001.
5043\end{thebibliography}
5044% \end{thebibliography}{99}
5045
5046\end{document}
Note: See TracBrowser for help on using the repository browser.