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

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