source: branches/anuga_1_2_1/anuga_core/documentation/user_manual/anuga_user_manual.tex @ 8083

Last change on this file since 8083 was 8083, checked in by habili, 13 years ago

Bug fixes, including correct use of starttime and duration.

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