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

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

Changed "\anuga" macro to insert just "ANUGA", also added line in introduction about v1.0.

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