source: documentation/user_manual/anuga_user_manual.tex @ 3087

Last change on this file since 3087 was 3087, checked in by howard, 18 years ago

Added material on geo_spatial data and alpha shapes to Appendix A

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 94.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
16
17\documentclass{manual}
18
19\usepackage{graphicx}
20\usepackage{datetime}
21
22\input{definitions}
23
24\title{\anuga User Manual}
25\author{Geoscience Australia and the Australian National University}
26
27% Please at least include a long-lived email address;
28% the rest is at your discretion.
29\authoraddress{Geoscience Australia \\
30  Email: \email{ole.nielsen@ga.gov.au}
31}
32
33%Draft date
34
35% update before release!
36% Use an explicit date so that reformatting
37% doesn't cause a new date to be used.  Setting
38% the date to \today can be used during draft
39% stages to make it easier to handle versions.
40
41
42\longdate       % Make date format long using datetime.sty
43%\settimeformat{xxivtime} % 24 hour Format
44\settimeformat{oclock} % Verbose
45\date{\today, \ \currenttime}
46%\hyphenation{set\_datadir}
47
48\ifhtml
49\date{\today} % latex2html does not know about datetime
50\fi
51
52
53
54
55
56\release{1.0}   % release version; this is used to define the
57                % \version macro
58
59\makeindex          % tell \index to actually write the .idx file
60\makemodindex       % If this contains a lot of module sections.
61
62\setcounter{tocdepth}{3}
63\setcounter{secnumdepth}{3}
64
65
66\begin{document}
67\maketitle
68
69
70% This makes the contents more accessible from the front page of the HTML.
71\ifhtml
72\chapter*{Front Matter\label{front}}
73\fi
74
75%Subversion keywords:
76%
77%$LastChangedDate: 2006-06-05 07:07:42 +0000 (Mon, 05 Jun 2006) $
78%$LastChangedRevision: 3087 $
79%$LastChangedBy: howard $
80
81\input{copyright}
82
83
84\begin{abstract}
85\label{def:anuga}
86
87\noindent \anuga\index{\anuga} is a hydrodynamic modelling tool that
88allows users to model realistic flow problems in complex geometries.
89Examples include dam breaks or the effects of natural hazards such
90as riverine flooding, storm surges and tsunami.
91
92The user must specify a study area represented by a mesh of triangular
93cells, the topography and bathymetry, frictional resistance, initial
94values for water level (called \emph{stage}\index{stage} within \anuga),
95boundary
96conditions and forces such as windstress or pressure gradients if
97applicable.
98
99\anuga tracks the evolution of water depth and horizontal momentum
100within each cell over time by solving the shallow water wave equation
101governing equation using a finite-volume method.
102
103\anuga cannot model details of breaking waves, flow under ceilings such
104as pipes, turbulence and vortices, vertical convection or viscous
105flows.
106
107\anuga also incorporates a mesh generator, called \code{pmesh}, that
108allows the user to set up the geometry of the problem interactively as
109well as tools for interpolation and surface fitting, and a number of
110auxiliary tools for visualising and interrogating the model output.
111
112Most \anuga components are written in the object-oriented programming
113language Python and most users will interact with \anuga by writing
114small Python programs based on the \anuga library
115functions. Computationally intensive components are written for
116efficiency in C routines working directly with the Numerical Python
117structures.
118
119
120\end{abstract}
121
122\tableofcontents
123
124
125\chapter{Introduction}
126
127
128\section{Purpose}
129
130The purpose of this user manual is to introduce the new user to
131the software, describe what it can do and give step-by-step
132instructions for setting up and running hydrodynamic simulations.
133
134\section{Scope}
135
136This manual covers only what is needed to operate the software after
137installation and configuration. It does not includes instructions
138for installing the software or detailed API documentation, both of
139which will be covered in separate publications and by documentation
140in the source code.
141
142\section{Audience}
143
144Readers are assumed to be familiar with the operating environment
145and have a general understanding of the problem background, as well
146as enough programming experience to adapt the code to different
147requirements, as described in this manual, and to understand the
148basic terminology of object-oriented programming.
149
150\pagebreak
151\chapter{Background}
152
153
154Modelling the effects on the built environment of natural hazards such
155as riverine flooding, storm surges and tsunami is critical for
156understanding their economic and social impact on our urban
157communities.  Geoscience Australia and the Australian National
158University are developing a hydrodynamic inundation modelling tool
159called \anuga to help simulate the impact of these hazards.
160
161The core of \anuga is the fluid dynamics module, called pyvolution,
162which is based on a finite-volume method for solving the shallow water
163wave equation.  The study area is represented by a mesh of triangular
164cells.  By solving the governing equation within each cell, water
165depth and horizontal momentum are tracked over time.
166
167A major capability of pyvolution is that it can model the process of
168wetting and drying as water enters and leaves an area.  This means
169that it is suitable for simulating water flow onto a beach or dry land
170and around structures such as buildings.  Pyvolution is also capable
171of modelling hydraulic jumps due to the ability of the finite-volume
172method to accommodate discontinuities in the solution.
173
174To set up a particular scenario the user specifies the geometry
175(bathymetry and topography), the initial water level, boundary
176conditions such as tide, and any forcing terms that may drive the
177system such as wind stress or atmospheric pressure gradients.
178Frictional resistance from the different terrains in the model is
179represented by predefined forcing terms.
180
181A mesh generator, called pmesh, allows the user to set up the geometry
182of the problem interactively and to identify boundary segments and
183regions using symbolic tags.  These tags may then be used to set the
184actual boundary conditions and attributes for different regions
185(e.g. the Manning friction coefficient) for each simulation.
186
187Most \anuga components are written in the object-oriented programming
188language Python.  Software written in Python can be produced quickly
189and can be readily adapted to changing requirements throughout its
190lifetime.  Computationally intensive components are written for
191efficiency in C routines working directly with the Numerical Python
192structures.  The animation tool developed for \anuga is based on
193OpenSceneGraph, an Open Source Software (OSS) component allowing high
194level interaction with sophisticated graphics primitives.
195
196See \cite{nielsen2005} and \cite{roberts1999} for more background on \anuga.
197%FIXME (Ole): Add bibliography
198
199\chapter{Getting Started}
200\label{ch:getstarted}
201
202This section is designed to assist the reader to get started with
203\anuga by working through simple examples. Two examples are discussed;
204the first is a simple but artificial example that is useful to illustrate
205many of the ideas, and the second is a more realistic example.
206
207\section{A Simple Example}
208
209\subsection{Overview}
210
211What follows is a discussion of the structure and operation of a script which we will call \file{runup.py},
212with just enough detail to allow
213the reader to appreciate what's involved in setting up a scenario
214like the one it depicts.
215
216This example carries out the solution of the shallow-water wave
217equation in the simple case of a configuration comprising a flat
218bed, sloping at a fixed angle in one direction and having a
219constant depth across each line in the perpendicular direction.
220
221The example demonstrates many of the basic ideas involved in
222setting up a more complex scenario. In the general case the user
223specifies the geometry (bathymetry and topography), the initial
224water level, boundary conditions such as tide, and any forcing
225terms that may drive the system such as wind stress or atmospheric
226pressure gradients. Frictional resistance from the different
227terrains in the model is represented by predefined forcing
228terms. The boundary is reflective on three sides and a time dependent wave on one side.
229
230The present example represents a simple scenario and does not
231include any forcing terms, nor is the data taken from a file as it
232would be in many typical cases.
233
234The conserved quantities involved in the
235problem are stage (absolute height of water surface),
236$x$-momentum and $y$-momentum. Other quantities
237involved in the computation are the friction and elevation.
238
239Water depth can be obtained through the equation
240
241\begin{tabular}{rcrcl}
242  \code{depth} &=& \code{stage} &$-$& \code{elevation}
243\end{tabular}
244
245
246%\emph{[More details of the problem background]}
247
248\subsection{Outline of the Program}
249
250In outline, \file{runup.py} performs the following steps:
251
252\begin{enumerate}
253
254   \item Sets up a triangular mesh.
255
256   \item Sets certain parameters governing the mode of
257operation of the model-specifying, for instance, where to store the model output.
258
259   \item Inputs various quantities describing physical measurements, such
260as the elevation, to be specified at each mesh point (vertex).
261
262   \item Sets up the boundary conditions.
263
264   \item Carries out the evolution of the model through a series of time
265steps and outputs the results, providing a results file that can
266be visualised.
267
268\end{enumerate}
269
270\subsection{The Code}
271
272%FIXME: we are using the \code function here.
273%This should be used wherever possible
274For reference we include below the complete code listing for
275\file{runup.py}. Subsequent paragraphs provide a
276`commentary' that describes each step of the program and explains it
277significance.
278
279\verbatiminput{examples/runup.py}
280%\verbatiminput{examples/bedslope.py}
281
282\subsection{Establishing the Mesh}
283
284The first task is to set up the triangular mesh to be used for the
285scenario. This is carried out through the statement:
286
287{\small \begin{verbatim}
288    points, vertices, boundary = rectangular(10, 10)
289\end{verbatim}}
290
291The function \function{rectangular} is imported from a module
292\module{mesh\_factory} defined elsewhere. (\anuga also contains
293several other schemes that can be used for setting up meshes, but we
294shall not discuss these now.) The above assignment sets up a $10
295\times 10$ rectangular mesh, triangulated in a specific way. In
296general, the assignment
297
298{\small \begin{verbatim}
299    points, vertices, boundary = rectangular(m, n)
300\end{verbatim}}
301
302returns:
303
304\begin{itemize}
305
306   \item a list \code{points} of length $N$, where $N = (m + 1)(n + 1)$,
307comprising the coordinates \code{(x, y)} of each of the $N$ mesh
308points,
309
310   \item a list \code{vertices} of length $2mn$ (each entry specifies the three
311vertices of one of the triangles used in the triangulation) , and
312
313   \item a dictionary \code{boundary}, used to tag the triangle edges on
314the boundaries. Each key corresponds to a triangle edge on one of
315the four boundaries and its value is one of \code{`left'},
316\code{`right'}, \code{`top'} and \code{`bottom'}, indicating
317which boundary the edge in question belongs to.
318
319\end{itemize}
320
321Since these three variables encapsulate the information needed to
322specify the grid, it may be helpful to consider how they look in a
323very simple case, not directly related to the example at hand.
324Figure \ref{fig:simplemesh} represents a very simple mesh comprising
325just 11 points and three triangles. (To avoid confusion, we should
326emphasise that this particular mesh is \emph{not} generated by
327\code{rectangular}---it is not even rectangular in nature. )
328
329
330\begin{figure}[h]
331  \begin{center}
332    \includegraphics[width=90mm, height=90mm]{triangularmesh.eps}
333  \end{center}
334
335  \caption{A simple mesh}
336  \label{fig:simplemesh}
337\end{figure}
338
339
340The variables \code{points}, \code{vertices} and \code{boundary}
341represent the data displayed in Figure \ref{fig:simplemesh} as
342follows. The list \code{points} stores the coordinates of the
343points, and may be displayed schematically as in Table \ref{tab:points}.
344
345
346\begin{table}
347  \begin{center}
348    \begin{tabular}[t]{|c|cc|} \hline
349      index & \code{x} & \code{y}\\  \hline
350      0 & 1 & 1\\
351      1 & 4 & 2\\
352      2 & 8 & 1\\
353      3 & 1 & 3\\
354      4 & 5 & 5\\
355      5 & 8 & 6\\
356      6 & 11 & 5\\
357      7 & 3 & 6\\
358      8 & 1 & 8\\
359      9 & 4 & 9\\
360      10 & 10 & 7\\  \hline
361    \end{tabular}
362  \end{center}
363
364  \caption{Point coordinates for mesh in
365    Figure \protect \ref{fig:simplemesh}}
366  \label{tab:points}
367\end{table}
368
369The list \code{vertices} specifies the triangles that make up the
370mesh. It does this by specifying, for each triangle, the indices
371(the numbers shown in the first column above) that correspond to the
372three points at its vertices, taken in an anti-clockwise order
373around the triangle. Thus, in the example shown in Figure
374\ref{fig:simplemesh}, the variable \code{vertices} contains the
375entries shown in Table \ref{tab:vertices}. The starting point is arbitrary so triangle $(0,1,3)$ is considered the same as $(1,3,0)$ and $(3,0,1)$.
376
377
378\begin{table}
379  \begin{center}
380    \begin{tabular}{|c|ccc|} \hline
381      index & \multicolumn{3}{c|}{\code{vertices}}\\ \hline
382      0 & 0 & 1 & 3\\
383      1 & 1 & 2 & 4\\
384      2 & 2 & 5 & 4\\
385      3 & 2 & 6 & 5\\
386      4 & 4 & 5 & 9\\
387      5 & 4 & 9 & 7\\
388      6 & 3 & 4 & 7\\
389      7 & 7 & 9 & 8\\
390      8 & 1 & 4 & 3\\
391      9 & 5 & 10 & 9\\  \hline
392    \end{tabular}
393  \end{center}
394
395  \caption{Vertices for mesh in
396    Figure \protect \ref{fig:simplemesh}}
397  \label{tab:vertices}
398\end{table}
399
400Finally, the variable \code{boundary} identifies the boundary
401triangles and associates a tag with each.
402
403
404
405
406\subsection{Initialising the Domain}
407
408These variables are then used to set up a data structure
409\code{domain}, through the assignment:
410
411{\small \begin{verbatim}
412    domain = Domain(points, vertices, boundary)
413\end{verbatim}}
414
415This uses a Python class \class{Domain}, imported from
416\module{shallow\_water}, which is an extension of a more generic
417class of the same name in the module \refmodule{pyvolution.domain}
418(page \pageref{mod:pyvolution.domain}),
419and inherits
420some methods from the generic class but has others specific to the
421shallow-water scenarios in which it is used. Specific options for
422domain are set at this point. One of them is to set the basename for
423the output file:
424
425{\small \begin{verbatim}
426    domain.set_name('bedslope')
427\end{verbatim}}
428
429
430\subsection{Initial Conditions}
431
432The next task is to specify a number of quantities that we wish to
433set for each mesh point. The class \class{Domain} has a method
434\method{set\_quantity}, used to specify these quantities. It is a
435particularly flexible method that allows the user to set quantities
436in a variety of ways---using constants, functions, numeric arrays or
437expressions involving other quantities, arbitrary data points with
438associated values, all of which can be passed as arguments. All
439quantities can be initialised using \method{set\_quantity}. For a
440conserved quantity (such as \code{stage, xmomentum, ymomentum}) this
441is called an \emph{initial condition}. However, other quantities
442that aren't updated by the equation are also assigned values using
443the same interface. The code in the present example demonstrates a
444number of forms in which we can invoke \method{set\_quantity}.
445
446
447\subsubsection{Elevation}
448
449The elevation, or height of the bed, is set using a function,
450defined through the statements below, which is specific to this
451example and specifies a particularly simple initial configuration
452for demonstration purposes:
453
454{\small \begin{verbatim}
455    def f(x,y):
456        return -x/2
457\end{verbatim}}
458
459This simply associates an elevation with each point \code{(x, y)} of
460the plane.  It specifies that the bed slopes linearly in the
461\code{x} direction, with slope $-\frac{1}{2}$,  and is constant in
462the \code{y} direction.  %[screen shot?]
463
464Once the function \function{f} is specified, the quantity
465\code{elevation} is assigned through the simple statement:
466
467{\small \begin{verbatim}
468    domain.set_quantity('elevation', f)
469\end{verbatim}}
470
471
472\subsubsection{Friction}
473
474The assignment of the friction quantity demonstrates another way we
475can use \method{set\_quantity} to set quantities---namely, assign
476them to a constant numerical value:
477
478{\small \begin{verbatim}
479    domain.set_quantity('friction', 0.1)
480\end{verbatim}}
481
482This just specifies that the Manning friction coefficient is set
483to 0.1 at every mesh point.
484
485\subsubsection{Stage}
486
487The stage (the height of the water surface) is related to the
488elevation and the depth at any time by the equation
489
490
491{\small \begin{verbatim}
492    stage = elevation + depth
493\end{verbatim}}
494
495
496For this example, we simply assign a constant value to \code{stage},
497using the statement
498
499{\small \begin{verbatim}
500    domain.set_quantity('stage', -.4)
501\end{verbatim}}
502
503which specifies that the surface level is set to a height of $-0.4$, i.e. 0.4 units
504below the zero level.
505
506(Although it is not necessary for this example, it may be useful to
507digress here and mention a variant to this requirement, which allows
508us to illustrate another way to use \method{set\_quantity}---namely,
509incorporating an expression involving other quantities. Suppose,
510instead of setting a constant value for the stage, we wished to
511specify a constant value for the \emph{depth}. For such a case we
512need to specify that \code{stage} is everywhere obtained by adding
513that value to the value already specified for \code{elevation}. We
514would do this by means of the statements:
515
516{\small \begin{verbatim}
517    h = 0.05 # Constant depth
518    domain.set_quantity('stage', expression = 'elevation + %f' %h)
519\end{verbatim}}
520
521That is, the value of \code{stage} is set to $\code{h} = 0.05$ plus
522the value of \code{elevation} already defined.
523
524The reader will probably appreciate that this capability to
525incorporate expressions into statements using \method{set\_quantity}
526greatly expands its power.)
527
528\subsection{Boundary Conditions}
529
530The boundary conditions are specified as follows:
531
532{\small \begin{verbatim}
533    Br = Reflective_boundary(domain)
534
535    Bt = Transmissive_boundary(domain)
536
537    Bd = Dirichlet_boundary([0.2,0.,0.])
538
539    Bw = Time_boundary(domain=domain,
540                f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
541\end{verbatim}}
542
543The effect of these statements is to set up a selection of different
544alternative boundary conditions and store them in variables that can be
545assigned as needed. Each boundary condition specifies the
546behaviour at a boundary in terms of the behaviour in neighbouring
547elements. The boundary conditions introduced here may be briefly described as
548follows:
549
550\begin{itemize}
551    \item \textbf{Reflective boundary} Returns same \code{stage} as
552      as present in its neighbour volume but momentum vector
553      reversed 180 degrees (reflected).
554      Specific to the shallow water equation as it works with the
555      momentum quantities assumed to be the second and third conserved
556      quantities. A reflective boundary condition models a solid wall.
557    \item \textbf{Transmissive boundary} Returns same conserved quantities as
558      those present in its neighbour volume. This is one way of modelling
559      outflow from a domain, but it should be used with caution if flow is
560      not steady state as replication of momentum at the boundary
561      may cause occasional spurious effects. If this occurs,
562      consider using e.g. a Dirichlet boundary condition.
563    \item \textbf{[Dirichlet boundary} Specifies a fixed value at the
564      boundary and assigns initial values to the $x$-momentum and $y$-momentum.
565    \item \textbf {Time boundary} Like a Dirichlet boundary but with behaviour
566      varying with time.
567\end{itemize}
568
569Once the boundary objects have been specified through the
570statements above, they can be applied through a statement of the
571form
572
573{\small \begin{verbatim}
574    domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
575\end{verbatim}}
576
577This statement stipulates that, in the current example, the right
578boundary varies with time, as defined by the lambda function, while the other
579boundaries are all reflective.
580
581The reader may wish to experiment by varying the choice of boundary
582types for one or more of the boundaries. In the case of \code{Bd}
583and \code{Bw}, the three arguments in each case represent the
584\code{stage}, $x$-momentum and $y$-momentum, respectively.
585
586{\small \begin{verbatim}
587    Bw = Time_boundary(domain=domain,
588                       f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])
589\end{verbatim}}
590
591
592
593\subsection{Evolution}
594
595The final statement \nopagebreak[3]
596{\small \begin{verbatim}
597    for t in domain.evolve(yieldstep = 0.1, duration = 4.0):
598        print domain.timestepping_statistics()
599\end{verbatim}}
600
601is the key step that causes the configuration of the domain to
602`evolve' in accordance with the model embodied in the code, over a
603series of steps indicated by the values of \code{yieldstep} and
604\code{duration}, which can be altered as required.  The value of
605\code{yieldstep} controls the time interval between successive model
606outputs.  Behind the scenes more time steps are generally taken.
607
608
609\subsection{Output}
610
611%Give details here of the form of the output and explain how it can
612%be used with swollen. Include screen shots.//
613
614The output is a NetCDF file with the extension \code{.sww}. It
615contains stage and momentum information and can be used with the
616\code{swollen} visualisation package to generate a visual display.
617See Section \ref{sec:file formats} (page \pageref{sec:file formats})
618for more on NetCDF and other file formats.
619
620
621\section{How to Run the Code}
622
623The code can be run in various ways:
624
625\begin{itemize}
626  \item{from a Windows command line} as in \code{python runup.py}
627  \item{within the Python IDLE environment}
628  \item{within emacs}
629  \item{from a Linux command line} as in \code{python runup.py}
630\end{itemize}
631
632
633\section{Exploring the Model Output}
634
635Figure \ref{fig:runupstart} shows the domain with water surface as
636specified by the initial condition, $t=0$. Figure
637\ref{fig:bedslope2} shows later snapshots for $t=2.3$ and $t=4$
638where the system has been evolved and the wave is encroaching on the
639previously dry bed.  All figures are screenshots from an interactive
640animation tool called Swollen which is part of \anuga. Swollen is
641described in more detail is Section \ref{sec:swollen}.
642
643
644
645\begin{figure}[hbt]
646
647  \centerline{\includegraphics[width=75mm, height=75mm]
648    {examples/runupstart.eps}}
649
650  \caption{Bedslope example viewed with Swollen}
651  \label{fig:runupstart}
652\end{figure}
653
654
655\begin{figure}[hbt]
656
657  \centerline{
658    \includegraphics[width=75mm, height=75mm]{examples/runupduring.eps}
659    \includegraphics[width=75mm, height=75mm]{examples/runupend.eps}
660   }
661
662  \caption{Bedslope example viewed with Swollen}
663  \label{fig:bedslope2}
664\end{figure}
665
666
667
668
669\clearpage
670
671%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
672
673\section{An Example with Real Data}
674
675The following discussion builds on the concepts introduced through
676the \file{runup.py} example and introduces a second
677example, \file{run\_sydney\_smf.py}, that follows the same basic
678outline, but incorporates more complex features and refers to a
679real-life scenario, rather than the artificial illustrative one used
680in \file{runup.py}. The domain of interest surrounds the
681Sydney region, and predominantly covers Sydney Harbour. A
682hypothetical tsunami wave is generated by a submarine mass failure
683situated on the edge of the continental shelf.
684
685\subsection{Overview}
686As in the case of \file{runup.py}, the actions carried
687out by the program can be organised according to this outline:
688
689\begin{enumerate}
690
691   \item Set up a triangular mesh.
692
693   \item Set certain parameters governing the mode of
694operation of the model---specifying, for instance, where to store the
695model output.
696
697   \item Input various quantities describing physical measurements, such
698as the elevation, to be specified at each mesh point (vertex).
699
700   \item Set up the boundary conditions.
701
702   \item Carry out the evolution of the model through a series of time
703steps and outputs the results, providing a results file that can be
704visualised.
705
706\end{enumerate}
707
708
709
710\subsection{The Code}
711
712Here is the code for \file{run\_sydney\_smf.py}:
713
714%\verbatiminput{"runsydneysmf.py"}
715\verbatiminput{examples/runsydney.py}
716
717In discussing the details of this example, we follow the outline
718given above, discussing each major step of the code in turn.
719
720\subsection{Establishing the Mesh}
721
722One obvious way that the present example differs from
723\file{runup.py} is in the use of a more complex method to
724create the mesh. Instead of imposing a mesh structure on a
725rectangular grid, the technique used for this example involves
726building mesh structures inside polygons specified by the user,
727using a mesh-generator referred to as \code{pmesh}.
728
729The following remarks may help the reader understand how
730\code{pmesh} is used.
731
732In its simplest form, \code{pmesh} creates the mesh within a single
733polygon whose vertices are at geographical locations specified by the
734user. The user specifies the \emph{resolution}---that is, the maximal
735area of a triangle used for triangulation---and mesh points are
736created inside the polygon through a random process. Figure
737\ref{fig:pentagon} shows a simple example of this, in which
738the triangulation is carried out within a pentagon.
739
740
741\begin{figure}[hbt]
742
743
744
745  \caption{Mesh points are created inside the polygon}
746  \label{fig:pentagon}
747\end{figure}
748
749Boundary tags are not restricted to \code{`left'}, \code{`right'},
750\code{`bottom'} and \code{`top'}, as in the case of
751\file{runup.py}. Instead the user specifies a list of
752tags appropriate to the configuration being modelled.
753
754While a mesh created inside a polygon offers more flexibility than
755one based on a rectangular grid, using \code{pmesh} in the limited
756form we have described so far still doesn't provide a way to adapt
757to geographic or other features in the landscape, whose presence may
758require us to vary the resolution in the neighbourhood of the
759features. To cope with this requirement, \code{pmesh} also allows
760the user to specify a number of \emph{interior polygons}, which are
761triangulated separately, each according to a separately specified
762resolution. See Figure \ref{fig:interior meshes}. It is also
763possible to specify one or more `holes'---that is, areas bounded by
764polygons in which no triangulation is required.
765
766\begin{figure}[hbt]
767
768
769
770  \caption{Interior meshes with individual resolution}
771  \label{fig:interior meshes}
772\end{figure}
773
774In its general form, \code{pmesh} takes for its input a bounding
775polygon and (optionally) a list of interior polygons. The user
776specifies resolutions, both for the bounding polygon and for each of
777the interior polygons. Given this data, \code{pmesh} first creates a
778triangular mesh inside the bounding polygon, using the specified
779resolution, and then creates a separate triangulation inside each of
780the interior polygons, using the resolution specified for that
781polygon.
782
783The function used to implement this process is
784\function{create\_mesh\_from\_regions}. Its arguments include the
785bounding polygon and its resolution, a list of boundary tags, and a
786list of pairs \code{[polygon, resolution]}, specifying the interior
787polygons and their resolutions.
788
789In practice, the details of the polygons used are read from a
790separate file \file{project.py}. The resulting mesh is output to a
791\emph{meshfile}\index{meshfile}. This term is used to describe a
792file of a specific format used to store the data specifying a mesh.
793(There are in fact two possible formats for such a file: it can
794either be a binary file, with extension \code{.msh}, or an ASCII
795file, with extension \code{.tsh}. In the present case, the binary
796file format \code{.msh} is used. See Section \ref{sec:file formats}
797(page \pageref{sec:file formats}) for more on file formats.)
798\code{pmesh} assigns a name to the file by appending the extension
799\code{.msh} to the name specified in the input file
800\file{project.py}. This name is stored in the variable
801\code{meshname}.
802
803The statements
804
805{\small \begin{verbatim}
806    interior_res = 5000%
807    interior_regions = [[project.harbour_polygon_2, interior_res],
808                    [project.botanybay_polygon_2, interior_res]]
809\end{verbatim}}
810
811are used to read in the specific polygons \code{project.harbour\_polygon\_2} and
812\code{botanybay\_polygon\_2} from \file{project.py} and assign a
813common resolution of 5000 to each. The statement
814
815{\small \begin{verbatim}
816    create_mesh_from_regions(project.diffpolygonall,%
817                         boundary_tags= {'bottom': [0],%
818                                         'right1': [1],%
819                                         'right0': [2],%
820                                         'right2': [3],%
821                                         'top': [4],%
822                                         'left1': [5],%
823                                         'left2': [6],%
824                                         'left3': [7]},%
825                         maximum_triangle_area=100000,%
826                         filename=meshname,%
827                         interior_regions=interior_regions)
828\end{verbatim}}
829
830is then used to create the mesh, taking the bounding polygon to be the polygon
831\code{diffpolygonall} specified in \file{project.py}. The
832argument \code{boundary\_tags} assigns a dictionary, whose keys are the
833names of the boundary tags used for the bounding polygon---\code{`bottom'},
834`right0', `right1', `right2', `top', `left1', `left2' and `left3'---
835and whose values identify the indices of the segments associated with each of these
836tags. (The value associated with each boundary tag is a one-element list.)
837
838
839\subsection{Initialising the Domain}
840
841As with \file{runup.py}, once we have created the mesh, the next
842step is to create the data structure \code{domain}. We did this for
843\file{runup.py} by inputting lists of points and triangles and
844specifying the boundary tags directly. However, in the present case,
845we use a method that works directly with the meshfile
846\code{meshname}, as follows:
847
848
849{\small \begin{verbatim}
850    domain = Domain(meshname, use_cache=True, verbose=True)
851\end{verbatim}}
852
853Providing a filename instead of the lists used in \file{runup.py}
854above causes \code{Domain} to convert a meshfile \code{meshname}
855into an instance of the data structure \code{domain}, allowing us to
856use methods like \method{set\_quantity} to set quantities and to
857apply other operations.
858
859%(In principle, the
860%second argument of \function{pmesh\_to\_domain\_instance} can be any
861%subclass of \class{Domain}, but for applications involving the
862%shallow-water wave equation, the second argument of
863%\function{pmesh\_to\_domain\_instance} can always be set simply to
864%\class{Domain}.)
865
866The following statements specify a basename and data directory, and
867identify quantities to be stored. For the first two, values are
868taken from \file{project.py}.
869
870{\small \begin{verbatim}
871    domain.set_name(project.basename)%
872    domain.set_datadir(project.outputdir)%
873    domain.set_quantities_to_be_stored(['stage', 'xmomentum',
874        'ymomentum'])
875\end{verbatim}}
876
877
878\subsection{Initial Conditions}
879Quantities for \file{runsydney.py} are set
880using similar methods to those in \file{runup.py}. However,
881in this case, many of the values are read from the auxiliary file
882\file{project.py} or, in the case of \code{elevation}, from an
883ancillary points file.
884
885
886
887\subsubsection{Stage}
888
889For the scenario we are modelling in this case, we use a callable
890object \code{tsunami\_source}, assigned by means of a function
891\function{slump\_tsunami}. This is similar to how we set elevation in
892\file{runup.py} using a function---however, in this case the
893function is both more complex and more interesting.
894
895The function returns the water displacement for all \code{x} and
896\code{y} in the domain. The water displacement is a double Gaussian
897function that depends on the characteristics of the slump (length,
898thickness, slope, etc), its location (origin) and the depth at that
899location.
900
901
902\subsubsection{Friction}
903
904We assign the friction exactly as we did for \file{runup.py}:
905
906{\small \begin{verbatim}
907    domain.set_quantity('friction', 0.03)
908\end{verbatim}}
909
910
911\subsubsection{Elevation}
912
913The elevation is specified by reading data from a file:
914
915{\small \begin{verbatim}
916    domain.set_quantity('elevation',
917                        filename = project.combineddemname + '.pts',
918                        use_cache = True,
919                        verbose = True)
920\end{verbatim}}
921
922However, before this step can be executed, some preliminary steps
923are needed to prepare the file from which the data is taken. Two
924source files are used for this data---their names are specified in
925the file \file{project.py}, in the variables \code{coarsedemname}
926and \code{finedemname}. They contain `coarse' and `fine' data,
927respectively---that is, data sampled at widely spaced points over a
928large region and data sampled at closely spaced points over a
929smaller subregion. The data in these files is combined through the
930statement
931
932{\small \begin{verbatim}
933combine_rectangular_points_files(project.finedemname + '.pts',
934                                 project.coarsedemname + '.pts',
935                                 project.combineddemname + '.pts')
936\end{verbatim}}
937
938The effect of this is simply to combine the datasets by eliminating
939any coarse data associated with points inside the smaller region
940common to both datasets. The name to be assigned to the resulting
941dataset is also derived from the name stored in the variable
942\code{combinedname} in the file \file{project.py}.
943
944\subsection{Boundary Conditions}
945
946Setting boundaries follows a similar pattern to the one used for
947\file{runup.py}, except that in this case we need to associate a
948boundary type with each of the
949boundary tag names introduced when we established the mesh. In place of the four
950boundary types introduced for \file{runup.py}, we use the reflective
951boundary for each of the
952eight tagged segments:
953
954{\small \begin{verbatim}
955    Br = Reflective_boundary(domain)
956    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
957                          'right2': Br, 'top': Br, 'left1': Br,
958                          'left2': Br, 'left3': Br} )
959\end{verbatim}}
960
961\subsection{Evolution}
962
963With the basics established, the running of the `evolve' step is
964very similar to the corresponding step in \file{runup.py}:
965
966{\small \begin{verbatim}
967    import time t0 = time.time()
968
969    for t in domain.evolve(yieldstep = 120, duration = 18000):
970        print domain.timestepping_statistics()
971        print domain.boundary_statistics(tags = 'bottom')
972
973    print 'That took %.2f seconds' %(time.time()
974\end{verbatim}}
975
976%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
977%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
978
979\chapter{\anuga Public Interface}
980\label{ch:interface}
981
982This chapter gives an overview of the features of \anuga available
983to the user at the public interface. These are grouped under the
984following headings:
985
986\begin{itemize}
987    \item Establishing the Mesh
988    \item Initialising the Domain
989    \item Specifying the Quantities
990    \item Initial Conditions
991    \item Boundary Conditions
992    \item Forcing Functions
993    \item Evolution
994\end{itemize}
995
996The listings are intended merely to give the reader an idea of what
997each feature is, where to find it and how it can be used---they do
998not give full specifications. For these the reader
999may consult the code. The code for every function or class contains
1000a documentation string, or `docstring', that specifies the precise
1001syntax for its use. This appears immediately after the line
1002introducing the code, between two sets of triple quotes.
1003
1004Each listing also describes the location of the module in which
1005the code for the feature being described can be found. All modules
1006are in the folder \file{inundation} or one of its subfolders, and the
1007location of each module is described relative to \file{inundation}. Rather
1008than using pathnames, whose syntax depends on the operating system,
1009we use the format adopted for importing the function or class for
1010use in Python code. For example, suppose we wish to specify that the
1011function \function{create\_mesh\_from\_regions} is in a module called
1012\module{mesh\_interface} in a subfolder of \module{inundation} called
1013\code{pmesh}. In Linux or Unix syntax, the pathname of the file
1014containing the function, relative to \file{inundation}, would be
1015
1016\begin{center}
1017%    \code{pmesh/mesh\_interface.py}
1018    \code{pmesh}$\slash$\code{mesh\_interface.py}
1019\end{center}
1020
1021while in Windows syntax it would be
1022
1023\begin{center}
1024    \code{pmesh}$\backslash$\code{mesh\_interface.py}
1025\end{center}
1026
1027Rather than using either of these forms, in this chapter we specify
1028the location simply as \code{pmesh.mesh_interface}, in keeping with
1029the usage in the Python statement for importing the function,
1030namely:
1031\begin{center}
1032    \code{from pmesh.mesh\_interface import create\_mesh\_from\_regions}
1033\end{center}
1034
1035Each listing details the full set of parameters for the class or
1036function; however, the description is generally limited to the most
1037important parameters and the reader is again referred to the code
1038for more details.
1039
1040The following parameters are common to many functions and classes
1041and are omitted from the descriptions given below:
1042
1043%\begin{center}
1044\begin{tabular}{ll}  %\hline
1045%\textbf{Name } & \textbf{Description}\\
1046%\hline
1047\emph{usecache} & Specifies whether caching is to be used\\
1048\emph{verbose} & If \code{True}, provides detailed terminal output
1049to the user\\  % \hline
1050\end{tabular}
1051%\end{center}
1052
1053\section{Mesh Generation}
1054\refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}\label{sec:meshgeneration}
1055
1056\begin{classdesc}  {Mesh}{userSegments=None,
1057                 userVertices=None,
1058                 holes=None,
1059                 regions=None,
1060                 geo_reference=None}
1061Module: \module{pmesh.mesh}
1062
1063A class used to store a representation of a two-dimensional
1064triangular mesh. The user may initialise the class by specifying
1065lists of segments and/or vertices, using the parameters
1066\code{userSegments} and \code{userVertices}, and may also specify
1067lists of regions and/or holes, using the parameters \code{regions}
1068and \code{holes}.
1069\end{classdesc}
1070
1071\begin{funcdesc}  {create\_mesh\_from\_regions}{bounding_polygon,
1072                             boundary_tags,
1073                             maximum_triangle_area,
1074                             filename=None,
1075                             interior_regions=None,
1076                             poly_geo_reference=None,
1077                             mesh_geo_reference=None,
1078                             minimum_triangle_angle=28.0}
1079Module: \module{pmesh.mesh\_interface}
1080
1081This function allows a user to initiate the automatic creation of a
1082mesh inside a specified polygon. Among the parameters that can be
1083set are the \emph{resolution} (maximal area for any triangle in the
1084mesh) and the minimal angle allowable in any triangle. The user can
1085specify a number of internal polygons within each of which a
1086separate mesh is to be created, generally with a smaller resolution.
1087Additionally, the user specifies a list of boundary tags, one for
1088each edge of the bounding polygon.
1089\end{funcdesc}
1090
1091
1092\begin{funcdesc} {add_region}{}
1093
1094\end{funcdesc}
1095
1096\begin{funcdesc} {add_hole}{}
1097
1098\end{funcdesc}
1099
1100
1101\begin{funcdesc}  {add_region_from_polygon}{self, polygon, tags=None,
1102                                max_triangle_area=None, geo_reference=None}
1103Module: \module{pmesh.mesh}  Class: \class{Mesh}
1104
1105% Translate following into layman's language
1106This method is used to add a region to a \class{Mesh} instance.  The
1107region is described by the polygon passed in.  Additionally, the
1108user specifies a list of boundary tags, one for each edge of the
1109bounding polygon.
1110
1111
1112\end{funcdesc}
1113
1114\begin{funcdesc}  {add_hole_from_polygon}{self, polygon, tags=None,
1115    geo_reference=None}
1116Module: \module{pmesh.mesh.Mesh}
1117
1118% Translate following into layman's language
1119This method is used to add a `hole'---that is, a region where the
1120triangular mesh will not be generated---to a \class{Mesh} instance.
1121The region is described by the polygon passed in.  Additionally, the
1122user specifies a list of boundary tags, one for each edge of the
1123bounding polygon.
1124\end{funcdesc}
1125
1126\begin{funcdesc}  {generate_mesh}{self,
1127                      maximum_triangle_area=None,
1128                      minimum_triangle_angle=28.0,
1129                      verbose=False}
1130Module: \module{pmesh.mesh.Mesh}
1131
1132% Translate following into layman's language
1133This method is used to generate the triangular mesh.  The  maximal
1134area of any triangle in the mesh can be specified, along with the
1135minimum angle in any triangle.
1136\end{funcdesc}
1137
1138
1139\begin{funcdesc}  {export_mesh_file}{self,ofile}
1140Module: \module{pmesh.mesh.Mesh}
1141
1142% Translate following into layman's language
1143This method is used to save the mesh to a file. \code{ofile} is the
1144name of the mesh file to be written, including the extension.  Use
1145the extension \code{.msh} for the file to be in NetCDF format and
1146\code{.tsh} for the file to be ASCII format.
1147\end{funcdesc}
1148
1149%%%%%%
1150\section{Initialising the Domain}
1151
1152%Include description of the class Domain and the module domain.
1153
1154%FIXME (Ole): This is also defined in a later chapter
1155%\declaremodule{standard}{pyvolution.domain}
1156
1157\begin{classdesc} {Domain} {source=None,
1158                 triangles=None,
1159                 boundary=None,
1160                 conserved_quantities=None,
1161                 other_quantities=None,
1162                 tagged_elements=None,
1163                 geo_reference=None,
1164                 use_inscribed_circle=False,
1165                 mesh_filename=None,
1166                 use_cache=False,
1167                 verbose=False,
1168                 full_send_dict=None,
1169                 ghost_recv_dict=None,
1170                 processor=0,
1171                 numproc=1}
1172Module: \refmodule{pyvolution.domain}
1173
1174This class is used to create an instance of a data structure used to
1175store and manipulate data associated with a mesh. The mesh is
1176specified either by assigning the name of a meshfile to
1177\code{source} or by
1178\end{classdesc}
1179
1180\begin{funcdesc}  {pmesh\_to\_domain\_instance}{file_name, DomainClass, use_cache = False, verbose = False}
1181Module: \module{pyvolution.pmesh2domain}
1182
1183Once the initial mesh file has been created, this function is
1184applied to convert it to a domain object---that is, to a member of
1185the special Python class \class{Domain} (or a subclass of
1186\class{Domain}), which provides access to properties and methods
1187that allow quantities to be set and other operations to be carried
1188out.
1189
1190\code{file\_name} is the name of the mesh file to be converted,
1191including the extension. \code{DomainClass} is the class to be
1192returned, which must be a subclass of \class{Domain} having the same
1193interface as \class{Domain}---in practice, it can usually be set
1194simply to \class{Domain}.
1195
1196This is now superseded by Domain(mesh_filename).
1197\end{funcdesc}
1198
1199
1200\subsection{Key Methods of Domain}
1201
1202\begin{funcdesc} {set\_name}{name}
1203    Module: \refmodule{pyvolution.domain}, page \pageref{mod:pyvolution.domain}  %\code{pyvolution.domain}
1204
1205    Assigns the name \code{name} to the domain.
1206\end{funcdesc}
1207
1208\begin{funcdesc} {get\_name}{}
1209    Module: \module{pyvolution.domain}
1210
1211    Returns the name assigned to the domain by \code{set_name}. If no name has been
1212    assigned, returns \code{`domain'}.
1213\end{funcdesc}
1214
1215\begin{funcdesc} {set\_datadir}{name}
1216    Module: \module{pyvolution.domain}
1217
1218    Specifies the directory used for data, assigning it to the pathname \code{name}. The default value, before
1219    \code{set\_datadir} has been run, is the value \code{default_datadir}
1220    specified in \code{config.py}.
1221
1222    Since different operating systems use different formats for specifying pathnames,
1223    it is necessary to specify path separators using the Python code \code{os.sep}, rather than
1224    the operating-specific ones such as `$\slash$' or `$\backslash$'.
1225    For this to work you will need to include the statement \code{import os}
1226    in your code, before the first appearance of \code{set\_datadir}.
1227
1228    For example, to set the data directory to a subdirectory
1229    \code{data} of the directory \code{project}, you could use
1230    the statements:
1231
1232    {\small \begin{verbatim}
1233        import os
1234        domain.set_datadir{'project' + os.sep + 'data'}
1235    \end{verbatim}}
1236\end{funcdesc}
1237
1238\begin{funcdesc} {get_datadir}{}
1239    Module: \module{pyvolution.domain}
1240
1241    Returns the data directory set by \code{set_datadir} or,
1242    if \code{set_datadir} has not
1243    been run, returns the value \code{default_datadir} specified in
1244    \code{config.py}.
1245\end{funcdesc}
1246
1247\begin{funcdesc} {set_time}{time=0.0}
1248    Module: \module{pyvolution.domain}
1249
1250    Sets the initial time, in seconds, for the simulation. The
1251    default is 0.0.
1252\end{funcdesc}
1253
1254\begin{funcdesc} {set_default_order}{n}
1255    Sets the default (spatial) order to the value specified by
1256    \code{n}, which must be either 1 or 2. (Assigning any other value
1257    to \code{n} will cause an error.)
1258\end{funcdesc}
1259
1260
1261%%%%%%
1262\section{Initial Conditions}
1263
1264\begin{funcdesc}{set\_quantity}{name,
1265    numeric = None,
1266    quantity = None,
1267    function = None,
1268    geospatial_data = None,
1269    filename = None,
1270    attribute_name = None,
1271    alpha = None,
1272    location = 'vertices',
1273    indices = None,
1274    verbose = False,
1275    use_cache = False}
1276  Module: \module{pyvolution.domain}
1277  (see also \module{pyvolution.quantity.set_values})
1278
1279This function is used to assign values to individual quantities for a
1280domain. It is very flexible and can be used with many data types: a
1281statement of the form \code{domain.set\_quantity(name, x)} can be used
1282to define a quantity having the name \code{name}, where the other
1283argument \code{x} can be any of the following:
1284
1285\begin{itemize}
1286\item a number, in which case all vertices in the mesh gets that for
1287the quantity in question.
1288\item a list of numbers or a Numeric array ordered the same way as the mesh vertices.
1289\item a function (e.g.\ see the samples introduced in Chapter 2)
1290\item an expression composed of other quantities and numbers, arrays, lists (for
1291example, a linear combination of quantities)
1292\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.
1293\item a geospatial dataset (See ?????). Optional argument attribute_name applies here as with files.
1294\end{itemize}
1295
1296
1297Exactly one of the arguments
1298  numeric, quantity, function, points, filename
1299must be present.
1300
1301
1302Set quantity will look at the type of the second argument (\code{numeric}) and
1303determine what action to take.
1304
1305Values can also be set using the appropriate keyword arguments.
1306If 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)}
1307are all equivalent.
1308
1309
1310Other optional arguments are
1311\begin{itemize}
1312\item \code{indices} which is a list of ids of triangles to which set_quantity should apply its assignment of values.
1313\item \code{location} determines which part of the triangles to assign to. Options are 'vertices' (default), 'edges', and 'centroids'.
1314\end{itemize}
1315
1316
1317\end{funcdesc}
1318
1319
1320
1321
1322
1323
1324
1325%%%
1326\anuga provides a number of predefined initial conditions to be used
1327with \code{set_quantity}.
1328
1329\begin{funcdesc}{tsunami_slump}{length, depth, slope, width=None, thickness=None,
1330                x0=0.0, y0=0.0, alpha=0.0,
1331                gravity=9.8, gamma=1.85,
1332                massco=1, dragco=1, frictionco=0, psi=0,
1333                dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
1334                domain=None,
1335                verbose=False}
1336This function returns a callable object representing an initial water
1337displacement generated by a submarine sediment failure. These failures can take the form of
1338a submarine slump or slide.
1339
1340The arguments include as a minimum, the slump or slide length, the water depth to the centre of sediment
1341mass, and the bathymetric slope. Other slump or slide parameters can be included if they are known.
1342\end{funcdesc}
1343
1344
1345%%%
1346\begin{funcdesc}{file_function}{filename,
1347                  domain = None,
1348                  quantities = None,
1349                  interpolation_points = None,
1350                  verbose = False,
1351                  use_cache = False}
1352Module: \module{pyvolution.util}
1353
1354Reads the time history of spatial data from NetCDF file and returns
1355a callable object. Returns interpolated values based on the input
1356file using the underlying \code{interpolation_function}.
1357
1358\code{quantities} is either the name of a single quantity to be
1359interpolated or a list of such quantity names. In the second case, the resulting
1360function will return a tuple of values---one for each quantity.
1361
1362\code{interpolation_points} is a list of absolute UTM coordinates
1363for points at which values are sought.
1364
1365To get access to the model time stored within the file function use
1366the method \code{f.get_time()}
1367\end{funcdesc}
1368
1369%%%
1370\begin{classdesc}{Interpolation\_function}{self,
1371                 time,
1372                 quantities,
1373                 quantity_names = None,
1374                 vertex_coordinates = None,
1375                 triangles = None,
1376                 interpolation_points = None,
1377                 verbose = False}
1378Module: \module{pyvolution.least\_squares}
1379
1380Given a time series, either as a sequence of numbers or defined at
1381the vertices of a triangular mesh (such as those stored in SWW
1382files), \code{Interpolation\_function} is used to create a callable
1383object that interpolates a value for an arbitrary time \code{t}
1384within the model limits and possibly a point \code{(x, y)} within a
1385mesh region.
1386
1387The actual time series at which data is available is specified by
1388means of an array \code{time} of monotonically increasing times. The
1389quantities containing the values to be interpolated are specified in
1390an array---or dictionary of arrays (used in conjunction with the
1391optional argument \code{quantity\_names}) --- called
1392\code{quantities}. The optional arguments \code{vertex_coordinates}
1393and \code{triangles} represent the spatial mesh associated with the
1394quantity arrays. If omitted the function created by
1395\code{Interpolation\_function} will be a function of \code{t} only.
1396
1397Since, in practice, values need to be computed at specified points,
1398the syntax allows the user to specify, once and for all, a list
1399\code{interpolation\_points} of points at which values are required.
1400In this case, the function may be called using the form \code{f(t,
1401id)}, where \code{id} is an index for the list
1402\code{interpolation\_points}.
1403
1404\end{classdesc}
1405
1406%%%
1407%\begin{funcdesc}{set\_region}{functions}
1408%[Low priority. Will be merged into set\_quantity]
1409
1410%Module:\module{pyvolution.domain}
1411%\end{funcdesc}
1412
1413
1414
1415%%%%%%
1416\section{Boundary Conditions}
1417
1418\anuga provides a large number of predefined boundary conditions,
1419represented by objects such as \code{Reflective\_boundary(domain)} and
1420\code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples
1421in Chapter 2. Alternatively, you may prefer to ``roll your own'',
1422following the method explained in Section \ref{sec:roll your own}.
1423
1424These boundary objects may be used with the function \code{set\_boundary} described below
1425to assign boundary conditions according to the tags used to label boundary segments.
1426
1427\begin{funcdesc}{set\_boundary}{boundary_map}
1428Module: \module{pyvolution.domain}
1429
1430This function allows you to assign a boundary object (corresponding to a
1431pre-defined or user-specified boundary condition) to every boundary segment that
1432has been assigned a particular tag.
1433
1434This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects
1435and whose keys are the symbolic tags.
1436
1437\end{funcdesc}
1438
1439\begin{funcdesc} {get_boundary_tags}{}
1440Module: \module{pyvolution.mesh}
1441
1442Returns a list of the available boundary tags.
1443\end{funcdesc}
1444
1445%%%
1446\subsection{Predefined boundary conditions}
1447
1448\begin{classdesc}{Reflective_boundary}{Boundary}
1449Module: \module{pyvolution.shallow\_water}
1450
1451Reflective boundary returns same conserved quantities as those present in
1452the neighbouring volume but reflected.
1453
1454This class is specific to the shallow water equation as it works with the
1455momentum quantities assumed to be the second and third conserved quantities.
1456\end{classdesc}
1457
1458%%%
1459\begin{classdesc}{Transmissive_boundary}{domain = None}
1460Module: \module{pyvolution.generic\_boundary\_conditions}
1461
1462A transmissive boundary returns the same conserved quantities as
1463those present in the neighbouring volume.
1464
1465The underlying domain must be specified when the boundary is instantiated.
1466\end{classdesc}
1467
1468%%%
1469\begin{classdesc}{Dirichlet_boundary}{conserved_quantities=None}
1470Module: \module{pyvolution.generic\_boundary\_conditions}
1471
1472A Dirichlet boundary returns constant values for the conserved
1473quantities.
1474\end{classdesc}
1475
1476%%%
1477\begin{classdesc}{Time_boundary}{domain = None, f = None}
1478Module: \module{pyvolution.generic\_boundary\_conditions}
1479
1480A time-dependent boundary returns values for the conserved
1481quantities as a function \code{f(t)} of time. The user must specify
1482the domain to get access to the model time.
1483\end{classdesc}
1484
1485%%%
1486\begin{classdesc}{File_boundary}{Boundary}
1487Module: \module{pyvolution.generic\_boundary\_conditions}
1488
1489The boundary values are obtained from a file and interpolated. The
1490file is assumed to contain a time series and possibly also spatial
1491information. The conserved quantities are given as a function of
1492time.
1493\end{classdesc}
1494
1495
1496\subsection{User-defined boundary conditions}
1497\label{sec:roll your own}
1498[How to roll your own]
1499
1500
1501
1502
1503
1504\section{Forcing Functions}
1505
1506\anuga provides a number of predefined forcing functions to be used with .....
1507
1508%\begin{itemize}
1509
1510
1511%  \item \indexedcode{}
1512%  [function, arguments]
1513
1514%  \item \indexedcode{}
1515
1516%\end{itemize}
1517
1518
1519
1520\section{Evolution}
1521
1522  \begin{funcdesc}{evolve}{yieldstep = None, finaltime = None, duration = None, skip_initial_step = False}
1523
1524  Module: \module{pyvolution.domain}
1525
1526  This function (a method of \class{domain}) is invoked once all the
1527  preliminaries have been completed, and causes the model to progress
1528  through successive steps in its evolution, storing results and
1529  outputting statistics whenever a user-specified period
1530  \code{yieldstep} is completed (generally during this period the
1531  model will evolve through several steps internally). The user
1532  specifies the total time period over which the evolution is to take
1533  place, by specifying values (in seconds) for either \code{duration}
1534  or \code{finaltime}, as well as the interval in seconds after which
1535  results are to be stored and statistics output.
1536
1537  You can include \method{evolve} in a statement of the type:
1538
1539  {\small \begin{verbatim}
1540      for t in domain.evolve(yieldstep, finaltime):
1541          <Do something with domain and t>
1542  \end{verbatim}}
1543
1544  \end{funcdesc}
1545
1546
1547
1548\subsection{Diagnostics}
1549
1550  \begin{funcdesc}{timestepping_statistics}{}
1551  Module: \module{pyvolution.domain}
1552
1553
1554  \end{funcdesc}
1555
1556
1557  \begin{funcdesc}{boundary\_statistics}{quantities = None, tags = None}
1558  Module: \module{pyvolution.domain}
1559
1560
1561  \end{funcdesc}
1562
1563
1564  \begin{funcdesc}{get_quantity}{name, location='vertices', indices = None}
1565  Module: \module{pyvolution.domain}
1566  Allow access to individual quantities and their methods
1567
1568  \end{funcdesc}
1569
1570
1571  \begin{funcdesc}{get_values}{location='vertices', indices = None}
1572  Module: \module{pyvolution.quantity}
1573
1574  Extract values for quantity as an array
1575
1576  \end{funcdesc}
1577
1578
1579  \begin{funcdesc}{get_integral}{}
1580  Module: \module{pyvolution.quantity}
1581
1582  Return computed integral over entire domain for this quantity
1583
1584  \end{funcdesc}
1585
1586
1587\section{Other}
1588
1589  \begin{funcdesc}{domain.create_quantity_from_expression}{???}
1590
1591  Handy for creating derived quantities on-the-fly.
1592  See \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
1593  \end{funcdesc}
1594
1595
1596  \begin{classdesc}{Geospatial_data}{data_points = None,
1597                 attributes = None,
1598                 geo_reference = None,
1599                 default_attribute_name = None,
1600                 file_name = None}
1601    Module: \module{geospatial_data.geo_spatial_data}
1602    Creates a georeferenced geospatial data object from either arrays or
1603    a file (pts or xya).
1604
1605    Objects of this class can be used with \method{set\_quantity}.
1606
1607    FIXME (Ole): Describe methods such as get_attributes() etc
1608  \end{classdesc}
1609
1610
1611
1612
1613
1614%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1615%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1616
1617\chapter{\anuga System Architecture}
1618
1619From pyvolution/documentation
1620
1621\section{File Formats}
1622\label{sec:file formats}
1623
1624\anuga makes use of a number of different file formats. The
1625following table lists all these formats, which are described in more
1626detail in the paragraphs below.
1627
1628\bigskip
1629
1630\begin{center}
1631
1632\begin{tabular}{|ll|}  \hline
1633
1634\textbf{Extension} & \textbf{Description} \\
1635\hline\hline
1636
1637\code{.sww} & NetCDF format for storing model output
1638\code{f(t,x,y)}\\
1639
1640\code{.tms} & NetCDF format for storing time series \code{f(t)}\\
1641
1642\code{.xya} & ASCII format for storing arbitrary points and
1643associated attributes\\
1644
1645\code{.pts} & NetCDF format for storing arbitrary points and
1646associated attributes\\
1647
1648\code{.asc} & ASCII format of regular DEMs as output from ArcView\\
1649
1650\code{.prj} & Associated ArcView file giving more metadata for
1651\code{.asc} format\\
1652
1653\code{.ers} & ERMapper header format of regular DEMs for ArcView\\
1654
1655\code{.dem} & NetCDF representation of regular DEM data\\
1656
1657\code{.tsh} & ASCII format for storing meshes and associated
1658boundary and region info\\
1659
1660\code{.msh} & NetCDF format for storing meshes and associated
1661boundary and region info\\
1662
1663\code{.nc} & Native ferret NetCDF format\\
1664
1665\code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
1666%\caption{File formats used by \anuga}
1667\end{tabular}
1668
1669
1670\end{center}
1671
1672The above table shows the file extensions used to identify the
1673formats of files. However, typically, in referring to a format we
1674capitalise the extension and omit the initial full stop---thus, we
1675refer, for example, to `SWW files' or `PRJ files'.
1676
1677\bigskip
1678
1679A typical dataflow can be described as follows:
1680
1681\subsection{Manually Created Files}
1682
1683\begin{tabular}{ll}
1684ASC, PRJ & Digital elevation models (gridded)\\
1685TSH & Triangular meshes (e.g. created from \code{pmesh})\\
1686NC & Model outputs for use as boundary conditions (e.g. from MOST)
1687\end{tabular}
1688
1689\subsection{Automatically Created Files}
1690
1691\begin{tabular}{ll}
1692ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert
1693DEMs to native \code{.pts} file\\
1694
1695NC $\rightarrow$ SWW & Convert MOST boundary files to
1696boundary \code{.sww}\\
1697
1698PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\
1699
1700TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using
1701Swollen\\
1702
1703TSH + Boundary SWW $\rightarrow$ SWW & Simulation using
1704\code{pyvolution}
1705\end{tabular}
1706
1707
1708
1709
1710\bigskip
1711
1712\subsection{SWW and TMS Formats}
1713
1714The SWW and TMS formats are both NetCDF formats, and are of key
1715importance for \anuga.
1716
1717An SWW file is used for storing \anuga output and therefore pertains
1718to a set of points and a set of times at which a model is evaluated.
1719It contains, in addition to dimension information, the following
1720variables:
1721
1722\begin{itemize}
1723    \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays
1724    \item \code{elevation}, a Numeric array storing bed-elevations
1725    \item \code{volumes}, a list specifying the points at the vertices of each of the
1726    triangles
1727    % Refer here to the example to be provided in describing the simple example
1728    \item \code{time}, a Numeric array containing times for model
1729    evaluation
1730\end{itemize}
1731
1732
1733The contents of an SWW file may be viewed using the visualisation
1734tool \code{swollen}, which creates an on-screen geometric
1735representation. See section \ref{sec:swollen} (page
1736\pageref{sec:swollen}) in Appendix \ref{ch:supportingtools} for more
1737on \code{swollen}.
1738
1739Alternatively, there are tools, such as \code{ncdump}, that allow
1740you to convert an NetCDF file into a readable format such as the
1741Class Definition Language (CDL). The following is an excerpt from a
1742CDL representation of the output file \file{bedslope.sww} generated
1743from running the simple example \file{runup.py} of
1744Chapter \ref{ch:getstarted}:
1745
1746\verbatiminput{examples/bedslopeexcerpt.cdl}
1747
1748The SWW format is used not only for output but also serves as input
1749for functions such as \function{file_boundary} and
1750\function{file_function}, described in Chapter \ref{ch:interface}.
1751
1752A TMS file is used to store time series data that is independent of
1753position.
1754
1755
1756\subsection{Meshfile Formats}
1757
1758A meshfile is a file that has a specific format suited to specifying
1759mesh data for \anuga. A meshfile can have one of two formats: it can
1760be either a TSH file, which is an ASCII file, or an MSH file, which
1761is a NetCDF file. A meshfile can be generated from the function
1762\function{create_mesh_from_regions} (see Section
1763\ref{sec:meshgeneration}) and used to initialise a domain.
1764
1765A meshfile describes the outline of the mesh---the vertices and line
1766segments that enclose the region in which the mesh is created---and
1767the triangular mesh itself, which is specified by listing the
1768triangles and their vertices, and the segments, which are those
1769sides of the triangles that are associated with boundary conditions.
1770
1771In addition, a meshfile may contain `holes' and/or `regions'. A hole
1772or region is defined by specifying a point and a number of segments
1773that enclose that point. A hole represents an area where no mesh is
1774to be created, while a region is a labelled area used for defining
1775properties of a mesh, such as friction values.
1776
1777A meshfile can also contain a georeference, which describes an
1778offset to be applied to $x$ and $y$ values---eg to the vertices.
1779
1780
1781\subsection{Formats for Storing Arbitrary Points and Attributes}
1782
1783An XYA file is used to store data representing arbitrary numerical
1784attributes associated with a set of points.
1785
1786The format for an XYA file is:
1787%\begin{verbatim}
1788
1789            first line:     \code{[attribute names]}\\
1790            other lines:  \code{x y [attributes]}\\
1791
1792            for example:\\
1793            \code{elevation, friction}\\
1794            \code{0.6, 0.7, 4.9, 0.3}\\
1795            \code{1.9, 2.8, 5, 0.3}\\
1796            \code{2.7, 2.4, 5.2, 0.3}
1797
1798        The first two columns are always implicitly assumed to be $x$, $y$ coordinates.
1799        Use the same delimiter for the attribute names and the data.
1800
1801        An XYA file can optionally end with lines of this type:
1802
1803            \code{\#geo reference}\\
1804            \code{56}\\
1805            \code{466600.0}\\
1806            \code{8644444.0}
1807
1808        Here the first number specifies the zone (in this case zone 56)  and other numbers specify the
1809        coordinates (in this case (466600.0, 8644444.0)) of the lower left corner.
1810
1811A PTS file is a NetCDF representation of the data held in an XYA
1812file. If the data is associated with a set of $N$ points, then the
1813data is stored using an $N \times 2$ Numeric array of float
1814variables for the points and an $N \times 1$ Numeric array for each
1815attribute.
1816
1817%\end{verbatim}
1818
1819\subsection{ArcView Formats}
1820
1821Files of the three formats ASC, PRJ and ERS are all associated with
1822data from ArcView.
1823
1824An ASC file is an ASCII representation of DEM output from ArcView.
1825It has the following format...
1826
1827A PRJ file is an ArcView file used in conjunction with an ASC file
1828to represent metadata for a DEM.
1829
1830
1831\subsection{DEM Format}
1832
1833A DEM file is a NetCDF representation of regular DEM data.
1834
1835
1836\subsection{Other Formats}
1837
1838
1839
1840
1841\subsection{Basic File Conversions}
1842
1843  \begin{funcdesc}{sww2dem}{basename_in, basename_out = None,
1844            quantity = None,
1845            timestep = None,
1846            reduction = None,
1847            cellsize = 10,
1848            NODATA_value = -9999,
1849            easting_min = None,
1850            easting_max = None,
1851            northing_min = None,
1852            northing_max = None,
1853            expand_search = False,
1854            verbose = False,
1855            origin = None,
1856            datum = 'WGS84',
1857            format = 'ers'}
1858  Module: \module{pyvolution.data\_manager}
1859
1860  Takes data from an SWW file and converts it to DEM format (ASC or
1861  ERS)
1862  \end{funcdesc}
1863
1864
1865  \begin{funcdesc}{dem2pts}{basename_in, basename_out=None,
1866            easting_min=None, easting_max=None,
1867            northing_min=None, northing_max=None,
1868            use_cache=False, verbose=False}
1869  Module: \module{pyvolution.data\_manager}
1870
1871  Takes DEM data (a NetCDF file representation of data from a regular Digital
1872  Elevation Model) and converts it to PTS format.
1873  \end{funcdesc}
1874
1875
1876
1877%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1878%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1879
1880\chapter{Basic \anuga Assumptions}
1881
1882(From pyvolution/documentation)
1883
1884
1885Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
1886If one wished to recreate scenarios prior to that date it must be done
1887using some relative time (e.g. 0).
1888
1889
1890All spatial data relates to the WGS84 datum (or GDA94) and has been
1891projected into UTM with false easting of 500000 and false northing of
18921000000 on the southern hemisphere (0 on the northern).
1893
1894It is assumed that all computations take place within one UTM zone.
1895
1896DEMs, meshes and boundary conditions can have different origins within
1897one UTM zone. However, the computation will use that of the mesh for
1898numerical stability.
1899
1900
1901%OLD
1902%The dataflow is: (See data_manager.py and from scenarios)
1903%
1904%
1905%Simulation scenarios
1906%--------------------%
1907%%
1908%
1909%Sub directories contain scrips and derived files for each simulation.
1910%The directory ../source_data contains large source files such as
1911%DEMs provided externally as well as MOST tsunami simulations to be used
1912%as boundary conditions.
1913%
1914%Manual steps are:
1915%  Creation of DEMs from argcview (.asc + .prj)
1916%  Creation of mesh from pmesh (.tsh)
1917%  Creation of tsunami simulations from MOST (.nc)
1918%%
1919%
1920%Typical scripted steps are%
1921%
1922%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
1923%                   native dem and pts formats%
1924%
1925%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
1926%                  as boundary condition%
1927%
1928%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
1929%                   smoothing. The outputs are tsh files with elevation data.%
1930%
1931%  run_simulation.py: Use the above together with various parameters to
1932%                     run inundation simulation.
1933
1934
1935%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1936%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1937
1938\appendix
1939
1940\chapter{Supporting Tools}
1941\label{ch:supportingtools}
1942
1943This section describes a number of supporting tools, supplied with \anuga, that offer a
1944variety of types of functionality and enhance the basic capabilities of \anuga.
1945
1946\section{caching}
1947
1948  The \code{cache} function is used to provide supervised caching of function results. A Python
1949  function call of the form
1950
1951      {\small \begin{verbatim}
1952      result = func(arg1,...,argn)
1953      \end{verbatim}}
1954
1955  can be replaced by
1956
1957      {\small \begin{verbatim}
1958      from caching import cache
1959      result = cache(func,(arg1,...,argn))
1960      \end{verbatim}}
1961
1962  which returns the same output but reuses cached
1963  results if the function has been computed previously in the same context.
1964  \code{result} and the arguments can be simple types, tuples, list, dictionaries or
1965  objects, but not unhashable types such as functions or open file objects.
1966  The function \code{func} may be a member function of an object or a module.
1967
1968  This type of caching is particularly useful for computationally intensive
1969  functions with few frequently used combinations of input arguments. Note that
1970  if the inputs or output are very large caching may not save time because
1971  disc access may dominate the execution time.
1972
1973  If the function definition changes after a result has been cached, this will be
1974  detected by examining the functions \code{bytecode (co_code, co_consts,
1975  func_defaults, co_argcount)} and the function will be recomputed.
1976
1977  Options are set by means of the function \code{set_option(key, value)},
1978  where \code{key} is a key associated with a
1979  Python dictionary \code{options}. This dictionary stores settings such as the name of
1980  the directory used, the maximum
1981  number of cached files allowed, and so on.
1982
1983  The \code{cache} function allows the user also to specify a list of dependent files. If any of these
1984  have been changed, the function is recomputed and the results stored again.
1985
1986  %Other features include support for compression and a capability to \ldots
1987
1988
1989   \textbf{USAGE:}
1990
1991    {\small \begin{verbatim}
1992    result = cache(func, args, kwargs, dependencies, cachedir, verbose,
1993                   compression, evaluate, test, return_filename)
1994    \end{verbatim}}
1995
1996
1997\section{swollen}
1998\label{sec:swollen}
1999 The output generated by \anuga may be viewed by
2000means of the visualisation tool \code{swollen}, which takes the
2001\code{sww} file output by \anuga and creates a visual representation
2002of the data. Examples may be seen in Figures \ref{fig:runupstart}
2003and \ref{fig:bedslope2}. To view an \code{sww} file with
2004\code{swollen} in the Windows environment, you can simply drag the
2005icon representing the file over an icon on the desktop for the
2006\code{swollen} executable file (or a shortcut to it). Alternatively,
2007you can operate \code{swollen} from the command line, in both
2008Windows and Linux environments.
2009
2010On successful operation, you will see an interactive moving-picture display. You can use keys and the mouse
2011to slow down, speed up or stop the display, change the viewing position or carry out a number of other
2012simple operations.
2013
2014The main keys operating the interactive screen are:\\
2015
2016\begin{center}
2017\begin{tabular}{|ll|}   \hline
2018
2019\code{w} & toggle wireframe\\
2020
2021space bar & start/stop\\
2022
2023up/down arrows & increase/decrease speed\\
2024
2025left/right arrows & direction in time \emph{(when running)}\\ & step through simulation \emph{(when stopped)}\\
2026
2027left mouse button & rotate\\
2028
2029middle mouse button & pan\\
2030
2031right mouse button & zoom\\  \hline
2032
2033\end{tabular}
2034\end{center}
2035
2036\vfill
2037
2038The following table describes how to operate swollen from the command line:
2039
2040Usage: \code{swollen [options] swwfile \ldots}\\  \nopagebreak
2041Options:\\  \nopagebreak
2042\begin{tabular}{ll}
2043  \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY_CENTER |}\\
2044                                    & \code{HEAD_MOUNTED_DISPLAY}\\
2045  \code{--rgba} & Request a RGBA colour buffer visual\\
2046  \code{--stencil} & Request a stencil buffer visual\\
2047  \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
2048                                    & overridden by environmental variable\\
2049  \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD_BUFFER | HORIZONTAL_SPLIT |}\\
2050                                    & \code{VERTICAL_SPLIT | LEFT_EYE | RIGHT_EYE |}\\
2051                                     & \code{ON | OFF} \\
2052  \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
2053  \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
2054  \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
2055  \code{-help} & Display this information\\
2056  \code{-hmax <float>} & Height above which transparency is set to
2057                                     \code{alphamax}\\
2058  \code{-hmin <float>} & Height below which transparency is set to
2059                                     zero\\
2060  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
2061                                     up, default is overhead)\\
2062  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
2063  \code{-movie <dirname>} & Save numbered images to named directory and
2064                                     quit\\
2065  \code{-nosky} & Omit background sky\\
2066  \code{-scale <float>} & Vertical scale factor\\
2067  \code{-texture <file>} & Image to use for bedslope topography\\
2068  \code{-tps <rate>} & Timesteps per second\\
2069  \code{-version} & Revision number and creation (not compile)
2070                                     date\\
2071\end{tabular}
2072
2073\section{utilities/polygons}
2074
2075  \declaremodule{standard}{utilities.polygon}
2076  \refmodindex{utilities.polygon}
2077
2078  \begin{classdesc}{Polygon_function}{regions, default = 0.0, geo_reference = None}
2079  Module: \code{utilities.polygon}
2080
2081  Creates a callable object that returns one of a specified list of values when
2082  evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the
2083  point belongs to. The parameter \code{regions} is a list of pairs
2084  \code{(P, v)}, where each \code{P} is a polygon and each \code{v}
2085  is either a constant value or a function of coordinates \code{x}
2086  and \code{y}, specifying the return value for a point inside \code{P}. The
2087  optional parameter \code{default} may be used to specify a value
2088  for a point not lying inside any of the specified polygons. When a
2089  point lies in more than one polygon, the return value is taken to
2090  be the value for whichever of these polygon appears later in the
2091  list.
2092  %FIXME (Howard): CAN x, y BE VECTORS?
2093
2094  \end{classdesc}
2095
2096  \begin{funcdesc}{read_polygon}{filename}
2097  Module: \code{utilities.polygon}
2098
2099  Reads the specified file and returns a polygon. Each
2100  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
2101  as coordinates of one vertex of the polygon.
2102  \end{funcdesc}
2103
2104  \begin{funcdesc}{populate_polygon}{polygon, number_of_points, seed = None, exclude = None}
2105  Module: \code{utilities.polygon}
2106
2107  Populates the interior of the specified polygon with the specified number of points,
2108  selected by means of a uniform distribution function.
2109  \end{funcdesc}
2110
2111  \begin{funcdesc}{point_in_polygon}{polygon, delta=1e-8}
2112  Module: \code{utilities.polygon}
2113
2114  Returns a point inside the specified polygon and close to the edge. The distance between
2115  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
2116  second argument \code{delta}, which is taken as $10^{-8}$ by default.
2117  \end{funcdesc}
2118
2119  \begin{funcdesc}{inside_polygon}{points, polygon, closed = True, verbose = False}
2120  Module: \code{utilities.polygon}
2121
2122  Used to test whether the members of a list of points
2123  are inside the specified polygon. Returns a Numeric
2124  array comprising the indices of the points in the list that lie inside the polygon.
2125  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
2126  Points on the edges of the polygon are regarded as inside if
2127  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
2128  \end{funcdesc}
2129
2130  \begin{funcdesc}{outside_polygon}{points, polygon, closed = True, verbose = False}
2131  Module: \code{utilities.polygon}
2132
2133  Exactly like \code{inside_polygon}, but with the words `inside' and `outside' interchanged.
2134  \end{funcdesc}
2135
2136  \begin{funcdesc}{is_inside_polygon}{point, polygon, closed=True, verbose=False}
2137  Module: \code{utilities.polygon}
2138
2139  Returns \code{True} if \code{point} is inside \code{polygon} or
2140  \code{False} otherwise. Points on the edges of the polygon are regarded as inside if
2141  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
2142  \end{funcdesc}
2143
2144  \begin{funcdesc}{is_outside_polygon}{point, polygon, closed=True, verbose=False}
2145  Module: \code{utilities.polygon}
2146
2147  Exactly like \code{is_outside_polygon}, but with the words `inside' and `outside' interchanged.
2148  \end{funcdesc}
2149
2150  \begin{funcdesc}{point_on_line}{x, y, x0, y0, x1, y1}
2151  Module: \code{utilities.polygon}
2152
2153  Returns \code{True} or \code{False}, depending on whether the point with coordinates
2154  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
2155  and \code{x1, y1} (extended if necessary at either end).
2156  \end{funcdesc}
2157
2158  \begin{funcdesc}{separate_points_by_polygon}{points, polygon, closed = True, verbose = False}
2159    \indexedcode{separate_points_by_polygon}
2160  Module: \code{utilities.polygon}
2161
2162  \end{funcdesc}
2163
2164  \begin{funcdesc}{polygon_area}{polygon}
2165  Module: \code{utilities.polygon}
2166
2167  Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html)
2168  \end{funcdesc}
2169
2170  \begin{funcdesc}{plot_polygons}{polygons, figname, verbose = False}
2171  Module: \code{utilities.polygon}
2172
2173  Plots each polygon contained in input polygon list, e.g. [poly1, poly2, poly3].
2174 Each polygon is closed for plotting purposes and subsequent plot saved to figname.
2175  Returns list containing the minimum and maximum of $x$ and $y$,
2176  i.e. [$x_{min}, x_{max}, y_{min}, y_{max}$].
2177  \end{funcdesc}
2178
2179\section{coordinate_transforms}
2180
2181\section{geospatial_data}
2182
2183This describes a class that represents arbitrary point data in UTM
2184coordinates along with named attribute values.
2185
2186%FIXME (Ole): This gives a LaTeX error
2187%\declaremodule{standard}{geospatial_data}
2188%\refmodindex{geospatial_data}
2189
2190
2191
2192\begin{classdesc}{Geospatial_data}
2193  {data_points = None,
2194    attributes = None,
2195    geo_reference = None,
2196    default_attribute_name = None,
2197    file_name = None}
2198Module: \code{geospatial_data}
2199
2200This class is used to store a set of data points and associated
2201attributes, allowing these to be manipulated by methods defined for
2202the class.
2203
2204The data points are specified either by reading them from a NetCDF
2205or XYA file, identified through the parameter \code{file_name}, or
2206by providing their \code{x}- and \code{y}-coordinates in metres,
2207either as a sequence of 2-tuples of floats or as an $M \times 2$
2208Numeric array of floats, where $M$ is the number of points.
2209Coordinates are interpreted relative to the origin specified by the
2210object \code{geo_reference}, which contains data indicating the UTM
2211zone, easting and northing. If \code{geo_reference} is not
2212specified, a default is used.
2213
2214Attributes are specified through the parameter \code{attributes},
2215set either to a list or array of length $M$ or to a dictionary whose
2216keys are the attribute names and whose values are lists or arrays of
2217length $M$. One of the attributes may be specified as the default
2218attribute, by assigning its name to \code{default_attribute_name}.
2219If no value is specified, the default attribute is taken to be the
2220first one.
2221\end{classdesc}
2222
2223
2224\begin{funcdesc}{import_points_file}{delimiter = None, verbose = False}
2225
2226\end{funcdesc}
2227
2228
2229\begin{funcdesc}{export_points_file}{ofile, absolute=False}
2230
2231\end{funcdesc}
2232
2233
2234\begin{funcdesc}{get_data_points}{absolute = True}
2235
2236\end{funcdesc}
2237
2238
2239\begin{funcdesc}{set_attributes}{attributes}
2240
2241\end{funcdesc}
2242
2243
2244\begin{funcdesc}{get_attributes}{attribute_name = None}
2245
2246\end{funcdesc}
2247
2248
2249\begin{funcdesc}{get_all_attributes}{}
2250
2251\end{funcdesc}
2252
2253
2254\begin{funcdesc}{set_default_attribute_name}{default_attribute_name}
2255
2256\end{funcdesc}
2257
2258
2259\begin{funcdesc}{set_geo_reference}{geo_reference}
2260
2261\end{funcdesc}
2262
2263
2264\begin{funcdesc}{__add__}
2265
2266\end{funcdesc}
2267
2268
2269\section{pmesh GUI}
2270
2271\section{alpha_shape}
2272An \emph{alpha shape} is a shape created from a set of points in the
2273plane, according to an algorithm. This shape is intended to capture
2274the intuitive notion of the `shape formed by the points'. For a
2275sufficiently simple set, the alpha shape reduces to the more precise
2276notion of the convex hull, but in general an alpha shape can be much
2277much more complex and may be far from convex.
2278
2279In the context of \anuga, an alpha shape is used to fit a polygonal
2280boundary around a set of points when fitting a mesh to a region. The
2281algorithm uses a parameter $\alpha$ that can be adjusted to make the
2282resultant shape resemble the shape suggested by intuition more
2283closely.
2284
2285The following paragraphs describe the class used to model an alpha
2286shape and some of the more important methods and attributes
2287associated with instances of this class.
2288
2289\begin{classdesc}{Alpha_Shape}{points, alpha = None}
2290Module: \code{alpha_shape}
2291
2292To instantiate this class the user supplies the points from which
2293the alpha shape is to be created (in the form of a list of 2-tuples
2294\code{[[x1, y1],[x2, y2]}\ldots\code{]}, assigned to the parameter
2295\code{points}) and, optionally, a value for the parameter $\alpha$.
2296The alpha shape is then computed and details of the required
2297boundary may be retrieved as attributes of the class instance.
2298\end{classdesc}
2299
2300
2301\begin{funcdesc}{alpha_shape_via_files}{point_file, boundary_file, alpha= None}
2302Module: \code{alpha_shape}
2303
2304This function reads points from the specified file
2305\code{point_file}, computes the associated alpha shape (either using
2306the specified value for \code{alpha} or, if no value is specified,
2307automatically setting it to a value judged to be optimal} and
2308outputs the boundary to a file named \code{boundary_file}. This
2309output file lists the coordinates \code{x, y} of each point in the
2310boundary, using one line per point.
2311\end{funcdesc}
2312
2313
2314\begin{funcdesc}{set_boundary_type}{self,raw_boundary=True,
2315                          remove_holes=False,
2316                          smooth_indents=False,
2317                          expand_pinch=False,
2318                          boundary_points_fraction=0.2}
2319Module: \code{alpha_shape}
2320
2321This function sets flags that govern the operation of the algorithm
2322that computes the boundary,
2323        raw_boundary    Return raw boundary i.e. the regular edges of the
2324                        alpha shape.
2325        remove_holes    filter to remove small holes
2326                        (small is defined by  boundary_points_fraction )
2327        smooth_indents  remove sharp triangular indents in boundary
2328        expand_pinch    test for pinch-off and correct
2329                           i.e. a boundary vertex with more than two edges.
2330\end{funcdesc}
2331
2332\section{utilities/numerical_tools}
2333
2334\begin{tabular}{|l l|}  \hline
2335\code{angle(v1, v2=None)} & Angle between two-dimensional vectors
2336\code{v1} and \code{v2}, in range $0$ to $2\pi$. If \code{v2} is
2337omitted it is taken as the unit vector in the $x$-direction.\\
2338
2339\code{normal_vector(v)} & Normal vector to \code{v}\\
2340
2341\code{mean(x)} & Mean value of a vector\\
2342
2343\code{cov(x, y=None)} & Covariance of vectors \code{x} and \code{y}.
2344If y is None, returns \code{cov(x, x)}\\
2345
2346\code{err(x, y=0, n=2, relative=True)} & Relative error of
2347$\parallel$\code{x}$-$\code{y}$\parallel$ to
2348$\parallel$\code{y}$\parallel$ (2-norm if \code{n} = 2 or Max norm
2349if \code{n} = None). If denominator evaluates to zero or if \code{y}
2350is omitted, absolute error is returned.\\
2351
2352\code{norm(x)} & 2-norm of \code{x}\\
2353
2354\code{corr(x, y=None)} & Correlation of \code{x} and \code{y}. If
2355\code{y} is \code{None} returns autocorrelation of \code{x}.\\
2356
2357\code{ensure_numeric(A, typecode = None)} & Ensure that sequence is
2358a Numeric array. If \code{A} is already a Numeric array it will be
2359returned unaltered. Otherwise, an attempt is made to convert it to a
2360Numeric array. This function is necessary as \code{array(A)} can
2361cause memory overflow.\\
2362
2363\code{histogram(a, bins, relative=False)} & Standard histogram. If
2364\code{relative} is \code{True}, values will be normalised against
2365the total and thus represent frequencies rather than counts.\\
2366
2367\code{create_bins(data, number_of_bins = None)} & Safely create bins
2368for use with histogram. If data contains only one point or is
2369constant, one bin will be created. If \code{number_of_bins} is
2370omitted, 10 bins will be created.\\  \hline
2371
2372\end{tabular}
2373
2374
2375\begin{itemize}
2376  \item \indexedcode{ensure_numeric}
2377  \item \indexedcode{mean}
2378  \item
2379\end{itemize}
2380
2381
2382\chapter{Modules available in \anuga}
2383
2384
2385\section{\module{pyvolution.general\_mesh} }
2386\declaremodule[pyvolution.generalmesh]{}{pyvolution.general\_mesh}
2387\label{mod:pyvolution.generalmesh}
2388
2389\section{\module{pyvolution.mesh} }
2390\declaremodule{}{pyvolution.mesh}
2391\label{mod:pyvolution.mesh}
2392
2393\section{\module{pyvolution.domain} --- Generic module for 2D triangular domains for finite-volume computations of conservation laws}
2394\declaremodule{}{pyvolution.domain}
2395\label{mod:pyvolution.domain}
2396
2397
2398\section{\module{pyvolution.quantity}}
2399\declaremodule{}{pyvolution.quantity}
2400\label{mod:pyvolution.quantity}
2401
2402
2403\section{\module{pyvolution.shallow\_water} --- 2D triangular domains for finite-volume computations of the shallow water wave equation. This module contains a specialisation of class Domain from module domain.py consisting of methods specific to the Shallow Water Wave Equation
2404}
2405\declaremodule[pyvolution.shallowwater]{}{pyvolution.shallow\_water}
2406\label{mod:pyvolution.shallowwater}
2407
2408
2409
2410
2411%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2412
2413\chapter{Frequently Asked Questions}
2414
2415
2416\section{General Questions}
2417
2418\subsubsection{What is \anuga?}
2419
2420\subsubsection{Why is it called \anuga?}
2421
2422\subsubsection{How do I obtain a copy of \anuga?}
2423
2424\subsubsection{What developments are expected for \anuga in the future?}
2425
2426\subsubsection{Are there any published articles about \anuga that I can reference?}
2427
2428\section{Modelling Questions}
2429
2430\subsubsection{Which type of problems are \anuga good for?}
2431
2432\subsubsection{Which type of problems are beyond the scope of \anuga?}
2433
2434\subsubsection{Can I start the simulation at an arbitrary time?}
2435Yes, using \code{domain.set_time()} you can specify an arbitrary starting time.
2436This is for example useful in conjunction with a file_boundary, which may start hours before anything hits the model boundary. By assigning a later time for the model to start, computational resources aren't wasted.
2437
2438\subsubsection{Can I change values for any quantity during the simulation?}
2439Yes, using \code{domain.set_quantity()} inside the domain.evolve loop you
2440can change values of any quantity. This is for example useful if you wish to
2441let the system settle for a while before assigning an initial condition. Another example would be changing the values for elevation to model e.g. erosion.
2442
2443\subsubsection{Can I change boundary conditions during the simulation?}
2444Not sure, but it would be nice :-)
2445
2446\subsubsection{Why does a file\_function return a list of numbers when evaluated?}
2447Currently, file\_function works by returning values for the conserved
2448quantities \code{stage}, \code{xmomentum} and \code{ymomentum} at a given point in time and space as a triplet. To access e.g.\ \code{stage} one must specify element 0 of the triplet returned by file\_function.
2449
2450\subsubsection{Which diagnostics are available to troubleshoot a simulation?}
2451
2452\subsubsection{How do I use a DEM in my simulation?}
2453You use \code{dem2pts} to convert your DEM to the required .pts format. This .pts file is then called When setting the elevation data to the mesh
2454in \code{domain.set_quantity}
2455
2456\subsubsection{What sort of DEM resolution should I use?}
2457Try and work with the "best" you have available. Onshore DEMs are typically available in 25m, 100m and 250m grids. Note, offshore data is often sparse,
2458or non-existant.
2459
2460\subsubsection{What sort of mesh resolution should I use?}
2461The 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,
2462if 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,
2463you need a fine mesh over regions where the DEM changes rapidly, and other areas of significant interest, such as the coast.
2464
2465\subsubsection{How often should I store the output?}
2466This 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
2467to look in detail at the evolution, then you will need to balance against your storage requirements and the duration of the simulation.
2468If 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
2469quantities on a mesh with approximately 300000 triangles on a 2 min interval for 5 hours will result in approximately 350Mb sww file.
2470
2471\subsection{Boundary Conditions}
2472
2473\subsubsection{How do I create a Dirichlet boundary condition?}
2474
2475\subsubsection{How do I know which boundary tags are available?}
2476The method \code{domain.get_boundary_tags()} will return a list of
2477available tags for use with \code{domain.set_boundary_condition()}.
2478
2479
2480
2481
2482
2483\chapter{Glossary}
2484
2485\begin{tabular}{|l p{10cm}| c|}  \hline
2486    \emph{Term} & \emph{Definition} & \emph{Page}\\  \hline
2487
2488    \indexedbold{\anuga} & Name of software (joint development between ANU and
2489    GA) & \pageref{def:anuga}\\
2490
2491    \indexedbold{domain} & The domain of a function is the set of all input values to the
2492    function.&\\
2493
2494    \indexedbold{Dirichlet boundary} & A Dirichlet boundary condition imposed on a differential equation
2495 which specifies the values the solution is to take on the boundary of the
2496 domain.&\\
2497
2498    \indexedbold{elevation} & refers to bathymetry and topography&\\
2499
2500    \indexedbold{bathymetry} & offshore elevation &\\
2501
2502    \indexedbold{topography} & onshore elevation &\\
2503
2504    \indexedbold{evolution} & integration of the shallow water wave equations
2505    over time &\\
2506
2507    \indexedbold{forcing term} & &\\
2508
2509    \indexedbold{IDLE} & Development environment shipped with
2510    Python &\\
2511
2512    \indexedbold{Manning friction coefficient} & &\\
2513
2514    \indexedbold{mesh} & Triangulation of domain &\\
2515
2516    \indexedbold{meshfile} & [generic word for either .tsh or
2517    .msh file] & \\
2518
2519    \indexedbold{points file} & [generic word for either .pts or
2520    .xya file] & \\
2521
2522    \indexedbold{grid} & evenly spaced mesh & \\
2523
2524    \indexedbold{NetCDF} & &\\
2525
2526    \indexedbold{conserved quantity} & conserved (stage, x and y
2527    momentum) & \\
2528
2529    \indexedbold{reflective boundary} & &\\
2530
2531    \indexedbold{stage} & &\\
2532
2533%    \indexedbold{try this}
2534
2535    \indexedbold{swollen} & visualisation tool &\\
2536
2537    \indexedbold{time boundary} & defined in the manual (flog from
2538    there)&\\
2539
2540    \indexedbold{transmissive boundary} & defined in the manual (flog from
2541    there)&\\
2542
2543    \indexedbold{xmomentum} & conserved quantity (note, two-dimensional SWW equations say only x and y and NOT
2544    z) &\\
2545
2546    \indexedbold{ymomentum}  & conserved quantity & \\
2547
2548
2549
2550    \indexedbold{resolution} &  The maximal area of a triangular cell in a
2551    mesh & \\
2552
2553    \indexedbold{polygon} & A sequence of points in the plane. (Arbitrary polygons can be created
2554    in this way.)
2555    \anuga represents a polygon in one of two ways. One way is to represent it as a
2556    list whose members are either Python tuples
2557    or Python lists of length 2. The unit square, for example, would be represented by the
2558    list
2559    [ [0,0], [1,0], [1,1], [0,1] ]. The alternative is to represent it as an
2560    $N \times 2$ Numeric array, where $N$ is the number of points.
2561
2562    NOTE: More can be read in the module utilities/polygon.py .... &
2563    \\  \hline
2564
2565    \end{tabular}
2566
2567    \begin{tabular}{|l p{10cm}|c|}   \hline
2568
2569    \emph{Term} & \emph{Definition} & \emph{Page}\\  \hline
2570    \indexedbold{easting} & A rectangular (x,y) coordinate measurement of distance east from a north-south reference line,
2571usually a meridian used as the axis of origin within a map zone or
2572projection. Easting is a UTM (Universal Transverse Mercator)
2573Coordinate. & \\
2574
2575    \indexedbold{northing} & A rectangular (x,y) coordinate measurement of distance north from a north-south reference line,
2576usually a meridian used as the axis of origin within a map zone or
2577projection. Northing is a UTM (Universal Transverse Mercator)
2578Coordinate. & \\
2579
2580
2581    \indexedbold{latitude} & The angular distance on a mericlear north and south of the equator, expressed in degrees and
2582    minutes. & \\
2583
2584    \indexedbold{longitude} & The angular distance east or west, between the meridian of a particular place on Earth and that of the
2585Prime Meridian (located in Greenwich, England) expressed in degrees
2586or time.& \\
2587
2588    \indexedbold{edge} & A triangular cell within the computational mesh can be depicted as a set of vertices joined by lines (the
2589    edges). & \\
2590
2591    \indexedbold{vertex} & A point at which edges meet. & \\
2592
2593    \indexedbold{finite volume} & The method evaluates the terms in the shallow water wave equation as fluxes at the surfaces of each
2594finite volume. Because the flux entering a given volume is identical
2595to that leaving the adjacent volume, these methods are conservative.
2596Another advantage of the finite volume method is that it is easily
2597formulated to allow for unstructured meshes. The method is used in
2598many computational fluid dynamics packages. & \\
2599
2600
2601    \indexedbold{flux} & the amount of flow through the volume per unit
2602    time & \\
2603
2604    \indexedbold{Digital Elevation Model (DEM)} & DEMs are digital files consisting of points of elevations,
2605sampled systematically at equally spaced intervals.& \\  \hline
2606
2607
2608\end{tabular}
2609
2610%The \code{\e appendix} markup need not be repeated for additional
2611%appendices.
2612
2613
2614%
2615%  The ugly "%begin{latexonly}" pseudo-environments are really just to
2616%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
2617%  not really valuable.
2618%
2619%  If you don't want the Module Index, you can remove all of this up
2620%  until the second \input line.
2621%
2622
2623%begin{latexonly}
2624%\renewcommand{\indexname}{Module Index}
2625%end{latexonly}
2626\input{mod\jobname.ind}        % Module Index
2627%
2628%begin{latexonly}
2629%\renewcommand{\indexname}{Index}
2630%end{latexonly}
2631\input{\jobname.ind}            % Index
2632
2633
2634
2635\end{document}
Note: See TracBrowser for help on using the repository browser.