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

Last change on this file since 4908 was 4908, checked in by nick, 16 years ago

updates

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