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

Last change on this file since 7552 was 7552, checked in by ole, 13 years ago

Added a little more explanation of the output from triangular_cross

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 207.8 KB
RevLine 
[6450]1% Complete documentation on the extended LaTeX markup used for Python
[7064]2% documentation is available in ''Documenting Python'', which is part
[6450]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;
[7064]15% 'modanuga_user_manual.ind' not found?
[6450]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}
[7086]31\usepackage[hang,small,bf]{caption}
[6450]32
33\input{definitions}
34
[7064]35%%%%%
36% Set penalties for widows, etc, very high
37%%%%%
38
39\widowpenalty=10000
40\clubpenalty=10000
41\raggedbottom
42
[6450]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
[7064]67  \date{\today} % latex2html does not know about datetime
[6450]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
[7064]72                % will be used.
73
[6450]74%\release{1.0}   % release version; this is used to define the
75%                % \version macro
76
77\makeindex          % tell \index to actually write the .idx file
78\makemodindex       % If this contains a lot of module sections.
79
80\setcounter{tocdepth}{3}
81\setcounter{secnumdepth}{3}
82
[7064]83%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84
[6450]85\begin{document}
86\maketitle
87
88% This makes the contents more accessible from the front page of the HTML.
89\ifhtml
[7064]90  \chapter*{Front Matter\label{front}}
[6450]91\fi
92
93%Subversion keywords:
94%
95%$LastChangedDate: 2009-10-29 05:18:30 +0000 (Thu, 29 Oct 2009) $
96%$LastChangedRevision: 7552 $
97%$LastChangedBy: ole $
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
[7064]118\anuga also incorporates a mesh generator that
[6450]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
[7064]127efficiency in C routines working directly with Python numpy structures.
[6450]128
129\end{abstract}
130
131\tableofcontents
132
[7064]133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[6450]134
135\chapter{Introduction}
136
137
138\section{Purpose}
139
140The purpose of this user manual is to introduce the new user to the
[7064]141inundation software system, describe what it can do and give step-by-step
[6450]142instructions for setting up and running hydrodynamic simulations.
[7064]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.
[6450]148
[7068]149This manual describes \anuga version 1.0. To check for later versions of this manual
150go to \url{https://datamining.anu.edu.au/anuga}.
[6450]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
[7071]160The latest installation instructions may be found at:
[7086]161\url{http://datamining.anu.edu.au/\~{}ole/anuga/user_manual/anuga_installation_guide.pdf}.
[7064]162
[6450]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,
[7064]174as well as enough programming experience to adapt the code to different
[6450]175requirements.
176
[7064]177\pagebreak
[6450]178
[7064]179%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[6450]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
[7138]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.}.
[6450]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
[7064]208drive the system such as rainfall, abstraction of water, wind stress or atmospheric pressure
[6450]209gradients. Gravity and frictional resistance from the different
210terrains in the model are represented by predefined forcing terms.
[7064]211See section \ref{sec:forcing terms} for details on forcing terms available in \anuga.
[6450]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
[7135]218(e.g.\ the Manning friction coefficient) for each simulation.
[6450]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
[7064]224efficiency in C routines working directly with Python numeric
[6450]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
[7064]234number of limitations that any potential user needs to be aware of. They are:
[6450]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
[7135]239  waves or 3D turbulence (e.g.\ vorticity).
240  %\item The surface is assumed to be open, e.g.\ \anuga cannot model
[6450]241  %flow under ceilings or in pipes
242  \item All spatial coordinates are assumed to be UTM (meters). As such,
[7064]243  \anuga is unsuitable for modelling flows in areas larger than one UTM zone
[6450]244  (6 degrees wide).
[7064]245  \item Fluid is assumed to be inviscid -- i.e.\ no kinematic viscosity included.
[6450]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.
[7064]250%\item Mesh resolutions near coastlines with steep gradients need to be...
[6450]251  \item Frictional resistance is implemented using Manning's formula, but
[7064]252  \anuga has not yet been fully validated in regard to bottom roughness.
[7134]253%\item \anuga contains no tsunami-genic functionality relating to earthquakes.
[6450]254\end{itemize}
255
[7064]256%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[6450]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
[7064]266
[6450]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
[7064]284system such as rainfall, abstraction of water, wind stress or atmospheric pressure gradients.
[6450]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
[7064]299Water depth can be obtained through the equation:
[6450]300
[7064]301\begin{verbatim}
302depth = stage - elevation
303\end{verbatim}
[6450]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
[7064]311         operation of the model, specifying, for instance,
312         where to store the model output.
[6450]313   \item Inputs various quantities describing physical measurements, such
[7064]314         as the elevation, to be specified at each mesh point (vertex).
[6450]315   \item Sets up the boundary conditions.
316   \item Carries out the evolution of the model through a series of time
[7064]317         steps and outputs the results, providing a results file that can
318         be viewed.
[6450]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
[7064]325'commentary' that describes each step of the program and explains it
[6450]326significance.
327
[7086]328\label{ref:runup_py_code}
[6450]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
[7064]336\begin{verbatim}
337points, vertices, boundary = rectangular_cross(10, 10)
338\end{verbatim}
[6450]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
[7064]34410$ rectangular mesh, triangulated in a regular way. The assignment:
[6450]345
[7064]346\begin{verbatim}
347points, vertices, boundary = rectangular_cross(m, n)
348\end{verbatim}
[6450]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
[7064]356         the boundary and associates each with one of the symbolic tags \code{'left'}, \code{'right'},
[7552]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.
[6450]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
[7064]373\begin{verbatim}
374domain = Domain(points, vertices, boundary)
375\end{verbatim}
[6450]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
[7064]382\begin{verbatim}
383domain.set_name('runup')
384domain.set_datadir('.')
385\end{verbatim}
[6450]386
[7086]387In addition, the following statement could be used to state that
[6450]388quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are
[7342]389to be stored at every timestep and \code{elevation} only once at
390the beginning of the simulation:
[6450]391
[7064]392\begin{verbatim}
[7342]393domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
[7064]394\end{verbatim}
[6450]395
[7342]396However, this is not necessary, as the above is the default behaviour.
[7086]397
[6450]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
[7064]404of ways -- using constants, functions, numeric arrays, expressions
[6450]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
[7064]416The elevation, or height of the bed, is set using a function
[6450]417defined through the statements below, which is specific to this
418example and specifies a particularly simple initial configuration
419for demonstration purposes:
420
[7064]421\begin{verbatim}
[7086]422def topography(x, y):
[7064]423    return -x/2
424\end{verbatim}
[6450]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
[7086]431Once the function \function{topography} is specified, the quantity
[6450]432\code{elevation} is assigned through the simple statement:
433
[7064]434\begin{verbatim}
[7086]435domain.set_quantity('elevation', topography)
[7064]436\end{verbatim}
[6450]437
438NOTE: If using function to set \code{elevation} it must be vector
[7086]439compatible. For example, using square root will not work.
[6450]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
[7064]445quantities -- namely, assign them to a constant numerical value:
[6450]446
[7064]447\begin{verbatim}
448domain.set_quantity('friction', 0.1)
449\end{verbatim}
[6450]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
[7064]457elevation and the depth at any time by the equation:
[6450]458
[7064]459\begin{verbatim}
460stage = elevation + depth
461\end{verbatim}
[6450]462
463For this example, we simply assign a constant value to \code{stage},
[7064]464using the statement:
[6450]465
[7064]466\begin{verbatim}
467domain.set_quantity('stage', -0.4)
468\end{verbatim}
[6450]469
470which specifies that the surface level is set to a height of $-0.4$,
[7135]471i.e.\ 0.4 units (metres) below the zero level.
[6450]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
[7064]475us to illustrate another way to use \method{set\_quantity} -- namely,
[6450]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
[7064]483\begin{verbatim}
484h = 0.05    # Constant depth
485domain.set_quantity('stage', expression='elevation + %f' % h)
486\end{verbatim}
[6450]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}
[7064]493greatly expands its power. See Section \ref{sec:initial conditions} for more
[6450]494details.
495
496\subsection{Boundary Conditions}\index{boundary conditions}
497
498The boundary conditions are specified as follows:
499
[7064]500\begin{verbatim}
501Br = Reflective_boundary(domain)
502Bt = Transmissive_boundary(domain)
503Bd = Dirichlet_boundary([0.2, 0.0, 0.0])
[7086]504Bw = Time_boundary(domain=domain,
505                   f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0])
[7064]506\end{verbatim}
[6450]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}
[7086]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).
[7064]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.
[6450]521    \item \textbf{Transmissive boundary}\label{def:transmissive boundary} 
[7064]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
[7134]527          eventually causing \anuga to crash. If this occurs,
[7135]528          consider using e.g.\ a Dirichlet boundary condition with a stage value
[7064]529          less than the elevation at the boundary.
[6450]530    \item \textbf{Dirichlet boundary}\label{def:dirichlet boundary} Specifies
[7064]531          constant values for stage, $x$-momentum and $y$-momentum at the boundary.
[6450]532    \item \textbf{Time boundary}\label{def:time boundary} Like a Dirichlet
[7064]533          boundary but with behaviour varying with time.
[6450]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
[7064]540the function \code{rectangular}. The example given in
[6450]541Section \ref{sec:realdataexample} illustrates another way of
542assigning the values, by means of the function
[7064]543\code{create_mesh_from_regions}.
[6450]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}
[7064]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.
[6450]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
[7064]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.)
[7134]578All segments in bounding polygon must be tagged. If a tag is not supplied, the default tag name 'exterior' will be assigned by \anuga.
[6450]579
580Using the boundary objects described above, we assign a boundary
[7064]581condition to each part of the boundary by means of a statement like:
[6450]582
[7064]583\begin{verbatim}
584domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
585\end{verbatim}
[6450]586
[7064]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:
[6450]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
[7064]601The command \code{set_boundary} stipulates that, in the current example, the right
[6450]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
[7064]610\begin{verbatim}
611Bw = Time_boundary(domain=domain, f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
612\end{verbatim}
[6450]613
[7064]614\subsection{Evolution}\index{evolution}
[6450]615
[7064]616The final statement:
[6450]617
[7064]618\begin{verbatim}
[7086]619for t in domain.evolve(yieldstep=0.1, duration=10.0):
[7064]620    print domain.timestepping_statistics()
621\end{verbatim}
[6450]622
[7064]623causes the configuration of the domain to 'evolve', over a series of
[6450]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
[7313]633\anuga viewer \code{anuga\_viewer} to generate a visual
[7064]634display (see Section \ref{sec:animate}). See Section \ref{sec:file formats}
[6450]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
[7064]663on the previously dry bed.
[6450]664
[7086]665\code{animate} is described in more detail in Section \ref{sec:animate}.
[6450]666
[7064]667\begin{figure}[htp]
[6450]668  \centerline{\includegraphics[width=75mm, height=75mm]
669    {graphics/bedslopestart.jpg}}
[7134]670  \caption{Runup example viewed with the \anuga viewer}
[6450]671  \label{fig:runupstart}
672\end{figure}
673
[7064]674\begin{figure}[htp]
[6450]675  \centerline{
[7064]676    \includegraphics[width=75mm, height=75mm]{graphics/bedslopeduring.jpg}
[6450]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
[7064]692more complex topographies. These examples build on the
[6450]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}
[7064]697
[6450]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
[7064]703         operation of the model -- specifying, for instance, where to store the
704         model output.
[6450]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
[7064]708         steps and output the results, providing a results file that can be
709         viewed.
[6450]710\end{enumerate}
711
712\subsection{The Code}
713
714Here is the code for the first version of the channel flow \file{channel1.py}:
[7086]715
[6450]716\verbatiminput{demos/channel1.py}
[7086]717
[6450]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
[7064]723In this example we use a similar simple structured triangular mesh as in \file{runup.py}
[6450]724for simplicity, but this time we will use a symmetric one and also
[7064]725change the physical extent of the domain. The assignment:
[6450]726
[7064]727\begin{verbatim}
728points, vertices, boundary = rectangular_cross(m, n, len1=length, len2=width)
729\end{verbatim}
730
[7134]731returns an \code{mxn} mesh similar to the one used in the previous example, except that now the
[6450]732extent in the x and y directions are given by the value of \code{length} and \code{width}
733respectively.
734
[7064]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:
[6450]738
[7064]739\begin{verbatim}
740length = 10.0
741width = 5.0
[6450]742dx = dy = 1           # Resolution: Length of subdivisions on both axes
743
744points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
745                                               len1=length, len2=width)
[7064]746\end{verbatim}
[6450]747
[7064]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
[7086]751The rest of this script is similar to the previous example on page \pageref{ref:runup_py_code}.
[7064]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}
[6450]755%  domain.set_quantity('stage', expression='elevation')
[7064]756%\end{verbatim}
[6450]757
[7064]758
[6450]759\section{Model Output}
760
761The following figure is a screenshot from the \anuga visualisation
762tool \code{animate} of output from this example.
[7064]763
764\begin{figure}[htp]
[6450]765  \centerline{\includegraphics[height=75mm]
766    {graphics/channel1.png}}%
[7134]767  \caption{Simple channel example viewed with the \anuga viewer.}
[6450]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}:
[7064]775
[6450]776\verbatiminput{demos/channel2.py}
[7064]777
[6450]778This example differs from the first version in that a constant outflow boundary condition has
[7064]779been defined:
780
781\begin{verbatim}
782Bo = Dirichlet_boundary([-5, 0, 0])    # Outflow
783\end{verbatim}
784
[6450]785and that it is applied to the right hand side boundary when the water level there exceeds 0m.
[7064]786
787\begin{verbatim}
[7086]788for t in domain.evolve(yieldstep=0.2, finaltime=40.0):
[6450]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})
[7064]794\end{verbatim}
795
[6450]796\label{sec:change boundary code}
[7086]797The \code{if} statement in the timestepping loop (\code{evolve}) gets the quantity
[7064]798\code{stage} and obtains the interpolated value at the point (10m,
[6450]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
[7064]801replaced by the outflow boundary using the method:
802
803\begin{verbatim}
804domain.set_boundary({'right': Bo})
805\end{verbatim}
806
[6450]807This type of dynamically varying boundary could for example be
[7064]808used to model the breakdown of a sluice door when water exceeds a certain level.
[6450]809
810\subsection{Output}
811
[7064]812The text output from this example looks like this:
813
814\begin{verbatim}
[6450]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...
[7064]823\end{verbatim}
[6450]824
825\subsection{Flow through more complex topograhies}
826
827Here is the code for the third version of the channel flow \file{channel3.py}:
[7064]828
[6450]829\verbatiminput{demos/channel3.py}
830
831This example differs from the first two versions in that the topography
832contains obstacles.
833
[7064]834This is accomplished here by defining the function \code{topography} as follows:
835
836\begin{verbatim}
[6450]837def topography(x,y):
[7064]838    """Complex topography defined by a function of vectors x and y."""
[6450]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
[7064]857\end{verbatim}
[6450]858
[7064]859In addition, changing the resolution to \code{dx = dy = 0.1} creates a finer mesh resolving the new features better.
[6450]860
[7064]861A screenshot of this model at time 15s is:
862\begin{figure}[htp]
[6450]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
[7064]870\section{An Example with Real Data}
[6450]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
[7064]896         operation of the model -- specifying, for instance, where to store the
897         model output.
[6450]898
899   \item Input various quantities describing physical measurements, such
[7064]900         as the elevation, to be specified at each mesh point (vertex).
[6450]901
902   \item Set up the boundary conditions.
903
904   \item Carry out the evolution of the model through a series of time
[7064]905         steps and output the results, providing a results file that can be
906         visualised.
[6450]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
[7086]927The mesh-generator creates the mesh within a single
[6450]928polygon whose vertices are at geographical locations specified by
[7064]929the user. The user specifies the \emph{resolution} -- that is, the
[7086]930maximal area of a triangle used for triangulation -- and a triangular
[6450]931mesh is created inside the polygon using a mesh generation engine.
[7134]932On any given platform, the same mesh will be returned each time the
933script is run.
[6450]934
[7064]935Boundary tags are not restricted to \code{'left'}, \code{'bottom'},
936\code{'right'} and \code{'top'}, as in the case of
[6450]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
[7064]945possible to specify one or more 'holes' -- that is, areas bounded by
[6450]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
[7086]955\function{create\_domain\_from\_regions} which creates a Domain object as
956well as a mesh file.  Its arguments include the
[6450]957bounding polygon and its resolution, a list of boundary tags, and a
[7086]958list of pairs \code{[polygon, resolution]} specifying the interior
[6450]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
[7064]968\pageref{sec:file formats}) for more on file formats.
[6450]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
[7064]986\clearpage
[6450]987
[7064]988\begin{figure}[htp]
989  \centerline{\includegraphics[scale=0.5]{graphics/cairns3.jpg}}
990  \caption{Landscape of the Cairns scenario.}
991  \label{fig:cairns3d}
[6450]992\end{figure}
[7064]993
[6450]994The following statements are used to read in the specific polygons
995from \code{project.cairns} and assign a defined resolution to
996each polygon.
997
[7064]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}
[6450]1009
1010Figure \ref{fig:cairnspolys}
1011illustrates the polygons used for the Cairns scenario.
1012
[7064]1013\clearpage
[6450]1014
[7064]1015\begin{figure}[htp]
[6450]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
[7064]1022The statement:
[6450]1023
[7064]1024\begin{verbatim}
[6450]1025remainder_res = 10000000
[7086]1026domain = 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)
[7064]1036\end{verbatim}
1037
[6450]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
[7064]1042polygon -- \code{'top'}, \code{'ocean\_east'}, \code{'bottom'}, and
1043\code{'onshore'} -- and whose values identify the indices of the
[6450]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.
[7086]1049Use the underlying mesh interface for such cases
1050(see Chapter \ref{sec:mesh interface}).
[6450]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
[7086]1060Since we used \code{create_domain_from_regions} to create the mesh file, we do not need to
[7134]1061create the domain explicitly, as the above function does both mesh and domain creation.
[6450]1062
1063The following statements specify a basename and data directory, and
[7134]1064sets a minimum storable height, which helps with visualisation and post-processing
1065if one wants to remove water less than 1cm deep (for instance).
[6450]1066
[7064]1067\begin{verbatim}
[7134]1068domain.set_name('cairns_' + project.scenario) # Name of SWW file
1069domain.set_datadir('.')                       # Store SWW output here
[7086]1070domain.set_minimum_storable_height(0.01)      # Store only depth > 1cm
[7064]1071\end{verbatim}
[6450]1072
[7064]1073\subsection{Initial Conditions}
[6450]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
[7134]1079auxiliary points file.
[6450]1080
1081\subsubsection{Stage}
1082
[7135]1083The stage is initially set to 0.0 (i.e.\ Mean Sea Level) by the following statements:
[6450]1084
[7086]1085\begin{verbatim}
1086tide = 0.0
1087domain.set_quantity('stage', tide)
1088\end{verbatim}
[6450]1089
[7134]1090It could also take the value of the highest astronomical tide.
1091
[7086]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
[6450]1107\subsubsection{Friction}
1108
1109We assign the friction exactly as we did for \file{runup.py}:
1110
[7064]1111\begin{verbatim}
1112domain.set_quantity('friction', 0.0)
1113\end{verbatim}
[6450]1114
1115\subsubsection{Elevation}
1116
[7134]1117The elevation is specified by reading data from a file with a name derived from
1118\code{project.demname} with the \code{.pts} extension:
[6450]1119
[7064]1120\begin{verbatim}
1121domain.set_quantity('elevation',
[7086]1122                    filename=project.demname + '.pts',
[7064]1123                    use_cache=True,
[7086]1124                    verbose=True,
1125                    alpha=0.1)
[7064]1126\end{verbatim}
[6450]1127
[7134]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
[6450]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
[7086]1144boundary for each of the tagged segments defined by \code{create_domain_from_regions}:
[6450]1145
[7064]1146\begin{verbatim}
[7086]1147Bd = Dirichlet_boundary([tide,0,0]) # Mean water level
1148Bs = 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 = 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})
[7064]1165\end{verbatim}
[6450]1166
[7086]1167Note that we use different boundary conditions depending on the \code{scenario}
1168defined in \file{project.py}.
1169
[7153]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
[6450]1173\subsection{Evolution}
1174
[7064]1175With the basics established, the running of the 'evolve' step is
[7086]1176very similar to the corresponding step in \file{runup.py}, except we have different \code{evolve}
1177loops for the two scenarios.
[6450]1178
[7086]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
[7153]1182every ten seconds:
[7086]1183
[7064]1184\begin{verbatim}
[7086]1185if project.scenario == 'slide':
[7134]1186    # Initial run without any event
[7086]1187    for t in domain.evolve(yieldstep=10, finaltime=60):
1188        print domain.timestepping_statistics()
1189        print domain.boundary_statistics(tags='ocean_east')
[6450]1190
[7134]1191    # Add slide to water surface
[7064]1192    if allclose(t, 60):
[7134]1193        domain.add_quantity('stage', tsunami_source)
[6450]1194
[7134]1195    # Continue propagating wave
[7064]1196    for t in domain.evolve(yieldstep=10, finaltime=5000,
[7086]1197                           skip_initial_step=True):
1198        print domain.timestepping_statistics()
1199        print domain.boundary_statistics(tags='ocean_east')
[6450]1200
[7086]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')
[7064]1212\end{verbatim}
1213
[6450]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
[7153]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
[6450]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.
[7064]1228As described earlier, the user may run \code{animate} to view a three-dimensional representation
[6450]1229of the simulation.
1230
1231The user may also be interested in a maximum inundation map. This simply shows the
[7064]1232maximum water depth over the domain and is achieved with the function \code{sww2dem}
1233described in Section \ref{sec:basicfileconversions}).
[6450]1234\file{ExportResults.py} demonstrates how this function can be used:
1235
1236\verbatiminput{demos/cairns/ExportResults.py}
1237
[7064]1238The script generates a maximum water depth ASCII grid at a defined
[6450]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
[7064]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}
[6450]1262\end{figure}
1263
[7064]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}
[6450]1272\end{figure}
1273
[7064]1274\clearpage
1275
[6450]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
[7064]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}
[6450]1285\end{figure}
1286
1287These locations
1288must be stored in either a .csv or .txt file. The corresponding .csv file for
[7064]1289the gauges shown in Figure \ref{fig:cairnsgauges} is \file{gauges.csv}:
[6450]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
[7064]1296\file{GetTimeseries.py} which generates the csv files for each point location. The CSV files
[6450]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
[7064]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}
[6450]1320\end{figure}
1321
[7064]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}
[6450]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}
[7064]1341%    \item Specifying the Quantities: Section \ref{sec:quantities}
[6450]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
[7064]1349each feature is, where to find it and how it can be used -- they do
[6450]1350not give full specifications; for these the reader
1351may consult the code. The code for every function or class contains
[7064]1352a documentation string, or 'docstring', that specifies the precise
[6450]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
[7064]1366containing the function, relative to \file{inundation}, would be:
[6450]1367
[7064]1368\begin{verbatim}
1369pmesh/mesh_interface.py
1370\end{verbatim}
1371
[6450]1372\label{sec:mesh interface}
[7064]1373while in Windows syntax it would be:
[6450]1374
[7064]1375\begin{verbatim}
1376pmesh\mesh_interface.py
1377\end{verbatim}
[6450]1378
1379Rather than using either of these forms, in this chapter we specify
[7064]1380the location simply as \code{pmesh.mesh_interface}, in keeping with
[6450]1381the usage in the Python statement for importing the function,
1382namely:
1383
[7064]1384\begin{verbatim}
1385from pmesh.mesh_interface import create_mesh_from_regions
1386\end{verbatim}
1387
[6450]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
[7064]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\\
[6450]1401\end{tabular}
1402
[7064]1403
[6450]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
[7064]1413\begin{figure}[htp]
[6450]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
[7064]1421\clearpage
[6450]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
[7064]1426points, and may be displayed schematically as in Table \ref{tab:points}.
[6450]1427
[7064]1428\begin{table}[htp]
[6450]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}
[7064]1445  \caption{Point coordinates for mesh in Figure \protect \ref{fig:simplemesh}}
[6450]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
[7064]1459\begin{table}[htp]
[6450]1460  \begin{center}
[7064]1461    \begin{tabular}{|c|ccc|}
1462      \hline
1463      index & \multicolumn{3}{c|}{\code{points}}\\
1464      \hline
[6450]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\\
[7064]1474      9 & 5 & 10 & 9\\
1475      \hline
[6450]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
[7064]1486% \refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}
1487\label{sec:meshgeneration}
[6450]1488
[7134]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}
[6450]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
[7134]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
[6450]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
[7134]1611\begin{classdesc}{Mesh}{userSegments=None,
1612                        userVertices=None,
1613                        holes=None,
1614                        regions=None,
1615                        geo_reference=None}
[6450]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
[7134]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.
[6450]1630\end{classdesc}
1631
1632\subsubsection{Key Methods of Class Mesh}
1633
[7134]1634\begin{methoddesc}{\emph{<mesh>}.add_hole}{x, y, geo_reference=None}
1635Module: \module{pmesh.mesh}
[6450]1636
[7134]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.
[6450]1642\end{methoddesc}
1643
[7134]1644\begin{methoddesc}{\emph{<mesh>}.add_hole_from_polygon}{polygon,
1645                                                        segment_tags=None,
1646                                                        geo_reference=None}
1647Module: \module{pmesh.mesh}
[6450]1648
[7064]1649This method is used to add a 'hole' within a region -- that is, to
[6450]1650define a interior region where the triangular mesh will not be
[7064]1651generated -- to a \class{Mesh} instance. The region boundary is described by
[6450]1652the polygon passed in.  Additionally, the user specifies a list of
1653boundary tags, one for each edge of the bounding polygon.
[7134]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.
[6450]1660\end{methoddesc}
1661
[7134]1662\begin{methoddesc}{\emph{<mesh>}.add_points_and_segments}{points,
1663                                                          segments=None,
1664                                                          segment_tags=None}
1665Module: \module{pmesh.mesh}
[6450]1666
[7134]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
[7135]1672of the line by its point index, e.g.\ use \code{segments = [[0,1],[1,2]]} to make a
[7134]1673polyline between points 0, 1 and 2.
1674
1675\code{segment_tags} may be used to optionally define a tag for each segment.
[6450]1676\end{methoddesc}
1677
[7134]1678\begin{methoddesc}{\emph{<mesh>}.add_region}{x,y, geo_reference=None, tag=None}
1679Module: \module{pmesh.mesh}
[6450]1680
[7134]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.
[6450]1691\end{methoddesc}
1692
[7134]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}
[6450]1699
[7134]1700This method adds a region to a
[6450]1701\class{Mesh} instance.  Regions are commonly used to describe an area
[7134]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.
[6450]1712\end{methoddesc}
1713
[7134]1714\begin{methoddesc}{\emph{<mesh>}.add_vertices}{point_data}
1715Module: \module{pmesh.mesh}
[6450]1716
[7134]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.
[6450]1721\end{methoddesc}
1722
[7134]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}
[6450]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.
[7134]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
[7135]1738boundary, i.e.\ the regular edges of the alpha shape.
[7134]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
[7135]1747(i.e.\ a boundary vertex with more than two edges).
[6450]1748\end{methoddesc}
1749
[7134]1750\begin{methoddesc}{\emph{<mesh>}.export_mesh_file}{ofile}
1751Module: \module{pmesh.mesh}
[6450]1752
[7134]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
[6450]1757\code{.tsh} for the file to be ASCII format.
1758\end{methoddesc}
1759
[7134]1760\begin{methoddesc}{\emph{<mesh>}.generate_mesh}{maximum_triangle_area="",
1761                                                minimum_triangle_angle=28.0,
1762                                                verbose=True}
1763Module: \module{pmesh.mesh}
[6450]1764
[7134]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.
[6450]1772\end{methoddesc}
1773
[7134]1774\begin{methoddesc}{\emph{<mesh>}.import_ungenerate_file}{ofile,
1775                                                         tag=None,
1776                                                         region_tag=None}
1777Module: \module{pmesh.mesh}
[6450]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
[7134]1781vertices and segments.
1782
1783\code{ofile} is the name of the polygon file.
1784
[6450]1785\code{tag} is the tag given to all the polygon's segments.
[7134]1786If \code{tag} is not supplied then the segment will not effect the water
1787flow, it will only effect the mesh generation.
1788
[6450]1789\code{region_tag} is the tag given to all the polygon's segments.  If
[7134]1790it is a string the tag will be assigned to all regions.  If it
[6450]1791is a list the first value in the list will be applied to the first
[7134]1792polygon etc.
[6450]1793
1794This function can be used to import building footprints.
1795\end{methoddesc}
1796
[7064]1797
[6450]1798\section{Initialising the Domain}\index{Initialising the Domain}
1799\label{sec:initialising the domain}
1800
[7134]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}
[6450]1818Module: \refmodule{abstract_2d_finite_volumes.domain}
1819
[7134]1820This class is used to create an instance of a structure used to
[6450]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
[7134]1829\begin{methoddesc}{\emph{<domain>}.set_name}{name}
1830Module: \refmodule{abstract_2d_finite_volumes.domain},
1831page \pageref{mod:domain}
[6450]1832
[7134]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'}.
[6450]1835\end{methoddesc}
1836
[7134]1837\begin{methoddesc}{\emph{<domain>}.get_name}{}
1838Module: \module{abstract_2d_finite_volumes.domain}
[6450]1839
[7134]1840Returns the name assigned to the domain by \code{set_name()}. If no name has been
1841assigned, returns \code{'domain'}.
[6450]1842\end{methoddesc}
1843
[7134]1844\begin{methoddesc}{\emph{<domain>}.set_datadir}{path}
1845Module: \module{abstract_2d_finite_volumes.domain}
[6450]1846
[7134]1847\code{path} specifies the path to the directory used to store SWW files.
[6450]1848
[7134]1849Before this method is used to set the SWW directory path, the assumed directory
1850path is \code{default_datadir} specified in \code{config.py}.
[6450]1851
[7134]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()}.
[6450]1857
[7134]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}
[7064]1863import os
1864domain.set_datadir{'project' + os.sep + 'data'}
[7134]1865\end{verbatim}
[6450]1866\end{methoddesc}
1867
[7134]1868\begin{methoddesc}{\emph{<domain>}.get_datadir}{}
1869Module: \module{abstract_2d_finite_volumes.domain}
[6450]1870
[7134]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}.
[6450]1875\end{methoddesc}
1876
[7134]1877\begin{methoddesc}{\emph{<domain>}.set_minimum_allowed_height}{minimum_allowed_height}
1878Module: \module{shallow_water.shallow_water_domain}
[6450]1879
[7134]1880Set the minimum depth (in metres) that will be recognised in
1881the numerical scheme (including limiters and flux computations)
[6450]1882
[7134]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,
[7135]1886e.g.\ for large scale simulations, the computation time can be
[7134]1887significantly reduced.
[6450]1888\end{methoddesc}
1889
[7134]1890\begin{methoddesc}{\emph{<domain>}.set_minimum_storable_height}{minimum_storable_height}
1891Module: \module{shallow_water.shallow_water_domain}
[6450]1892
[7134]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.
[6450]1898\end{methoddesc}
1899
[7134]1900\begin{methoddesc}{\emph{<domain>}.set_maximum_allowed_speed}{maximum_allowed_speed}
1901Module: \module{shallow_water.shallow_water_domain}
[6450]1902
[7134]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.
[6450]1909\end{methoddesc}
1910
[7134]1911\begin{methoddesc}{\emph{<domain>}.set_time}{time=0.0}
1912Module: \module{abstract_2d_finite_volumes.domain}
[6450]1913
[7134]1914\code{time} sets the initial time, in seconds, for the simulation. The
1915default is 0.0.
[6450]1916\end{methoddesc}
1917
[7134]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.)
[6450]1924\end{methoddesc}
1925
[7134]1926\begin{methoddesc}{\emph{<domain>}.set_store_vertices_uniquely}{flag, reduction=None}
1927Module: \module{shallow_water.shallow_water_domain}
1928
[6450]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
[7312]1933\code{flag} may be \code{True} (meaning allow surface to be discontinuous) or \code{False} (meaning smooth vertex values).
[7134]1934
[7312]1935\code{reduction} defines the smoothing operation if \code{flag} is \code{False}.  If not
[7134]1936supplied, \code{reduction} is assumed to be \code{mean}.
1937
1938Triangles stored in the SWW file can be discontinuous reflecting
[6450]1939the internal representation of the finite-volume scheme
[7134]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
[6450]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
[7134]1948Storing one way or the other is controlled in \anuga through the method
1949\code{<domain>.store_vertices_uniquely()}. Options are:
[6450]1950\begin{itemize}
[7134]1951  \item \code{<domain>.store_vertices_uniquely(True)}: Allow discontinuities in the SWW file
1952  \item \code{<domain>.store_vertices_uniquely(False)}: (Default).
[6450]1953  Average values
[7134]1954  to ensure continuity in SWW file. The latter also makes for smaller
1955  SWW files.
[6450]1956\end{itemize}
1957
[7135]1958Note that when model data in the SWW file are averaged (i.e.\ not stored uniquely),
[7134]1959then there will most likely be a small discrepancy between values extracted from the SWW
[6634]1960file and the same data stored in the model domain. This must be borne in mind when comparing
[7134]1961data from the SWW files with that of the model internally.
[6450]1962\end{methoddesc}
1963
[7342]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
[6450]1984% Structural methods
[7134]1985\begin{methoddesc}{\emph{<domain>}.get_nodes}{absolute=False}
1986Module: \module{abstract_2d_finite_volumes.domain}
[6450]1987
[7134]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
[7135]1990they were provided in the constructor i.e.\ without any duplication.
[6450]1991
[7134]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.
[6450]1995\end{methoddesc}
1996
[7134]1997\begin{methoddesc}{\emph{<domain>}.get_vertex_coordinates}{absolute=False}
1998Module: \module{abstract_2d_finite_volumes.domain}
[6450]1999
[7134]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.
[6450]2004
[7134]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.
[6450]2008\end{methoddesc}
2009
[7134]2010\begin{methoddesc}{\emph{<domain>}.get_centroid_coordinates}{absolute=False}
2011Module: \module{abstract_2d_finite_volumes.domain}
[6450]2012
[7134]2013Return centroid coordinates for all triangles as an \code{Mx2} array.
[6450]2014   
[7134]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.
[6450]2018\end{methoddesc}
2019
[7134]2020\begin{methoddesc}{\emph{<domain>}.get_triangles}{indices=None}
2021Module: \module{abstract_2d_finite_volumes.domain}
[6450]2022
[7134]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()}.
[6450]2027
[7134]2028\code{indices}, if specified, is the set of triangle \code{id}s of interest.
[6450]2029\end{methoddesc}
2030
[7134]2031\begin{methoddesc}{\emph{<domain>}.get_disconnected_triangles}{}
2032Module: \module{abstract_2d_finite_volumes.domain}
[6450]2033
[7134]2034Get the domain mesh based on nodes obtained from \code{get_vertex_coordinates()}.
[6450]2035
[7134]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.
[6450]2040
[7134]2041This provides a mesh where no triangles share nodes
2042(hence the name disconnected triangles) and different
2043nodes may have the same coordinates.
[6450]2044
[7134]2045This version of the mesh is useful for storing meshes with
[7135]2046discontinuities at each node and is e.g.\ used for storing
[7134]2047data in SWW files.
[6450]2048
[7134]2049The triangles created will have the format:
2050
2051\begin{verbatim}
[7064]2052[[0,1,2],
2053 [3,4,5],
2054 [6,7,8],
2055 ...
2056 [3*M-3 3*M-2 3*M-1]]
[7134]2057\end{verbatim}
[6450]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
[7134]2071that \code{set_quantity()} does not necessarily have to be used in the initial
[6450]2072condition setting; it can be used at any time throughout the simulation.
2073
[7134]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()})
[6450]2089
[7134]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:
[6450]2095
[7134]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
[7135]2100  \item a function (e.g.\ see the samples introduced in Chapter 2)
[7134]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}
[6450]2110
[7134]2111Exactly one of the arguments \code{numeric}, \code{quantity}, \code{function},
2112\code{geospatial_data} and \code{filename} must be present.
[6450]2113
[7134]2114\code{set_quantity()} will look at the type of the \code{numeric} and
2115determine what action to take.
[6450]2116
[7134]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.
[6450]2120
[7134]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}
[6450]2133
[7134]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.
[6450]2136\end{methoddesc}
2137
[7134]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()})
[6450]2153
[7134]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)}.
[6450]2157
[7134]2158A typical use of this function is to add structures to an existing elevation model:
[6450]2159
[7134]2160\begin{verbatim} 
[7064]2161# Create digital elevation model from points file
2162domain.set_quantity('elevation', filename='elevation_file.pts, verbose=True)
[6450]2163
[7064]2164# Add buildings from file
2165building_polygons, building_heights = csv2building_polygons(building_file)
[6450]2166
[7064]2167B = []
2168for key in building_polygons:
2169    poly = building_polygons[key]
2170    elev = building_heights[key]
2171    B.append((poly, elev))
[6450]2172
[7064]2173domain.add_quantity('elevation', Polygon_function(B, default=0.0))
[7134]2174\end{verbatim}
[6450]2175\end{methoddesc}
2176
[7134]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})
[6450]2180
[7134]2181This function is used to assign values to individual quantities given
2182a regional tag.   It is similar to \code{set_quantity()}.
[6450]2183
[7134]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'.
[6450]2190
[7134]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}
[6450]2196
2197\begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None,
[7134]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}
[6450]2205
[7134]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.
[6450]2209
[7134]2210\code{length} is the length of the slide or slump.
[6450]2211
[7134]2212\code{depth} is the water depth to the centre of the sediment mass.
[6450]2213
[7134]2214\code{slope} is the bathymetric slope.
[6450]2215
[7134]2216Other slump or slide parameters can be included if they are known.
2217\end{funcdesc}
[6450]2218
[7134]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}
[6450]2229
[7134]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}.
[6450]2233
[7282]2234\code{filename} is the name of the input file. 
2235This would be either an SWW, TMS or STS file.
[6450]2236
[7134]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.
[7283]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.
[6450]2244
[7134]2245\code{interpolation_points} is a list of absolute coordinates or a
2246geospatial object for points at which values are sought.
[6450]2247
[7134]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.
[6450]2251
[7134]2252The model time stored within the file function can be accessed using
2253the method \code{\emph{<callable_object>}.get_time()}
[6450]2254
[7134]2255The underlying algorithm used is as follows:\\
[7135]2256Given a time series (i.e.\ a series of values associated with
[7134]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.
[6450]2265
[7134]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.
[6450]2275
[7134]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}
[6450]2282
[7134]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}
[6450]2293
[7135]2294Given a time series (i.e.\ a series of values associated with
[7134]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.
[6450]2301
[7134]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.
[6450]2311
[7134]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}
[6450]2318
[7134]2319
[6450]2320\section{Boundary Conditions}\index{boundary conditions}
2321\label{sec:boundary conditions}
2322
2323\anuga provides a large number of predefined boundary conditions,
[7134]2324represented by objects such as \code{Reflective_boundary(domain)} and
2325\code{Dirichlet_boundary([0.2, 0.0, 0.0])}, described in the examples
[7064]2326in Chapter 2. Alternatively, you may prefer to ''roll your own'',
[6450]2327following the method explained in Section \ref{sec:roll your own}.
2328
[7134]2329These boundary objects may be used with the function \code{set_boundary} described below
[6450]2330to assign boundary conditions according to the tags used to label boundary segments.
2331
[7134]2332\begin{methoddesc}{\emph{<domain>}.set_boundary}{boundary_map}
2333Module: \module{abstract_2d_finite_volumes.domain}
[6450]2334
[7134]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.
[6450]2338
[7134]2339\code{boundary_map} is a dictionary of boundary objects keyed by symbolic tags.
[6450]2340\end{methoddesc}
2341
[7134]2342\begin{methoddesc} {\emph{<domain>}.get_boundary_tags}{}
2343Module: \module{abstract\_2d\_finite\_volumes.domain}
[6450]2344
[7134]2345Returns a list of the available boundary tags.
[6450]2346\end{methoddesc}
2347
2348\subsection{Predefined boundary conditions}
2349
[7134]2350\begin{classdesc}{Reflective_boundary}{domain=None}
2351Module: \module{shallow_water}
[6450]2352
[7134]2353Reflective boundary returns same conserved quantities as those present in
2354the neighbouring volume but reflected.
[6450]2355
[7134]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.
[6450]2358\end{classdesc}
2359
[7134]2360\begin{classdesc}{Transmissive_boundary}{domain=None}
[7064]2361  \label{pg: transmissive boundary}
[7134]2362Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
[6450]2363
[7134]2364A transmissive boundary returns the same conserved quantities as
2365those present in the neighbouring volume.
[6450]2366
[7134]2367The underlying domain must be specified when the boundary is instantiated.
[6450]2368\end{classdesc}
2369
[7134]2370\begin{classdesc}{Dirichlet_boundary}{conserved_quantities=None}
2371Module: \module{abstract_2d_finite_volumes.generic_boundary_conditions}
[6450]2372
[7134]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.
[6450]2378\end{classdesc}
2379
[7134]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}
[6450]2385
[7134]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.
[6450]2389
[7134]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.
[6450]2394\end{classdesc}
2395
[7134]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}
[6450]2405
[7134]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.
[6450]2410
[7134]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.
[6450]2415\end{classdesc}
2416
[7134]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}
[6450]2427
[7134]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}.
[6450]2431
[7134]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.
[6450]2437
[7134]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.
[6450]2441
[7134]2442See \code{File_boundary} for further details.
[6450]2443\end{classdesc}
2444
[7134]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}
[6450]2449
[7134]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.
[6450]2453
[7134]2454This type of boundary is useful when stage is known at the boundary as a
2455function of time, but momenta (or speeds) aren't.
[6450]2456
[7134]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.
[6450]2459
[7134]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.
[6450]2464\end{classdesc}
2465
[7134]2466\begin{classdesc}{Transmissive_stage_zero_momentum_boundary}{domain=None}
2467Module: \module{shallow_water}
2468\label{pg: transmissive stage zero momentum boundary}
[6450]2469
[7134]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
[6450]2473
[7134]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.
[6450]2478
[7134]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.
[6450]2481
[7134]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}).
[6450]2485\end{classdesc}
2486
[7134]2487\begin{classdesc}{Dirichlet_discharge_boundary}{domain=None, stage0=None, wh0=None}
2488Module: \module{shallow_water.shallow_water_domain}
[6450]2489
[7134]2490\code{stage0} sets the stage.
2491
2492\code{wh0} sets momentum in the inward normal direction.
[6450]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
[7134]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
[6450]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
[7134]2507particular values or flags to be used by the \code{evaluate()} method.
[6450]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.
[7064]2516Gravity and friction are always calculated using the elevation and friction quantities,
2517but the user may additionally add forcing terms to the list
[7134]2518\code{domain.forcing_terms} and have them affect the model.
[6450]2519
[7134]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}
[6450]2530
[7134]2531This is a general class to modify any quantity according to a given rate of change.
[7135]2532Other specific forcing terms are based on this class but it can be used by itself as well (e.g.\ to modify momentum).
[6450]2533
[7134]2534\code{domain} is the domain being evolved.
[6450]2535
[7134]2536\code{quantity_name} is the name of the quantity that will be affected by this forcing term.
[6450]2537
[7134]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
[7135]2541before a forcing term is applied (i.e.\ simulation has started).
2542The default value is 0.0 -- i.e.\ no forcing.
[7134]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}
[7064]2554P = [[x0, y0], [x1, y0], [x1, y1], [x0, y1]]    # Square polygon
[6450]2555
[7064]2556xmom = General_forcing(domain, 'xmomentum', polygon=P)
2557ymom = General_forcing(domain, 'ymomentum', polygon=P)
[6450]2558
[7064]2559xmom.rate = f
2560ymom.rate = g
2561
2562domain.forcing_terms.append(xmom)
2563domain.forcing_terms.append(ymom)
[7134]2564\end{verbatim}
[7064]2565
[7134]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}
[6450]2570
[7134]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}
[6450]2578
[7134]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.
[6450]2581
[7134]2582Inflow is based on the \code{General_forcing} class so the functionality is similar.
[6450]2583
[7134]2584\code{domain} is the domain being evolved.
[7064]2585
[7134]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.
[6450]2590
[7134]2591\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
[6450]2592
[7134]2593\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
[7064]2594
[7134]2595Example:
2596
2597\begin{verbatim}
2598hydrograph = Inflow(center=(320, 300), radius=10,
2599                    rate=file_function('QPMF_Rot_Sub13.tms'))
2600
[7064]2601domain.forcing_terms.append(hydrograph)
[7134]2602\end{verbatim}
[7064]2603
[7134]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}
[6450]2606
[7134]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}
[6450]2615
[7134]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.
[6450]2618
[7134]2619Rainfall is based on the \code{General_forcing} class so the functionality is similar.
[6450]2620
[7134]2621\code{domain} is the domain being evolved.
[7064]2622
[7134]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).
[6450]2630
[7134]2631\code{center} and \code{ radius} optionally restrict forcing to a circle with given center and radius.
[6450]2632
[7134]2633\code{polygon} optionally restricts forcing to an area enclosed by the given polygon.
2634
2635Example:
2636
2637\begin{verbatim}
[7064]2638catchmentrainfall = Rainfall(rate=file_function('Q100_2hr_Rain.tms'))
2639domain.forcing_terms.append(catchmentrainfall)
[7134]2640\end{verbatim}
[6450]2641
[7134]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}
[6450]2644
[7134]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}
[6450]2669
[7134]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.
[6450]2672
[7134]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.
[6450]2676
[7189]2677The \code{Culvert_flow} class takes as input:
[7134]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}
[6450]2702
[7134]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:
[6450]2705
[7134]2706\begin{verbatim}
[7064]2707from anuga.culvert_flows.culvert_class import Culvert_flow
2708from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
[6450]2709
[7064]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],
[7134]2715                        width=1.20,
2716                        height=0.75,
[7064]2717                        culvert_routine=boyd_generalised_culvert_model,
2718                        number_of_barrels=1,
2719                        verbose=True)
[6450]2720
[7064]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)
[6450]2731
[7064]2732domain.forcing_terms.append(culvert1)
2733domain.forcing_terms.append(culvert2)
[7134]2734\end{verbatim}
2735\end{classdesc}
[6450]2736
2737
2738\section{Evolution}\index{evolution}
2739\label{sec:evolution}
2740
[7134]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}
[6450]2746
[7134]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.
[6450]2755
[7134]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.
[6450]2760
[7134]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}
[7064]2778for t in domain.evolve(yieldstep, finaltime):
2779    <Do something with domain and t>
[7134]2780\end{verbatim}
[7064]2781\end{methoddesc}
[6450]2782
2783\subsection{Diagnostics}
2784\label{sec:diagnostics}
2785
[7134]2786\begin{methoddesc}{\emph{<domain>}.statistics}{}
2787Module: \module{abstract\_2d\_finite\_volumes.domain}
[6450]2788
[7134]2789Outputs statistics about the mesh within the \code{Domain}.
2790\end{methoddesc}
[6450]2791
[7134]2792\begin{methoddesc}{\emph{<domain>}.timestepping_statistics}{track_speeds=False, triangle_id=None}
2793Module: \module{abstract_2d_finite_volumes.domain}
[6450]2794
[7134]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)}
[6450]2797
[7134]2798Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and
2799the number of first-order steps, respectively.
[6450]2800
[7134]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.
[6450]2806
[7134]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}
[6450]2810
[7134]2811\begin{methoddesc}{\emph{<domain>}.boundary_statistics}{quantities=None,
2812                                                      tags=None}
2813Module: \module{abstract_2d_finite_volumes.domain}
[6450]2814
[7134]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}
[7064]2826Boundary values at time 0.5000:
[6450]2827    top:
2828        stage in [ -0.25821218,  -0.02499998]
2829    bottom:
2830        stage in [ -0.27098821,  -0.02499974]
[7134]2831\end{verbatim}
2832\end{methoddesc}
[6450]2833
[7134]2834\begin{methoddesc}{Q = \emph{<domain>}.get_quantity}{name,
2835                                               location='vertices',
2836                                               indices=None}
2837Module: \module{abstract_2d_finite_volumes.domain}
[6450]2838
[7134]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}.
[6450]2841
[7134]2842\code{name} is the name of the quantity to retrieve.
[6450]2843
[7134]2844\code{location} is
[6450]2845
[7134]2846\code{indices} is
2847\end{methoddesc}
[6450]2848
[7134]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}
[6450]2853
[7134]2854Selects quantities and derived quantities for which extrema attained at internal timesteps
2855will be collected.
[6450]2856
[7134]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}
[6450]2863
[7134]2864\code{polygon} can be used to monitor only triangles inside the polygon. Otherwise
2865all triangles will be included.
[6450]2866
[7134]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}
[7064]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)
[7134]2895\end{verbatim}
[6450]2896
[7134]2897The quantities (and derived quantities) listed here must be selected at model
2898initialisation time using the method \code{domain.set_quantities_to_be_monitored()}.
[6450]2899
[7134]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.
[6450]2904
[7134]2905These values are also stored in the SWW file for post-processing.
2906\end{methoddesc}
[6450]2907
[7134]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}
[6450]2915
[7134]2916Extract values for quantity as a numeric array.
[6450]2917
[7134]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.
[7064]2922   
[7134]2923\code{location} specifies where values are to be stored.
2924Permissible options are 'vertices', 'edges', 'centroids' or 'unique vertices'.
[6450]2925
[7134]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}.
[7064]2928
[7134]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}.
[6450]2935
[7134]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.
[6450]2941
[7134]2942The values will be stored in elements following their internal ordering.
2943\end{methoddesc}
[6450]2944
[7134]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}
[6450]2959
[7134]2960Assign values to a quantity object.
[6450]2961
[7134]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.
[6450]2967
[7134]2968\code{location} is ??
[6450]2969
[7134]2970\code{indices} is ??
[6450]2971
[7134]2972The method \code{set_values()} is always called by \code{set_quantity()}
2973behind the scenes.
2974\end{methoddesc}
[6450]2975
[7134]2976\begin{methoddesc}{\emph{<quantity>}.get_integral}{}
2977Module: \module{abstract_2d_finite_volumes.quantity}
[6450]2978
[7134]2979Return the computed integral over the entire domain for the quantity.
2980\end{methoddesc}
[6450]2981
[7134]2982\begin{methoddesc}{\emph{<quantity>}.get_maximum_value}{indices=None}
2983Module: \module{abstract_2d_finite_volumes.quantity}
[6450]2984
[7134]2985Return the maximum value of a quantity (on centroids).
[6450]2986
[7134]2987\code{indices} is the optional set of element \code{id}s that
2988the operation applies to.
[6450]2989
[7134]2990We do not seek the maximum at vertices as each vertex can
2991have multiple values -- one for each triangle sharing it.
2992\end{methoddesc}
[6450]2993
[7134]2994\begin{methoddesc}{\emph{<quantity>}.get_maximum_location}{indices=None}
2995Module: \module{abstract_2d_finite_volumes.quantity}
[6450]2996
[7134]2997Return the location of the maximum value of a quantity (on centroids).
[6450]2998
[7134]2999\code{indices} is the optional set of element \code{id}s that
3000the operation applies to.
[6450]3001
[7134]3002We do not seek the maximum at vertices as each vertex can
3003have multiple values -- one for each triangle sharing it.
[6450]3004
[7134]3005If there are multiple cells with the same maximum value, the
3006first cell encountered in the triangle array is returned.
3007\end{methoddesc}
[6450]3008
[7134]3009\begin{methoddesc}{\emph{<domain>}.get_wet_elements}{indices=None}
3010Module: \module{shallow_water.shallow_water_domain}
[6450]3011
[7134]3012Returns the indices for elements where h $>$ minimum_allowed_height
[6450]3013
[7134]3014\code{indices} is the optional set of element \code{id}s that
3015the operation applies to.
3016\end{methoddesc}
[6450]3017
[7134]3018\begin{methoddesc}{\emph{<domain>}.get_maximum_inundation_elevation}{indices=None}
3019Module: \module{shallow_water.shallow_water_domain}
[6450]3020
[7134]3021Return highest elevation where h $>$ 0.
[6450]3022
[7134]3023\code{indices} is the optional set of element \code{id}s that
3024the operation applies to.
[6450]3025
[7134]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
[6450]3047\section{Queries of SWW model output files}
[7134]3048After a model has been run, it is often useful to extract various information from the SWW
[6450]3049output file (see Section \ref{sec:sww format}). This is typically more convenient than using the
[7064]3050diagnostics described in Section \ref{sec:diagnostics} which rely on the model running -- something
[7134]3051that can be very time consuming. The SWW files are easy and quick to read and offer information
[6450]3052about the model results such as runup heights, time histories of selected quantities,
3053flow through cross sections and much more.
3054
[7134]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}
[6450]3060
[7134]3061Return the highest elevation where depth is positive ($h > 0$).
[6450]3062
[7134]3063\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
[6450]3064
[7134]3065\code{polygon} restricts the query to the points that lie within the polygon.
[6450]3066
[7134]3067\code {time_interval} restricts the query to within the time interval.
[6450]3068
[7134]3069Example to find maximum runup elevation:
[6450]3070
[7134]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".
[6450]3077\end{funcdesc}
3078
[7134]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}
[6450]3084
[7134]3085Return location (x,y) of the highest elevation where depth is positive ($h > 0$).
[6450]3086
[7134]3087\code{filename} is the path to a NetCDF SWW file containing \anuga model output.
[6450]3088
[7134]3089\code{polygon} restricts the query to the points that lie within the polygon.
[6450]3090
[7134]3091\code {time_interval} restricts the query to within the time interval.
[6450]3092
[7134]3093Example to find maximum runup location:
[6450]3094
[7134]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".
[6450]3102\end{funcdesc}
3103
[7134]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}
[6450]3120
[7134]3121Read a set of SWW files and plot the time series for the prescribed quantities
3122at defined gauge locations and prescribed time range.
[6450]3123
[7134]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!
[6450]3129\end{funcdesc}
3130
[7134]3131\begin{funcdesc}{(time, Q) = get_flow_through_cross_section}{filename, polyline, verbose=False}
3132Module: \module{shallow_water.data_manager}
[6450]3133
[7134]3134Obtain flow ($m^3/s$) perpendicular to specified cross section.
[6450]3135
[7134]3136\code{filename} is the path to the SWW file.
[6450]3137
[7134]3138\code{polyline} is the representation of the desired cross section -- it may contain
3139multiple sections allowing for complex shapes. Assumes absolute UTM coordinates.
[6450]3140
[7134]3141Returns a tuple \code{time, Q} where:
[6450]3142
[7134]3143\code{time} is all the stored times in the SWW file.
[6450]3144
[7134]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}
[7064]3154cross_section = [[x, 0], [x, width]]
[7134]3155time, Q = get_flow_through_cross_section(filename, cross_section)
3156\end{verbatim}
[6450]3157\end{funcdesc}
3158
3159
3160\section{Other}
[7134]3161\begin{methoddesc}{quantity = \emph{<domain>}.create_quantity_from_expression}{string}
3162Module: \module{abstract_2d_finite_volumes.domain}
[6450]3163
[7134]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}
[7064]3173Depth = domain.create_quantity_from_expression('stage-elevation')
[6450]3174
[7064]3175exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
3176Absolute_momentum = domain.create_quantity_from_expression(exp)
[7134]3177\end{verbatim}
[6450]3178
[7134]3179%See also \file{Analytical_solution_circular_hydraulic_jump.py} for an example.
3180\end{methoddesc}
[6450]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
[7064]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)}\\
[7134]3198  \code{.sts} & NetCDF format for storing model ouput \code{f(t,x,y)} without mesh information\\
[7064]3199  \code{.tms} & NetCDF format for storing time series \code{f(t)}\\
[7134]3200  \code{.csv/.txt} & ASCII format for storing arbitrary points and associated attributes\\
[7064]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
[6450]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
[7064]3215capitalise the extension and omit the initial full stop -- thus, we
[7134]3216refer to 'SWW files' or 'PRJ files', not 'sww files' or '.prj files'.
[6450]3217
3218\bigskip
3219
3220A typical dataflow can be described as follows:
3221
[7134]3222SOMETHING MISSING HERE!?
3223
[6450]3224\subsection{Manually Created Files}
3225
3226\begin{tabular}{ll}
3227ASC, PRJ & Digital elevation models (gridded)\\
[7135]3228NC & Model outputs for use as boundary conditions (e.g.\ from MOST)
[6450]3229\end{tabular}
3230
3231\subsection{Automatically Created Files}
3232
3233\begin{tabular}{ll}
[7064]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
[6450]3240\end{tabular}
3241
3242\bigskip
3243
3244\subsection{SWW, STS and TMS Formats}
3245\label{sec:sww format}
[7134]3246The SWW, STS and TMS formats are all NetCDF formats and are of key importance for \anuga.
[6450]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}
[7064]3254  \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
[7134]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
[6450]3257    % Refer here to the example to be provided in describing the simple example
[7134]3258  \item \code{time}: a numeric array containing times for model evaluation
[6450]3259\end{itemize}
3260
[7134]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}.
[6450]3264
3265Alternatively, there are tools, such as \code{ncdump}, that allow
[7134]3266you to convert a NetCDF file into a readable format such as the
[6450]3267Class Definition Language (CDL). The following is an excerpt from a
3268CDL representation of the output file \file{runup.sww} generated
[7134]3269from running the simple example \file{runup.py} of Chapter \ref{ch:getstarted}:
[6450]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
[7134]3278An STS file is used for storing a set of points and associated times.
[6450]3279It contains, in addition to dimension information, the following
3280variables:
3281\begin{itemize}
[7064]3282  \item \code{x} and \code{y}: coordinates of the points, represented as numeric arrays
[7134]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})
[7064]3285  \item \code{elevation}: a numeric array storing bed-elevations
[6450]3286    % Refer here to the example to be provided in describing the simple example
[7064]3287  \item \code{time}: a numeric array containing times for model evaluation
[6450]3288\end{itemize}
3289
[7064]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.
[6450]3295
[7064]3296A TMS file is used to store time series data that is independent of position.
[6450]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
[7134]3304from the function \function{create_mesh_from_regions()} (see
[7064]3305Section \ref{sec:meshgeneration}) and be used to initialise a domain.
[6450]3306
[7064]3307A mesh file can define the outline of the mesh -- the vertices and
[6450]3308line segments that enclose the region in which the mesh is
[7064]3309created -- and the triangular mesh itself, which is specified by
[6450]3310listing the triangles and their vertices, and the segments, which
3311are those sides of the triangles that are associated with boundary
3312conditions.
3313
[7064]3314In addition, a mesh file may contain 'holes' and/or 'regions'. A
[6450]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
[7135]3321offset to be applied to $x$ and $y$ values -- e.g.\ to the vertices.
[6450]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:\\
[7064]3329 \\
3330first line:   \code{[column names]}\\
3331other lines:  \code{[x value], [y value], [attributes]}\\
[6450]3332
[7064]3333for example:
[6450]3334
[7064]3335\begin{verbatim}
[7134]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
[7064]3340\end{verbatim}
[6450]3341
[7064]3342The delimiter is a comma. The first two columns are assumed to
[7134]3343be $x$ and $y$ coordinates.
[6450]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
[7064]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
[6450]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
[7134]3383\begin{funcdesc}{sww2dem}{(basename_in,
3384            basename_out=None,
[7064]3385            quantity=None,
3386            timestep=None,
3387            reduction=None,
3388            cellsize=10,
3389            number_of_decimal_places=None,
3390            NODATA_value=-9999,
3391            easting_min=None,
3392            easting_max=None,
3393            northing_min=None,
3394            northing_max=None,
3395            verbose=False,
3396            origin=None,
3397            datum='WGS84',
[7134]3398            format='ers',
3399            block_size=None}
3400Module: \module{shallow_water.data_manager}
[6450]3401
[7134]3402Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or
3403ERS) of a desired grid size \code{cellsize} in metres. The user can select how
3404many decimal places the output data is represented with by using \code{number_of_decimal_places},
3405with the default being 3.
3406
3407The $easting$ and $northing$ values are used if the user wishes to determine the output
3408within a specified rectangular area. The \code{reduction} input refers to a function
[7135]3409to reduce the quantities over all time step of the SWW file, e.g.\ maximum.
[7064]3410\end{funcdesc}
[6450]3411
[7064]3412\begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
[6450]3413            easting_min=None, easting_max=None,
3414            northing_min=None, northing_max=None,
3415            use_cache=False, verbose=False}
3416  Module: \module{shallow\_water.data\_manager}
3417
3418  Takes DEM data (a NetCDF file representation of data from a regular Digital
3419  Elevation Model) and converts it to PTS format.
[7064]3420\end{funcdesc}
[6450]3421
[7064]3422\begin{funcdesc}{urs2sts}{basename_in, basename_out=None,
[6450]3423            weights=None, verbose=False,
3424            origin=None,mean_stage=0.0,
3425            zscale=1.0, ordering_filename=None}
3426  Module: \module{shallow\_water.data\_manager}
3427
[7064]3428  Takes URS data (timeseries data in mux2 format) and converts it to STS format.
3429  The optional filename \code{ordering\_filename} specifies the permutation of indices
3430  of points to select along with their longitudes and latitudes. This permutation will also be
3431  stored in the STS file. If absent, all points are taken and the permutation will be trivial,
[7135]3432  i.e.\ $0, 1, \ldots, N-1$, where $N$ is the total number of points stored. 
[7064]3433\end{funcdesc}
[6450]3434
3435\begin{funcdesc}{csv2building\_polygons}{file\_name, floor\_height=3}
[7064]3436  Module: \module{shallow\_water.data\_manager}
[6450]3437
[7064]3438  Convert CSV files of the form:
[6450]3439
[7064]3440  \begin{verbatim} 
3441easting,northing,id,floors
3442422664.22,870785.46,2,0
3443422672.48,870780.14,2,0
3444422668.17,870772.62,2,0
3445422660.35,870777.17,2,0
3446422664.22,870785.46,2,0
3447422661.30,871215.06,3,1
3448422667.50,871215.70,3,1
3449422668.30,871204.86,3,1
3450422662.21,871204.33,3,1
3451422661.30,871215.06,3,1
3452  \end{verbatim}
[6450]3453
[7064]3454  to a dictionary of polygons with \code{id} as key.
3455  The associated number of \code{floors} are converted to m above MSL and
3456  returned as a separate dictionary also keyed by \code{id}.
[6450]3457   
[7064]3458  Optional parameter \code{floor_height} is the height of each building story.
[6450]3459
[7135]3460  These can e.g.\ be converted to a \code{Polygon_function} for use with \code{add_quantity}
[7064]3461  as shown on page \pageref{add quantity}.
[6450]3462\end{funcdesc}
3463
3464%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3465
3466\chapter{\anuga mathematical background}
3467\label{cd:mathematical background}
3468
[7064]3469
[6450]3470\section{Introduction}
3471
3472This chapter outlines the mathematics underpinning \anuga.
3473
3474
3475\section{Model}
3476\label{sec:model}
3477
3478The shallow water wave equations are a system of differential
3479conservation equations which describe the flow of a thin layer of
3480fluid over terrain. The form of the equations are:
3481\[
3482\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
3483x}+\frac{\partial \GG}{\partial y}=\SSS
3484\]
3485where $\UU=\left[ {{\begin{array}{*{20}c}
3486 h & {uh} & {vh} \\
3487\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
3488$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
3489entering the system are bed elevation $z$ and stage (absolute water
3490level) $w$, where the relation $w = z + h$ holds true at all times.
3491The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
3492by
3493\[
3494\EE=\left[ {{\begin{array}{*{20}c}
3495 {uh} \hfill \\
3496 {u^2h+gh^2/2} \hfill \\
3497 {uvh} \hfill \\
3498\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
3499 {vh} \hfill \\
3500 {vuh} \hfill \\
3501 {v^2h+gh^2/2} \hfill \\
3502\end{array} }} \right]
3503\]
3504and the source term (which includes gravity and friction) is given
3505by
3506\[
3507\SSS=\left[ {{\begin{array}{*{20}c}
3508 0 \hfill \\
3509 -{gh(z_{x} + S_{fx} )} \hfill \\
3510 -{gh(z_{y} + S_{fy} )} \hfill \\
3511\end{array} }} \right]
3512\]
3513where $S_f$ is the bed friction. The friction term is modelled using
3514Manning's resistance law
3515\[
3516S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
3517=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
3518\]
3519in which $\eta$ is the Manning resistance coefficient.
3520The model does not currently include consideration of kinematic viscosity or dispersion.
3521
3522As demonstrated in our papers, \cite{ZR1999,nielsen2005} these
3523equations and their implementation in \anuga provide a reliable
3524model of general flows associated with inundation such as dam breaks
3525and tsunamis.
3526
[7064]3527
[6450]3528\section{Finite Volume Method}
3529\label{sec:fvm}
3530
3531We use a finite-volume method for solving the shallow water wave
3532equations \cite{ZR1999}. The study area is represented by a mesh of
3533triangular cells as in Figure~\ref{fig:mesh} in which the conserved
3534quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
3535in each volume are to be determined. The size of the triangles may
3536be varied within the mesh to allow greater resolution in regions of
3537particular interest.
3538
[7064]3539\begin{figure}[htp] \begin{center}
3540  \includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-five}
3541  \caption{Triangular mesh used in our finite volume method. Conserved
3542           quantities $h$, $uh$ and $vh$ are associated with the centroid of
3543           each triangular cell.}
3544  \label{fig:mesh}
3545\end{center} \end{figure}
[6450]3546
3547The equations constituting the finite-volume method are obtained by
3548integrating the differential conservation equations over each
3549triangular cell of the mesh. Introducing some notation we use $i$ to
3550refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
3551set of indices referring to the cells neighbouring the $i$th cell.
3552Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
3553the length of the edge between the $i$th and $j$th cells.
3554