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

Last change on this file since 6045 was 6045, checked in by ole, 15 years ago

Implemented, tested and documunted new boundary condition:

Transmissive_stage_zero_momentum_boundary

Renamed two others according to style guide but retained backward compatibility

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