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

Last change on this file since 6450 was 6450, checked in by jgriffin, 15 years ago

Added description of new ability for Time_boundary to call a default_boundary if the time exceeds the length of the input time series.

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