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

Last change on this file since 7064 was 7064, checked in by rwilson, 11 years ago

Fiddling with layout of user guide.

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