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

Last change on this file since 3925 was 3925, checked in by sexton, 18 years ago

updates

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