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

Last change on this file since 5475 was 5475, checked in by ole, 11 years ago

FAQ update: get_revision_number()

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