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

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

Added the Patong Beach full-scale validation information.

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