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

Last change on this file since 5555 was 5555, checked in by jakeman, 11 years ago

John added information on STS format to ANUGA manual

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