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

Last change on this file since 4265 was 4265, checked in by ole, 17 years ago

Move mathematical destcription from linuxconf paper and limiting.tex into user manual.
Have not tried to compile it though.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 140.6 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-02-16 06:26:04 +0000 (Fri, 16 Feb 2007) $
83%$LastChangedRevision: 4265 $
84%$LastChangedBy: ole $
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  \end{funcdesc}
2259
2260
2261  \begin{funcdesc}{boundary\_statistics}{quantities = None, tags = None}
2262  Module: \module{abstract\_2d\_finite\_volumes.domain}
2263
2264  Returns a string of the following type when \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}:
2265
2266  {\small \begin{verbatim}
2267 Boundary values at time 0.5000:
2268    top:
2269        stage in [ -0.25821218,  -0.02499998]
2270    bottom:
2271        stage in [ -0.27098821,  -0.02499974]
2272  \end{verbatim}}
2273
2274  \end{funcdesc}
2275
2276
2277  \begin{funcdesc}{get\_quantity}{name, location='vertices', indices = None}
2278  Module: \module{abstract\_2d\_finite\_volumes.domain}
2279
2280  Allow access to individual quantities and their methods
2281
2282  \end{funcdesc}
2283
2284
2285  \begin{funcdesc}{get\_values}{location='vertices', indices = None}
2286  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2287
2288  Extract values for quantity as an array
2289
2290  \end{funcdesc}
2291
2292
2293  \begin{funcdesc}{get\_integral}{}
2294  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2295
2296  Return computed integral over entire domain for this quantity
2297
2298  \end{funcdesc}
2299
2300
2301
2302
2303  \begin{funcdesc}{get\_maximum\_value}{indices = None}
2304  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2305
2306  Return maximum value of quantity (on centroids)
2307
2308  Optional argument indices is the set of element ids that
2309  the operation applies to. If omitted all elements are considered.
2310
2311  We do not seek the maximum at vertices as each vertex can
2312  have multiple values - one for each triangle sharing it.
2313  \end{funcdesc}
2314
2315
2316
2317  \begin{funcdesc}{get\_maximum\_location}{indices = None}
2318  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2319
2320  Return location of maximum value of quantity (on centroids)
2321
2322  Optional argument indices is the set of element ids that
2323  the operation applies to.
2324
2325  We do not seek the maximum at vertices as each vertex can
2326  have multiple values - one for each triangle sharing it.
2327
2328  If there are multiple cells with same maximum value, the
2329  first cell encountered in the triangle array is returned.
2330  \end{funcdesc}
2331
2332
2333
2334  \begin{funcdesc}{get\_wet\_elements}{indices=None}
2335  Module: \module{shallow\_water.shallow\_water\_domain}
2336
2337  Return indices for elements where h $>$ minimum_allowed_height
2338  Optional argument indices is the set of element ids that the operation applies to.
2339  \end{funcdesc}
2340
2341
2342  \begin{funcdesc}{get\_maximum\_inundation\_elevation}{indices=None}
2343  Module: \module{shallow\_water.shallow\_water\_domain}
2344
2345  Return highest elevation where h $>$ 0.\\
2346  Optional argument indices is the set of element ids that the operation applies to.\\
2347
2348  Example to find maximum runup elevation:\\
2349     z = domain.get_maximum_inundation_elevation()
2350  \end{funcdesc}
2351
2352
2353  \begin{funcdesc}{get\_maximum\_inundation\_location}{indices=None}
2354  Module: \module{shallow\_water.shallow\_water\_domain}
2355
2356  Return location (x,y) of highest elevation where h $>$ 0.\\
2357  Optional argument indices is the set of element ids that the operation applies to.\\
2358
2359  Example to find maximum runup location:\\
2360     x, y = domain.get_maximum_inundation_location()
2361  \end{funcdesc}
2362
2363
2364
2365\section{Other}
2366
2367  \begin{funcdesc}{domain.create\_quantity\_from\_expression}{string}
2368
2369  Handy for creating derived quantities on-the-fly as for example
2370  \begin{verbatim}
2371  Depth = domain.create_quantity_from_expression('stage-elevation')
2372
2373  exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5')
2374  Absolute_momentum = domain.create_quantity_from_expression(exp)
2375  \end{verbatim}
2376
2377  %See also \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
2378  \end{funcdesc}
2379
2380
2381
2382
2383
2384%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2385%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2386
2387\chapter{\anuga System Architecture}
2388
2389
2390\section{File Formats}
2391\label{sec:file formats}
2392
2393\anuga makes use of a number of different file formats. The
2394following table lists all these formats, which are described in more
2395detail in the paragraphs below.
2396
2397\bigskip
2398
2399\begin{center}
2400
2401\begin{tabular}{|ll|}  \hline
2402
2403\textbf{Extension} & \textbf{Description} \\
2404\hline\hline
2405
2406\code{.sww} & NetCDF format for storing model output
2407\code{f(t,x,y)}\\
2408
2409\code{.tms} & NetCDF format for storing time series \code{f(t)}\\
2410
2411\code{.xya} & ASCII format for storing arbitrary points and
2412associated attributes\\
2413
2414\code{.pts} & NetCDF format for storing arbitrary points and
2415associated attributes\\
2416
2417\code{.asc} & ASCII format of regular DEMs as output from ArcView\\
2418
2419\code{.prj} & Associated ArcView file giving more metadata for
2420\code{.asc} format\\
2421
2422\code{.ers} & ERMapper header format of regular DEMs for ArcView\\
2423
2424\code{.dem} & NetCDF representation of regular DEM data\\
2425
2426\code{.tsh} & ASCII format for storing meshes and associated
2427boundary and region info\\
2428
2429\code{.msh} & NetCDF format for storing meshes and associated
2430boundary and region info\\
2431
2432\code{.nc} & Native ferret NetCDF format\\
2433
2434\code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
2435%\caption{File formats used by \anuga}
2436\end{tabular}
2437
2438
2439\end{center}
2440
2441The above table shows the file extensions used to identify the
2442formats of files. However, typically, in referring to a format we
2443capitalise the extension and omit the initial full stop---thus, we
2444refer, for example, to `SWW files' or `PRJ files'.
2445
2446\bigskip
2447
2448A typical dataflow can be described as follows:
2449
2450\subsection{Manually Created Files}
2451
2452\begin{tabular}{ll}
2453ASC, PRJ & Digital elevation models (gridded)\\
2454NC & Model outputs for use as boundary conditions (e.g. from MOST)
2455\end{tabular}
2456
2457\subsection{Automatically Created Files}
2458
2459\begin{tabular}{ll}
2460ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert
2461DEMs to native \code{.pts} file\\
2462
2463NC $\rightarrow$ SWW & Convert MOST boundary files to
2464boundary \code{.sww}\\
2465
2466PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\
2467
2468TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using
2469\code{animate}\\
2470
2471TSH + Boundary SWW $\rightarrow$ SWW & Simulation using
2472\code{\anuga}\\
2473
2474Polygonal mesh outline $\rightarrow$ & TSH or MSH
2475\end{tabular}
2476
2477
2478
2479
2480\bigskip
2481
2482\subsection{SWW and TMS Formats}
2483
2484The SWW and TMS formats are both NetCDF formats, and are of key
2485importance for \anuga.
2486
2487An SWW file is used for storing \anuga output and therefore pertains
2488to a set of points and a set of times at which a model is evaluated.
2489It contains, in addition to dimension information, the following
2490variables:
2491
2492\begin{itemize}
2493    \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays
2494    \item \code{elevation}, a Numeric array storing bed-elevations
2495    \item \code{volumes}, a list specifying the points at the vertices of each of the
2496    triangles
2497    % Refer here to the example to be provided in describing the simple example
2498    \item \code{time}, a Numeric array containing times for model
2499    evaluation
2500\end{itemize}
2501
2502
2503The contents of an SWW file may be viewed using the anuga viewer
2504\code{animate}, which creates an on-screen geometric
2505representation. See section \ref{sec:animate} (page
2506\pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more
2507on \code{animate}.
2508
2509Alternatively, there are tools, such as \code{ncdump}, that allow
2510you to convert an NetCDF file into a readable format such as the
2511Class Definition Language (CDL). The following is an excerpt from a
2512CDL representation of the output file \file{runup.sww} generated
2513from running the simple example \file{runup.py} of
2514Chapter \ref{ch:getstarted}:
2515
2516\verbatiminput{examples/bedslopeexcerpt.cdl}
2517
2518The SWW format is used not only for output but also serves as input
2519for functions such as \function{file\_boundary} and
2520\function{file\_function}, described in Chapter \ref{ch:interface}.
2521
2522A TMS file is used to store time series data that is independent of
2523position.
2524
2525
2526\subsection{Mesh File Formats}
2527
2528A mesh file is a file that has a specific format suited to
2529triangular meshes and their outlines. A mesh file can have one of
2530two formats: it can be either a TSH file, which is an ASCII file, or
2531an MSH file, which is a NetCDF file. A mesh file can be generated
2532from the function \function{create\_mesh\_from\_regions} (see
2533Section \ref{sec:meshgeneration}) and used to initialise a domain.
2534
2535A mesh file can define the outline of the mesh---the vertices and
2536line segments that enclose the region in which the mesh is
2537created---and the triangular mesh itself, which is specified by
2538listing the triangles and their vertices, and the segments, which
2539are those sides of the triangles that are associated with boundary
2540conditions.
2541
2542In addition, a mesh file may contain `holes' and/or `regions'. A
2543hole represents an area where no mesh is to be created, while a
2544region is a labelled area used for defining properties of a mesh,
2545such as friction values.  A hole or region is specified by a point
2546and bounded by a number of segments that enclose that point.
2547
2548A mesh file can also contain a georeference, which describes an
2549offset to be applied to $x$ and $y$ values---eg to the vertices.
2550
2551
2552\subsection{Formats for Storing Arbitrary Points and Attributes}
2553
2554An XYA file is used to store data representing arbitrary numerical
2555attributes associated with a set of points.
2556
2557The format for an XYA file is:\\
2558%\begin{verbatim}
2559
2560            first line:     \code{[attribute names]}\\
2561            other lines:  \code{x y [attributes]}\\
2562
2563            for example:\\
2564            \code{elevation, friction}\\
2565            \code{0.6, 0.7, 4.9, 0.3}\\
2566            \code{1.9, 2.8, 5, 0.3}\\
2567            \code{2.7, 2.4, 5.2, 0.3}
2568
2569        The first two columns are always implicitly assumed to be $x$, $y$ coordinates.
2570        Use the same delimiter for the attribute names and the data.
2571
2572        An XYA file can optionally end with a description of the georeference:
2573
2574            \code{\#geo reference}\\
2575            \code{56}\\
2576            \code{466600.0}\\
2577            \code{8644444.0}
2578
2579        Here the first number specifies the UTM zone (in this case zone 56)  and other numbers specify the
2580        easting and northing
2581        coordinates (in this case (466600.0, 8644444.0)) of the lower left corner.
2582
2583A PTS file is a NetCDF representation of the data held in an XYA
2584file. If the data is associated with a set of $N$ points, then the
2585data is stored using an $N \times 2$ Numeric array of float
2586variables for the points and an $N \times 1$ Numeric array for each
2587attribute.
2588
2589%\end{verbatim}
2590
2591\subsection{ArcView Formats}
2592
2593Files of the three formats ASC, PRJ and ERS are all associated with
2594data from ArcView.
2595
2596An ASC file is an ASCII representation of DEM output from ArcView.
2597It contains a header with the following format:
2598
2599\begin{tabular}{l l}
2600\code{ncols}      &   \code{753}\\
2601\code{nrows}      &   \code{766}\\
2602\code{xllcorner}  &   \code{314036.58727982}\\
2603\code{yllcorner}  & \code{6224951.2960092}\\
2604\code{cellsize}   & \code{100}\\
2605\code{NODATA_value} & \code{-9999}
2606\end{tabular}
2607
2608The remainder of the file contains the elevation data for each grid point
2609in the grid defined by the above information.
2610
2611A PRJ file is an ArcView file used in conjunction with an ASC file
2612to represent metadata for a DEM.
2613
2614
2615\subsection{DEM Format}
2616
2617A DEM file is a NetCDF representation of regular DEM data.
2618
2619
2620\subsection{Other Formats}
2621
2622
2623
2624
2625\subsection{Basic File Conversions}
2626\label{sec:basicfileconversions}
2627
2628  \begin{funcdesc}{sww2dem}{basename_in, basename_out = None,
2629            quantity = None,
2630            timestep = None,
2631            reduction = None,
2632            cellsize = 10,
2633            NODATA_value = -9999,
2634            easting_min = None,
2635            easting_max = None,
2636            northing_min = None,
2637            northing_max = None,
2638            expand_search = False,
2639            verbose = False,
2640            origin = None,
2641            datum = 'WGS84',
2642            format = 'ers'}
2643  Module: \module{shallow\_water.data\_manager}
2644
2645  Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
2646  ERS) of a desired grid size \code{cellsize} in metres.
2647  The easting and northing values are used if the user wished to clip the output
2648  file to a specified rectangular area. The \code{reduction} input refers to a function
2649  to reduce the quantities over all time step of the SWW file, example, maximum.
2650  \end{funcdesc}
2651
2652
2653  \begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
2654            easting_min=None, easting_max=None,
2655            northing_min=None, northing_max=None,
2656            use_cache=False, verbose=False}
2657  Module: \module{shallow\_water.data\_manager}
2658
2659  Takes DEM data (a NetCDF file representation of data from a regular Digital
2660  Elevation Model) and converts it to PTS format.
2661  \end{funcdesc}
2662
2663
2664 
2665%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2666%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2667
2668\chapter{\anuga mathematical background}
2669\label{cd:mathematical background}
2670
2671\section{Introduction}
2672
2673This chapter outlines the mathematics underpinning \anuga.
2674 
2675 
2676
2677\section{Model}
2678\label{sec:model}
2679
2680The shallow water wave equations are a system of differential
2681conservation equations which describe the flow of a thin layer of
2682fluid over terrain. The form of the equations are:
2683\[
2684\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
2685x}+\frac{\partial \GG}{\partial y}=\SSS
2686\]
2687where $\UU=\left[ {{\begin{array}{*{20}c}
2688 h & {uh} & {vh} \\
2689\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
2690$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
2691entering the system are bed elevation $z$ and stage (absolute water
2692level) $w$, where the relation $w = z + h$ holds true at all times.
2693The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
2694by
2695\[
2696\EE=\left[ {{\begin{array}{*{20}c}
2697 {uh} \hfill \\
2698 {u^2h+gh^2/2} \hfill \\
2699 {uvh} \hfill \\
2700\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
2701 {vh} \hfill \\
2702 {vuh} \hfill \\
2703 {v^2h+gh^2/2} \hfill \\
2704\end{array} }} \right]
2705\]
2706and the source term (which includes gravity and friction) is given
2707by
2708\[
2709\SSS=\left[ {{\begin{array}{*{20}c}
2710 0 \hfill \\
2711 -{gh(z_{x} + S_{fx} )} \hfill \\
2712 -{gh(z_{y} + S_{fy} )} \hfill \\
2713\end{array} }} \right]
2714\]
2715where $S_f$ is the bed friction. The friction term is modelled using
2716Manning's resistance law
2717\[
2718S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
2719=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
2720\]
2721in which $\eta$ is the Manning resistance coefficient.
2722
2723As demonstrated in our papers, \cite{modsim2005,Rob99l} these
2724equations provide an excellent model of flows associated with
2725inundation such as dam breaks and tsunamis.
2726
2727\section{Finite Volume Method}
2728\label{sec:fvm}
2729
2730We use a finite-volume method for solving the shallow water wave
2731equations \cite{Rob99l}. The study area is represented by a mesh of
2732triangular cells as in Figure~\ref{fig:mesh} in which the conserved
2733quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
2734in each volume are to be determined. The size of the triangles may
2735be varied within the mesh to allow greater resolution in regions of
2736particular interest.
2737
2738\begin{figure}
2739\begin{center}
2740\includegraphics[width=5.0cm,keepaspectratio=true]{step-five}
2741\caption{Triangular mesh used in our finite volume method. Conserved
2742quantities $h$, $uh$ and $vh$ are associated with the centroid of
2743each triangular cell.} \label{fig:mesh}
2744\end{center}
2745\end{figure}
2746
2747The equations constituting the finite-volume method are obtained by
2748integrating the differential conservation equations over each
2749triangular cell of the mesh. Introducing some notation we use $i$ to
2750refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
2751set of indices referring to the cells neighbouring the $i$th cell.
2752Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
2753the length of the edge between the $i$th and $j$th cells.
2754
2755By applying the divergence theorem we obtain for each volume an
2756equation which describes the rate of change of the average of the
2757conserved quantities within each cell, in terms of the fluxes across
2758the edges of the cells and the effect of the source terms. In
2759particular, rate equations associated with each cell have the form
2760$$
2761 \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
2762$$
2763where
2764\begin{itemize}
2765\item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
2766\item $\SSS_i$ is the source term associated with the $i$th cell,
2767and
2768\item $\HH_{ij}$ is the outward normal flux of
2769material across the \textit{ij}th edge.
2770\end{itemize}
2771
2772
2773%\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
2774%cells
2775%\item $m_{ij}$ is the midpoint of
2776%the \textit{ij}th edge,
2777%\item
2778%$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
2779%normal along the \textit{ij}th edge, and The
2780
2781The flux $\HH_{ij}$ is evaluated using a numerical flux function
2782$\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
2783water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
2784$$
2785H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
2786$$
2787
2788Then
2789$$
2790\HH_{ij}  = \HH(\UU_i(m_{ij}),
2791\UU_j(m_{ij}); \mathbf{n}_{ij})
2792$$
2793where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
2794$\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
2795\textit{ij}th edge. The function $\UU_i(x)$ for $x \in
2796T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
2797neighbouring  cells.
2798
2799We use a second order reconstruction to produce a piece-wise linear
2800function construction of the conserved quantities for  all $x \in
2801T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This
2802function is allowed to be discontinuous across the edges of the
2803cells, but the slope of this function is limited to avoid
2804artificially introduced oscillations.
2805
2806Godunov's method (see \cite{Toro-92}) involves calculating the
2807numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
2808solving the corresponding one dimensional Riemann problem normal to
2809the edge. We use the central-upwind scheme of \cite{KurNP2001} to
2810calculate an approximation of the flux across each edge.
2811
2812\begin{figure}
2813\begin{center}
2814\includegraphics[width=5.0cm,keepaspectratio=true]{step-reconstruct}
2815\caption{From the values of the conserved quantities at the centroid
2816of the cell and its neighbouring cells, a discontinuous piecewise
2817linear reconstruction of the conserved quantities is obtained.}
2818\label{fig:mesh:reconstruct}
2819\end{center}
2820\end{figure}
2821
2822In the computations presented in this paper we use an explicit Euler
2823time stepping method with variable timestepping adapted to the
2824observed CFL condition.
2825
2826
2827\section{Flux limiting}
2828
2829The shallow water equations are solved numerically using a
2830finite volume method on unstructured triangular grid.
2831The upwind central scheme due to Kurganov and Petrova is used as an
2832approximate Riemann solver for the computation of inviscid flux functions.
2833This makes it possible to handle discontinuous solutions.
2834
2835To alleviate the problems associated with numerical instabilities due to
2836small water depths near a wet/dry boundary we employ a new flux limiter that
2837ensures that unphysical fluxes are never encounted.
2838
2839
2840Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction,
2841$w$ the absolute water level (stage) and
2842$z$ the bed elevation. The latter are assumed to be relative to the
2843same height datum.
2844The conserved quantities tracked by ANUGA are momentum in the
2845$x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$)
2846and depth ($h = w-z$).
2847
2848The flux calculation requires access to the velocity vector $(u, v)$ 
2849where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively.
2850In the presence of very small water depths, these calculations become
2851numerically unreliable and will typically cause unphysical speeds.
2852
2853We have employed a flux limiter which replaces the calculations above with
2854the limited approximations.
2855\begin{equation}
2856  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h},
2857\end{equation}
2858where $h_0$ is a regularisation parameter that controls the minimal
2859magnitude of the denominator. Taking the limits we have for $\hat{u}$
2860\[
2861  \lim_{h \rightarrow 0} \hat{u} = 
2862  \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0
2863\]
2864and
2865\[
2866  \lim_{h \rightarrow \infty} \hat{u} = 
2867  \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u
2868\]
2869with similar results for $\hat{v}$.
2870
2871The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator)
2872\[
2873  1 - h_0/h^2 = 0
2874\]
2875or
2876\[
2877  h_0 = h^2
2878\]
2879
2880
2881ANUGA has a global parameter $H_0$ that controls the minimal depth which
2882is considered in the various equations. This parameter is typically set to
2883$10^{-3}$. Setting
2884\[
2885  h_0 = H_0^2
2886\]
2887provides a reasonable balance between accurracy and stability. In fact,
2888setting $h=H_0$ will scale the predicted speed by a factor of $0.5$:
2889\[
2890  \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0}
2891\]
2892In general, for multiples of the minimal depth $N H_0$ one obtains
2893\[
2894  \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} = 
2895  \frac{\mu}{H_0 (1 + 1/N^2)}
2896\]
2897which converges quadratically to the true value with the multiple N.
2898
2899
2900%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.
2901
2902
2903
2904
2905
2906\section{Slope limiting}
2907A 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.
2908
2909However 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.
2910
2911
2912Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
2913let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.
2914Define the minimal depth across all vertices as $\hmin$ as
2915\[
2916  \hmin = \min_i h_i
2917\]
2918
2919Let $\tilde{w_i}$ be the stage obtained from a gradient limiter
2920limiting on stage only. The corresponding depth is then defined as
2921\[
2922  \tilde{h_i} = \tilde{w_i} - z_i
2923\]
2924We would use this limiter in deep water which we will define (somewhat boldly)
2925as
2926\[
2927  \hmin \ge \epsilon
2928\]
2929
2930
2931Similarly, let $\bar{w_i}$ be the stage obtained from a gradient
2932limiter limiting on depth respecting the bed slope.
2933The corresponding depth is defined as
2934\[
2935  \bar{h_i} = \bar{w_i} - z_i
2936\]
2937
2938
2939We introduce the concept of a balanced stage $w_i$ which is obtained as
2940the linear combination
2941
2942\[
2943  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
2944\]
2945or
2946\[
2947  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
2948\]
2949where $\alpha \in [0, 1]$.
2950
2951Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope
2952is ignored we have immediately that
2953\[
2954  \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0
2955\]
2956%where the maximal bed elevation range $dz$ is defined as
2957%\[
2958%  dz = \max_i |z_i - z|
2959%\]
2960
2961If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that
2962no negative depths occur. Formally, we will require that
2963\[
2964  \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i
2965\]
2966or
2967\begin{equation} 
2968  \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i
2969  \label{eq:limiter bound}
2970\end{equation} 
2971
2972There are two cases:
2973\begin{enumerate} 
2974  \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage)
2975  vertex is at least as far away from the bed than the shallow water
2976  (limited using depth). In this case we won't need any contribution from
2977  $\bar{h_i}$ and can accept any $alpha$.
2978 
2979  E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to
2980  \[
2981    \tilde{h_i} > \epsilon 
2982  \]
2983  whereas $\alpha=0$ yields
2984  \[
2985    \bar{h_i} > \epsilon       
2986  \]
2987  all well and good.
2988  \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is
2989  closer to the bed than the shallow water vertex or even below the bed.
2990  In this case we need to find an $alpha$ that will ensure a positive depth.   
2991  Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one
2992  obtains the bound
2993  \[
2994    \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i
2995  \] 
2996\end{enumerate} 
2997
2998Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one
2999arrives at the definition
3000\[
3001  \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}}
3002\]
3003which will guarantee that no vertex 'cuts' through the bed. Finally, should
3004$\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting
3005$alpha=0$ and similarly capping $\alpha$ at 1 just in case.
3006
3007%Furthermore,
3008%dropping the $\epsilon$ ensures that alpha is always positive and also 
3009%provides a numerical safety {??)
3010 
3011 
3012
3013
3014
3015%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3016%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3017
3018\chapter{Basic \anuga Assumptions}
3019
3020
3021Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
3022If one wished to recreate scenarios prior to that date it must be done
3023using some relative time (e.g. 0).
3024
3025
3026All spatial data relates to the WGS84 datum (or GDA94) and has been
3027projected into UTM with false easting of 500000 and false northing of
30281000000 on the southern hemisphere (0 on the northern).
3029
3030It is assumed that all computations take place within one UTM zone.
3031
3032DEMs, meshes and boundary conditions can have different origins within
3033one UTM zone. However, the computation will use that of the mesh for
3034numerical stability.
3035
3036When generating a mesh it is assumed that polygons do not cross.
3037Having polygons tht cross can cause the mesh generation to fail or bad
3038meshes being produced.
3039
3040
3041%OLD
3042%The dataflow is: (See data_manager.py and from scenarios)
3043%
3044%
3045%Simulation scenarios
3046%--------------------%
3047%%
3048%
3049%Sub directories contain scrips and derived files for each simulation.
3050%The directory ../source_data contains large source files such as
3051%DEMs provided externally as well as MOST tsunami simulations to be used
3052%as boundary conditions.
3053%
3054%Manual steps are:
3055%  Creation of DEMs from argcview (.asc + .prj)
3056%  Creation of mesh from pmesh (.tsh)
3057%  Creation of tsunami simulations from MOST (.nc)
3058%%
3059%
3060%Typical scripted steps are%
3061%
3062%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
3063%                   native dem and pts formats%
3064%
3065%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
3066%                  as boundary condition%
3067%
3068%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
3069%                   smoothing. The outputs are tsh files with elevation data.%
3070%
3071%  run_simulation.py: Use the above together with various parameters to
3072%                     run inundation simulation.
3073
3074
3075%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3076%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3077
3078\appendix
3079
3080\chapter{Supporting Tools}
3081\label{ch:supportingtools}
3082
3083This section describes a number of supporting tools, supplied with \anuga, that offer a
3084variety of types of functionality and enhance the basic capabilities of \anuga.
3085
3086\section{caching}
3087\label{sec:caching}
3088
3089The \code{cache} function is used to provide supervised caching of function
3090results. A Python function call of the form
3091
3092      {\small \begin{verbatim}
3093      result = func(arg1,...,argn)
3094      \end{verbatim}}
3095
3096  can be replaced by
3097
3098      {\small \begin{verbatim}
3099      from caching import cache
3100      result = cache(func,(arg1,...,argn))
3101      \end{verbatim}}
3102
3103  which returns the same output but reuses cached
3104  results if the function has been computed previously in the same context.
3105  \code{result} and the arguments can be simple types, tuples, list, dictionaries or
3106  objects, but not unhashable types such as functions or open file objects.
3107  The function \code{func} may be a member function of an object or a module.
3108
3109  This type of caching is particularly useful for computationally intensive
3110  functions with few frequently used combinations of input arguments. Note that
3111  if the inputs or output are very large caching may not save time because
3112  disc access may dominate the execution time.
3113
3114  If the function definition changes after a result has been cached, this will be
3115  detected by examining the functions \code{bytecode (co\_code, co\_consts,
3116  func\_defaults, co\_argcount)} and the function will be recomputed.
3117  However, caching will not detect changes in modules used by \code{func}.
3118  In this case cache must be cleared manually.
3119
3120  Options are set by means of the function \code{set\_option(key, value)},
3121  where \code{key} is a key associated with a
3122  Python dictionary \code{options}. This dictionary stores settings such as the name of
3123  the directory used, the maximum
3124  number of cached files allowed, and so on.
3125
3126  The \code{cache} function allows the user also to specify a list of dependent files. If any of these
3127  have been changed, the function is recomputed and the results stored again.
3128
3129  %Other features include support for compression and a capability to \ldots
3130
3131
3132   \textbf{USAGE:} \nopagebreak
3133
3134    {\small \begin{verbatim}
3135    result = cache(func, args, kwargs, dependencies, cachedir, verbose,
3136                   compression, evaluate, test, return_filename)
3137    \end{verbatim}}
3138
3139
3140\section{ANUGA viewer - animate}
3141\label{sec:animate}
3142 The output generated by \anuga may be viewed by
3143means of the visualisation tool \code{animate}, which takes the
3144\code{SWW} file output by \anuga and creates a visual representation
3145of the data. Examples may be seen in Figures \ref{fig:runupstart}
3146and \ref{fig:runup2}. To view an \code{SWW} file with
3147\code{animate} in the Windows environment, you can simply drag the
3148icon representing the file over an icon on the desktop for the
3149\code{animate} executable file (or a shortcut to it), or set up a
3150file association to make files with the extension \code{.sww} open
3151with \code{animate}. Alternatively, you can operate \code{animate}
3152from the command line, in both Windows and Linux environments.
3153
3154On successful operation, you will see an interactive moving-picture
3155display. You can use keys and the mouse to slow down, speed up or
3156stop the display, change the viewing position or carry out a number
3157of other simple operations. Help is also displayed when you press
3158the \code{h} key.
3159
3160The main keys operating the interactive screen are:\\
3161
3162\begin{center}
3163\begin{tabular}{|ll|}   \hline
3164
3165\code{w} & toggle wireframe \\
3166
3167space bar & start/stop\\
3168
3169up/down arrows & increase/decrease speed\\
3170
3171left/right arrows & direction in time \emph{(when running)}\\
3172& step through simulation \emph{(when stopped)}\\
3173
3174left mouse button & rotate\\
3175
3176middle mouse button & pan\\
3177
3178right mouse button & zoom\\  \hline
3179
3180\end{tabular}
3181\end{center}
3182
3183\vfill
3184
3185The following table describes how to operate animate from the command line:
3186
3187Usage: \code{animate [options] swwfile \ldots}\\  \nopagebreak
3188Options:\\  \nopagebreak
3189\begin{tabular}{ll}
3190  \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\
3191                                    & \code{HEAD\_MOUNTED\_DISPLAY}\\
3192  \code{--rgba} & Request a RGBA colour buffer visual\\
3193  \code{--stencil} & Request a stencil buffer visual\\
3194  \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
3195                                    & overridden by environmental variable\\
3196  \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD\_BUFFER | HORIZONTAL\_SPLIT |}\\
3197                                    & \code{VERTICAL\_SPLIT | LEFT\_EYE | RIGHT\_EYE |}\\
3198                                     & \code{ON | OFF} \\
3199  \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
3200  \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
3201\end{tabular}
3202
3203\begin{tabular}{ll}
3204  \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
3205  \code{-help} & Display this information\\
3206  \code{-hmax <float>} & Height above which transparency is set to
3207                                     \code{alphamax}\\
3208\end{tabular}
3209
3210\begin{tabular}{ll}
3211
3212  \code{-hmin <float>} & Height below which transparency is set to
3213                                     zero\\
3214\end{tabular}
3215
3216\begin{tabular}{ll}
3217  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
3218                                     up, default is overhead)\\
3219\end{tabular}
3220
3221\begin{tabular}{ll}
3222  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
3223
3224\end{tabular}
3225
3226\begin{tabular}{ll}
3227  \code{-movie <dirname>} & Save numbered images to named directory and
3228                                     quit\\
3229
3230  \code{-nosky} & Omit background sky\\
3231
3232
3233  \code{-scale <float>} & Vertical scale factor\\
3234  \code{-texture <file>} & Image to use for bedslope topography\\
3235  \code{-tps <rate>} & Timesteps per second\\
3236  \code{-version} & Revision number and creation (not compile)
3237                                     date\\
3238\end{tabular}
3239
3240\section{utilities/polygons}
3241
3242  \declaremodule{standard}{utilities.polygon}
3243  \refmodindex{utilities.polygon}
3244
3245  \begin{classdesc}{Polygon\_function}{regions, default = 0.0, geo_reference = None}
3246  Module: \code{utilities.polygon}
3247
3248  Creates a callable object that returns one of a specified list of values when
3249  evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the
3250  point belongs to. The parameter \code{regions} is a list of pairs
3251  \code{(P, v)}, where each \code{P} is a polygon and each \code{v}
3252  is either a constant value or a function of coordinates \code{x}
3253  and \code{y}, specifying the return value for a point inside \code{P}. The
3254  optional parameter \code{default} may be used to specify a value
3255  for a point not lying inside any of the specified polygons. When a
3256  point lies in more than one polygon, the return value is taken to
3257  be the value for whichever of these polygon appears later in the
3258  list.
3259  %FIXME (Howard): CAN x, y BE VECTORS?
3260
3261  \end{classdesc}
3262
3263  \begin{funcdesc}{read\_polygon}{filename}
3264  Module: \code{utilities.polygon}
3265
3266  Reads the specified file and returns a polygon. Each
3267  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
3268  as coordinates of one vertex of the polygon.
3269  \end{funcdesc}
3270
3271  \begin{funcdesc}{populate\_polygon}{polygon, number_of_points, seed = None, exclude = None}
3272  Module: \code{utilities.polygon}
3273
3274  Populates the interior of the specified polygon with the specified number of points,
3275  selected by means of a uniform distribution function.
3276  \end{funcdesc}
3277
3278  \begin{funcdesc}{point\_in\_polygon}{polygon, delta=1e-8}
3279  Module: \code{utilities.polygon}
3280
3281  Returns a point inside the specified polygon and close to the edge. The distance between
3282  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
3283  second argument \code{delta}, which is taken as $10^{-8}$ by default.
3284  \end{funcdesc}
3285
3286  \begin{funcdesc}{inside\_polygon}{points, polygon, closed = True, verbose = False}
3287  Module: \code{utilities.polygon}
3288
3289  Used to test whether the members of a list of points
3290  are inside the specified polygon. Returns a Numeric
3291  array comprising the indices of the points in the list that lie inside the polygon.
3292  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
3293  Points on the edges of the polygon are regarded as inside if
3294  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
3295  \end{funcdesc}
3296
3297  \begin{funcdesc}{outside\_polygon}{points, polygon, closed = True, verbose = False}
3298  Module: \code{utilities.polygon}
3299
3300  Exactly like \code{inside\_polygon}, but with the words `inside' and `outside' interchanged.
3301  \end{funcdesc}
3302
3303  \begin{funcdesc}{is\_inside\_polygon}{point, polygon, closed=True, verbose=False}
3304  Module: \code{utilities.polygon}
3305
3306  Returns \code{True} if \code{point} is inside \code{polygon} or
3307  \code{False} otherwise. Points on the edges of the polygon are regarded as inside if
3308  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
3309  \end{funcdesc}
3310
3311  \begin{funcdesc}{is\_outside\_polygon}{point, polygon, closed=True, verbose=False}
3312  Module: \code{utilities.polygon}
3313
3314  Exactly like \code{is\_outside\_polygon}, but with the words `inside' and `outside' interchanged.
3315  \end{funcdesc}
3316
3317  \begin{funcdesc}{point\_on\_line}{x, y, x0, y0, x1, y1}
3318  Module: \code{utilities.polygon}
3319
3320  Returns \code{True} or \code{False}, depending on whether the point with coordinates
3321  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
3322  and \code{x1, y1} (extended if necessary at either end).
3323  \end{funcdesc}
3324
3325  \begin{funcdesc}{separate\_points\_by\_polygon}{points, polygon, closed = True, verbose = False}
3326    \indexedcode{separate\_points\_by\_polygon}
3327  Module: \code{utilities.polygon}
3328
3329  \end{funcdesc}
3330
3331  \begin{funcdesc}{polygon\_area}{polygon}
3332  Module: \code{utilities.polygon}
3333
3334  Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html)
3335  \end{funcdesc}
3336
3337  \begin{funcdesc}{plot\_polygons}{polygons, figname, verbose = False}
3338  Module: \code{utilities.polygon}
3339
3340  Plots each polygon contained in input polygon list, e.g.
3341 \code{polygons = [poly1, poly2, poly3]} where \code{poly1 = [[x11,y11],[x12,y12],[x13,y13]]}
3342 etc.  Each polygon is closed for plotting purposes and subsequent plot saved to \code{figname}.
3343  Returns list containing the minimum and maximum of \code{x} and \code{y},
3344  i.e. \code{[x_{min}, x_{max}, y_{min}, y_{max}]}.
3345  \end{funcdesc}
3346
3347\section{coordinate\_transforms}
3348
3349\section{geospatial\_data}
3350\label{sec:geospatial}
3351
3352This describes a class that represents arbitrary point data in UTM
3353coordinates along with named attribute values.
3354
3355%FIXME (Ole): This gives a LaTeX error
3356%\declaremodule{standard}{geospatial_data}
3357%\refmodindex{geospatial_data}
3358
3359
3360
3361\begin{classdesc}{Geospatial\_data}
3362  {data_points = None,
3363    attributes = None,
3364    geo_reference = None,
3365    default_attribute_name = None,
3366    file_name = None}
3367Module: \code{geospatial\_data}
3368
3369This class is used to store a set of data points and associated
3370attributes, allowing these to be manipulated by methods defined for
3371the class.
3372
3373The data points are specified either by reading them from a NetCDF
3374or XYA file, identified through the parameter \code{file\_name}, or
3375by providing their \code{x}- and \code{y}-coordinates in metres,
3376either as a sequence of 2-tuples of floats or as an $M \times 2$
3377Numeric array of floats, where $M$ is the number of points.
3378Coordinates are interpreted relative to the origin specified by the
3379object \code{geo\_reference}, which contains data indicating the UTM
3380zone, easting and northing. If \code{geo\_reference} is not
3381specified, a default is used.
3382
3383Attributes are specified through the parameter \code{attributes},
3384set either to a list or array of length $M$ or to a dictionary whose
3385keys are the attribute names and whose values are lists or arrays of
3386length $M$. One of the attributes may be specified as the default
3387attribute, by assigning its name to \code{default\_attribute\_name}.
3388If no value is specified, the default attribute is taken to be the
3389first one.
3390\end{classdesc}
3391
3392
3393\begin{methoddesc}{import\_points\_file}{delimiter = None, verbose = False}
3394
3395\end{methoddesc}
3396
3397
3398\begin{methoddesc}{export\_points\_file}{ofile, absolute=False}
3399
3400\end{methoddesc}
3401
3402
3403\begin{methoddesc}{get\_data\_points}{absolute = True}
3404
3405\end{methoddesc}
3406
3407
3408\begin{methoddesc}{set\_attributes}{attributes}
3409
3410\end{methoddesc}
3411
3412
3413\begin{methoddesc}{get\_attributes}{attribute_name = None}
3414
3415\end{methoddesc}
3416
3417
3418\begin{methoddesc}{get\_all\_attributes}{}
3419
3420\end{methoddesc}
3421
3422
3423\begin{methoddesc}{set\_default\_attribute\_name}{default_attribute_name}
3424
3425\end{methoddesc}
3426
3427
3428\begin{methoddesc}{set\_geo\_reference}{geo_reference}
3429
3430\end{methoddesc}
3431
3432
3433\begin{methoddesc}{add}{}
3434
3435\end{methoddesc}
3436
3437
3438\begin{methoddesc}{clip}{}
3439Clip geospatial data by a polygon
3440
3441Inputs are \code{polygon} which is either a list of points, an Nx2 array or
3442a Geospatial data object and \code{closed}(optional) which determines
3443whether points on boundary should be regarded as belonging to the polygon
3444(\code{closed=True}) or not (\code{closed=False}).
3445Default is \code{closed=True}.
3446
3447Returns new Geospatial data object representing points
3448inside specified polygon.
3449\end{methoddesc}
3450
3451
3452\section{pmesh GUI}
3453 \emph{Pmesh}
3454allows the user to set up the mesh of the problem interactively.
3455It can be used to build the outline of a mesh or to visualise a mesh
3456automatically generated.
3457
3458Pmesh will let the user select various modes. The current allowable
3459modes are vertex, segment, hole or region.  The mode describes what
3460sort of object is added or selected in response to mouse clicks.  When
3461changing modes any prior selected objects become deselected.
3462
3463In general the left mouse button will add an object and the right
3464mouse button will select an object.  A selected object can de deleted
3465by pressing the the middle mouse button (scroll bar).
3466
3467\section{alpha\_shape}
3468\emph{Alpha shapes} are used to generate close-fitting boundaries
3469around sets of points. The alpha shape algorithm produces a shape
3470that approximates to the `shape formed by the points'---or the shape
3471that would be seen by viewing the points from a coarse enough
3472resolution. For the simplest types of point sets, the alpha shape
3473reduces to the more precise notion of the convex hull. However, for
3474many sets of points the convex hull does not provide a close fit and
3475the alpha shape usually fits more closely to the original point set,
3476offering a better approximation to the shape being sought.
3477
3478In \anuga, an alpha shape is used to generate a polygonal boundary
3479around a set of points before mesh generation. The algorithm uses a
3480parameter $\alpha$ that can be adjusted to make the resultant shape
3481resemble the shape suggested by intuition more closely. An alpha
3482shape can serve as an initial boundary approximation that the user
3483can adjust as needed.
3484
3485The following paragraphs describe the class used to model an alpha
3486shape and some of the important methods and attributes associated
3487with instances of this class.
3488
3489\begin{classdesc}{Alpha\_Shape}{points, alpha = None}
3490Module: \code{alpha\_shape}
3491
3492To instantiate this class the user supplies the points from which
3493the alpha shape is to be created (in the form of a list of 2-tuples
3494\code{[[x1, y1],[x2, y2]}\ldots\code{]}, assigned to the parameter
3495\code{points}) and, optionally, a value for the parameter
3496\code{alpha}. The alpha shape is then computed and the user can then
3497retrieve details of the boundary through the attributes defined for
3498the class.
3499\end{classdesc}
3500
3501
3502\begin{funcdesc}{alpha\_shape\_via\_files}{point_file, boundary_file, alpha= None}
3503Module: \code{alpha\_shape}
3504
3505This function reads points from the specified point file
3506\code{point\_file}, computes the associated alpha shape (either
3507using the specified value for \code{alpha} or, if no value is
3508specified, automatically setting it to an optimal value) and outputs
3509the boundary to a file named \code{boundary\_file}. This output file
3510lists the coordinates \code{x, y} of each point in the boundary,
3511using one line per point.
3512\end{funcdesc}
3513
3514
3515\begin{methoddesc}{set\_boundary\_type}{self,raw_boundary=True,
3516                          remove_holes=False,
3517                          smooth_indents=False,
3518                          expand_pinch=False,
3519                          boundary_points_fraction=0.2}
3520Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3521
3522This function sets flags that govern the operation of the algorithm
3523that computes the boundary, as follows:
3524
3525\code{raw\_boundary = True} returns raw boundary, i.e. the regular edges of the
3526                alpha shape.\\
3527\code{remove\_holes = True} removes small holes (`small' is defined by
3528\code{boundary\_points\_fraction})\\
3529\code{smooth\_indents = True} removes sharp triangular indents in
3530boundary\\
3531\code{expand\_pinch = True} tests for pinch-off and
3532corrects---preventing a boundary vertex from having more than two edges.
3533\end{methoddesc}
3534
3535
3536\begin{methoddesc}{get\_boundary}{}
3537Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3538
3539Returns a list of tuples representing the boundary of the alpha
3540shape. Each tuple represents a segment in the boundary by providing
3541the indices of its two endpoints.
3542\end{methoddesc}
3543
3544
3545\begin{methoddesc}{write\_boundary}{file_name}
3546Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3547
3548Writes the list of 2-tuples returned by \code{get\_boundary} to the
3549file \code{file\_name}, using one line per tuple.
3550\end{methoddesc}
3551
3552\section{Numerical Tools}
3553
3554The following table describes some useful numerical functions that
3555may be found in the module \module{utilities.numerical\_tools}:
3556
3557\begin{tabular}{|p{8cm} p{8cm}|}  \hline
3558\code{angle(v1, v2=None)} & Angle between two-dimensional vectors
3559\code{v1} and \code{v2}, or between \code{v1} and the $x$-axis if
3560\code{v2} is \code{None}. Value is in range $0$ to $2\pi$. \\
3561
3562\code{normal\_vector(v)} & Normal vector to \code{v}.\\
3563
3564\code{mean(x)} & Mean value of a vector \code{x}.\\
3565
3566\code{cov(x, y=None)} & Covariance of vectors \code{x} and \code{y}.
3567If \code{y} is \code{None}, returns \code{cov(x, x)}.\\
3568
3569\code{err(x, y=0, n=2, relative=True)} & Relative error of
3570$\parallel$\code{x}$-$\code{y}$\parallel$ to
3571$\parallel$\code{y}$\parallel$ (2-norm if \code{n} = 2 or Max norm
3572if \code{n} = \code{None}). If denominator evaluates to zero or if
3573\code{y}
3574is omitted or if \code{relative = False}, absolute error is returned.\\
3575
3576\code{norm(x)} & 2-norm of \code{x}.\\
3577
3578\code{corr(x, y=None)} & Correlation of \code{x} and \code{y}. If
3579\code{y} is \code{None} returns autocorrelation of \code{x}.\\
3580
3581\code{ensure\_numeric(A, typecode = None)} & Returns a Numeric array
3582for any sequence \code{A}. If \code{A} is already a Numeric array it
3583will be returned unaltered. Otherwise, an attempt is made to convert
3584it to a Numeric array. (Needed because \code{array(A)} can
3585cause memory overflow.)\\
3586
3587\code{histogram(a, bins, relative=False)} & Standard histogram. If
3588\code{relative} is \code{True}, values will be normalised against
3589the total and thus represent frequencies rather than counts.\\
3590
3591\code{create\_bins(data, number\_of\_bins = None)} & Safely create
3592bins for use with histogram. If \code{data} contains only one point
3593or is constant, one bin will be created. If \code{number\_of\_bins}
3594is omitted, 10 bins will be created.\\  \hline
3595
3596\end{tabular}
3597
3598
3599\chapter{Modules available in \anuga}
3600
3601
3602\section{\module{abstract\_2d\_finite\_volumes.general\_mesh} }
3603\declaremodule[generalmesh]{}{general\_mesh}
3604\label{mod:generalmesh}
3605
3606\section{\module{abstract\_2d\_finite\_volumes.neighbour\_mesh} }
3607\declaremodule[neighbourmesh]{}{neighbour\_mesh}
3608\label{mod:neighbourmesh}
3609
3610\section{\module{abstract\_2d\_finite\_volumes.domain} --- Generic module for 2D triangular domains for finite-volume computations of conservation laws}
3611\declaremodule{}{domain}
3612\label{mod:domain}
3613
3614
3615\section{\module{abstract\_2d\_finite\_volumes.quantity}}
3616\declaremodule{}{quantity}
3617\label{mod:quantity}
3618
3619\begin{verbatim}
3620Class Quantity - Implements values at each triangular element
3621
3622To create:
3623
3624   Quantity(domain, vertex_values)
3625
3626   domain: Associated domain structure. Required.
3627
3628   vertex_values: N x 3 array of values at each vertex for each element.
3629                  Default None
3630
3631   If vertex_values are None Create array of zeros compatible with domain.
3632   Otherwise check that it is compatible with dimenions of domain.
3633   Otherwise raise an exception
3634
3635\end{verbatim}
3636
3637
3638
3639
3640\section{\module{shallow\_water} --- 2D triangular domains for finite-volume
3641computations of the shallow water wave equation. This module contains a specialisation
3642of class Domain from module domain.py consisting of methods specific to the Shallow Water
3643Wave Equation
3644}
3645\declaremodule[shallowwater]{}{shallow\_water}
3646\label{mod:shallowwater}
3647
3648
3649
3650
3651%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3652
3653\chapter{Frequently Asked Questions}
3654
3655
3656\section{General Questions}
3657
3658\subsubsection{What is \anuga?}
3659It is a software package suitable for simulating 2D water flows in
3660complex geometries.
3661
3662\subsubsection{Why is it called \anuga?}
3663The software was developed in collaboration between the
3664Australian National University (ANU) and Geoscience Australia (GA).
3665
3666\subsubsection{How do I obtain a copy of \anuga?}
3667See \url{https://datamining.anu.edu.au/anuga} for all things ANUGA.
3668
3669%\subsubsection{What developments are expected for \anuga in the future?}
3670%This
3671
3672\subsubsection{Are there any published articles about \anuga that I can reference?}
3673See \url{https://datamining.anu.edu.au/anuga} for links.
3674
3675
3676\section{Modelling Questions}
3677
3678\subsubsection{Which type of problems are \anuga good for?}
3679General 2D waterflows in complex geometries such as
3680dam breaks, flows amoung structurs, coastal inundation etc.
3681
3682\subsubsection{Which type of problems are beyond the scope of \anuga?}
3683See Chapter \ref{ch:limitations}.
3684
3685\subsubsection{Can I start the simulation at an arbitrary time?}
3686Yes, using \code{domain.set\_time()} you can specify an arbitrary
3687starting time. This is for example useful in conjunction with a
3688file\_boundary, which may start hours before anything hits the model
3689boundary. By assigning a later time for the model to start,
3690computational resources aren't wasted.
3691
3692\subsubsection{Can I change values for any quantity during the simulation?}
3693Yes, using \code{domain.set\_quantity()} inside the domain.evolve
3694loop you can change values of any quantity. This is for example
3695useful if you wish to let the system settle for a while before
3696assigning an initial condition. Another example would be changing
3697the values for elevation to model e.g. erosion.
3698
3699\subsubsection{Can I change boundary conditions during the simulation?}
3700Yes - see example on page \pageref{sec:change boundary code} in Section
3701\ref{sec:change boundary}.
3702
3703\subsubsection{How do I access model time during the simulation?}
3704The variable \code{t} in the evolve for loop is the model time.
3705For 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})
3706one would write something like
3707{\small \begin{verbatim}
3708    for t in domain.evolve(yieldstep = 0.2, duration = 40.0):
3709
3710        if Numeric.allclose(t, 15):
3711            print 'Changing boundary to outflow'
3712            domain.set_boundary({'right': Bo})
3713
3714\end{verbatim}}
3715The 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()}.
3716
3717
3718\subsubsection{Why does a file\_function return a list of numbers when evaluated?}
3719Currently, file\_function works by returning values for the conserved
3720quantities \code{stage}, \code{xmomentum} and \code{ymomentum} at a given point in time
3721and space as a triplet. To access e.g.\ \code{stage} one must specify element 0 of the
3722triplet returned by file\_function.
3723
3724\subsubsection{Which diagnostics are available to troubleshoot a simulation?}
3725
3726\subsubsection{How do I use a DEM in my simulation?}
3727You use \code{dem2pts} to convert your DEM to the required .pts format. This .pts file is then called
3728when setting the elevation data to the mesh in \code{domain.set_quantity}
3729
3730\subsubsection{What sort of DEM resolution should I use?}
3731Try and work with the \emph{best} you have available. Onshore DEMs
3732are typically available in 25m, 100m and 250m grids. Note, offshore
3733data is often sparse, or non-existent.
3734
3735\subsubsection{What sort of mesh resolution should I use?}
3736The 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,
3737if 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,
3738you need a fine mesh over regions where the DEM changes rapidly, and other areas of significant interest, such as the coast.
3739If 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.
3740
3741
3742\subsubsection{How do I tag interior polygons?}
3743At the moment create_mesh_from_regions does not allow interior
3744polygons with symbolic tags. If tags are needed, the interior
3745polygons must be created subsequently. For example, given a filename
3746of polygons representing solid walls (in Arc Ungenerate format) can
3747be tagged as such using the code snippet:
3748\begin{verbatim}
3749  # Create mesh outline with tags
3750  mesh = create_mesh_from_regions(bounding_polygon,
3751                                  boundary_tags=boundary_tags)
3752  # Add buildings outlines with tags set to 'wall'. This would typically
3753  # bind to a Reflective boundary
3754  mesh.import_ungenerate_file(buildings_filename, tag='wall')
3755
3756  # Generate and write mesh to file
3757  mesh.generate_mesh(maximum_triangle_area=max_area)
3758  mesh.export_mesh_file(mesh_filename)
3759\end{verbatim}
3760
3761Note that a mesh object is returned from \code{create_mesh_from_regions}
3762when file name is omitted.
3763
3764\subsubsection{How often should I store the output?}
3765This 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
3766to look in detail at the evolution, then you will need to balance your storage requirements and the duration of the simulation.
3767If 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
3768quantities on a mesh with approximately 300000 triangles on a 2 min interval for 5 hours will result in approximately 350Mb SWW file
3769(as for the \file{run\_sydney\_smf.py} example).
3770
3771\subsection{Boundary Conditions}
3772
3773\subsubsection{How do I create a Dirichlet boundary condition?}
3774
3775A Dirichlet boundary condition sets a constant value for the
3776conserved quantities at the boundaries. A list containing
3777the constant values for stage, xmomentum and ymomentum is constructed
3778and used in the function call, e.g. \code{Dirichlet_boundary([0.2,0.,0.])}
3779
3780\subsubsection{How do I know which boundary tags are available?}
3781The method \code{domain.get\_boundary\_tags()} will return a list of
3782available tags for use with
3783\code{domain.set\_boundary\_condition()}.
3784
3785
3786
3787
3788
3789\chapter{Glossary}
3790
3791\begin{tabular}{|lp{10cm}|c|}  \hline
3792%\begin{tabular}{|llll|}  \hline
3793    \emph{Term} & \emph{Definition} & \emph{Page}\\  \hline
3794
3795    \indexedbold{\anuga} & Name of software (joint development between ANU and
3796    GA) & \pageref{def:anuga}\\
3797
3798    \indexedbold{bathymetry} & offshore elevation &\\
3799
3800    \indexedbold{conserved quantity} & conserved (stage, x and y
3801    momentum) & \\
3802
3803%    \indexedbold{domain} & The domain of a function is the set of all input values to the
3804%    function.&\\
3805
3806    \indexedbold{Digital Elevation Model (DEM)} & DEMs are digital files consisting of points of elevations,
3807sampled systematically at equally spaced intervals.& \\
3808
3809    \indexedbold{Dirichlet boundary} & A boundary condition imposed on a differential equation
3810 that specifies the values the solution is to take on the boundary of the
3811 domain. & \pageref{def:dirichlet boundary}\\
3812
3813    \indexedbold{edge} & A triangular cell within the computational mesh can be depicted
3814    as a set of vertices joined by lines (the edges). & \\
3815
3816    \indexedbold{elevation} & refers to bathymetry and topography &\\
3817
3818    \indexedbold{evolution} & integration of the shallow water wave equations
3819    over time &\\
3820
3821    \indexedbold{finite volume method} & The method evaluates the terms in the shallow water
3822    wave equation as fluxes at the surfaces of each finite volume. Because the
3823    flux entering a given volume is identical to that leaving the adjacent volume,
3824    these methods are conservative. Another advantage of the finite volume method is
3825    that it is easily formulated to allow for unstructured meshes. The method is used
3826    in many computational fluid dynamics packages. & \\
3827
3828    \indexedbold{forcing term} & &\\
3829
3830    \indexedbold{flux} & the amount of flow through the volume per unit
3831    time & \\
3832
3833    \indexedbold{grid} & Evenly spaced mesh & \\
3834
3835    \indexedbold{latitude} & The angular distance on a mericlear north and south of the
3836    equator, expressed in degrees and minutes. & \\
3837
3838    \indexedbold{longitude} & The angular distance east or west, between the meridian
3839    of a particular place on Earth and that of the Prime Meridian (located in Greenwich,
3840    England) expressed in degrees or time.& \\
3841
3842    \indexedbold{Manning friction coefficient} & &\\
3843
3844    \indexedbold{mesh} & Triangulation of domain &\\
3845
3846    \indexedbold{mesh file} & A TSH or MSH file & \pageref{def:mesh file}\\
3847
3848    \indexedbold{NetCDF} & &\\
3849
3850    \indexedbold{node} & A point at which edges meet & \\
3851
3852    \indexedbold{northing} & A rectangular (x,y) coordinate measurement of distance
3853    north from a north-south reference line, usually a meridian used as the axis of
3854    origin within a map zone or projection. Northing is a UTM (Universal Transverse
3855    Mercator) coordinate. & \\
3856
3857    \indexedbold{points file} & A PTS or XYA file & \\  \hline
3858
3859    \end{tabular}
3860
3861    \begin{tabular}{|lp{10cm}|c|}  \hline
3862
3863    \indexedbold{polygon} & A sequence of points in the plane. \anuga represents a polygon
3864    either as a list consisting of Python tuples or lists of length 2 or as an $N \times 2$
3865    Numeric array, where $N$ is the number of points.
3866
3867    The unit square, for example, would be represented either as
3868    \code{[ [0,0], [1,0], [1,1], [0,1] ]} or as \code{array( [0,0], [1,0], [1,1],
3869    [0,1] )}.
3870
3871    NOTE: For details refer to the module \module{utilities/polygon.py}. &
3872    \\     \indexedbold{resolution} &  The maximal area of a triangular cell in a
3873    mesh & \\
3874
3875
3876    \indexedbold{reflective boundary} & Models a solid wall. Returns same conserved
3877    quantities as those present in the neighbouring volume but reflected. Specific to the
3878    shallow water equation as it works with the momentum quantities assumed to be the
3879    second and third conserved quantities. & \pageref{def:reflective boundary}\\
3880
3881    \indexedbold{stage} & &\\
3882
3883%    \indexedbold{try this}
3884
3885    \indexedbold{animate} & visualisation tool used with \anuga &
3886    \pageref{sec:animate}\\
3887
3888    \indexedbold{time boundary} & Returns values for the conserved
3889quantities as a function of time. The user must specify
3890the domain to get access to the model time. & \pageref{def:time boundary}\\
3891
3892    \indexedbold{topography} & onshore elevation &\\
3893
3894    \indexedbold{transmissive boundary} & & \pageref{def:transmissive boundary}\\
3895
3896    \indexedbold{vertex} & A point at which edges meet. & \\
3897
3898    \indexedbold{xmomentum} & conserved quantity (note, two-dimensional SWW equations say
3899    only \code{x} and \code{y} and NOT \code{z}) &\\
3900
3901    \indexedbold{ymomentum}  & conserved quantity & \\  \hline
3902
3903    \end{tabular}
3904
3905
3906%The \code{\e appendix} markup need not be repeated for additional
3907%appendices.
3908
3909
3910%
3911%  The ugly "%begin{latexonly}" pseudo-environments are really just to
3912%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
3913%  not really valuable.
3914%
3915%  If you don't want the Module Index, you can remove all of this up
3916%  until the second \input line.
3917%
3918
3919%begin{latexonly}
3920%\renewcommand{\indexname}{Module Index}
3921%end{latexonly}
3922\input{mod\jobname.ind}        % Module Index
3923%
3924%begin{latexonly}
3925%\renewcommand{\indexname}{Index}
3926%end{latexonly}
3927\input{\jobname.ind}            % Index
3928
3929
3930
3931\begin{thebibliography}{99}
3932\bibitem[nielsen2005]{nielsen2005}
3933{\it Hydrodynamic modelling of coastal inundation}.
3934Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman.
3935In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
3936Modelling and Simulation. Modelling and Simulation Society of Australia and
3937New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.\\
3938http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
3939
3940\bibitem[grid250]{grid250}
3941Australian Bathymetry and Topography Grid, June 2005.
3942Webster, M.A. and Petkovic, P.
3943Geoscience Australia Record 2005/12. ISBN: 1-920871-46-2.\\
3944http://www.ga.gov.au/meta/ANZCW0703008022.html
3945
3946
3947\end{thebibliography}{99}
3948
3949\end{document}
Note: See TracBrowser for help on using the repository browser.