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

Last change on this file since 7312 was 7312, checked in by ole, 14 years ago

Clarified uniquely stored vertices

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