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

Last change on this file since 7134 was 7134, checked in by rwilson, 15 years ago

Major wide-spread changes to the user manual.

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