source: trunk/anuga_core/documentation/anuga_user_manual.tex @ 7935

Last change on this file since 7935 was 7848, checked in by steve, 14 years ago

Changed the logging levels in log.py so that the information about openning the
file ./anuga.log is now only an info log as opposed to critical log. I.e. by default
doesn't write to the console.

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