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

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

more updates

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