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

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

Warning about polygon points close together causing very small triangles.

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