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

Last change on this file since 5673 was 5673, checked in by duncan, 16 years ago

update to add_points_and_segments

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