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

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

Run time - flow through cross section and some refactoring + clean up

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