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

Last change on this file since 3756 was 3756, checked in by ole, 18 years ago

Updates to channel examples

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