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

Last change on this file since 5778 was 5765, checked in by ole, 16 years ago

Updated ANUGA_flyer draftImproved caution regarding Transmissive_boundary following talks with Will Power and Biljana Lukovic

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