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

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

Added channel examples to user manual

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