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

Last change on this file since 4674 was 4674, checked in by ole, 17 years ago

Typo

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