source: trunk/anuga_core/documentation/user_manual/anuga_user_manual.tex @ 7804

Last change on this file since 7804 was 7804, checked in by hudson, 12 years ago

Fixed up failing tests, updated user guide with new API (first few chapters only).

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