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

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

Updated installation guide

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