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

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

Manual update

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