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

Last change on this file since 5130 was 5130, checked in by ole, 16 years ago

Struggled with getting ~ to show up in hyperlink

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 157.5 KB
Line 
1% Complete documentation on the extended LaTeX markup used for Python
2% documentation is available in ``Documenting Python'', which is part
3% of the standard documentation for Python.  It may be found online
4% at:
5%
6%     http://www.python.org/doc/current/doc/doc.html
7
8
9%labels
10%Sections and subsections \label{sec: }
11%Chapters \label{ch: }
12%Equations \label{eq: }
13%Figures \label{fig: }
14
15% Is latex failing with;
16% `modanuga_user_manual.ind' not found?
17% try this command-line
18%   makeindex modanuga_user_manual.idx
19% To produce the modanuga_user_manual.ind file.
20
21
22\documentclass{manual}
23
24\usepackage{graphicx}
25\usepackage{datetime}
26
27\input{definitions}
28
29\title{\anuga User Manual}
30\author{Geoscience Australia and the Australian National University}
31
32% Please at least include a long-lived email address;
33% the rest is at your discretion.
34\authoraddress{Geoscience Australia \\
35  Email: \email{ole.nielsen@ga.gov.au}
36}
37
38%Draft date
39
40% update before release!
41% Use an explicit date so that reformatting
42% doesn't cause a new date to be used.  Setting
43% the date to \today can be used during draft
44% stages to make it easier to handle versions.
45
46
47\longdate       % Make date format long using datetime.sty
48%\settimeformat{xxivtime} % 24 hour Format
49\settimeformat{oclock} % Verbose
50\date{\today, \ \currenttime}
51%\hyphenation{set\_datadir}
52
53\ifhtml
54\date{\today} % latex2html does not know about datetime
55\fi
56
57
58
59
60\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-03-05 07:20:05 +0000 (Wed, 05 Mar 2008) $
86%$LastChangedRevision: 5130 $
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}
1707Module: \module{pmesh.mesh},  Class: \class{Mesh}
1708
1709This method is used to import a polygon file in the ungenerate
1710format, which is used by arcGIS. The polygons from the file are converted to
1711vertices and segments. \code{ofile} is the name of the polygon file.
1712\code{tag} is the tag given to all the polygon's segments.
1713
1714This function can be used to import building footprints.
1715\end{methoddesc}
1716
1717%%%%%%
1718\section{Initialising the Domain}
1719
1720%Include description of the class Domain and the module domain.
1721
1722%FIXME (Ole): This is also defined in a later chapter
1723%\declaremodule{standard}{...domain}
1724
1725\begin{classdesc} {Domain} {source=None,
1726                 triangles=None,
1727                 boundary=None,
1728                 conserved_quantities=None,
1729                 other_quantities=None,
1730                 tagged_elements=None,
1731                 use_inscribed_circle=False,
1732                 mesh_filename=None,
1733                 use_cache=False,
1734                 verbose=False,
1735                 full_send_dict=None,
1736                 ghost_recv_dict=None,
1737                 processor=0,
1738                 numproc=1}
1739Module: \refmodule{abstract_2d_finite_volumes.domain}
1740
1741This class is used to create an instance of a data structure used to
1742store and manipulate data associated with a mesh. The mesh is
1743specified either by assigning the name of a mesh file to
1744\code{source} or by specifying the points, triangle and boundary of the
1745mesh.
1746\end{classdesc}
1747
1748\subsection{Key Methods of Domain}
1749
1750\begin{methoddesc} {set\_name}{name}
1751    Module: \refmodule{abstract\_2d\_finite\_volumes.domain},
1752    page \pageref{mod:domain}
1753
1754    Assigns the name \code{name} to the domain.
1755\end{methoddesc}
1756
1757\begin{methoddesc} {get\_name}{}
1758    Module: \module{abstract\_2d\_finite\_volumes.domain}
1759
1760    Returns the name assigned to the domain by \code{set\_name}. If no name has been
1761    assigned, returns \code{`domain'}.
1762\end{methoddesc}
1763
1764\begin{methoddesc} {set\_datadir}{name}
1765    Module: \module{abstract\_2d\_finite\_volumes.domain}
1766
1767    Specifies the directory used for SWW files, assigning it to the
1768    pathname \code{name}. The default value, before
1769    \code{set\_datadir} has been run, is the value \code{default\_datadir}
1770    specified in \code{config.py}.
1771
1772    Since different operating systems use different formats for specifying pathnames,
1773    it is necessary to specify path separators using the Python code \code{os.sep}, rather than
1774    the operating-specific ones such as `$\slash$' or `$\backslash$'.
1775    For this to work you will need to include the statement \code{import os}
1776    in your code, before the first appearance of \code{set\_datadir}.
1777
1778    For example, to set the data directory to a subdirectory
1779    \code{data} of the directory \code{project}, you could use
1780    the statements:
1781
1782    {\small \begin{verbatim}
1783        import os
1784        domain.set_datadir{'project' + os.sep + 'data'}
1785    \end{verbatim}}
1786\end{methoddesc}
1787
1788\begin{methoddesc} {get\_datadir}{}
1789    Module: \module{abstract\_2d\_finite\_volumes.domain}
1790
1791    Returns the data directory set by \code{set\_datadir} or,
1792    if \code{set\_datadir} has not
1793    been run, returns the value \code{default\_datadir} specified in
1794    \code{config.py}.
1795\end{methoddesc}
1796
1797
1798\begin{methoddesc} {set\_minimum_allowed_height}{}
1799    Module: \module{shallow\_water.shallow\_water\_domain}
1800
1801    Set the minimum depth (in meters) that will be recognised in
1802    the numerical scheme (including limiters and flux computations)
1803
1804    Default value is $10^{-3}$ m, but by setting this to a greater value,
1805    e.g.\ for large scale simulations, the computation time can be
1806    significantly reduced.
1807\end{methoddesc}
1808
1809
1810\begin{methoddesc} {set\_minimum_storable_height}{}
1811    Module: \module{shallow\_water.shallow\_water\_domain}
1812
1813    Sets the minimum depth that will be recognised when writing
1814    to an sww file. This is useful for removing thin water layers
1815    that seems to be caused by friction creep.
1816\end{methoddesc}
1817
1818
1819\begin{methoddesc} {set\_maximum_allowed_speed}{}
1820    Module: \module{shallow\_water.shallow\_water\_domain}
1821
1822    Set the maximum particle speed that is allowed in water
1823    shallower than minimum_allowed_height. This is useful for
1824    controlling speeds in very thin layers of water and at the same time
1825    allow some movement avoiding pooling of water.
1826\end{methoddesc}
1827
1828
1829\begin{methoddesc} {set\_time}{time=0.0}
1830    Module: \module{abstract\_2d\_finite\_volumes.domain}
1831
1832    Sets the initial time, in seconds, for the simulation. The
1833    default is 0.0.
1834\end{methoddesc}
1835
1836\begin{methoddesc} {set\_default\_order}{n}
1837    Sets the default (spatial) order to the value specified by
1838    \code{n}, which must be either 1 or 2. (Assigning any other value
1839    to \code{n} will cause an error.)
1840\end{methoddesc}
1841
1842
1843\begin{methoddesc} {set\_store\_vertices\_uniquely}{flag}
1844Decide whether vertex values should be stored uniquely as
1845computed in the model or whether they should be reduced to one
1846value per vertex using averaging.
1847
1848Triangles stored in the sww file can be discontinuous reflecting
1849the internal representation of the finite-volume scheme
1850(this is a feature allowing for arbitrary steepness).
1851However, for visual purposes and also for use with \code{Field\_boundary}
1852(and \code{File\_boundary}) it is often desirable to store triangles
1853with values at each vertex point as the average of the potentially
1854discontinuous numbers found at vertices of different triangles sharing the
1855same vertex location.
1856
1857Storing one way or the other is controlled in ANUGA through the method
1858\code{domain.store\_vertices\_uniquely}. Options are
1859\begin{itemize}
1860  \item \code{domain.store\_vertices\_uniquely(True)}: Allow discontinuities in the sww file
1861  \item \code{domain.store\_vertices\_uniquely(False)}: (Default).
1862  Average values
1863  to ensure continuity in sww file. The latter also makes for smaller
1864  sww files.
1865\end{itemize}
1866
1867\end{methoddesc}
1868
1869
1870% Structural methods
1871\begin{methoddesc}{get\_nodes}{absolute=False}
1872    Return x,y coordinates of all nodes in mesh.
1873
1874    The nodes are ordered in an Nx2 array where N is the number of nodes.
1875    This is the same format they were provided in the constructor
1876    i.e. without any duplication.
1877
1878    Boolean keyword argument absolute determines whether coordinates
1879    are to be made absolute by taking georeference into account
1880    Default is False as many parts of ANUGA expects relative coordinates.
1881\end{methoddesc}
1882
1883
1884\begin{methoddesc}{get\_vertex_coordinates}{absolute=False}
1885
1886    Return vertex coordinates for all triangles.
1887
1888    Return all vertex coordinates for all triangles as a 3*M x 2 array
1889    where the jth vertex of the ith triangle is located in row 3*i+j and
1890    M the number of triangles in the mesh.
1891
1892    Boolean keyword argument absolute determines whether coordinates
1893    are to be made absolute by taking georeference into account
1894    Default is False as many parts of ANUGA expects relative coordinates.
1895\end{methoddesc}
1896
1897
1898\begin{methoddesc}{get\_triangles}{indices=None}
1899
1900        Return Mx3 integer array where M is the number of triangles.
1901        Each row corresponds to one triangle and the three entries are
1902        indices into the mesh nodes which can be obtained using the method
1903        get\_nodes()
1904
1905        Optional argument, indices is the set of triangle ids of interest.
1906\end{methoddesc}
1907
1908\begin{methoddesc}{get\_disconnected\_triangles}{}
1909
1910Get mesh based on nodes obtained from get_vertex_coordinates.
1911
1912        Return array Mx3 array of integers where each row corresponds to
1913        a triangle. A triangle is a triplet of indices into
1914        point coordinates obtained from get_vertex_coordinates and each
1915        index appears only once.\\
1916
1917        This provides a mesh where no triangles share nodes
1918        (hence the name disconnected triangles) and different
1919        nodes may have the same coordinates.\\
1920
1921        This version of the mesh is useful for storing meshes with
1922        discontinuities at each node and is e.g. used for storing
1923        data in sww files.\\
1924
1925        The triangles created will have the format
1926
1927    {\small \begin{verbatim}
1928        [[0,1,2],
1929         [3,4,5],
1930         [6,7,8],
1931         ...
1932         [3*M-3 3*M-2 3*M-1]]
1933     \end{verbatim}}
1934\end{methoddesc}
1935
1936
1937
1938%%%%%%
1939\section{Initial Conditions}
1940\label{sec:Initial Conditions}
1941In standard usage of partial differential equations, initial conditions
1942refers to the values associated to the system variables (the conserved
1943quantities here) for \code{time = 0}. In setting up a scenario script
1944as described in Sections \ref{sec:simpleexample} and \ref{sec:realdataexample},
1945\code{set_quantity} is used to define the initial conditions of variables
1946other than the conserved quantities, such as friction. Here, we use the terminology
1947of initial conditions to refer to initial values for variables which need
1948prescription to solve the shallow water wave equation. Further, it must be noted
1949that \code{set_quantity} does not necessarily have to be used in the initial
1950condition setting; it can be used at any time throughout the simulation.
1951
1952\begin{methoddesc}{set\_quantity}{name,
1953    numeric = None,
1954    quantity = None,
1955    function = None,
1956    geospatial_data = None,
1957    filename = None,
1958    attribute_name = None,
1959    alpha = None,
1960    location = 'vertices',
1961    indices = None,
1962    verbose = False,
1963    use_cache = False}
1964  Module: \module{abstract\_2d\_finite\_volumes.domain}
1965  (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values})
1966
1967This function is used to assign values to individual quantities for a
1968domain. It is very flexible and can be used with many data types: a
1969statement of the form \code{domain.set\_quantity(name, x)} can be used
1970to define a quantity having the name \code{name}, where the other
1971argument \code{x} can be any of the following:
1972
1973\begin{itemize}
1974\item a number, in which case all vertices in the mesh gets that for
1975the quantity in question.
1976\item a list of numbers or a Numeric array ordered the same way as the mesh vertices.
1977\item a function (e.g.\ see the samples introduced in Chapter 2)
1978\item an expression composed of other quantities and numbers, arrays, lists (for
1979example, a linear combination of quantities, such as
1980\code{domain.set\_quantity('stage','elevation'+x))}
1981\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.
1982\item a geospatial dataset (See Section \ref{sec:geospatial}).
1983Optional argument attribute\_name applies here as with files.
1984\end{itemize}
1985
1986
1987Exactly one of the arguments
1988  numeric, quantity, function, points, filename
1989must be present.
1990
1991
1992Set quantity will look at the type of the second argument (\code{numeric}) and
1993determine what action to take.
1994
1995Values can also be set using the appropriate keyword arguments.
1996If 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)}
1997are all equivalent.
1998
1999
2000Other optional arguments are
2001\begin{itemize}
2002\item \code{indices} which is a list of ids of triangles to which set\_quantity should apply its assignment of values.
2003\item \code{location} determines which part of the triangles to assign
2004  to. Options are 'vertices' (default), 'edges', 'unique vertices', and 'centroids'.
2005\end{itemize}
2006
2007%%%
2008\anuga provides a number of predefined initial conditions to be used
2009with \code{set\_quantity}. See for example callable object
2010\code{slump\_tsunami} below.
2011
2012\end{methoddesc}
2013
2014
2015
2016
2017\begin{funcdesc}{set_region}{tag, quantity, X, location='vertices'}
2018  Module: \module{abstract\_2d\_finite\_volumes.domain}
2019
2020  (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values})
2021
2022This function is used to assign values to individual quantities given
2023a regional tag.   It is similar to \code{set\_quantity}.
2024For example, if in the mesh-generator a regional tag of 'ditch' was
2025used, set\_region can be used to set elevation of this region to
2026-10m. X is the constant or function to be applied to the quantity,
2027over the tagged region.  Location describes how the values will be
2028applied.  Options are 'vertices' (default), 'edges', 'unique
2029vertices', and 'centroids'.
2030
2031This method can also be called with a list of region objects.  This is
2032useful for adding quantities in regions, and having one quantity
2033value based on another quantity. See  \module{abstract\_2d\_finite\_volumes.region} for
2034more details.
2035\end{funcdesc}
2036
2037
2038
2039
2040\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
2041                x0=0.0, y0=0.0, alpha=0.0,
2042                gravity=9.8, gamma=1.85,
2043                massco=1, dragco=1, frictionco=0, psi=0,
2044                dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
2045                domain=None,
2046                verbose=False}
2047Module: \module{shallow\_water.smf}
2048
2049This function returns a callable object representing an initial water
2050displacement generated by a submarine sediment failure. These failures can take the form of
2051a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead.
2052
2053The arguments include as a minimum, the slump or slide length, the water depth to the centre of sediment
2054mass, and the bathymetric slope. Other slump or slide parameters can be included if they are known.
2055\end{funcdesc}
2056
2057
2058%%%
2059\begin{funcdesc}{file\_function}{filename,
2060    domain = None,
2061    quantities = None,
2062    interpolation_points = None,
2063    verbose = False,
2064    use_cache = False}
2065Module: \module{abstract\_2d\_finite\_volumes.util}
2066
2067Reads the time history of spatial data for
2068specified interpolation points from a NetCDF file (\code{filename})
2069and returns
2070a callable object. \code{filename} could be a \code{sww} file.
2071Returns interpolated values based on the input
2072file using the underlying \code{interpolation\_function}.
2073
2074\code{quantities} is either the name of a single quantity to be
2075interpolated or a list of such quantity names. In the second case, the resulting
2076function will return a tuple of values---one for each quantity.
2077
2078\code{interpolation\_points} is a list of absolute coordinates or a
2079geospatial object
2080for points at which values are sought.
2081
2082The model time stored within the file function can be accessed using
2083the method \code{f.get\_time()}
2084
2085
2086The underlying algorithm used is as follows:\\
2087Given a time series (i.e.\ a series of values associated with
2088different times), whose values are either just numbers or a set of
2089 numbers defined at the vertices of a triangular mesh (such as those
2090 stored in SWW files), \code{Interpolation\_function} is used to
2091 create a callable object that interpolates a value for an arbitrary
2092 time \code{t} within the model limits and possibly a point \code{(x,
2093 y)} within a mesh region.
2094
2095 The actual time series at which data is available is specified by
2096 means of an array \code{time} of monotonically increasing times. The
2097 quantities containing the values to be interpolated are specified in
2098 an array---or dictionary of arrays (used in conjunction with the
2099 optional argument \code{quantity\_names}) --- called
2100 \code{quantities}. The optional arguments \code{vertex\_coordinates}
2101 and \code{triangles} represent the spatial mesh associated with the
2102 quantity arrays. If omitted the function created by
2103 \code{Interpolation\_function} will be a function of \code{t} only.
2104
2105 Since, in practice, values need to be computed at specified points,
2106 the syntax allows the user to specify, once and for all, a list
2107 \code{interpolation\_points} of points at which values are required.
2108 In this case, the function may be called using the form \code{f(t,
2109 id)}, where \code{id} is an index for the list
2110 \code{interpolation\_points}.
2111
2112
2113\end{funcdesc}
2114
2115%%%
2116%% \begin{classdesc}{Interpolation\_function}{self,
2117%%     time,
2118%%     quantities,
2119%%     quantity_names = None,
2120%%     vertex_coordinates = None,
2121%%     triangles = None,
2122%%     interpolation_points = None,
2123%%     verbose = False}
2124%% Module: \module{abstract\_2d\_finite\_volumes.least\_squares}
2125
2126%% Given a time series (i.e.\ a series of values associated with
2127%% different times), whose values are either just numbers or a set of
2128%% numbers defined at the vertices of a triangular mesh (such as those
2129%% stored in SWW files), \code{Interpolation\_function} is used to
2130%% create a callable object that interpolates a value for an arbitrary
2131%% time \code{t} within the model limits and possibly a point \code{(x,
2132%% y)} within a mesh region.
2133
2134%% The actual time series at which data is available is specified by
2135%% means of an array \code{time} of monotonically increasing times. The
2136%% quantities containing the values to be interpolated are specified in
2137%% an array---or dictionary of arrays (used in conjunction with the
2138%% optional argument \code{quantity\_names}) --- called
2139%% \code{quantities}. The optional arguments \code{vertex\_coordinates}
2140%% and \code{triangles} represent the spatial mesh associated with the
2141%% quantity arrays. If omitted the function created by
2142%% \code{Interpolation\_function} will be a function of \code{t} only.
2143
2144%% Since, in practice, values need to be computed at specified points,
2145%% the syntax allows the user to specify, once and for all, a list
2146%% \code{interpolation\_points} of points at which values are required.
2147%% In this case, the function may be called using the form \code{f(t,
2148%% id)}, where \code{id} is an index for the list
2149%% \code{interpolation\_points}.
2150
2151%% \end{classdesc}
2152
2153%%%
2154%\begin{funcdesc}{set\_region}{functions}
2155%[Low priority. Will be merged into set\_quantity]
2156
2157%Module:\module{abstract\_2d\_finite\_volumes.domain}
2158%\end{funcdesc}
2159
2160
2161
2162%%%%%%
2163\section{Boundary Conditions}\index{boundary conditions}
2164
2165\anuga provides a large number of predefined boundary conditions,
2166represented by objects such as \code{Reflective\_boundary(domain)} and
2167\code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples
2168in Chapter 2. Alternatively, you may prefer to ``roll your own'',
2169following the method explained in Section \ref{sec:roll your own}.
2170
2171These boundary objects may be used with the function \code{set\_boundary} described below
2172to assign boundary conditions according to the tags used to label boundary segments.
2173
2174\begin{methoddesc}{set\_boundary}{boundary_map}
2175Module: \module{abstract\_2d\_finite\_volumes.domain}
2176
2177This function allows you to assign a boundary object (corresponding to a
2178pre-defined or user-specified boundary condition) to every boundary segment that
2179has been assigned a particular tag.
2180
2181This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects
2182and whose keys are the symbolic tags.
2183
2184\end{methoddesc}
2185
2186\begin{methoddesc} {get\_boundary\_tags}{}
2187Module: \module{abstract\_2d\_finite\_volumes.domain}
2188
2189Returns a list of the available boundary tags.
2190\end{methoddesc}
2191
2192%%%
2193\subsection{Predefined boundary conditions}
2194
2195\begin{classdesc}{Reflective\_boundary}{Boundary}
2196Module: \module{shallow\_water}
2197
2198Reflective boundary returns same conserved quantities as those present in
2199the neighbouring volume but reflected.
2200
2201This class is specific to the shallow water equation as it works with the
2202momentum quantities assumed to be the second and third conserved quantities.
2203\end{classdesc}
2204
2205%%%
2206\begin{classdesc}{Transmissive\_boundary}{domain = None}
2207Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
2208
2209A transmissive boundary returns the same conserved quantities as
2210those present in the neighbouring volume.
2211
2212The underlying domain must be specified when the boundary is instantiated.
2213\end{classdesc}
2214
2215%%%
2216\begin{classdesc}{Dirichlet\_boundary}{conserved_quantities=None}
2217Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
2218
2219A Dirichlet boundary returns constant values for each of conserved
2220quantities. In the example of \code{Dirichlet\_boundary([0.2, 0.0, 0.0])},
2221the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and
2222\code{ymomentum} at the boundary are set to 0.0. The list must contain
2223a value for each conserved quantity.
2224\end{classdesc}
2225
2226%%%
2227\begin{classdesc}{Time\_boundary}{domain = None, f = None}
2228Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
2229
2230A time-dependent boundary returns values for the conserved
2231quantities as a function \code{f(t)} of time. The user must specify
2232the domain to get access to the model time.
2233\end{classdesc}
2234
2235%%%
2236\begin{classdesc}{File\_boundary}{Boundary}
2237Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions}
2238
2239This method may be used if the user wishes to apply a SWW file or
2240a time series file to a boundary segment or segments.
2241The boundary values are obtained from a file and interpolated to the
2242appropriate segments for each conserved quantity.
2243\end{classdesc}
2244
2245
2246
2247%%%
2248\begin{classdesc}{Transmissive\_Momentum\_Set\_Stage\_boundary}{Boundary}
2249Module: \module{shallow\_water}
2250
2251This boundary returns same momentum conserved quantities as
2252those present in its neighbour volume but sets stage as in a Time\_boundary.
2253The underlying domain must be specified when boundary is instantiated
2254
2255This type of boundary is useful when stage is known at the boundary as a
2256function of time, but momenta (or speeds) aren't.
2257
2258This class is specific to the shallow water equation as it works with the
2259momentum quantities assumed to be the second and third conserved quantities.
2260\end{classdesc}
2261
2262
2263\begin{classdesc}{Dirichlet\_Discharge\_boundary}{Boundary}
2264Module: \module{shallow\_water}
2265
2266Sets stage (stage0)
2267Sets momentum (wh0) in the inward normal direction.
2268\end{classdesc}
2269
2270
2271
2272\subsection{User-defined boundary conditions}
2273\label{sec:roll your own}
2274
2275All boundary classes must inherit from the generic boundary class
2276\code{Boundary} and have a method called \code{evaluate} which must
2277take as inputs \code{self, vol\_id, edge\_id} where self refers to the
2278object itself and vol\_id and edge\_id are integers referring to
2279particular edges. The method must return a list of three floating point
2280numbers representing values for \code{stage},
2281\code{xmomentum} and \code{ymomentum}, respectively.
2282
2283The constructor of a particular boundary class may be used to specify
2284particular values or flags to be used by the \code{evaluate} method.
2285Please refer to the source code for the existing boundary conditions
2286for examples of how to implement boundary conditions.
2287
2288
2289
2290%\section{Forcing Functions}
2291%
2292%\anuga provides a number of predefined forcing functions to be used with .....
2293
2294
2295
2296
2297\section{Evolution}\index{evolution}
2298
2299  \begin{methoddesc}{evolve}{yieldstep = None, finaltime = None, duration = None, skip_initial_step = False}
2300
2301  Module: \module{abstract\_2d\_finite\_volumes.domain}
2302
2303  This function (a method of \class{domain}) is invoked once all the
2304  preliminaries have been completed, and causes the model to progress
2305  through successive steps in its evolution, storing results and
2306  outputting statistics whenever a user-specified period
2307  \code{yieldstep} is completed (generally during this period the
2308  model will evolve through several steps internally
2309  as the method forces the water speed to be calculated
2310  on successive new cells). The user
2311  specifies the total time period over which the evolution is to take
2312  place, by specifying values (in seconds) for either \code{duration}
2313  or \code{finaltime}, as well as the interval in seconds after which
2314  results are to be stored and statistics output.
2315
2316  You can include \method{evolve} in a statement of the type:
2317
2318  {\small \begin{verbatim}
2319      for t in domain.evolve(yieldstep, finaltime):
2320          <Do something with domain and t>
2321  \end{verbatim}}
2322
2323  \end{methoddesc}
2324
2325
2326
2327\subsection{Diagnostics}
2328\label{sec:diagnostics}
2329
2330
2331  \begin{funcdesc}{statistics}{}
2332  Module: \module{abstract\_2d\_finite\_volumes.domain}
2333
2334  \end{funcdesc}
2335
2336  \begin{funcdesc}{timestepping\_statistics}{}
2337  Module: \module{abstract\_2d\_finite\_volumes.domain}
2338
2339  Returns a string of the following type for each
2340  timestep:
2341
2342  \code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12
2343  (12)}
2344
2345  Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
2346  the number of first-order steps, respectively.\\
2347
2348  The optional keyword argument \code{track_speeds=True} will
2349  generate a histogram of speeds generated by each triangle. The
2350  speeds relate to the size of the timesteps used by ANUGA and
2351  this diagnostics may help pinpoint problem areas where excessive speeds
2352  are generated.
2353
2354  \end{funcdesc}
2355
2356
2357  \begin{funcdesc}{boundary\_statistics}{quantities = None, tags = None}
2358  Module: \module{abstract\_2d\_finite\_volumes.domain}
2359
2360  Returns a string of the following type when \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}:
2361
2362  {\small \begin{verbatim}
2363 Boundary values at time 0.5000:
2364    top:
2365        stage in [ -0.25821218,  -0.02499998]
2366    bottom:
2367        stage in [ -0.27098821,  -0.02499974]
2368  \end{verbatim}}
2369
2370  \end{funcdesc}
2371
2372
2373  \begin{funcdesc}{get\_quantity}{name, location='vertices', indices = None}
2374  Module: \module{abstract\_2d\_finite\_volumes.domain}
2375
2376  Allow access to individual quantities and their methods
2377
2378  \end{funcdesc}
2379
2380
2381  \begin{funcdesc}{set\_quantities\_to\_be\_monitored}{}
2382  Module: \module{abstract\_2d\_finite\_volumes.domain}
2383
2384  Selects quantities and derived quantities for which extrema attained at internal timesteps
2385  will be collected.
2386
2387  Information can be tracked in the evolve loop by printing \code{quantity\_statistics} and
2388  collected data will be stored in the sww file.
2389
2390  Optional parameters \code{polygon} and \code{time\_interval} may be specified to restrict the
2391  extremum computation.
2392  \end{funcdesc}
2393
2394  \begin{funcdesc}{quantity\_statistics}{}
2395  Module: \module{abstract\_2d\_finite\_volumes.domain}
2396
2397  Reports on extrema attained by selected quantities.
2398
2399  Returns a string of the following type for each
2400  timestep:
2401
2402  \begin{verbatim}
2403  Monitored quantities at time 1.0000:
2404    stage-elevation:
2405      values since time = 0.00 in [0.00000000, 0.30000000]
2406      minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333)
2407      maximum attained at time = 0.00000000, location = (0.83333333, 0.16666667)
2408    ymomentum:
2409      values since time = 0.00 in [0.00000000, 0.06241221]
2410      minimum attained at time = 0.00000000, location = (0.33333333, 0.16666667)
2411      maximum attained at time = 0.22472667, location = (0.83333333, 0.66666667)
2412    xmomentum:
2413      values since time = 0.00 in [-0.06062178, 0.47886313]
2414      minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333)
2415      maximum attained at time = 0.35103646, location = (0.83333333, 0.16666667)
2416  \end{verbatim}
2417
2418  The quantities (and derived quantities) listed here must be selected at model
2419  initialisation using the method \code{domain.set_quantities_to_be_monitored}.\\
2420
2421  The optional keyword argument \code{precision='\%.4f'} will
2422  determine the precision used for floating point values in the output.
2423  This diagnostics helps track extrema attained by the selected quantities
2424  at every internal timestep.
2425
2426  These values are also stored in the sww file for post processing.
2427
2428  \end{funcdesc}
2429
2430
2431
2432  \begin{funcdesc}{get\_values}{location='vertices', indices = None}
2433  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2434
2435  Extract values for quantity as an array
2436
2437  \end{funcdesc}
2438
2439
2440  \begin{funcdesc}{get\_integral}{}
2441  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2442
2443  Return computed integral over entire domain for this quantity
2444
2445  \end{funcdesc}
2446
2447
2448
2449
2450  \begin{funcdesc}{get\_maximum\_value}{indices = None}
2451  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2452
2453  Return maximum value of quantity (on centroids)
2454
2455  Optional argument indices is the set of element ids that
2456  the operation applies to. If omitted all elements are considered.
2457
2458  We do not seek the maximum at vertices as each vertex can
2459  have multiple values - one for each triangle sharing it.
2460  \end{funcdesc}
2461
2462
2463
2464  \begin{funcdesc}{get\_maximum\_location}{indices = None}
2465  Module: \module{abstract\_2d\_finite\_volumes.quantity}
2466
2467  Return location of maximum value of quantity (on centroids)
2468
2469  Optional argument indices is the set of element ids that
2470  the operation applies to.
2471
2472  We do not seek the maximum at vertices as each vertex can
2473  have multiple values - one for each triangle sharing it.
2474
2475  If there are multiple cells with same maximum value, the
2476  first cell encountered in the triangle array is returned.
2477  \end{funcdesc}
2478
2479
2480
2481  \begin{funcdesc}{get\_wet\_elements}{indices=None}
2482  Module: \module{shallow\_water.shallow\_water\_domain}
2483
2484  Return indices for elements where h $>$ minimum_allowed_height
2485  Optional argument indices is the set of element ids that the operation applies to.
2486  \end{funcdesc}
2487
2488
2489  \begin{funcdesc}{get\_maximum\_inundation\_elevation}{indices=None}
2490  Module: \module{shallow\_water.shallow\_water\_domain}
2491
2492  Return highest elevation where h $>$ 0.\\
2493  Optional argument indices is the set of element ids that the operation applies to.\\
2494
2495  Example to find maximum runup elevation:\\
2496     z = domain.get_maximum_inundation_elevation()
2497  \end{funcdesc}
2498
2499
2500  \begin{funcdesc}{get\_maximum\_inundation\_location}{indices=None}
2501  Module: \module{shallow\_water.shallow\_water\_domain}
2502
2503  Return location (x,y) of highest elevation where h $>$ 0.\\
2504  Optional argument indices is the set of element ids that the operation applies to.\\
2505
2506  Example to find maximum runup location:\\
2507     x, y = domain.get_maximum_inundation_location()
2508  \end{funcdesc}
2509
2510
2511\section{Queries of SWW model output files}
2512After a model has been run, it is often useful to extract various information from the sww
2513output file (see Section \ref{sec:sww format}). This is typically more convenient than using the
2514diagnostics described in Section \ref{sec:diagnostics} which rely on the model running - something
2515that can be very time consuming. The sww files are easy and quick to read and offer much information
2516about the model results such as runup heights, time histories of selected quantities,
2517flow through cross sections and much more.
2518
2519\begin{funcdesc}{get\_maximum\_inundation\_elevation}{filename, polygon=None,
2520    time_interval=None, verbose=False}
2521  Module: \module{shallow\_water.data\_manager}
2522
2523  Return highest elevation where depth is positive ($h > 0$)
2524
2525  Example to find maximum runup elevation:\\
2526  max_runup = get_maximum_inundation_elevation(filename,
2527  polygon=None,
2528  time_interval=None,
2529  verbose=False)
2530
2531
2532  filename is a NetCDF sww file containing ANUGA model output.
2533  Optional arguments polygon and time_interval restricts the maximum runup calculation
2534  to a points that lie within the specified polygon and time interval.
2535
2536  If no inundation is found within polygon and time_interval the return value
2537  is None signifying "No Runup" or "Everything is dry".
2538
2539  See doc string for general function get_maximum_inundation_data for details.
2540\end{funcdesc}
2541
2542
2543\begin{funcdesc}{get\_maximum\_inundation\_location}{filename, polygon=None,
2544    time_interval=None, verbose=False}
2545  Module: \module{shallow\_water.data\_manager}
2546
2547  Return location (x,y) of highest elevation where depth is positive ($h > 0$)
2548
2549  Example to find maximum runup location:\\
2550  max_runup_location = get_maximum_inundation_location(filename,
2551  polygon=None,
2552  time_interval=None,
2553  verbose=False)
2554
2555
2556  filename is a NetCDF sww file containing ANUGA model output.
2557  Optional arguments polygon and time_interval restricts the maximum runup calculation
2558  to a points that lie within the specified polygon and time interval.
2559
2560  If no inundation is found within polygon and time_interval the return value
2561  is None signifying "No Runup" or "Everything is dry".
2562
2563  See doc string for general function get_maximum_inundation_data for details.
2564\end{funcdesc}
2565
2566
2567\begin{funcdesc}{sww2timeseries}{swwfiles, gauge_filename, production_dirs, report = None, reportname = None,
2568plot_quantity = None, generate_fig = False, surface = None, time_min = None, time_max = None, time_thinning = 1,
2569time_unit = None, title_on = None, use_cache = False, verbose = False}
2570
2571  Module: \module{anuga.abstract\_2d\_finite\_volumes.util}
2572
2573  Return csv files for the location in the \code{gauge_filename} and can also return plots of them
2574
2575  See doc string for general function sww2timeseries for details.
2576
2577\end{funcdesc}
2578
2579
2580
2581\section{Other}
2582
2583  \begin{funcdesc}{domain.create\_quantity\_from\_expression}{string}
2584
2585  Handy for creating derived quantities on-the-fly as for example
2586  \begin{verbatim}
2587  Depth = domain.create_quantity_from_expression('stage-elevation')
2588
2589  exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5')
2590  Absolute_momentum = domain.create_quantity_from_expression(exp)
2591  \end{verbatim}
2592
2593  %See also \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
2594  \end{funcdesc}
2595
2596
2597
2598
2599
2600%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2601%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2602
2603\chapter{\anuga System Architecture}
2604
2605
2606\section{File Formats}
2607\label{sec:file formats}
2608
2609\anuga makes use of a number of different file formats. The
2610following table lists all these formats, which are described in more
2611detail in the paragraphs below.
2612
2613\bigskip
2614
2615\begin{center}
2616
2617\begin{tabular}{|ll|}  \hline
2618
2619\textbf{Extension} & \textbf{Description} \\
2620\hline\hline
2621
2622\code{.sww} & NetCDF format for storing model output
2623\code{f(t,x,y)}\\
2624
2625\code{.tms} & NetCDF format for storing time series \code{f(t)}\\
2626
2627\code{.csv/.txt} & ASCII format called points csv for storing
2628arbitrary points and associated attributes\\
2629
2630\code{.pts} & NetCDF format for storing arbitrary points and
2631associated attributes\\
2632
2633\code{.asc} & ASCII format of regular DEMs as output from ArcView\\
2634
2635\code{.prj} & Associated ArcView file giving more metadata for
2636\code{.asc} format\\
2637
2638\code{.ers} & ERMapper header format of regular DEMs for ArcView\\
2639
2640\code{.dem} & NetCDF representation of regular DEM data\\
2641
2642\code{.tsh} & ASCII format for storing meshes and associated
2643boundary and region info\\
2644
2645\code{.msh} & NetCDF format for storing meshes and associated
2646boundary and region info\\
2647
2648\code{.nc} & Native ferret NetCDF format\\
2649
2650\code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
2651%\caption{File formats used by \anuga}
2652\end{tabular}
2653
2654
2655\end{center}
2656
2657The above table shows the file extensions used to identify the
2658formats of files. However, typically, in referring to a format we
2659capitalise the extension and omit the initial full stop---thus, we
2660refer, for example, to `SWW files' or `PRJ files'.
2661
2662\bigskip
2663
2664A typical dataflow can be described as follows:
2665
2666\subsection{Manually Created Files}
2667
2668\begin{tabular}{ll}
2669ASC, PRJ & Digital elevation models (gridded)\\
2670NC & Model outputs for use as boundary conditions (e.g. from MOST)
2671\end{tabular}
2672
2673\subsection{Automatically Created Files}
2674
2675\begin{tabular}{ll}
2676ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert
2677DEMs to native \code{.pts} file\\
2678
2679NC $\rightarrow$ SWW & Convert MOST boundary files to
2680boundary \code{.sww}\\
2681
2682PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\
2683
2684TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using
2685\code{animate}\\
2686
2687TSH + Boundary SWW $\rightarrow$ SWW & Simulation using
2688\code{\anuga}\\
2689
2690Polygonal mesh outline $\rightarrow$ & TSH or MSH
2691\end{tabular}
2692
2693
2694
2695
2696\bigskip
2697
2698\subsection{SWW and TMS Formats}
2699\label{sec:sww format}
2700
2701The SWW and TMS formats are both NetCDF formats, and are of key
2702importance for \anuga.
2703
2704An SWW file is used for storing \anuga output and therefore pertains
2705to a set of points and a set of times at which a model is evaluated.
2706It contains, in addition to dimension information, the following
2707variables:
2708
2709\begin{itemize}
2710    \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays
2711    \item \code{elevation}, a Numeric array storing bed-elevations
2712    \item \code{volumes}, a list specifying the points at the vertices of each of the
2713    triangles
2714    % Refer here to the example to be provided in describing the simple example
2715    \item \code{time}, a Numeric array containing times for model
2716    evaluation
2717\end{itemize}
2718
2719
2720The contents of an SWW file may be viewed using the anuga viewer
2721\code{animate}, which creates an on-screen geometric
2722representation. See section \ref{sec:animate} (page
2723\pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more
2724on \code{animate}.
2725
2726Alternatively, there are tools, such as \code{ncdump}, that allow
2727you to convert an NetCDF file into a readable format such as the
2728Class Definition Language (CDL). The following is an excerpt from a
2729CDL representation of the output file \file{runup.sww} generated
2730from running the simple example \file{runup.py} of
2731Chapter \ref{ch:getstarted}:
2732
2733\verbatiminput{examples/bedslopeexcerpt.cdl}
2734
2735The SWW format is used not only for output but also serves as input
2736for functions such as \function{file\_boundary} and
2737\function{file\_function}, described in Chapter \ref{ch:interface}.
2738
2739A TMS file is used to store time series data that is independent of
2740position.
2741
2742
2743\subsection{Mesh File Formats}
2744
2745A mesh file is a file that has a specific format suited to
2746triangular meshes and their outlines. A mesh file can have one of
2747two formats: it can be either a TSH file, which is an ASCII file, or
2748an MSH file, which is a NetCDF file. A mesh file can be generated
2749from the function \function{create\_mesh\_from\_regions} (see
2750Section \ref{sec:meshgeneration}) and used to initialise a domain.
2751
2752A mesh file can define the outline of the mesh---the vertices and
2753line segments that enclose the region in which the mesh is
2754created---and the triangular mesh itself, which is specified by
2755listing the triangles and their vertices, and the segments, which
2756are those sides of the triangles that are associated with boundary
2757conditions.
2758
2759In addition, a mesh file may contain `holes' and/or `regions'. A
2760hole represents an area where no mesh is to be created, while a
2761region is a labelled area used for defining properties of a mesh,
2762such as friction values.  A hole or region is specified by a point
2763and bounded by a number of segments that enclose that point.
2764
2765A mesh file can also contain a georeference, which describes an
2766offset to be applied to $x$ and $y$ values---eg to the vertices.
2767
2768
2769\subsection{Formats for Storing Arbitrary Points and Attributes}
2770
2771
2772A CSV/TXT file is used to store data representing
2773arbitrary numerical attributes associated with a set of points.
2774
2775The format for an CSV/TXT file is:\\
2776%\begin{verbatim}
2777
2778            first line:     \code{[column names]}\\
2779            other lines:  \code{[x value], [y value], [attributes]}\\
2780
2781            for example:\\
2782            \code{x, y, elevation, friction}\\
2783            \code{0.6, 0.7, 4.9, 0.3}\\
2784            \code{1.9, 2.8, 5, 0.3}\\
2785            \code{2.7, 2.4, 5.2, 0.3}
2786
2787        The delimiter is a comma. The first two columns are assumed to
2788        be x, y coordinates.
2789       
2790
2791A PTS file is a NetCDF representation of the data held in an points CSV
2792file. If the data is associated with a set of $N$ points, then the
2793data is stored using an $N \times 2$ Numeric array of float
2794variables for the points and an $N \times 1$ Numeric array for each
2795attribute.
2796
2797%\end{verbatim}
2798
2799\subsection{ArcView Formats}
2800
2801Files of the three formats ASC, PRJ and ERS are all associated with
2802data from ArcView.
2803
2804An ASC file is an ASCII representation of DEM output from ArcView.
2805It contains a header with the following format:
2806
2807\begin{tabular}{l l}
2808\code{ncols}      &   \code{753}\\
2809\code{nrows}      &   \code{766}\\
2810\code{xllcorner}  &   \code{314036.58727982}\\
2811\code{yllcorner}  & \code{6224951.2960092}\\
2812\code{cellsize}   & \code{100}\\
2813\code{NODATA_value} & \code{-9999}
2814\end{tabular}
2815
2816The remainder of the file contains the elevation data for each grid point
2817in the grid defined by the above information.
2818
2819A PRJ file is an ArcView file used in conjunction with an ASC file
2820to represent metadata for a DEM.
2821
2822
2823\subsection{DEM Format}
2824
2825A DEM file is a NetCDF representation of regular DEM data.
2826
2827
2828\subsection{Other Formats}
2829
2830
2831
2832
2833\subsection{Basic File Conversions}
2834\label{sec:basicfileconversions}
2835
2836  \begin{funcdesc}{sww2dem}{basename_in, basename_out = None,
2837            quantity = None,
2838            timestep = None,
2839            reduction = None,
2840            cellsize = 10,
2841            NODATA_value = -9999,
2842            easting_min = None,
2843            easting_max = None,
2844            northing_min = None,
2845            northing_max = None,
2846            expand_search = False,
2847            verbose = False,
2848            origin = None,
2849            datum = 'WGS84',
2850            format = 'ers'}
2851  Module: \module{shallow\_water.data\_manager}
2852
2853  Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
2854  ERS) of a desired grid size \code{cellsize} in metres.
2855  The easting and northing values are used if the user wished to clip the output
2856  file to a specified rectangular area. The \code{reduction} input refers to a function
2857  to reduce the quantities over all time step of the SWW file, example, maximum.
2858  \end{funcdesc}
2859
2860
2861  \begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
2862            easting_min=None, easting_max=None,
2863            northing_min=None, northing_max=None,
2864            use_cache=False, verbose=False}
2865  Module: \module{shallow\_water.data\_manager}
2866
2867  Takes DEM data (a NetCDF file representation of data from a regular Digital
2868  Elevation Model) and converts it to PTS format.
2869  \end{funcdesc}
2870
2871
2872
2873%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2874%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2875
2876\chapter{\anuga mathematical background}
2877\label{cd:mathematical background}
2878
2879\section{Introduction}
2880
2881This chapter outlines the mathematics underpinning \anuga.
2882
2883
2884
2885\section{Model}
2886\label{sec:model}
2887
2888The shallow water wave equations are a system of differential
2889conservation equations which describe the flow of a thin layer of
2890fluid over terrain. The form of the equations are:
2891\[
2892\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
2893x}+\frac{\partial \GG}{\partial y}=\SSS
2894\]
2895where $\UU=\left[ {{\begin{array}{*{20}c}
2896 h & {uh} & {vh} \\
2897\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
2898$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
2899entering the system are bed elevation $z$ and stage (absolute water
2900level) $w$, where the relation $w = z + h$ holds true at all times.
2901The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
2902by
2903\[
2904\EE=\left[ {{\begin{array}{*{20}c}
2905 {uh} \hfill \\
2906 {u^2h+gh^2/2} \hfill \\
2907 {uvh} \hfill \\
2908\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
2909 {vh} \hfill \\
2910 {vuh} \hfill \\
2911 {v^2h+gh^2/2} \hfill \\
2912\end{array} }} \right]
2913\]
2914and the source term (which includes gravity and friction) is given
2915by
2916\[
2917\SSS=\left[ {{\begin{array}{*{20}c}
2918 0 \hfill \\
2919 -{gh(z_{x} + S_{fx} )} \hfill \\
2920 -{gh(z_{y} + S_{fy} )} \hfill \\
2921\end{array} }} \right]
2922\]
2923where $S_f$ is the bed friction. The friction term is modelled using
2924Manning's resistance law
2925\[
2926S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
2927=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
2928\]
2929in which $\eta$ is the Manning resistance coefficient.
2930
2931As demonstrated in our papers, \cite{ZR1999,nielsen2005} these
2932equations provide an excellent model of flows associated with
2933inundation such as dam breaks and tsunamis.
2934
2935\section{Finite Volume Method}
2936\label{sec:fvm}
2937
2938We use a finite-volume method for solving the shallow water wave
2939equations \cite{ZR1999}. The study area is represented by a mesh of
2940triangular cells as in Figure~\ref{fig:mesh} in which the conserved
2941quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
2942in each volume are to be determined. The size of the triangles may
2943be varied within the mesh to allow greater resolution in regions of
2944particular interest.
2945
2946\begin{figure}
2947\begin{center}
2948\includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-five}
2949\caption{Triangular mesh used in our finite volume method. Conserved
2950quantities $h$, $uh$ and $vh$ are associated with the centroid of
2951each triangular cell.} \label{fig:mesh}
2952\end{center}
2953\end{figure}
2954
2955The equations constituting the finite-volume method are obtained by
2956integrating the differential conservation equations over each
2957triangular cell of the mesh. Introducing some notation we use $i$ to
2958refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
2959set of indices referring to the cells neighbouring the $i$th cell.
2960Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
2961the length of the edge between the $i$th and $j$th cells.
2962
2963By applying the divergence theorem we obtain for each volume an
2964equation which describes the rate of change of the average of the
2965conserved quantities within each cell, in terms of the fluxes across
2966the edges of the cells and the effect of the source terms. In
2967particular, rate equations associated with each cell have the form
2968$$
2969 \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
2970$$
2971where
2972\begin{itemize}
2973\item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
2974\item $\SSS_i$ is the source term associated with the $i$th cell,
2975and
2976\item $\HH_{ij}$ is the outward normal flux of
2977material across the \textit{ij}th edge.
2978\end{itemize}
2979
2980
2981%\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
2982%cells
2983%\item $m_{ij}$ is the midpoint of
2984%the \textit{ij}th edge,
2985%\item
2986%$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
2987%normal along the \textit{ij}th edge, and The
2988
2989The flux $\HH_{ij}$ is evaluated using a numerical flux function
2990$\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
2991water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
2992$$
2993H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
2994$$
2995
2996Then
2997$$
2998\HH_{ij}  = \HH(\UU_i(m_{ij}),
2999\UU_j(m_{ij}); \mathbf{n}_{ij})
3000$$
3001where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
3002$\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
3003\textit{ij}th edge. The function $\UU_i(x)$ for $x \in
3004T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
3005neighbouring  cells.
3006
3007We use a second order reconstruction to produce a piece-wise linear
3008function construction of the conserved quantities for  all $x \in
3009T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This
3010function is allowed to be discontinuous across the edges of the
3011cells, but the slope of this function is limited to avoid
3012artificially introduced oscillations.
3013
3014Godunov's method (see \cite{Toro1992}) involves calculating the
3015numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
3016solving the corresponding one dimensional Riemann problem normal to
3017the edge. We use the central-upwind scheme of \cite{KurNP2001} to
3018calculate an approximation of the flux across each edge.
3019
3020\begin{figure}
3021\begin{center}
3022\includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-reconstruct}
3023\caption{From the values of the conserved quantities at the centroid
3024of the cell and its neighbouring cells, a discontinuous piecewise
3025linear reconstruction of the conserved quantities is obtained.}
3026\label{fig:mesh:reconstruct}
3027\end{center}
3028\end{figure}
3029
3030In the computations presented in this paper we use an explicit Euler
3031time stepping method with variable timestepping adapted to the
3032observed CFL condition.
3033
3034
3035\section{Flux limiting}
3036
3037The shallow water equations are solved numerically using a
3038finite volume method on unstructured triangular grid.
3039The upwind central scheme due to Kurganov and Petrova is used as an
3040approximate Riemann solver for the computation of inviscid flux functions.
3041This makes it possible to handle discontinuous solutions.
3042
3043To alleviate the problems associated with numerical instabilities due to
3044small water depths near a wet/dry boundary we employ a new flux limiter that
3045ensures that unphysical fluxes are never encounted.
3046
3047
3048Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction,
3049$w$ the absolute water level (stage) and
3050$z$ the bed elevation. The latter are assumed to be relative to the
3051same height datum.
3052The conserved quantities tracked by ANUGA are momentum in the
3053$x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$)
3054and depth ($h = w-z$).
3055
3056The flux calculation requires access to the velocity vector $(u, v)$
3057where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively.
3058In the presence of very small water depths, these calculations become
3059numerically unreliable and will typically cause unphysical speeds.
3060
3061We have employed a flux limiter which replaces the calculations above with
3062the limited approximations.
3063\begin{equation}
3064  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h},
3065\end{equation}
3066where $h_0$ is a regularisation parameter that controls the minimal
3067magnitude of the denominator. Taking the limits we have for $\hat{u}$
3068\[
3069  \lim_{h \rightarrow 0} \hat{u} =
3070  \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0
3071\]
3072and
3073\[
3074  \lim_{h \rightarrow \infty} \hat{u} =
3075  \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u
3076\]
3077with similar results for $\hat{v}$.
3078
3079The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator)
3080\[
3081  1 - h_0/h^2 = 0
3082\]
3083or
3084\[
3085  h_0 = h^2
3086\]
3087
3088
3089ANUGA has a global parameter $H_0$ that controls the minimal depth which
3090is considered in the various equations. This parameter is typically set to
3091$10^{-3}$. Setting
3092\[
3093  h_0 = H_0^2
3094\]
3095provides a reasonable balance between accurracy and stability. In fact,
3096setting $h=H_0$ will scale the predicted speed by a factor of $0.5$:
3097\[
3098  \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0}
3099\]
3100In general, for multiples of the minimal depth $N H_0$ one obtains
3101\[
3102  \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} =
3103  \frac{\mu}{H_0 (1 + 1/N^2)}
3104\]
3105which converges quadratically to the true value with the multiple N.
3106
3107
3108%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.
3109
3110
3111
3112
3113
3114\section{Slope limiting}
3115A 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.
3116
3117However 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.
3118
3119
3120Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
3121let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.
3122Define the minimal depth across all vertices as $\hmin$ as
3123\[
3124  \hmin = \min_i h_i
3125\]
3126
3127Let $\tilde{w_i}$ be the stage obtained from a gradient limiter
3128limiting on stage only. The corresponding depth is then defined as
3129\[
3130  \tilde{h_i} = \tilde{w_i} - z_i
3131\]
3132We would use this limiter in deep water which we will define (somewhat boldly)
3133as
3134\[
3135  \hmin \ge \epsilon
3136\]
3137
3138
3139Similarly, let $\bar{w_i}$ be the stage obtained from a gradient
3140limiter limiting on depth respecting the bed slope.
3141The corresponding depth is defined as
3142\[
3143  \bar{h_i} = \bar{w_i} - z_i
3144\]
3145
3146
3147We introduce the concept of a balanced stage $w_i$ which is obtained as
3148the linear combination
3149
3150\[
3151  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
3152\]
3153or
3154\[
3155  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
3156\]
3157where $\alpha \in [0, 1]$.
3158
3159Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope
3160is ignored we have immediately that
3161\[
3162  \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0
3163\]
3164%where the maximal bed elevation range $dz$ is defined as
3165%\[
3166%  dz = \max_i |z_i - z|
3167%\]
3168
3169If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that
3170no negative depths occur. Formally, we will require that
3171\[
3172  \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i
3173\]
3174or
3175\begin{equation}
3176  \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i
3177  \label{eq:limiter bound}
3178\end{equation}
3179
3180There are two cases:
3181\begin{enumerate}
3182  \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage)
3183  vertex is at least as far away from the bed than the shallow water
3184  (limited using depth). In this case we won't need any contribution from
3185  $\bar{h_i}$ and can accept any $\alpha$.
3186
3187  E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to
3188  \[
3189    \tilde{h_i} > \epsilon
3190  \]
3191  whereas $\alpha=0$ yields
3192  \[
3193    \bar{h_i} > \epsilon
3194  \]
3195  all well and good.
3196  \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is
3197  closer to the bed than the shallow water vertex or even below the bed.
3198  In this case we need to find an $\alpha$ that will ensure a positive depth.
3199  Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one
3200  obtains the bound
3201  \[
3202    \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i
3203  \]
3204\end{enumerate}
3205
3206Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one
3207arrives at the definition
3208\[
3209  \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}}
3210\]
3211which will guarantee that no vertex 'cuts' through the bed. Finally, should
3212$\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting
3213$\alpha=0$ and similarly capping $\alpha$ at 1 just in case.
3214
3215%Furthermore,
3216%dropping the $\epsilon$ ensures that alpha is always positive and also
3217%provides a numerical safety {??)
3218
3219
3220
3221
3222
3223%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3224%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3225
3226\chapter{Basic \anuga Assumptions}
3227
3228
3229Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
3230If one wished to recreate scenarios prior to that date it must be done
3231using some relative time (e.g. 0).
3232
3233
3234All spatial data relates to the WGS84 datum (or GDA94) and has been
3235projected into UTM with false easting of 500000 and false northing of
32361000000 on the southern hemisphere (0 on the northern).
3237
3238It is assumed that all computations take place within one UTM zone and
3239all locations must consequently be specified in Cartesian coordinates
3240(eastings, northings) or (x,y) where the unit is metres.
3241
3242DEMs, meshes and boundary conditions can have different origins within
3243one UTM zone. However, the computation will use that of the mesh for
3244numerical stability.
3245
3246When generating a mesh it is assumed that polygons do not cross.
3247Having polygons tht cross can cause the mesh generation to fail or bad
3248meshes being produced.
3249
3250
3251%OLD
3252%The dataflow is: (See data_manager.py and from scenarios)
3253%
3254%
3255%Simulation scenarios
3256%--------------------%
3257%%
3258%
3259%Sub directories contain scrips and derived files for each simulation.
3260%The directory ../source_data contains large source files such as
3261%DEMs provided externally as well as MOST tsunami simulations to be used
3262%as boundary conditions.
3263%
3264%Manual steps are:
3265%  Creation of DEMs from argcview (.asc + .prj)
3266%  Creation of mesh from pmesh (.tsh)
3267%  Creation of tsunami simulations from MOST (.nc)
3268%%
3269%
3270%Typical scripted steps are%
3271%
3272%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
3273%                   native dem and pts formats%
3274%
3275%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
3276%                  as boundary condition%
3277%
3278%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
3279%                   smoothing. The outputs are tsh files with elevation data.%
3280%
3281%  run_simulation.py: Use the above together with various parameters to
3282%                     run inundation simulation.
3283
3284
3285%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3286%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3287
3288\appendix
3289
3290\chapter{Supporting Tools}
3291\label{ch:supportingtools}
3292
3293This section describes a number of supporting tools, supplied with \anuga, that offer a
3294variety of types of functionality and enhance the basic capabilities of \anuga.
3295
3296\section{caching}
3297\label{sec:caching}
3298
3299The \code{cache} function is used to provide supervised caching of function
3300results. A Python function call of the form
3301
3302      {\small \begin{verbatim}
3303      result = func(arg1,...,argn)
3304      \end{verbatim}}
3305
3306  can be replaced by
3307
3308      {\small \begin{verbatim}
3309      from caching import cache
3310      result = cache(func,(arg1,...,argn))
3311      \end{verbatim}}
3312
3313  which returns the same output but reuses cached
3314  results if the function has been computed previously in the same context.
3315  \code{result} and the arguments can be simple types, tuples, list, dictionaries or
3316  objects, but not unhashable types such as functions or open file objects.
3317  The function \code{func} may be a member function of an object or a module.
3318
3319  This type of caching is particularly useful for computationally intensive
3320  functions with few frequently used combinations of input arguments. Note that
3321  if the inputs or output are very large caching may not save time because
3322  disc access may dominate the execution time.
3323
3324  If the function definition changes after a result has been cached, this will be
3325  detected by examining the functions \code{bytecode (co\_code, co\_consts,
3326  func\_defaults, co\_argcount)} and the function will be recomputed.
3327  However, caching will not detect changes in modules used by \code{func}.
3328  In this case cache must be cleared manually.
3329
3330  Options are set by means of the function \code{set\_option(key, value)},
3331  where \code{key} is a key associated with a
3332  Python dictionary \code{options}. This dictionary stores settings such as the name of
3333  the directory used, the maximum
3334  number of cached files allowed, and so on.
3335
3336  The \code{cache} function allows the user also to specify a list of dependent files. If any of these
3337  have been changed, the function is recomputed and the results stored again.
3338
3339  %Other features include support for compression and a capability to \ldots
3340
3341
3342   \textbf{USAGE:} \nopagebreak
3343
3344    {\small \begin{verbatim}
3345    result = cache(func, args, kwargs, dependencies, cachedir, verbose,
3346                   compression, evaluate, test, return_filename)
3347    \end{verbatim}}
3348
3349
3350\section{ANUGA viewer - animate}
3351\label{sec:animate}
3352 The output generated by \anuga may be viewed by
3353means of the visualisation tool \code{animate}, which takes the
3354\code{SWW} file output by \anuga and creates a visual representation
3355of the data. Examples may be seen in Figures \ref{fig:runupstart}
3356and \ref{fig:runup2}. To view an \code{SWW} file with
3357\code{animate} in the Windows environment, you can simply drag the
3358icon representing the file over an icon on the desktop for the
3359\code{animate} executable file (or a shortcut to it), or set up a
3360file association to make files with the extension \code{.sww} open
3361with \code{animate}. Alternatively, you can operate \code{animate}
3362from the command line, in both Windows and Linux environments.
3363
3364On successful operation, you will see an interactive moving-picture
3365display. You can use keys and the mouse to slow down, speed up or
3366stop the display, change the viewing position or carry out a number
3367of other simple operations. Help is also displayed when you press
3368the \code{h} key.
3369
3370The main keys operating the interactive screen are:\\
3371
3372\begin{center}
3373\begin{tabular}{|ll|}   \hline
3374
3375\code{w} & toggle wireframe \\
3376
3377space bar & start/stop\\
3378
3379up/down arrows & increase/decrease speed\\
3380
3381left/right arrows & direction in time \emph{(when running)}\\
3382& step through simulation \emph{(when stopped)}\\
3383
3384left mouse button & rotate\\
3385
3386middle mouse button & pan\\
3387
3388right mouse button & zoom\\  \hline
3389
3390\end{tabular}
3391\end{center}
3392
3393\vfill
3394
3395The following table describes how to operate animate from the command line:
3396
3397Usage: \code{animate [options] swwfile \ldots}\\  \nopagebreak
3398Options:\\  \nopagebreak
3399\begin{tabular}{ll}
3400  \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\
3401                                    & \code{HEAD\_MOUNTED\_DISPLAY}\\
3402  \code{--rgba} & Request a RGBA colour buffer visual\\
3403  \code{--stencil} & Request a stencil buffer visual\\
3404  \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
3405                                    & overridden by environmental variable\\
3406  \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD\_BUFFER | HORIZONTAL\_SPLIT |}\\
3407                                    & \code{VERTICAL\_SPLIT | LEFT\_EYE | RIGHT\_EYE |}\\
3408                                     & \code{ON | OFF} \\
3409  \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
3410  \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
3411\end{tabular}
3412
3413\begin{tabular}{ll}
3414  \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
3415  \code{-help} & Display this information\\
3416  \code{-hmax <float>} & Height above which transparency is set to
3417                                     \code{alphamax}\\
3418\end{tabular}
3419
3420\begin{tabular}{ll}
3421
3422  \code{-hmin <float>} & Height below which transparency is set to
3423                                     zero\\
3424\end{tabular}
3425
3426\begin{tabular}{ll}
3427  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
3428                                     up, default is overhead)\\
3429\end{tabular}
3430
3431\begin{tabular}{ll}
3432  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
3433
3434\end{tabular}
3435
3436\begin{tabular}{ll}
3437  \code{-movie <dirname>} & Save numbered images to named directory and
3438                                     quit\\
3439
3440  \code{-nosky} & Omit background sky\\
3441
3442
3443  \code{-scale <float>} & Vertical scale factor\\
3444  \code{-texture <file>} & Image to use for bedslope topography\\
3445  \code{-tps <rate>} & Timesteps per second\\
3446  \code{-version} & Revision number and creation (not compile)
3447                                     date\\
3448\end{tabular}
3449
3450\section{utilities/polygons}
3451
3452  \declaremodule{standard}{utilities.polygon}
3453  \refmodindex{utilities.polygon}
3454
3455  \begin{classdesc}{Polygon\_function}{regions, default=0.0, geo_reference=None}
3456  Module: \code{utilities.polygon}
3457
3458  Creates a callable object that returns one of a specified list of values when
3459  evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the
3460  point belongs to. The parameter \code{regions} is a list of pairs
3461  \code{(P, v)}, where each \code{P} is a polygon and each \code{v}
3462  is either a constant value or a function of coordinates \code{x}
3463  and \code{y}, specifying the return value for a point inside \code{P}. The
3464  optional parameter \code{default} may be used to specify a value
3465  (or a function)
3466  for a point not lying inside any of the specified polygons. When a
3467  point lies in more than one polygon, the return value is taken to
3468  be the value for whichever of these polygon appears later in the
3469  list.
3470  %FIXME (Howard): CAN x, y BE VECTORS?
3471  The optional parameter geo\_reference refers to the status of points
3472  that are passed into the function. Typically they will be relative to
3473  some origin. In ANUGA, a typical call will take the form:
3474  {\small \begin{verbatim}
3475     set_quantity('elevation',
3476                  Polygon_function([(P1, v1), (P2, v2)],
3477                                   default=v3,
3478                                   geo_reference=domain.geo_reference))
3479  \end{verbatim}}
3480 
3481
3482  \end{classdesc}
3483
3484  \begin{funcdesc}{read\_polygon}{filename}
3485  Module: \code{utilities.polygon}
3486
3487  Reads the specified file and returns a polygon. Each
3488  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
3489  as coordinates of one vertex of the polygon.
3490  \end{funcdesc}
3491
3492  \begin{funcdesc}{populate\_polygon}{polygon, number_of_points, seed = None, exclude = None}
3493  Module: \code{utilities.polygon}
3494
3495  Populates the interior of the specified polygon with the specified number of points,
3496  selected by means of a uniform distribution function.
3497  \end{funcdesc}
3498
3499  \begin{funcdesc}{point\_in\_polygon}{polygon, delta=1e-8}
3500  Module: \code{utilities.polygon}
3501
3502  Returns a point inside the specified polygon and close to the edge. The distance between
3503  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
3504  second argument \code{delta}, which is taken as $10^{-8}$ by default.
3505  \end{funcdesc}
3506
3507  \begin{funcdesc}{inside\_polygon}{points, polygon, closed = True, verbose = False}
3508  Module: \code{utilities.polygon}
3509
3510  Used to test whether the members of a list of points
3511  are inside the specified polygon. Returns a Numeric
3512  array comprising the indices of the points in the list that lie inside the polygon.
3513  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
3514  Points on the edges of the polygon are regarded as inside if
3515  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
3516  \end{funcdesc}
3517
3518  \begin{funcdesc}{outside\_polygon}{points, polygon, closed = True, verbose = False}
3519  Module: \code{utilities.polygon}
3520
3521  Exactly like \code{inside\_polygon}, but with the words `inside' and `outside' interchanged.
3522  \end{funcdesc}
3523
3524  \begin{funcdesc}{is\_inside\_polygon}{point, polygon, closed=True, verbose=False}
3525  Module: \code{utilities.polygon}
3526
3527  Returns \code{True} if \code{point} is inside \code{polygon} or
3528  \code{False} otherwise. Points on the edges of the polygon are regarded as inside if
3529  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
3530  \end{funcdesc}
3531
3532  \begin{funcdesc}{is\_outside\_polygon}{point, polygon, closed=True, verbose=False}
3533  Module: \code{utilities.polygon}
3534
3535  Exactly like \code{is\_outside\_polygon}, but with the words `inside' and `outside' interchanged.
3536  \end{funcdesc}
3537
3538  \begin{funcdesc}{point\_on\_line}{x, y, x0, y0, x1, y1}
3539  Module: \code{utilities.polygon}
3540
3541  Returns \code{True} or \code{False}, depending on whether the point with coordinates
3542  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
3543  and \code{x1, y1} (extended if necessary at either end).
3544  \end{funcdesc}
3545
3546  \begin{funcdesc}{separate\_points\_by\_polygon}{points, polygon, closed = True, verbose = False}
3547    \indexedcode{separate\_points\_by\_polygon}
3548  Module: \code{utilities.polygon}
3549
3550  \end{funcdesc}
3551
3552  \begin{funcdesc}{polygon\_area}{polygon}
3553  Module: \code{utilities.polygon}
3554
3555  Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html)
3556  \end{funcdesc}
3557
3558  \begin{funcdesc}{plot\_polygons}{polygons, figname, verbose = False}
3559  Module: \code{utilities.polygon}
3560
3561  Plots each polygon contained in input polygon list, e.g.
3562 \code{polygons = [poly1, poly2, poly3]} where \code{poly1 = [[x11,y11],[x12,y12],[x13,y13]]}
3563 etc.  Each polygon is closed for plotting purposes and subsequent plot saved to \code{figname}.
3564  Returns list containing the minimum and maximum of \code{x} and \code{y},
3565  i.e. \code{[x_{min}, x_{max}, y_{min}, y_{max}]}.
3566  \end{funcdesc}
3567
3568\section{coordinate\_transforms}
3569
3570\section{geospatial\_data}
3571\label{sec:geospatial}
3572
3573This describes a class that represents arbitrary point data in UTM
3574coordinates along with named attribute values.
3575
3576%FIXME (Ole): This gives a LaTeX error
3577%\declaremodule{standard}{geospatial_data}
3578%\refmodindex{geospatial_data}
3579
3580
3581
3582\begin{classdesc}{Geospatial\_data}
3583  {data_points = None,
3584    attributes = None,
3585    geo_reference = None,
3586    default_attribute_name = None,
3587    file_name = None}
3588Module: \code{geospatial\_data}
3589
3590This class is used to store a set of data points and associated
3591attributes, allowing these to be manipulated by methods defined for
3592the class.
3593
3594The data points are specified either by reading them from a NetCDF
3595or CSV file, identified through the parameter \code{file\_name}, or
3596by providing their \code{x}- and \code{y}-coordinates in metres,
3597either as a sequence of 2-tuples of floats or as an $M \times 2$
3598Numeric array of floats, where $M$ is the number of points.
3599Coordinates are interpreted relative to the origin specified by the
3600object \code{geo\_reference}, which contains data indicating the UTM
3601zone, easting and northing. If \code{geo\_reference} is not
3602specified, a default is used.
3603
3604Attributes are specified through the parameter \code{attributes},
3605set either to a list or array of length $M$ or to a dictionary whose
3606keys are the attribute names and whose values are lists or arrays of
3607length $M$. One of the attributes may be specified as the default
3608attribute, by assigning its name to \code{default\_attribute\_name}.
3609If no value is specified, the default attribute is taken to be the
3610first one.
3611\end{classdesc}
3612
3613
3614\begin{methoddesc}{import\_points\_file}{delimiter = None, verbose = False}
3615
3616\end{methoddesc}
3617
3618
3619\begin{methoddesc}{export\_points\_file}{ofile, absolute=False}
3620
3621\end{methoddesc}
3622
3623
3624\begin{methoddesc}{get\_data\_points}{absolute = True, as\_lat\_long =
3625    False}
3626    If \code{as\_lat\_long} is\code{True} the point information
3627    returned will be in Latitudes and Longitudes.
3628
3629\end{methoddesc}
3630
3631
3632\begin{methoddesc}{set\_attributes}{attributes}
3633
3634\end{methoddesc}
3635
3636
3637\begin{methoddesc}{get\_attributes}{attribute_name = None}
3638
3639\end{methoddesc}
3640
3641
3642\begin{methoddesc}{get\_all\_attributes}{}
3643
3644\end{methoddesc}
3645
3646
3647\begin{methoddesc}{set\_default\_attribute\_name}{default_attribute_name}
3648
3649\end{methoddesc}
3650
3651
3652\begin{methoddesc}{set\_geo\_reference}{geo_reference}
3653
3654\end{methoddesc}
3655
3656
3657\begin{methoddesc}{add}{}
3658
3659\end{methoddesc}
3660
3661
3662\begin{methoddesc}{clip}{}
3663Clip geospatial data by a polygon
3664
3665Inputs are \code{polygon} which is either a list of points, an Nx2 array or
3666a Geospatial data object and \code{closed}(optional) which determines
3667whether points on boundary should be regarded as belonging to the polygon
3668(\code{closed=True}) or not (\code{closed=False}).
3669Default is \code{closed=True}.
3670
3671Returns new Geospatial data object representing points
3672inside specified polygon.
3673\end{methoddesc}
3674
3675
3676\begin{methoddesc}{clip_outside}{}
3677Clip geospatial data by a polygon
3678
3679Inputs are \code{polygon} which is either a list of points, an Nx2 array or
3680a Geospatial data object and \code{closed}(optional) which determines
3681whether points on boundary should be regarded as belonging to the polygon
3682(\code{closed=True}) or not (\code{closed=False}).
3683Default is \code{closed=True}.
3684
3685Returns new Geospatial data object representing points
3686\emph{out}side specified polygon.
3687\end{methoddesc}
3688
3689
3690\section{Graphical Mesh Generator GUI}
3691The program \code{graphical\_mesh\_generator.py} in the pmesh module
3692allows the user to set up the mesh of the problem interactively.
3693It can be used to build the outline of a mesh or to visualise a mesh
3694automatically generated.
3695
3696Graphical Mesh Generator will let the user select various modes. The
3697current allowable modes are vertex, segment, hole or region.  The mode
3698describes what sort of object is added or selected in response to
3699mouse clicks.  When changing modes any prior selected objects become
3700deselected.
3701
3702In general the left mouse button will add an object and the right
3703mouse button will select an object.  A selected object can de deleted
3704by pressing the the middle mouse button (scroll bar).
3705
3706\section{alpha\_shape}
3707\emph{Alpha shapes} are used to generate close-fitting boundaries
3708around sets of points. The alpha shape algorithm produces a shape
3709that approximates to the `shape formed by the points'---or the shape
3710that would be seen by viewing the points from a coarse enough
3711resolution. For the simplest types of point sets, the alpha shape
3712reduces to the more precise notion of the convex hull. However, for
3713many sets of points the convex hull does not provide a close fit and
3714the alpha shape usually fits more closely to the original point set,
3715offering a better approximation to the shape being sought.
3716
3717In \anuga, an alpha shape is used to generate a polygonal boundary
3718around a set of points before mesh generation. The algorithm uses a
3719parameter $\alpha$ that can be adjusted to make the resultant shape
3720resemble the shape suggested by intuition more closely. An alpha
3721shape can serve as an initial boundary approximation that the user
3722can adjust as needed.
3723
3724The following paragraphs describe the class used to model an alpha
3725shape and some of the important methods and attributes associated
3726with instances of this class.
3727
3728\begin{classdesc}{Alpha\_Shape}{points, alpha = None}
3729Module: \code{alpha\_shape}
3730
3731To instantiate this class the user supplies the points from which
3732the alpha shape is to be created (in the form of a list of 2-tuples
3733\code{[[x1, y1],[x2, y2]}\ldots\code{]}, assigned to the parameter
3734\code{points}) and, optionally, a value for the parameter
3735\code{alpha}. The alpha shape is then computed and the user can then
3736retrieve details of the boundary through the attributes defined for
3737the class.
3738\end{classdesc}
3739
3740
3741\begin{funcdesc}{alpha\_shape\_via\_files}{point_file, boundary_file, alpha= None}
3742Module: \code{alpha\_shape}
3743
3744This function reads points from the specified point file
3745\code{point\_file}, computes the associated alpha shape (either
3746using the specified value for \code{alpha} or, if no value is
3747specified, automatically setting it to an optimal value) and outputs
3748the boundary to a file named \code{boundary\_file}. This output file
3749lists the coordinates \code{x, y} of each point in the boundary,
3750using one line per point.
3751\end{funcdesc}
3752
3753
3754\begin{methoddesc}{set\_boundary\_type}{self,raw_boundary=True,
3755                          remove_holes=False,
3756                          smooth_indents=False,
3757                          expand_pinch=False,
3758                          boundary_points_fraction=0.2}
3759Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3760
3761This function sets flags that govern the operation of the algorithm
3762that computes the boundary, as follows:
3763
3764\code{raw\_boundary = True} returns raw boundary, i.e. the regular edges of the
3765                alpha shape.\\
3766\code{remove\_holes = True} removes small holes (`small' is defined by
3767\code{boundary\_points\_fraction})\\
3768\code{smooth\_indents = True} removes sharp triangular indents in
3769boundary\\
3770\code{expand\_pinch = True} tests for pinch-off and
3771corrects---preventing a boundary vertex from having more than two edges.
3772\end{methoddesc}
3773
3774
3775\begin{methoddesc}{get\_boundary}{}
3776Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3777
3778Returns a list of tuples representing the boundary of the alpha
3779shape. Each tuple represents a segment in the boundary by providing
3780the indices of its two endpoints.
3781\end{methoddesc}
3782
3783
3784\begin{methoddesc}{write\_boundary}{file_name}
3785Module: \code{alpha\_shape},  Class: \class{Alpha\_Shape}
3786
3787Writes the list of 2-tuples returned by \code{get\_boundary} to the
3788file \code{file\_name}, using one line per tuple.
3789\end{methoddesc}
3790
3791\section{Numerical Tools}
3792
3793The following table describes some useful numerical functions that
3794may be found in the module \module{utilities.numerical\_tools}:
3795
3796\begin{tabular}{|p{8cm} p{8cm}|}  \hline
3797\code{angle(v1, v2=None)} & Angle between two-dimensional vectors
3798\code{v1} and \code{v2}, or between \code{v1} and the $x$-axis if
3799\code{v2} is \code{None}. Value is in range $0$ to $2\pi$. \\
3800
3801\code{normal\_vector(v)} & Normal vector to \code{v}.\\
3802
3803\code{mean(x)} & Mean value of a vector \code{x}.\\
3804
3805\code{cov(x, y=None)} & Covariance of vectors \code{x} and \code{y}.
3806If \code{y} is \code{None}, returns \code{cov(x, x)}.\\
3807
3808\code{err(x, y=0, n=2, relative=True)} & Relative error of
3809$\parallel$\code{x}$-$\code{y}$\parallel$ to
3810$\parallel$\code{y}$\parallel$ (2-norm if \code{n} = 2 or Max norm
3811if \code{n} = \code{None}). If denominator evaluates to zero or if
3812\code{y}
3813is omitted or if \code{relative = False}, absolute error is returned.\\
3814
3815\code{norm(x)} & 2-norm of \code{x}.\\
3816
3817\code{corr(x, y=None)} & Correlation of \code{x} and \code{y}. If
3818\code{y} is \code{None} returns autocorrelation of \code{x}.\\
3819
3820\code{ensure\_numeric(A, typecode = None)} & Returns a Numeric array
3821for any sequence \code{A}. If \code{A} is already a Numeric array it
3822will be returned unaltered. Otherwise, an attempt is made to convert
3823it to a Numeric array. (Needed because \code{array(A)} can
3824cause memory overflow.)\\
3825
3826\code{histogram(a, bins, relative=False)} & Standard histogram. If
3827\code{relative} is \code{True}, values will be normalised against
3828the total and thus represent frequencies rather than counts.\\
3829
3830\code{create\_bins(data, number\_of\_bins = None)} & Safely create
3831bins for use with histogram. If \code{data} contains only one point
3832or is constant, one bin will be created. If \code{number\_of\_bins}
3833is omitted, 10 bins will be created.\\  \hline
3834
3835\end{tabular}
3836
3837
3838\chapter{Modules available in \anuga}
3839
3840
3841\section{\module{abstract\_2d\_finite\_volumes.general\_mesh} }
3842\declaremodule[generalmesh]{}{general\_mesh}
3843\label{mod:generalmesh}
3844
3845\section{\module{abstract\_2d\_finite\_volumes.neighbour\_mesh} }
3846\declaremodule[neighbourmesh]{}{neighbour\_mesh}
3847\label{mod:neighbourmesh}
3848
3849\section{\module{abstract\_2d\_finite\_volumes.domain} --- Generic module for 2D triangular domains for finite-volume computations of conservation laws}
3850\declaremodule{}{domain}
3851\label{mod:domain}
3852
3853
3854\section{\module{abstract\_2d\_finite\_volumes.quantity}}
3855\declaremodule{}{quantity}
3856\label{mod:quantity}
3857
3858\begin{verbatim}
3859Class Quantity - Implements values at each triangular element
3860
3861To create:
3862
3863   Quantity(domain, vertex_values)
3864
3865   domain: Associated domain structure. Required.
3866
3867   vertex_values: N x 3 array of values at each vertex for each element.
3868                  Default None
3869
3870   If vertex_values are None Create array of zeros compatible with domain.
3871   Otherwise check that it is compatible with dimenions of domain.
3872   Otherwise raise an exception
3873
3874\end{verbatim}
3875
3876
3877
3878
3879\section{\module{shallow\_water} --- 2D triangular domains for finite-volume
3880computations of the shallow water wave equation. This module contains a specialisation
3881of class Domain from module domain.py consisting of methods specific to the Shallow Water
3882Wave Equation
3883}
3884\declaremodule[shallowwater]{}{shallow\_water}
3885\label{mod:shallowwater}
3886
3887
3888
3889
3890%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3891
3892\chapter{Frequently Asked Questions}
3893
3894
3895\section{General Questions}
3896
3897\subsubsection{What is \anuga?}
3898It is a software package suitable for simulating 2D water flows in
3899complex geometries.
3900
3901\subsubsection{Why is it called \anuga?}
3902The software was developed in collaboration between the
3903Australian National University (ANU) and Geoscience Australia (GA).
3904
3905\subsubsection{How do I obtain a copy of \anuga?}
3906See \url{https://datamining.anu.edu.au/anuga} for all things ANUGA.
3907
3908%\subsubsection{What developments are expected for \anuga in the future?}
3909%This
3910
3911\subsubsection{Are there any published articles about \anuga that I can reference?}
3912See \url{https://datamining.anu.edu.au/anuga} for links.
3913
3914
3915\section{Modelling Questions}
3916
3917\subsubsection{Which type of problems are \anuga good for?}
3918General 2D waterflows in complex geometries such as
3919dam breaks, flows amoung structurs, coastal inundation etc.
3920
3921\subsubsection{Which type of problems are beyond the scope of \anuga?}
3922See Chapter \ref{ch:limitations}.
3923
3924\subsubsection{Can I start the simulation at an arbitrary time?}
3925Yes, using \code{domain.set\_time()} you can specify an arbitrary
3926starting time. This is for example useful in conjunction with a
3927file\_boundary, which may start hours before anything hits the model
3928boundary. By assigning a later time for the model to start,
3929computational resources aren't wasted.
3930
3931\subsubsection{Can I change values for any quantity during the simulation?}
3932Yes, using \code{domain.set\_quantity()} inside the domain.evolve
3933loop you can change values of any quantity. This is for example
3934useful if you wish to let the system settle for a while before
3935assigning an initial condition. Another example would be changing
3936the values for elevation to model e.g. erosion.
3937
3938\subsubsection{Can I change boundary conditions during the simulation?}
3939Yes - see example on page \pageref{sec:change boundary code} in Section
3940\ref{sec:change boundary}.
3941
3942\subsubsection{How do I access model time during the simulation?}
3943The variable \code{t} in the evolve for loop is the model time.
3944For 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})
3945one would write something like
3946{\small \begin{verbatim}
3947    for t in domain.evolve(yieldstep = 0.2, duration = 40.0):
3948
3949        if Numeric.allclose(t, 15):
3950            print 'Changing boundary to outflow'
3951            domain.set_boundary({'right': Bo})
3952
3953\end{verbatim}}
3954The 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()}.
3955
3956
3957\subsubsection{Why does a file\_function return a list of numbers when evaluated?}
3958Currently, file\_function works by returning values for the conserved
3959quantities \code{stage}, \code{xmomentum} and \code{ymomentum} at a given point in time
3960and space as a triplet. To access e.g.\ \code{stage} one must specify element 0 of the
3961triplet returned by file\_function.
3962
3963\subsubsection{Which diagnostics are available to troubleshoot a simulation?}
3964
3965\subsubsection{How do I use a DEM in my simulation?}
3966You use \code{dem2pts} to convert your DEM to the required .pts format. This .pts file is then called
3967when setting the elevation data to the mesh in \code{domain.set_quantity}
3968
3969\subsubsection{What sort of DEM resolution should I use?}
3970Try and work with the \emph{best} you have available. Onshore DEMs
3971are typically available in 25m, 100m and 250m grids. Note, offshore
3972data is often sparse, or non-existent.
3973
3974\subsubsection{What sort of mesh resolution should I use?}
3975The 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,
3976if 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,
3977you need a fine mesh over regions where the DEM changes rapidly, and other areas of significant interest, such as the coast.
3978If 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.
3979
3980
3981\subsubsection{How do I tag interior polygons?}
3982At the moment create_mesh_from_regions does not allow interior
3983polygons with symbolic tags. If tags are needed, the interior
3984polygons must be created subsequently. For example, given a filename
3985of polygons representing solid walls (in Arc Ungenerate format) can
3986be tagged as such using the code snippet:
3987\begin{verbatim}
3988  # Create mesh outline with tags
3989  mesh = create_mesh_from_regions(bounding_polygon,
3990                                  boundary_tags=boundary_tags)
3991  # Add buildings outlines with tags set to 'wall'. This would typically
3992  # bind to a Reflective boundary
3993  mesh.import_ungenerate_file(buildings_filename, tag='wall')
3994
3995  # Generate and write mesh to file
3996  mesh.generate_mesh(maximum_triangle_area=max_area)
3997  mesh.export_mesh_file(mesh_filename)
3998\end{verbatim}
3999
4000Note that a mesh object is returned from \code{create_mesh_from_regions}
4001when file name is omitted.
4002
4003\subsubsection{How often should I store the output?}
4004This 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
4005to look in detail at the evolution, then you will need to balance your storage requirements and the duration of the simulation.
4006If 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
4007quantities on a mesh with approximately 300000 triangles on a 2 min interval for 5 hours will result in approximately 350Mb SWW file
4008(as for the \file{run\_sydney\_smf.py} example).
4009
4010\subsubsection{How can I set the friction in different areas in the domain?}
4011The model area will typically be estimating the water height and momentum over varying
4012topographies which will have different friction values. One way of assigning
4013different friction values is to create polygons (say \code{poly1, poly2 and poly3}) describing each
4014area and then set the corresponding friction values in the following way
4015
4016\code{domain.set_quantity('friction',Polygon_function([(poly1,f1),(poly2,f2),
4017(poly3,f3))]))}
4018
4019The values of \code{f1,f2} and \code{f3} could be constant or functions
4020as determined by the user.
4021
4022\subsubsection{How can I combine data sets?}
4023
4024A user may have access to a range of different resolution DEMs and raw data points (such
4025as beach profiles, spot heights, single or multi-beam data etc) and will need
4026to combine them to create an overall elevation data set.
4027
4028If there are multiple DEMs, say of 10m and 25m resolution, then the technique is similar to
4029that defined in the Cairns example described earlier, that is
4030
4031{\small \begin{verbatim}
4032convert_dem_from_ascii2netcdf(10m_dem_name, use_cache=True, verbose=True)
4033convert_dem_from_ascii2netcdf(25m_dem_name, use_cache=True, verbose=True)
4034\end{verbatim}}
4035followed by
4036{\small \begin{verbatim}
4037dem2pts(10m_dem_name, use_cache=True, verbose=True)
4038dem2pts(25m_dem_name, use_cache=True, verbose=True)
4039\end{verbatim}}
4040These data sets can now be combined by
4041{\small \begin{verbatim}
4042from anuga.geospatial_data.geospatial_data import *
4043G1 = Geospatial_data(file_name = 10m_dem_name + '.pts')
4044G2 = Geospatial_data(file_name = 25m_dem_name + '.pts')
4045G = G1 + G2
4046G.export_points_file(combined_dem_name + ‘.pts’)
4047\end{verbatim}}
4048this is the basic way of combining data sets, however, the user will need to
4049assess the boundaries of each data set and whether they overlap. For example, consider
4050if the 10m DEM is describing by \code{poly1} and the 25m DEM is described by \code{poly2}
4051with \code{poly1} completely enclosed in \code{poly2} as shown in Figure \ref{fig:polydata}
4052\begin{figure}[hbt]
4053  \centerline{\includegraphics{graphics/polyanddata.jpg}}
4054  \caption{Polygons describing the extent of the 10m and 25m DEM.}
4055  \label{fig:polydata}
4056\end{figure}
4057To combine the data sets, the geospatial addition is updated to
4058{\small \begin{verbatim}
4059G = G1 + G2.clip_outside(Geospatial_data(poly1))
4060\end{verbatim}}
4061For this example, we assume that \code{poly2} is the domain, otherwise an additional dataset
4062would be required for the remainder of the domain.
4063
4064This technique can be expanded to handle point data sets as well. In the case
4065of a bathymetry data set available in text format in an \code{.csv} file, then
4066the geospatial addition is updated to
4067{\small \begin{verbatim}
4068G3 = Geospatial_data(file_name = bathy_data_name + '.csv')
4069G = G1 + G2.clip_outside(Geospatial_data(poly1)) + G3
4070\end{verbatim}}
4071The \code{.csv} file has the data stored as \code{x,y,elevation} with the text \code{elevation}
4072on the first line.
4073
4074The coastline could be included
4075as part of the clipping polygon to separate the offshore and onshore datasets if required.
4076Assume that \code{poly1} crosses the coastline
4077In this case, two new polygons could be created out of \code{poly1} which uses the coastline
4078as the divider. As shown in Figure \ref{fig:polycoast}, \code{poly3} describes the
4079onshore data and \code{poly4} describes the offshore data.
4080\begin{figure}[hbt]
4081  \centerline{\includegraphics{graphics/polyanddata2.jpg}}
4082  \caption{Inclusion of new polygons separating the 10m DEM area into an
4083  onshore (poly3) and offshore (poly4) data set.}
4084  \label{fig:polycoast}
4085\end{figure}
4086Let's include the bathymetry
4087data described above, so to combine the datasets in this case,
4088{\small \begin{verbatim}
4089G = G1.clip(Geospatial_data(poly3)) + G2.clip_outside(Geospatial_data(poly1)) + G3
4090\end{verbatim}}
4091
4092Finally, to fit the elevation data to the mesh, the script is adjusted in this way
4093{\small \begin{verbatim}
4094    domain.set_quantity('elevation',
4095                        filename = combined_dem_name + '.pts',
4096                        use_cache = True,
4097                        verbose = True)
4098\end{verbatim}}
4099\subsection{Boundary Conditions}
4100
4101\subsubsection{How do I create a Dirichlet boundary condition?}
4102
4103A Dirichlet boundary condition sets a constant value for the
4104conserved quantities at the boundaries. A list containing
4105the constant values for stage, xmomentum and ymomentum is constructed
4106and used in the function call, e.g. \code{Dirichlet_boundary([0.2,0.,0.])}
4107
4108\subsubsection{How do I know which boundary tags are available?}
4109The method \code{domain.get\_boundary\_tags()} will return a list of
4110available tags for use with
4111\code{domain.set\_boundary\_condition()}.
4112
4113\subsubsection{What is the difference between file_boundary and field_boundary?}
4114The only difference is field_boundary will allow you to change the level of the stage height when you read in the boundary condition.
4115This is very useful when running different tide heights in the same area as you need only to convert
4116one boundary condition to a SWW file, ideally for tide height of 0 m (saving disk space). Then you can
4117use field_boundary to read this SWW file and change the stage height (tide) on the fly depending on the scenario.
4118
4119
4120
4121
4122\subsection{Analysing Results}
4123
4124\subsubsection{How do I easily plot "tide gauges" timeseries graphs from a SWW file?}
4125
4126There is two ways to do this.
4127
41281) Create csv files from the sww file using \code{anuga.abstract_2d_finite_volumes.util sww2csv_gauges}
4129and then use \code{anuga.abstract_2d_finite_volumes.util csv2timeseries_graphs} to
4130create the plots. This code is newer, has unit tests and might be easier to use. Read doc strings for more information and
4131review section 4.7 of this manual.
4132
4133Or
4134
41352) Use \code{anuga.abstract_2d_finite_volumes.util sww2timeseries} to do the whole thing
4136however this doesn't have a much control on the file name and plots. Plus there is no unit tests yet.
4137
4138Read the doc string for more information.
4139
4140\chapter{Glossary}
4141
4142\begin{tabular}{|lp{10cm}|c|}  \hline
4143%\begin{tabular}{|llll|}  \hline
4144    \emph{Term} & \emph{Definition} & \emph{Page}\\  \hline
4145
4146    \indexedbold{\anuga} & Name of software (joint development between ANU and
4147    GA) & \pageref{def:anuga}\\
4148
4149    \indexedbold{bathymetry} & offshore elevation &\\
4150
4151    \indexedbold{conserved quantity} & conserved (stage, x and y
4152    momentum) & \\
4153
4154%    \indexedbold{domain} & The domain of a function is the set of all input values to the
4155%    function.&\\
4156
4157    \indexedbold{Digital Elevation Model (DEM)} & DEMs are digital files consisting of points of elevations,
4158sampled systematically at equally spaced intervals.& \\
4159
4160    \indexedbold{Dirichlet boundary} & A boundary condition imposed on a differential equation
4161 that specifies the values the solution is to take on the boundary of the
4162 domain. & \pageref{def:dirichlet boundary}\\
4163
4164    \indexedbold{edge} & A triangular cell within the computational mesh can be depicted
4165    as a set of vertices joined by lines (the edges). & \\
4166
4167    \indexedbold{elevation} & refers to bathymetry and topography &\\
4168
4169    \indexedbold{evolution} & integration of the shallow water wave equations
4170    over time &\\
4171
4172    \indexedbold{finite volume method} & The method evaluates the terms in the shallow water
4173    wave equation as fluxes at the surfaces of each finite volume. Because the
4174    flux entering a given volume is identical to that leaving the adjacent volume,
4175    these methods are conservative. Another advantage of the finite volume method is
4176    that it is easily formulated to allow for unstructured meshes. The method is used
4177    in many computational fluid dynamics packages. & \\
4178
4179    \indexedbold{forcing term} & &\\
4180
4181    \indexedbold{flux} & the amount of flow through the volume per unit
4182    time & \\
4183
4184    \indexedbold{grid} & Evenly spaced mesh & \\
4185
4186    \indexedbold{latitude} & The angular distance on a mericlear north and south of the
4187    equator, expressed in degrees and minutes. & \\
4188
4189    \indexedbold{longitude} & The angular distance east or west, between the meridian
4190    of a particular place on Earth and that of the Prime Meridian (located in Greenwich,
4191    England) expressed in degrees or time.& \\
4192
4193    \indexedbold{Manning friction coefficient} & &\\
4194
4195    \indexedbold{mesh} & Triangulation of domain &\\
4196
4197    \indexedbold{mesh file} & A TSH or MSH file & \pageref{def:mesh file}\\
4198
4199    \indexedbold{NetCDF} & &\\
4200
4201    \indexedbold{node} & A point at which edges meet & \\
4202
4203    \indexedbold{northing} & A rectangular (x,y) coordinate measurement of distance
4204    north from a north-south reference line, usually a meridian used as the axis of
4205    origin within a map zone or projection. Northing is a UTM (Universal Transverse
4206    Mercator) coordinate. & \\
4207
4208    \indexedbold{points file} & A PTS or CSV file & \\  \hline
4209
4210    \end{tabular}
4211
4212    \begin{tabular}{|lp{10cm}|c|}  \hline
4213
4214    \indexedbold{polygon} & A sequence of points in the plane. \anuga represents a polygon
4215    either as a list consisting of Python tuples or lists of length 2 or as an $N \times 2$
4216    Numeric array, where $N$ is the number of points.
4217
4218    The unit square, for example, would be represented either as
4219    \code{[ [0,0], [1,0], [1,1], [0,1] ]} or as \code{array( [0,0], [1,0], [1,1],
4220    [0,1] )}.
4221
4222    NOTE: For details refer to the module \module{utilities/polygon.py}. &
4223    \\     \indexedbold{resolution} &  The maximal area of a triangular cell in a
4224    mesh & \\
4225
4226
4227    \indexedbold{reflective boundary} & Models a solid wall. Returns same conserved
4228    quantities as those present in the neighbouring volume but reflected. Specific to the
4229    shallow water equation as it works with the momentum quantities assumed to be the
4230    second and third conserved quantities. & \pageref{def:reflective boundary}\\
4231
4232    \indexedbold{stage} & &\\
4233
4234%    \indexedbold{try this}
4235
4236    \indexedbold{animate} & visualisation tool used with \anuga &
4237    \pageref{sec:animate}\\
4238
4239    \indexedbold{time boundary} & Returns values for the conserved
4240quantities as a function of time. The user must specify
4241the domain to get access to the model time. & \pageref{def:time boundary}\\
4242
4243    \indexedbold{topography} & onshore elevation &\\
4244
4245    \indexedbold{transmissive boundary} & & \pageref{def:transmissive boundary}\\
4246
4247    \indexedbold{vertex} & A point at which edges meet. & \\
4248
4249    \indexedbold{xmomentum} & conserved quantity (note, two-dimensional SWW equations say
4250    only \code{x} and \code{y} and NOT \code{z}) &\\
4251
4252    \indexedbold{ymomentum}  & conserved quantity & \\  \hline
4253
4254    \end{tabular}
4255
4256
4257%The \code{\e appendix} markup need not be repeated for additional
4258%appendices.
4259
4260
4261%
4262%  The ugly "%begin{latexonly}" pseudo-environments are really just to
4263%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
4264%  not really valuable.
4265%
4266%  If you don't want the Module Index, you can remove all of this up
4267%  until the second \input line.
4268%
4269
4270%begin{latexonly}
4271%\renewcommand{\indexname}{Module Index}
4272%end{latexonly}
4273\input{mod\jobname.ind}        % Module Index
4274%
4275%begin{latexonly}
4276%\renewcommand{\indexname}{Index}
4277%end{latexonly}
4278\input{\jobname.ind}            % Index
4279
4280
4281
4282\begin{thebibliography}{99}
4283\bibitem[nielsen2005]{nielsen2005}
4284{\it Hydrodynamic modelling of coastal inundation}.
4285Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman.
4286In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
4287Modelling and Simulation. Modelling and Simulation Society of Australia and
4288New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.\\
4289http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
4290
4291\bibitem[grid250]{grid250}
4292Australian Bathymetry and Topography Grid, June 2005.
4293Webster, M.A. and Petkovic, P.
4294Geoscience Australia Record 2005/12. ISBN: 1-920871-46-2.\\
4295http://www.ga.gov.au/meta/ANZCW0703008022.html
4296
4297\bibitem[ZR1999]{ZR1999}
4298\newblock {Catastrophic Collapse of Water Supply Reservoirs in Urban Areas}.
4299\newblock C.~Zoppou and S.~Roberts.
4300\newblock {\em ASCE J. Hydraulic Engineering}, 125(7):686--695, 1999.
4301
4302\bibitem[Toro1999]{Toro1992}
4303\newblock Riemann problems and the waf method for solving the two-dimensional
4304  shallow water equations.
4305\newblock E.~F. Toro.
4306\newblock {\em Philosophical Transactions of the Royal Society, Series A},
4307  338:43--68, 1992.
4308
4309\bibitem{KurNP2001}
4310\newblock Semidiscrete central-upwind schemes for hyperbolic conservation laws
4311  and hamilton-jacobi equations.
4312\newblock A.~Kurganov, S.~Noelle, and G.~Petrova.
4313\newblock {\em SIAM Journal of Scientific Computing}, 23(3):707--740, 2001.
4314\end{thebibliography}{99}
4315
4316\end{document}
Note: See TracBrowser for help on using the repository browser.