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

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

Typo

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