source: documentation/user_manual/anuga_user_manual.tex @ 2637

Last change on this file since 2637 was 2637, checked in by steve, 18 years ago
  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 61.5 KB
RevLine 
[2329]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
[2439]9%labels
10%Sections and subsections \label{sec: }
11%Chapters \label{ch: }
12%Equations \label{eq: }
13%Figures \label{fig: }
14
15
[2384]16
[2329]17\documentclass{manual}
[2274]18
[2499]19\usepackage{graphicx}
20\input{definitions}
21
[2384]22\title{\anuga User Manual}
[2284]23\author{Howard Silcock, Ole Nielsen, Duncan Gray, Jane Sexton}
[2274]24
[2329]25% Please at least include a long-lived email address;
26% the rest is at your discretion.
27\authoraddress{Geoscience Australia \\
28  Email: \email{ole.nielsen@ga.gov.au}
29}
30
[2377]31%Draft date
[2491]32\date{\today}   % update before release!
[2422]33                % Use an explicit date so that reformatting
34                % doesn't cause a new date to be used.  Setting
35                % the date to \today can be used during draft
36                % stages to make it easier to handle versions.
[2329]37
[2491]38\release{1.0}   % release version; this is used to define the
[2422]39                % \version macro
[2329]40
[2422]41\makeindex          % tell \index to actually write the .idx file
42%\makemodindex          % If this contains a lot of module sections.
[2329]43
[2467]44\setcounter{tocdepth}{3}
[2529]45\setcounter{secnumdepth}{3}
[2274]46\begin{document}
47\maketitle
48
[2329]49% This makes the contents more accessible from the front page of the HTML.
50\ifhtml
51\chapter*{Front Matter\label{front}}
52\fi
[2285]53
54%Subversion keywords:
55%
[2387]56%$LastChangedDate: 2006-03-30 09:20:39 +0000 (Thu, 30 Mar 2006) $
57%$LastChangedRevision: 2637 $
58%$LastChangedBy: steve $
[2285]59
[2383]60\input{copyright}
[2274]61
[2383]62
[2329]63\begin{abstract}
64
[2467]65\noindent \anuga\index{\anuga} is a hydrodynamic modelling tool that
66allows users to model realistic flow problems in complex geometries.
67Examples include dam breaks or the effects of natural hazards such
68as riverine flooding, storm surges and tsunami.
[2274]69
[2283]70The user must specify a study area represented by a mesh of triangular
71cells, the topography and bathymetry, frictional resistance, initial
[2422]72values for water level (called \emph{stage}\index{stage} within \anuga),
[2329]73boundary
[2283]74conditions and forces such as windstress or pressure gradients if
75applicable.
76
[2384]77\anuga tracks the evolution of water depth and horizontal momentum
[2283]78within each cell over time by solving the shallow water wave equation
79governing equation using a finite-volume method.
80
[2384]81\anuga cannot model details of breaking waves, flow under ceilings such
[2283]82as pipes, turbulence and vortices, vertical convection or viscous
83flows.
84
[2384]85\anuga also incorporates a mesh generator, called \code{pmesh}, that
[2283]86allows the user to set up the geometry of the problem interactively as
87well as tools for interpolation and surface fitting, and a number of
88auxiliary tools for visualising and interrogating the model output.
89
[2384]90Most \anuga components are written in the object-oriented programming
[2434]91language Python and most users will interact with \anuga by writing
[2384]92small Python programs based on the \anuga library
[2283]93functions. Computationally intensive components are written for
94efficiency in C routines working directly with the Numerical Python
95structures.
96
97
[2329]98\end{abstract}
[2283]99
[2329]100\tableofcontents
[2283]101
[2274]102
[2329]103\chapter{Introduction}
104
105
106\section{Purpose}
107
[2274]108The purpose of this user manual is to introduce the new user to
109the software, describe what it can do and give step-by-step
[2529]110instructions for setting up and running hydrodynamic simulations.
[2274]111
[2329]112\section{Scope}
[2274]113
114This manual covers only what is needed to operate the software
[2529]115after installation and configuration. It does not includes instructions for
[2274]116installing the software or detailed API documentation, both of
117which will be covered in separate publications.
118
[2329]119\section{Audience}
[2274]120
[2283]121Readers are assumed to be familiar with the operating environment
[2274]122and have a general understanding of the problem background, as
123well as enough programming experience to adapt the code to
124different requirements, as described in this manual,  and to
125understand the basic terminology of object-oriented programming.
126
[2329]127\section{Structure of This Manual}
[2274]128
129This manual is structured as follows:
130
131\begin{itemize}
[2384]132  \item Background (What \anuga does)
[2283]133  \item A \emph{Getting Started} section
[2434]134  \item A detailed description of the public interface
135  \item \anuga 's overall architecture, components and file formats
136  \item Assumptions
[2274]137\end{itemize}
138
139
[2329]140\pagebreak
141\chapter{Getting Started}
[2274]142
143This section is designed to assist the reader to get started with
[2434]144\anuga by working through simple examples. Two examples are discussed;
145the first is a simple but artificial example that is useful to illustrate
[2529]146many of the ideas, and the second is a more realistic example.
[2434]147
[2455]148\section{A Simple Example}
[2434]149
[2455]150\subsection{Overview}
151
[2529]152What follows is a discussion of the structure and operation of the file
153\code{bedslopephysical.py}, with just enough detail to allow the reader
[2274]154to appreciate what's involved in setting up a scenario like the
155one it depicts.
156
157This example carries out the solution of the shallow-water wave
158equation in the simple case of a configuration comprising a flat
159bed, sloping at a fixed angle in one direction and having a
[2283]160constant depth across each line in the perpendicular direction.
161
[2274]162The example demonstrates many of the basic ideas involved in
163setting up a more complex scenario. In the general case the user
164specifies the geometry (bathymetry and topography), the initial
165water level, boundary conditions such as tide, and any forcing
166terms that may drive the system such as wind stress or atmospheric
167pressure gradients. Frictional resistance from the different
168terrains in the model is represented by predefined forcing
[2283]169terms. The boundary is reflective on three sides and a time dependent wave on one side.
170
[2529]171The present example represents a simple scenario and does not
172include any forcing terms, nor is the data taken from a file as it
[2637]173would be in many typical cases.
[2274]174
[2529]175The conserved quantities involved in the
176problem are water depth, $x$-momentum and $y$-momentum. Other quantities
[2637]177involved in the computation are the friction, stage (absolute height of water surface)
[2529]178and elevation, the last two being related to
179the depth through the equation
180
[2637]181\begin{tabular}{rcrcl}
[2574]182  \code{stage} &=& \code{elevation} &+& \code{depth}
[2637]183\end{tabular}
[2529]184
[2574]185
[2274]186%\emph{[More details of the problem background]}
187
[2455]188\subsection{Outline of the Program}
[2274]189
[2529]190In outline, \code{bedslopephysical.py} performs the following steps:
[2274]191
192\begin{enumerate}
193
194   \item Sets up a triangular mesh.
195
196   \item Sets certain parameters governing the mode of
[2283]197operation of the model-specifying, for instance, where to store the model output.
[2274]198
199   \item Inputs various quantities describing physical measurements, such
[2283]200as the elevation, to be specified at each mesh point (vertex).
[2274]201
202   \item Sets up the boundary conditions.
203
204   \item Carries out the evolution of the model through a series of time
205steps and outputs the results, providing a results file that can
206be visualised.
207
208\end{enumerate}
209
[2455]210\subsection{The Code}
[2274]211
[2422]212%FIXME: we are using the \code function here.
[2434]213%This should be used wherever possible
[2274]214For reference we include below the complete code listing for
[2529]215\code{bedslopephysical.py}. Subsequent paragraphs provide a `commentary'
[2283]216that describes each step of the program and explains it significance.
[2274]217
[2600]218%\verbatiminput{examples/bedslopephysical.py}
219\verbatiminput{examples/bedslope.py}
[2283]220
[2455]221\subsection{Establishing the Mesh}
[2274]222
223The first task is to set up the triangular mesh to be used for the
224scenario. This is carried out through the statement:
225
[2277]226{\small \begin{verbatim}
[2422]227    points, vertices, boundary = rectangular(10, 10)
[2277]228\end{verbatim}}
[2274]229
[2358]230The function \code{rectangular} is imported from a module
[2463]231\code{mesh\_factory} defined elsewhere. (\anuga also contains
232several other schemes that can be used for setting up meshes, but we
233shall not discuss these now.) The above assignment sets up a $10
234\times 10$ rectangular mesh, triangulated in a specific way. In
235general, the assignment
[2274]236
[2277]237{\small \begin{verbatim}
[2274]238    points, vertices, boundary = rectangular(m, n)
[2277]239\end{verbatim}}
[2274]240
241returns:
242
243\begin{itemize}
244
[2463]245   \item a list \code{points} of length $N$, where $N = (m + 1)(n + 1)$,
246comprising the coordinates $(x, y)$ of each of the $N$ mesh points,
[2274]247
[2463]248   \item a list \code{vertices} of length $2mn$ (each entry specifies the three
249vertices of one of the triangles used in the triangulation) , and
[2274]250
[2358]251   \item a dictionary \code{boundary}, used to tag the triangle edges on
[2274]252the boundaries. Each key corresponds to a triangle edge on one of
[2358]253the four boundaries and its value is one of \code{`left'},
254\code{`right'}, \code{`top'} and \code{`bottom'}, indicating
[2274]255which boundary the edge in question belongs to.
256
257\end{itemize}
258
259
[2455]260\subsection{Initialising the Domain}
[2274]261
262These variables are then used to set up a data structure
[2358]263\code{domain}, through the assignment:
[2274]264
[2277]265{\small \begin{verbatim}
[2274]266    domain = Domain(points, vertices, boundary)
[2277]267\end{verbatim}}
[2274]268
[2358]269This uses a Python class \code{Domain}, imported from
270\code{shallow\_water}, which is an extension of a more generic
271class of the same name in the module \code{domain}, and inherits
[2274]272some methods from the generic class but has others specific to the
[2422]273shallow-water scenarios in which it is used. Specific options for domain
[2434]274are set at this point. One of them is to set the basename for the output file:
[2274]275
[2463]276{\small \begin{verbatim}
[2289]277    domain.set_name('bedslope')
[2277]278\end{verbatim}}
[2274]279
280
[2455]281\subsection{Specifying the Quantities}
[2274]282
[2283]283The next task is to specify a number of quantities that we wish to set
[2358]284for each mesh point. The class \code{Domain} has a method
285\code{set\_quantity}, used to specify these quantities. It is a
[2283]286particularly flexible method that allows the user to set quantities in
287a variety of ways---using constants, functions, numeric arrays or
288expressions involving other quantities, arbitrary data points with
289associated values, all of which can be passed as arguments. All
[2358]290quantities can be initialised using \code{set\_quantity}. For
291conserved quantities (\code{stage, xmomentum, ymomentum}) this is
[2283]292called the \emph{initial condition}, for other quantities that aren't
293updated by the equation, the same interface is used to assign their
294values. The code in the present example demonstrates a number of forms
[2358]295in which we can invoke \code{set\_quantity}.
[2274]296
297
[2455]298\subsubsection{Elevation}
[2274]299
[2463]300The elevation, or height of the bed, is set using a function,
301defined through the statements below, which is specific to this
302example and specifies a particularly simple initial configuration
303for demonstration purposes:
[2274]304
[2277]305{\small \begin{verbatim}
[2274]306    def f(x,y):
307        return -x/2
[2277]308\end{verbatim}}
[2274]309
[2485]310This simply associates an elevation with each point \code{(x, y)} of
311the plane.  It specifies that the bed slopes linearly in the
312\code{x} direction, with slope $-\frac{1}{2}$,  and is constant in
313the \code{y} direction.  %[screen shot?]
314
315Once the function \code{f} is specified, the quantity
[2358]316\code{elevation} is assigned through the simple statement:
[2274]317
[2277]318{\small \begin{verbatim}
[2274]319    domain.set_quantity('elevation', f)
[2277]320\end{verbatim}}
[2274]321
322
[2455]323\subsubsection{Friction}
[2274]324
325The assignment of the friction quantity demonstrates another way
[2358]326we can use \code{set\_quantity} to set quantities---namely,
[2274]327assign them to a constant numerical value:
328
[2277]329{\small \begin{verbatim}
[2274]330    domain.set_quantity('friction', 0.1)
[2277]331\end{verbatim}}
[2274]332
333This just specifies that the Manning friction coefficient is set
334to 0.1 at every mesh point.
335
[2463]336\subsubsection{Stage}
[2274]337
[2463]338The stage (the height of the water surface) is related to the
[2637]339elevation and the depth at any time by the equation
[2274]340
[2574]341
342{\small \begin{verbatim}
343    stage = elevation + depth
[2637]344\end{verbatim}}
[2574]345
[2637]346
[2529]347For this example, we simply assign a constant value to \code{stage},
348using the statement
349
[2277]350{\small \begin{verbatim}
[2529]351    domain.set_quantity('stage', -.4)
352\end{verbatim}}
353
[2637]354which specifies that the surface level is set to a height of $-0.4$, i.e. 0.4 units
[2529]355below the zero level.
356
357(Although it is not necessary for this example, it may be useful to digress here
358and mention a variant to this requirement, which allows us to illustrate
359another way to use \code{set\_quantity}---namely, incorporating an expression
360involving other quantities. Suppose, instead of setting a constant value
361for the stage, we wished
[2637]362to specify a constant value for the \emph{depth}. For such a case we
[2529]363need to specify that \code{stage} is
364everywhere obtained by adding that value to the value already
365specified for \code{elevation}. We would do this by means of the statements:
366
367{\small \begin{verbatim}
[2274]368    h = 0.05 \# Constant depth
369    domain.set_quantity('stage', expression = 'elevation + %f' %h)
[2277]370\end{verbatim}}
[2274]371
[2637]372That is, the value of \code{stage} is set to $\code{h} = 0.05$ plus the
[2529]373value of \code{elevation} already defined.
[2274]374
[2637]375The reader will probably appreciate that this capability to incorporate
[2529]376expressions into statements using \code{set\_quantity} greatly expands
377its power.)
378
[2455]379\subsubsection{Boundary Conditions}
[2274]380
381The boundary conditions are specified as follows:
382
[2277]383{\small \begin{verbatim}
[2274]384    Br = Reflective_boundary(domain)
385
386    Bt = Transmissive_boundary(domain)
387
388    Bd = Dirichlet_boundary([0.2,0.,0.])
389
390    Bw = Time_boundary(domain=domain,
391                f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
[2277]392\end{verbatim}}
[2274]393
394The effect of these statements is to set up four alternative
395boundary conditions and store them in variables that can be
396assigned as needed. Each boundary condition specifies the
397behaviour at a boundary in terms of the behaviour in neighbouring
[2529]398elements. The boundary conditions introduced here may be briefly described as
[2274]399follows:
400
401\begin{description}
[2358]402    \item[Reflective boundary] Returns same \code{stage} as
[2283]403    as present in its neighbour volume but momentum vector reversed 180 degrees (reflected).
404    Specific to the shallow water equation as it works with the
405    momentum quantities assumed to be the second and third conserved
406    quantities.
[2450]407
[2274]408    \item[Transmissive boundary]Returns same conserved quantities as
409    those present in its neighbour volume.
[2450]410
[2274]411    \item[Dirichlet boundary]Specifies a fixed value at the
[2637]412boundary and assigns initial values to the $x$-momentum and $y$-momentum.
[2434]413
[2274]414    \item[Time boundary.]A Dirichlet boundary whose behaviour varies with time.
415\end{description}
416
417Once the four boundary types have been specified through the
418statements above, they can be applied through a statement of the
419form
420
[2277]421{\small \begin{verbatim}
422    domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
423\end{verbatim}}
[2274]424
425This statement stipulates that, in the current example, the left
426boundary is fixed, with an elevation of 0.2, while the other
427boundaries are all reflective.
428
[2529]429The reader may wish to experiment by varying the choice of boundary types
[2637]430for one or more of the boundaries. In the case of \code{Bd} and \code{Bw},
431the three arguments in each case represent the
[2274]432
[2529]433{\small \begin{verbatim}
434    Bw = Time_boundary(domain=domain,
435                f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
436\end{verbatim}}
437
438
439
[2455]440\subsection{Evolution}
[2274]441
442The final statement \nopagebreak[3]
[2277]443{\small \begin{verbatim}
[2497]444    for t in domain.evolve(yieldstep = 0.1, duration = 4.0):
445        print domain.timestepping_statistics()
[2277]446\end{verbatim}}
[2274]447
448is the key step that causes the configuration of the domain to
449`evolve' in accordance with the model embodied in the code, over a
[2358]450series of steps indicated by the values of \code{yieldstep} and
[2497]451\code{duration}, which can be altered as required.
[2450]452The value of \code{yieldstep} controls the time interval between successive model outputs.
[2434]453Behind the scenes more time steps are generally taken.
[2274]454
455
[2455]456\subsection{Output}
[2274]457
458%Give details here of the form of the output and explain how it can
459%be used with swollen. Include screen shots.//
460
[2358]461The output is a NetCDF file with the extension \code{.sww}. It
[2274]462contains stage and momentum information and can be used with the
[2358]463\code{swollen} visualisation package to generate a visual display.
[2274]464
465
[2329]466\section{How to Run the Code}
[2274]467
468The code can be run in various ways:
469
470\begin{itemize}
[2529]471  \item{from a Windows command line} as in \code{python bedslopephysical.py}
[2463]472  \item{within the Python IDLE environment}
473  \item{within emacs}
[2529]474  \item{from a Linux command line} as in \code{python bedslopephysical.py}
[2274]475\end{itemize}
476
[2499]477
478\section{Exploring the model output}
479
[2529]480Figure \ref{fig:bedslopestart} shows the \\
[2499]481Figure \ref{fig:bedslope2} shows later snapshots...
482
483
484
[2637]485\begin{figure}[hbt]
[2499]486
[2600]487%  \centerline{ \includegraphics[width=75mm, height=75mm]{examples/bedslopestart.eps}}
[2637]488
[2574]489  \caption{Bedslope example viewed with Swollen}
490  \label{fig:bedslopestart}
[2637]491\end{figure}
[2542]492
[2499]493
[2637]494\begin{figure}[hbt]
[2499]495
[2637]496  \centerline{
[2600]497%    \includegraphics[width=75mm, height=75mm]{examples/bedslopeduring.eps}
498%    \includegraphics[width=75mm, height=75mm]{examples/bedslopeend.eps}
[2637]499   }
500
[2574]501  \caption{Bedslope example viewed with Swollen}
502  \label{fig:bedslope2}
[2637]503\end{figure}
[2499]504
[2542]505
[2499]506
507
508\clearpage
509
[2463]510%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2274]511
[2434]512\section{An Example with Real Data}
[2274]513
[2450]514The following discussion builds on the concepts introduced through
[2529]515the \code{bedslopephysical.py} example and introduces a second example,
[2463]516\code{run\_sydney\_smf.py}, that follows the same basic outline, but
[2457]517incorporates more complex features and refers to a real-life
518scenario, rather than the artificial illustrative one used in
[2529]519\code{bedslopephysical.py}. The domain of interest surrounds the Sydney region,
[2637]520and predominantly covers Sydney Harbour. A hypothetical tsunami wave is
[2529]521generated by a submarine mass failure situated on the edge of the
522continental shelf.
[2274]523
[2455]524\subsection{Overview}
[2529]525As in the case of \code{bedslopephysical.py}, the actions carried out by the
[2455]526program can be organised according to this outline:
527
528\begin{enumerate}
529
530   \item Set up a triangular mesh.
531
532   \item Set certain parameters governing the mode of
[2529]533operation of the model---specifying, for instance, where to store the
[2455]534model output.
535
536   \item Input various quantities describing physical measurements, such
537as the elevation, to be specified at each mesh point (vertex).
538
539   \item Set up the boundary conditions.
540
541   \item Carry out the evolution of the model through a series of time
542steps and outputs the results, providing a results file that can be
543visualised.
544
545\end{enumerate}
546
547
548
549\subsection{The Code}
550
[2463]551Here is the code for \code{run\_sydney\_smf.py}:
[2357]552
[2475]553%\verbatiminput{"runsydneysmf.py"}
[2482]554\verbatiminput{examples/runsydney.py}
[2357]555
[2463]556In discussing the details of this example, we follow the outline
557given above, discussing each major step of the code in turn.
[2450]558
[2455]559\subsection{Establishing the Mesh}
[2450]560
[2463]561One obvious way that the present example differs from
[2529]562\code{bedslopephysical.py} is in the use of a more complex method to create
[2463]563the mesh. Instead of imposing a mesh structure on a rectangular
564grid, the technique used for this example involves building mesh
565structures inside polygons specified by the user, using a
566mesh-generator referred to as \code{pmesh}.
[2450]567
[2463]568The following remarks may help the reader understand how
569\code{pmesh} is used.
[2450]570
[2463]571In its simplest form, \code{pmesh} creates the mesh within a single
[2497]572polygon whose vertices are at geographical locations specified by the
573user. The user specifies the \emph{resolution}---that is, the maximal
574area of a triangle used for triangulation---and mesh points are
575created inside the polygon through a random process. Figure
[2637]576\ref{fig:pentagon} shows a simple example of this, in which
[2497]577the triangulation is carried out within a pentagon.
[2357]578
[2497]579
[2637]580\begin{figure}[hbt]
[2497]581
582
[2637]583
[2574]584  \caption{Mesh points are created inside the polygon}
[2637]585  \label{fig:pentagon}
586\end{figure}
[2542]587
[2457]588Boundary tags are not restricted to \code{`left'}, \code{`right'},
589\code{`bottom'} and \code{`top'}, as in the case of
[2529]590\code{bedslopephysical.py}. Instead the user specifies a list of tags
[2457]591appropriate to the configuration being modelled.
[2357]592
[2457]593While a mesh created inside a polygon offers more flexibility than
[2463]594one based on a rectangular grid, using \code{pmesh} in the limited
595form we have described so far still doesn't provide a way to adapt
596to geographic or other features in the landscape, whose presence may
[2455]597require us to vary the resolution in the neighbourhood of the
[2467]598features. To cope with this requirement, \code{pmesh} also allows
599the user to specify a number of \emph{interior polygons}, which are
600triangulated separately, each according to a separately specified
[2497]601resolution. See Figure \ref{fig:interior meshes}.
[2358]602
[2637]603\begin{figure}[hbt]
[2497]604
605
[2637]606
[2574]607  \caption{Interior meshes with individual resolution}
[2637]608  \label{fig:interior meshes}
609\end{figure}
[2542]610
[2463]611In its general form, \code{pmesh} takes for its input a bounding
612polygon and (optionally) a list of interior polygons. The user
613specifies resolutions, both for the bounding polygon and for each of
614the interior polygons. Given this data, \code{pmesh} first creates a
615triangular mesh inside the bounding polygon, using the specified
616resolution, and then creates a separate triangulation inside each of
617the interior polygons, using the resolution specified for that
618polygon.
[2449]619
[2457]620The function used to implement this process is
[2463]621\code{create\_mesh\_from\_regions}. Its arguments include the
622bounding polygon and its resolution, a list of boundary tags, and a
623list of pairs \code{[polygon, resolution]}, specifying the interior
624polygons and their resolutions.
[2450]625
[2455]626In practice, the details of the polygons used are read from a
[2457]627separate file \code{project.py}. The resulting mesh is output to a
[2463]628\emph{meshfile}\index{meshfile}. This term is used to describe a
629file of a specific format used to store the data specifying a mesh.
630(There are in fact two possible formats for such a file: it can
631either be a binary file, with extension \code{.msh}, or an ASCII
632file, with extension \code{.tsh}. In the present case, the binary
633file format \code{.msh} is used. See Section \ref{sec:file formats}
[2467]634(p. \pageref{sec:file formats}) for more on file formats.)
635\code{pmesh} assigns a name to the file by appending the extension
636\code{.msh} to the name specified in the input file
637\code{project.py}. This name is stored in the variable
[2457]638\code{meshname}.
[2450]639
[2637]640The statements
[2450]641
[2529]642{\small \begin{verbatim}
643    interior_res = 5000%
644    interior_regions = [[project.harbour_polygon_2, interior_res],
645                    [project.botanybay_polygon_2, interior_res]]
646\end{verbatim}}
647
648are used to read in the specific polygons \code{project.harbour\_polygon\_2} and
649\code{botanybay\_polygon\_2} from \code{project.py} and assign a
[2637]650common resolution of 5000 to each. The statement
[2529]651
652{\small \begin{verbatim}
653    create_mesh_from_regions(project.diffpolygonall,%
[2637]654                         boundary_tags= {'bottom': [0],%
[2529]655                                         'right1': [1],%
656                                         'right0': [2],%
657                                         'right2': [3],%
658                                         'top': [4],%
659                                         'left1': [5],%
660                                         'left2': [6],%
661                                         'left3': [7]},%
662                         maximum_triangle_area=100000,%
[2637]663                         filename=meshname,%
[2529]664                         interior_regions=interior_regions)
665\end{verbatim}}
666
667is then used to create the mesh, taking the bounding polygon to be the polygon
[2637]668\code{diffpolygonall} specified in \code{project.py}. The
669argument \code{boundary\_tags} assigns a dictionary, whose keys are the
670names of the boundary tags used for the bounding polygon---\code{`bottom'},
671`right0', `right1', `right2', `top', `left1', `left2' and `left3'---
[2529]672and whose values identify the indices of the segments associated with each of these
673tags. (The value associated with each boundary tag is a one-element list.)
674
675
[2455]676\subsection{Initialising the Domain}
[2450]677
[2529]678As with \code{bedslopephysical.py}, once we have created the mesh, the next
[2463]679step is to create the data structure \code{domain}. We did this for
[2529]680\code{bedslopephysical.py} by inputting lists of points and triangles and
[2463]681specifying the boundary tags directly. However, in the present case,
682we use a method that works directly with the meshfile
683\code{meshname}, as follows:
[2450]684
[2463]685{\small \begin{verbatim}
[2529]686    domain = pmesh_to_domain_instance(meshname,
687    Domain, use_cache = True, verbose = True)
[2455]688\end{verbatim}}
689
[2463]690The function \code{pmesh\_to\_domain\_instance} converts a meshfile
[2457]691\code{meshname} into an instance of the data structure
692\code{domain}, allowing us to use methods like \code{set\_quantity}
[2463]693to set quantities and to apply other operations. (In principle, the
694second argument of \code{pmesh\_to\_domain\_instance} can be any
695subclass of \code{Domain}, but for applications involving the
696shallow-water wave equation, the second argument of
697\code{pmesh\_to\_domain\_instance} can always be set simply to
698\code{Domain}.)
[2457]699
[2529]700The following statements specify a basename and data directory, and
701identify quantities to be stored. For the first two, values are
702taken from \code{project.py}.
703
[2637]704{\small \begin{verbatim}
[2529]705    domain.set_name(project.basename)%
706    domain.set_datadir(project.outputdir)%
707    domain.set_quantities_to_be_stored(['stage', 'xmomentum',
708        'ymomentum'])
709\end{verbatim}}
710
711
[2455]712\subsection{Specifying the Quantities}
[2529]713Quantities for \code{run\_sydney\_smf.py} are set
714using similar methods to those in \code{bedslopephysical.py}. However,
[2463]715in this case, many of the values are read from the auxiliary file
716\code{project.py} or, in the case of \code{elevation}, from an
717ancillary points file.
[2455]718
[2457]719
720
[2463]721\subsubsection{Stage}
[2455]722
[2463]723For the scenario we are modelling in this case, we use a callable
724object \code{tsunami\_source}, assigned by means of a function
725\code{slump\_tsunami}. This is similar to how we set elevation in
[2529]726\code{bedslopephysical.py} using a function---however, in this case the
[2637]727function is both more complex and more interesting.
[2463]728
[2637]729The function returns the water displacement for all \code{x}
730and \code{y} in the domain. The water displacement is a ?? function that depends
[2529]731on the characteristics of the slump (length, thickness, slope, etc), its
[2637]732location (origin) and the depth at that location.
[2455]733
734
[2463]735\subsubsection{Friction}
736
[2529]737We assign the friction exactly as we did for \code{bedslopephysical.py}:
[2463]738
739{\small \begin{verbatim}
[2529]740    domain.set_quantity('friction', 0.03)
[2463]741\end{verbatim}}
742
743
744\subsubsection{Elevation}
745
[2529]746The elevation is specified by reading data from a file:
[2463]747
[2637]748{\small \begin{verbatim}
749    domain.set_quantity('elevation',
750                    filename = project.combineddemname + '.pts',
751                    use_cache = True,
752                    verbose = True)
[2455]753\end{verbatim}}
754
[2463]755However, before this step can be executed, some preliminary steps
756are needed to prepare the file from which the data is taken. Two
757source files are used for this data---their names are specified in
758the file \code{project.py}, in the variables \code{coarsedemname}
759and \code{finedemname}. They contain `coarse' and `fine' data,
760respectively---that is, data sampled at widely spaced points over a
761large region and data sampled at closely spaced points over a
762smaller subregion. The data in these files is combined through the
763statement
764
765{\small \begin{verbatim}
766combine_rectangular_points_files(project.finedemname + '.pts',
767                                 project.coarsedemname + '.pts',
768                                 project.combineddemname + '.pts')
769\end{verbatim}}
770
771The effect of this is simply to combine the datasets by eliminating
772any coarse data associated with points inside the smaller region
773common to both datasets. The name to be assigned to the resulting
774dataset is also derived from the name stored in the variable
775\code{combinedname} in the file \code{project.py}.
776
[2455]777\subsection{Boundary Conditions}
778
[2637]779Setting boundaries follows a similar pattern to the one used for
780\code{bedslopephysical.py}, except that in this case we need to associate a
[2529]781boundary type with each of the
[2637]782boundary tag names introduced when we established the mesh. In place of the four
783boundary types introduced for \code{bedslopephysical.py}, we use the reflective
[2529]784boundary for each of the
785eight tagged segments:
[2463]786
787{\small \begin{verbatim}
[2529]788    Br = Reflective_boundary(domain)
789    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
790                          'right2': Br, 'top': Br, 'left1': Br,
791                          'left2': Br, 'left3': Br} )
[2455]792\end{verbatim}}
793
[2463]794\subsection{Evolution}
795
796With the basics established, the running of the `evolve' step is
[2529]797very similar to the corresponding step in \code{bedslopephysical.py}:
[2463]798
799{\small \begin{verbatim}
[2529]800    import time t0 = time.time()
[2637]801
[2529]802    for t in domain.evolve(yieldstep = 120, duration = 18000):
803        print domain.timestepping_statistics()
804        print domain.boundary_statistics(tags = 'bottom')
[2637]805
[2529]806    print 'That took %.2f seconds' %(time.time()
[2463]807\end{verbatim}}
808
[2467]809%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
810%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
811
[2434]812\chapter{\anuga Public Interface}
[2358]813
[2474]814This chapter gives an overview of the features of \anuga available
815to the user at the public interface. These are grouped under the
816following headings:
[2358]817
[2474]818\begin{itemize}
[2485]819    \item Establishing the Mesh
820    \item Initialising the Domain
821    \item Specifying the Quantities
822    \item Initial Conditions
823    \item Boundary Conditions
824    \item Forcing Functions
[2497]825    \item Evolution
[2474]826\end{itemize}
827
[2485]828The listings are intended merely to give the reader an idea of what
829each feature is, where to find it and how it can be used---they do
[2559]830not give full specifications. For these the reader
[2485]831may consult the code. The code for every function or class contains
832a documentation string, or `docstring', that specifies the precise
833syntax for its use. This appears immediately after the line
834introducing the code, between two sets of triple quotes.
[2474]835
[2529]836Each listing also describes the location of the module in which
[2485]837the code for the feature being described can be found. All modules
[2529]838are in the folder \code{inundation} or one of its subfolders, and the
839location of each module is described relative to \code{inundation}. Rather
[2485]840than using pathnames, whose syntax depends on the operating system,
841we use the format adopted for importing the function or class for
842use in Python code. For example, suppose we wish to specify that the
843function \code{create\_mesh\_from\_regions} is in a module called
844\code{mesh\_interface} in a subfolder of \code{inundation} called
845\code{pmesh}. In Linux or Unix syntax, the pathname of the file
846containing the function, relative to \code{inundation}, would be
[2358]847
[2485]848\begin{center}
[2529]849    \code{pmesh/mesh\_interface.py}
[2485]850\end{center}
[2358]851
[2485]852while in Windows syntax it would be
[2450]853
[2485]854\begin{center}
[2529]855    \code{pmesh}$\backslash$\code{mesh\_interface.py}
[2485]856\end{center}
[2450]857
[2485]858Rather than using either of these forms, in this chapter we specify
859the location simply as \code{pmesh.mesh_interface}, in keeping with
860the usage in the Python statement for importing the function,
861namely:
862\begin{center}
[2529]863    \code{from pmesh.mesh\_interface import create\_mesh\_from\_regions}
[2485]864\end{center}
[2450]865
866
[2485]867\section{Mesh Generation}
[2439]868
[2467]869\begin{funcdesc}  {create\_mesh\_from\_regions}{bounding_polygon,
870                             boundary_tags,
871                             maximum_triangle_area,
872                             filename=None,
873                             interior_regions=None,
874                             poly_geo_reference=None,
875                             mesh_geo_reference=None,
876                             minimum_triangle_angle=28.0}
[2485]877Module: \code{pmesh.mesh\_interface}
878
[2529]879% Translate following into layman's language
[2637]880This function is used to create a triangular mesh suitable for use with
881\anuga, within a specified region. The region is specified as the interior of a polygon
882(the \emph{bounding polygon}). The user specifies the bounding polygon and the
883\emph{resolution}---that is, maximal area of any triangle in the mesh. There is
[2559]884also an option to specify a number of internal polygons within each of which a
885separate mesh is created, generally with a smaller resolution. Additionally,
[2637]886the user specifies a list of boundary tags, one for each edge of the bounding
887polygon.
[2467]888\end{funcdesc}
889
[2485]890%%%%%%
891\section{Initialising Domain}
[2467]892
[2485]893\begin{funcdesc}  {pmesh\_to\_domain\_instance}{file_name, DomainClass, use_cache = False, verbose = False}
894Module: \code{pyvolution.pmesh2domain}
[2467]895
[2637]896Once the initial mesh file has been created, this function is applied
[2559]897to convert it to a domain object---that is, to a member of
898the special Python class Domain (or a subclass of Domain), which provides access to properties and
899methods that allow quantities to be set and other operations to be carried out.
[2467]900
[2485]901\code{file\_name} is the name of the mesh file to be converted,
902including the extension. \code{DomainClass} is the class to be
[2529]903returned, which must be a subclass of \code{Domain} having the same
[2485]904interface as \code{Domain}---in practice, it can usually be set
905simply to \code{Domain}.
[2467]906\end{funcdesc}
907
[2497]908
909\subsection{Setters and getters of class Domain}
910
[2628]911\begin{funcdesc} {set\_name}{name}
[2631]912    Module: \code{pyvolution.domain}
[2637]913
[2628]914    Assigns the name \code{name} to the domain
[2497]915\end{funcdesc}
916
[2628]917\begin{funcdesc} {get\_name}{}
918    Module: \code{pyvolution.domain}
[2637]919
[2628]920    Returns the name assigned to the domain by \code{set_name}. If no name has been
921    assigned, returns `domain'.
[2497]922\end{funcdesc}
923
[2628]924\begin{funcdesc} {set\_datadir}{name}
925    Module: \code{pyvolution.domain}
[2637]926
[2628]927    Sets the directory used for data to the value \code{name}. The default value, before
[2637]928    \code{set\_datadir} is run, is the value \code{default_datadir}
[2631]929    specified in \code{config.py}.
[2497]930\end{funcdesc}
931
[2628]932\begin{funcdesc} {get_datadir}{}
933    Module: \code{pyvolution.domain}
[2637]934
[2628]935    Returns the data directory set by \code{set\_datadir} or, if \code{set\_datadir} has not
[2637]936    been run, returns the value \code{default_datadir} specified in
[2631]937    \code{config.py}.
[2497]938\end{funcdesc}
939
940\begin{funcdesc} {set_time}{??}
941\end{funcdesc}
942
943\begin{funcdesc} {set_default_order}{??}
944\end{funcdesc}
945
946
[2485]947%%%%%%
948\section{Setting Quantities}
[2467]949
[2637]950\begin{funcdesc}{set\_quantity}{name,
951    numeric = None,
952    quantity = None,
953    function = None,
954    geospatial_data = None,
[2631]955    filename = None,
[2637]956    attribute_name = None,
957    alpha = None,
958    location = 'vertices',
959    indices = None,
[2631]960    verbose = False,
[2637]961    use_cache = False}
[2631]962  Module: \code{pyvolution.domain}
963  (see also \code{pyvolution.quantity.set_values})
[2467]964
[2631]965This function is used to assign values to individual quantities for a
966domain. It is very flexible and can be used with many data types: a
967statement of the form \code{domain.set\_quantity(name, x)} can be used
968to define a quantity having the name \code{name}, where the other
969argument \code{x} can be any of the following:
[2422]970
[2559]971\begin{itemize}
[2631]972\item a number, in which case all vertices in the mesh gets that for
973the quantity in question.
[2574]974\item a list of numbers or a Numeric array ordered the same way as the mesh vertices.
[2637]975\item a function (e.g.\ see the samples introduced in Chapter 2)
[2559]976\item an expression composed of other quantities and numbers, arrays, lists (for
977example, a linear combination of quantities)
[2637]978\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.
[2631]979\item a geospatial dataset (See ?????). Optional argument attribute_name applies here as with files.
[2559]980\end{itemize}
[2450]981
982
[2485]983Exactly one of the arguments
984  numeric, quantity, function, points, filename
985must be present.
[2631]986
987
[2637]988Set quantity will look at the type of the second argument (\code{numeric}) and
989determine what action to take.
[2631]990
991Values can also be set using the appropriate keyword arguments.
[2637]992If 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)}
993are all equivalent.
[2631]994
995
996Other optional arguments are
997\begin{itemize}
998\item \code{indices} which is a list of ids of triangles to which set_quantity should apply its assignment of values.
999\item \code{location} determines which part of the triangles to assign to. Options are 'vertices' (default), 'edges', and 'centroids'.
[2637]1000\end{itemize}
[2631]1001
1002
1003a number, in which case all vertices in the mesh gets that for
1004the quantity in question.
1005\item a list of numbers or a Numeric array ordered the same way as the mesh vertices.
1006
1007
[2485]1008\end{funcdesc}
[2434]1009
[2631]1010
1011
1012
1013
1014
1015
[2485]1016%%%
1017\anuga provides a number of predefined initial conditions to be used
1018with \code{set_quantity}.
[2434]1019
[2637]1020\begin{funcdesc}{tsunami_slump}{length, depth, slope, width=None, thickness=None,
1021                x0=0.0, y0=0.0, alpha=0.0,
1022                gravity=9.8, gamma=1.85,
1023                massco=1, dragco=1, frictionco=0, psi=0,
1024                dx=None, kappa=3.0, kappad=0.8, zsmall=0.01,
[2485]1025                domain=None,
1026                verbose=False}
1027This function returns a callable object representing an initial water
1028displacement generated by a submarine sediment slide.
[2450]1029
[2485]1030The arguments include the downslope slide length, the water depth to the slide centre of mass,
1031and the bathymetric slope.
1032\end{funcdesc}
[2450]1033
1034
[2485]1035%%%
1036\begin{funcdesc}{file_function}{filename,
1037                  domain = None,
1038                  quantities = None,
1039                  interpolation_points = None,
1040                  verbose = False,
1041                  use_cache = False}
1042Module: \code{pyvolution.util}
[2450]1043
[2485]1044Reads the time history of spatial data from NetCDF file and returns
1045a callable object. Returns interpolated values based on the input
1046file using the underlying \code{interpolation\_function}.
[2450]1047
[2485]1048\code{quantities} is either the name of a single quantity to be
[2559]1049interpolated or a list of such quantity names. In the second case, the resulting
[2485]1050function will return a tuple of values---one for each quantity.
[2450]1051
[2485]1052\code{interpolation\_points} is a list of absolute UTM coordinates
1053for points at which values are sought.
1054\end{funcdesc}
[2450]1055
[2485]1056%%%
[2626]1057\begin{classdesc}{Interpolation\_function}{self,
[2474]1058                 time,
1059                 quantities,
1060                 quantity_names = None,
1061                 vertex_coordinates = None,
1062                 triangles = None,
1063                 interpolation_points = None,
1064                 verbose = False}
[2485]1065Module: \code{pyvolution.least\_squares}
1066
[2637]1067Given a time series, either as a sequence of numbers or
[2631]1068defined at the vertices of a triangular mesh (such
1069as those stored in \code{sww} files), \code{Interpolation\_function}
1070is used to create a callable object that interpolates a value for
[2637]1071an arbitrary time \code{t} within the model limits and possibly a
1072point \code{(x, y)} within a mesh region.
[2631]1073
[2637]1074The actual time series at which data is available is specified by means
[2631]1075of an array \code{time} of monotonically increasing times.
[2637]1076The quantities containing the values to be interpolated are specified in
1077an array---or dictionary of arrays (used in conjunction with the optional argument
[2631]1078\code{quantitity\_names}) --- called \code{quantities}.
[2637]1079The optional arguments \code{vertex_coordinates} and \code{triangles} represent
[2631]1080the spatial mesh associated with the quantity arrays. If omitted the function created
1081by \code{Interpolation\_function} will be a function of \code{t} only.
1082
1083
[2637]1084Since, in practice, values need to be computed at specified
[2626]1085points, the syntax allows the user to specify, once and for all, a list
1086\code{interpolation\_points} of points at which values are required. In this case,
1087the function may be called using the form \code{f(t, id)}, where \code{id} is an
1088index identifying a member of \code{interpolation\_points}.
[2474]1089
[2631]1090
1091
1092
1093
[2474]1094\end{classdesc}
[2434]1095
[2485]1096%%%
1097\begin{funcdesc}{set\_region}{functions}
[2474]1098[Low priority. Will be merged into set\_quantity]
[2450]1099
[2485]1100Module:\code{pyvolution.domain}
[2474]1101\end{funcdesc}
[2450]1102
[2358]1103
1104
[2485]1105%%%%%%
[2434]1106\section{Boundary Conditions}
[2358]1107
[2631]1108\anuga provides a large number of predefined boundary conditions,
1109represented by objects such as \code{Reflective\_boundary(domain)} and
1110\code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples
1111in Chapter 2. Alternatively, you may prefer to ``roll your own'',
1112following the method explained in Section \ref{sec:roll your own}.
[2358]1113
[2627]1114These boundary objects may be used with the function \code{set\_boundary} described below
1115to assign boundary conditions according to the tags used to label boundary segments.
1116
[2485]1117\begin{funcdesc}{set\_boundary}{boundary_map}
1118Module: \code{pyvolution.domain}
[2486]1119
[2627]1120This function allows you to assign a boundary object (corresponding to a
1121pre-defined or user-specified boundary condition) to every boundary segment that
[2637]1122has been assigned a particular tag.
[2486]1123
[2627]1124This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects
1125and whose keys are the symbolic tags.
[2486]1126
[2485]1127\end{funcdesc}
[2358]1128
[2631]1129\begin{funcdesc} {get_boundary_tags}{}
1130Module: \code{pyvolution.mesh}
[2497]1131\end{funcdesc}
1132
[2485]1133%%%
1134\subsection{Predefined boundary conditions}
[2422]1135
[2485]1136\begin{classdesc}{Reflective_boundary}{Boundary}
1137Module: \code{pyvolution.shallow\_water}
[2422]1138
[2485]1139Reflective boundary returns same conserved quantities as those present in
1140the neighbouring volume but reflected.
[2422]1141
[2485]1142This class is specific to the shallow water equation as it works with the
1143momentum quantities assumed to be the second and third conserved quantities.
1144\end{classdesc}
[2422]1145
[2485]1146%%%
1147\begin{classdesc}{Transmissive_boundary}{domain = None}
1148Module: \code{pyvolution.generic\_boundary\_conditions}
[2422]1149
[2485]1150A transmissive boundary returns the same conserved quantities as
1151those present in the neighbouring volume.
[2422]1152
[2485]1153The underlying domain must be specified when the boundary is instantiated.
1154\end{classdesc}
[2358]1155
[2485]1156%%%
1157\begin{classdesc}{Dirichlet_boundary}{conserved_quantities=None}
1158Module: \code{pyvolution.generic\_boundary\_conditions}
[2422]1159
[2485]1160A Dirichlet boundary returns constant values for the conserved
1161quantities.
1162\end{classdesc}
[2422]1163
[2485]1164%%%
1165\begin{classdesc}{Time_boundary}{domain = None, f = None}
1166Module: \code{pyvolution.generic\_boundary\_conditions}
[2358]1167
[2485]1168A time-dependent boundary returns values for the conserved
1169quantities as a function \code{f(t)} of time. The user must specify
1170the domain to get access to the model time.
1171\end{classdesc}
[2358]1172
[2485]1173%%%
1174\begin{classdesc}{File_boundary}{Boundary}
1175Module: \code{pyvolution.generic\_boundary\_conditions}
[2358]1176
[2485]1177The boundary values are obtained from a file and interpolated. The
1178file is assumed to contain a time series and possibly also spatial
1179information. The conserved quantities are given as a function of
1180time.
1181\end{classdesc}
[2358]1182
[2422]1183
[2485]1184\subsection{User-defined boundary conditions}
[2631]1185\label{sec:roll your own}
[2485]1186[How to roll your own]
[2358]1187
1188
1189
[2474]1190
[2485]1191
[2434]1192\section{Forcing Functions}
[2358]1193
[2434]1194\anuga provides a number of predefined forcing functions to be used with .....
[2358]1195
[2474]1196%\begin{itemize}
[2358]1197
1198
[2474]1199%  \item \indexedcode{}
1200%  [function, arguments]
[2422]1201
[2474]1202%  \item \indexedcode{}
[2358]1203
[2474]1204%\end{itemize}
[2358]1205
1206
[2485]1207
[2497]1208\section{Evolution}
[2485]1209
[2600]1210  \begin{funcdesc}{evolve}{yieldstep = None, finaltime = None, duration = None, skip_initial_step = False}
[2637]1211
[2600]1212  Module: pyvolution.domain
[2637]1213
1214  This function (a method of \code{domain}) is invoked once all the preliminaries have been
1215  completed, and causes the model to progress through successive steps in its evolution, storing results and
[2626]1216  outputting statistics whenever a user-specified period \code{yieldstep} is completed (generally
1217  during this period the model will evolve through several steps internally). The user specifies
1218  the total time period over which the evolution is to take place, by specifying values (in seconds) for either
[2637]1219  \code{duration} or \code{finaltime}, as well as the interval in seconds after which results are to be
1220  stored and statistics output.
1221
[2600]1222  You can include \code{evolve} in a statement of the type:
[2637]1223
[2600]1224  {\small \begin{verbatim}
1225      for t in domain.evolve(yieldstep, finaltime):
1226          <Do something with domain and t>
1227  \end{verbatim}}
[2485]1228
1229  \end{funcdesc}
1230
1231
1232
[2497]1233\subsection{Diagnostics}
1234
1235  \begin{funcdesc}{timestepping_statistics}{???}
1236
[2485]1237  \end{funcdesc}
1238
[2497]1239
1240  \begin{funcdesc}{boundary\_statistics}{???}
1241
1242  \end{funcdesc}
[2637]1243
1244
[2497]1245  \begin{funcdesc}{get_quantity}{???}
[2637]1246  Module: \code{pyvolution.domain}
[2497]1247  Allow access to individual quantities and their methods
[2637]1248
1249  \end{funcdesc}
1250
1251
[2497]1252  \begin{funcdesc}{get_values}{???}
[2637]1253  Module: \code{pyvolution.quantity}
1254
[2497]1255  Extract values for quantity as an array
[2637]1256
1257  \end{funcdesc}
1258
1259
[2497]1260  \begin{funcdesc}{get_integral}{???}
[2637]1261  Module: \code{pyvolution.quantity}
1262
[2497]1263  Return computed integral over entire domain for this quantity
[2524]1264
[2637]1265  \end{funcdesc}
1266
1267
1268\section{Other}
1269
1270  \begin{funcdesc}{domain.create_quantity_from_expression}{???}
1271
1272  Handy for creating derived quantities on-the-fly.
[2529]1273  See \code{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use.
[2637]1274  \end{funcdesc}
[2497]1275
[2474]1276%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1277%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2358]1278
[2434]1279\chapter{\anuga System Architecture}
[2370]1280
[2413]1281From pyvolution/documentation
[2370]1282
[2434]1283\section{File Formats}
[2439]1284\label{sec:file formats}
[2370]1285
[2467]1286\anuga makes use of a number of different file formats. The
1287following table lists all these formats, which are described in more
1288detail in the paragraphs below.
[2439]1289
[2475]1290\bigskip
1291
[2485]1292\begin{center}
[2475]1293
[2467]1294\begin{tabular}{|ll|}  \hline
[2475]1295
1296\textbf{Extension} & \textbf{Description} \\
1297\hline\hline
1298
[2467]1299\code{.sww} & NetCDF format for storing model output
1300\code{f(t,x,y)}\\
[2439]1301
[2467]1302\code{.tms} & NetCDF format for storing time series \code{f(t)}\\
[2439]1303
[2467]1304\code{.xya} & ASCII format for storing arbitrary points and
1305associated attributes\\
[2463]1306
[2467]1307\code{.pts} & NetCDF format for storing arbitrary points and
1308associated attributes\\
[2463]1309
[2467]1310\code{.asc} & ASCII format of regular DEMs as output from ArcView\\
[2463]1311
[2467]1312\code{.prj} & Associated ArcView file giving more metadata for
1313\code{.asc} format\\
[2463]1314
[2467]1315\code{.ers} & ERMapper header format of regular DEMs for ArcView\\
[2463]1316
[2467]1317\code{.dem} & NetCDF representation of regular DEM data\\
[2463]1318
[2467]1319\code{.tsh} & ASCII format for storing meshes and associated
1320boundary and region info\\
[2463]1321
[2467]1322\code{.msh} & NetCDF format for storing meshes and associated
1323boundary and region info\\
[2463]1324
[2467]1325\code{.nc} & Native ferret NetCDF format\\
[2463]1326
[2467]1327\code{.geo} & Houdinis ASCII geometry format (?) \\  \par \hline
1328%\caption{File formats used by \anuga}
1329\end{tabular}
[2463]1330
1331
[2485]1332\end{center}
[2475]1333
1334\bigskip
1335
[2467]1336A typical dataflow can be described as follows:
[2463]1337
[2467]1338\subsection{Manually Created Files}
[2463]1339
[2467]1340\begin{tabular}{ll}
1341ASC, PRJ & Digital elevation models (gridded)\\
1342TSH & Triangular
1343meshes (e.g. created from \code{pmesh})\\
1344NC & Model outputs for use as boundary conditions (e.g. from MOST)
1345\end{tabular}
[2463]1346
[2467]1347\subsection{Automatically Created Files}
[2463]1348
[2467]1349\begin{tabular}{ll}
1350ASC, PRJ  $\rightarrow$  DEM  $\rightarrow$  PTS & Convert
1351DEMs to native \code{.pts} file\\
[2463]1352
[2467]1353NC $\rightarrow$ SWW & Convert MOST boundary files to
1354boundary \code{.sww}\\
[2463]1355
[2467]1356PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\
1357
1358TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using
1359Swollen\\
1360
1361TSH + Boundary SWW $\rightarrow$ SWW & Simulation using
1362\code{pyvolution}
1363\end{tabular}
1364
[2497]1365
1366\subsection{Basic file conversions}
1367
1368  \begin{funcdesc}{sww2dem}{???}
[2637]1369  Module: \code{pyvolution.data\_manager}
[2497]1370
[2637]1371
1372  \end{funcdesc}
1373
1374
[2497]1375  \begin{funcdesc}{dem2pts}{???}
[2637]1376  Module: \code{pyvolution.data_manager}
[2497]1377
[2637]1378
1379  \end{funcdesc}
1380
[2463]1381%\[
1382%  \left[
1383%  \begin{array}{ccr}
1384%  2 & 4 & 4\\
1385%  1 & 1 & 1
1386%  \end{array}
1387%  \right]
1388%\]
1389
1390%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1391%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1392
[2434]1393\chapter{Basic \anuga Assumptions}
[2370]1394
[2413]1395(From pyvolution/documentation)
[2370]1396
[2413]1397
1398Physical model time cannot be earlier than 1 Jan 1970 00:00:00.
1399If one wished to recreate scenarios prior to that date it must be done
1400using some relative time (e.g. 0).
1401
1402
1403All spatial data relates to the WGS84 datum (or GDA94) and has been
1404projected into UTM with false easting of 500000 and false northing of
14051000000 on the southern hemisphere (0 on the northern).
1406
1407It is assumed that all computations take place within one UTM zone.
1408
1409DEMs, meshes and boundary conditions can have different origins within
[2422]1410one UTM zone. However, the computation will use that of the mesh for
[2413]1411numerical stability.
1412
1413
1414%OLD
1415%The dataflow is: (See data_manager.py and from scenarios)
1416%
1417%
1418%Simulation scenarios
1419%--------------------%
1420%%
1421%
1422%Sub directories contain scrips and derived files for each simulation.
1423%The directory ../source_data contains large source files such as
1424%DEMs provided externally as well as MOST tsunami simulations to be used
1425%as boundary conditions.
1426%
1427%Manual steps are:
1428%  Creation of DEMs from argcview (.asc + .prj)
1429%  Creation of mesh from pmesh (.tsh)
1430%  Creation of tsunami simulations from MOST (.nc)
1431%%
1432%
1433%Typical scripted steps are%
1434%
1435%  prepare_dem.py:  Convert asc and prj files supplied by arcview to
1436%                   native dem and pts formats%
1437%
1438%  prepare_pts.py: Convert netcdf output from MOST to an sww file suitable
1439%                  as boundary condition%
1440%
1441%  prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares
1442%                   smoothing. The outputs are tsh files with elevation data.%
1443%
1444%  run_simulation.py: Use the above together with various parameters to
[2434]1445%                     run inundation simulation.
[2413]1446
1447
[2463]1448%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1449%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2413]1450
[2329]1451\appendix
[2358]1452
[2434]1453\chapter{Supporting Tools}
[2358]1454
[2559]1455This section describes a number of supporting tools, supplied with \anuga, that offer a
1456variety of types of functionality and enhance the basic capabilities of \anuga.
[2358]1457
[2434]1458\section{caching}
[2358]1459
[2450]1460  The \code{cache} function is used to provide supervised caching of function results. A Python
[2434]1461  function call of the form
[2358]1462
[2463]1463      {\small \begin{verbatim}
[2434]1464      result = func(arg1,...,argn)
1465      \end{verbatim}}
[2358]1466
[2434]1467  can be replaced by
[2358]1468
[2463]1469      {\small \begin{verbatim}
[2434]1470      from caching import cache
1471      result = cache(func,(arg1,...,argn))
1472      \end{verbatim}}
[2450]1473
[2434]1474  which returns the same output but reuses cached
1475  results if the function has been computed previously in the same context.
1476  \code{result} and the arguments can be simple types, tuples, list, dictionaries or
[2450]1477  objects, but not unhashable types such as functions or open file objects.
[2434]1478  The function \code{func} may be a member function of an object or a module.
[2358]1479
[2434]1480  This type of caching is particularly useful for computationally intensive
1481  functions with few frequently used combinations of input arguments. Note that
[2559]1482  if the inputs or output are very large caching may not save time because
[2434]1483  disc access may dominate the execution time.
1484
[2559]1485  If the function definition changes after a result has been cached, this will be
[2434]1486  detected by examining the functions \code{bytecode (co_code, co_consts,
[2559]1487  func_defualts, co_argcount)} and the function will be recomputed.
[2450]1488
[2637]1489  Options are set by means of the function \code{set_option(key, value)},
[2559]1490  where \code{key} is a key associated with a
[2637]1491  Python dictionary \code{options}. This dictionary stores settings such as the name of
[2559]1492  the directory used, the maximum
[2434]1493  number of cached files allowed, and so on.
[2450]1494
[2434]1495  The \code{cache} function allows the user also to specify a list of dependent files. If any of these
1496  have been changed, the function is recomputed and the results stored again.
[2450]1497
[2559]1498  %Other features include support for compression and a capability to \ldots
[2434]1499
[2450]1500
[2434]1501   \textbf{USAGE:}
[2450]1502
[2463]1503    {\small \begin{verbatim}
[2434]1504    result = cache(func, args, kwargs, dependencies, cachedir, verbose,
[2559]1505                   compression, evaluate, test, return_filename)
[2434]1506    \end{verbatim}}
1507
[2637]1508
[2434]1509\section{swollen}
[2637]1510The output generated by \anuga may be viewed by means of the visualisation tool \code{swollen},
1511which takes the \code{sww} file output by \anuga and creates a visual representation of the data.
1512Examples may be seen in Figures \ref{fig:bedslopestart} and \ref{fig:bedslope2}.
1513To view an \code{sww} file with \code{swollen} in the
1514Windows environment, you can simply drag the icon representing the file over an icon on the desktop
[2559]1515for the \code{swollen} executable file (or a shortcut to it). Alternatively, you can operate \code{swollen}
1516from the command line, in both Windows and Linux environments.
[2434]1517
[2559]1518On successful operation, you will see an interactive moving-picture display. You can use keys and the mouse
[2637]1519to slow down, speed up or stop the display, change the viewing position or carry out a number of other
[2559]1520simple operations.
[2637]1521
[2434]1522The main keys operating the interactive screen are:\\
1523
1524\begin{tabular}{|ll|}   \hline
1525
1526\code{w} & toggle wireframe\\
1527
1528space bar & start/stop\\
1529
1530up/down arrows & increase/decrease speed\\
1531
1532left/right arrows & direction in time \emph{(when running)}\\ & step through simulation \emph{(when stopped)}\\
1533
1534left mouse button & rotate\\
1535
1536middle mouse button & pan\\
1537
1538right mouse button & zoom\\  \hline
1539
1540\end{tabular}
1541
1542\vfill
1543
1544The following table describes how to operate swollen from the command line:
1545
1546Usage: \code{swollen [options] swwfile \ldots}\\  \nopagebreak
1547Options:\\  \nopagebreak
1548\begin{tabular}{ll}
1549  \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY_CENTER |}\\
1550                                    & \code{HEAD_MOUNTED_DISPLAY}\\
1551  \code{--rgba} & Request a RGBA colour buffer visual\\
1552  \code{--stencil} & Request a stencil buffer visual\\
1553  \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\
1554                                    & overridden by environmental variable\\
1555  \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD_BUFFER | HORIZONTAL_SPLIT |}\\
1556                                    & \code{VERTICAL_SPLIT | LEFT_EYE | RIGHT_EYE |}\\
1557                                     & \code{ON | OFF} \\
1558  \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\
1559  \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\
1560  \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\
1561  \code{-help} & Display this information\\
1562  \code{-hmax <float>} & Height above which transparency is set to
1563                                     \code{alphamax}\\
1564  \code{-hmin <float>} & Height below which transparency is set to
1565                                     zero\\
1566  \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is
1567                                     up, default is overhead)\\
1568  \code{-loop}  & Repeated (looped) playback of \code{.swm} files\\
1569  \code{-movie <dirname>} & Save numbered images to named directory and
1570                                     quit\\
1571  \code{-nosky} & Omit background sky\\
1572  \code{-scale <float>} & Vertical scale factor\\
1573  \code{-texture <file>} & Image to use for bedslope topography\\
1574  \code{-tps <rate>} & Timesteps per second\\
1575  \code{-version} & Revision number and creation (not compile)
1576                                     date\\
1577\end{tabular}
1578
[2574]1579\section{utilities/polygons}
[2358]1580
[2560]1581  \begin{classdesc}{Polygon_function}{regions, default = 0.0, geo_reference = None}
[2637]1582  Module: \code{utilities.polygon}
1583
1584
[2560]1585  \end{classdesc}
[2637]1586
[2560]1587  \begin{funcdesc}{read_polygon}{filename}
[2637]1588  Module: \code{utilities.polygon}
1589
[2560]1590  Reads the specified file and returns a polygon. Each
1591  line of the file must contain exactly two numbers, separated by a comma, which are interpreted
1592  as coordinates of one vertex of the polygon.
1593  \end{funcdesc}
[2637]1594
[2560]1595  \begin{funcdesc}{populate_polygon}{polygon, number_of_points, seed = None, exclude = None}
[2637]1596  Module: \code{utilities.polygon}
1597
[2560]1598  Populates the interior of the specified polygon with the specified number of points,
1599  selected by means of a uniform distribution function.
1600  \end{funcdesc}
[2637]1601
[2560]1602  \begin{funcdesc}{point_in_polygon}{polygon, delta=1e-8}
[2637]1603  Module: \code{utilities.polygon}
1604
[2560]1605  Returns a point inside the specified polygon and close to the edge. The distance between
1606  the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the
1607  second argument \code{delta}, which is taken as $10^{-8}$ by default.
1608  \end{funcdesc}
[2637]1609
[2560]1610  \begin{funcdesc}{inside_polygon}{points, polygon, closed = True, verbose = False}
[2637]1611  Module: \code{utilities.polygon}
1612
[2560]1613  Used to test whether a single point---or the members of a list of points---
1614  are inside the specified polygon. If the first argument is a single point,
1615  returns \code{True} if the point is inside the polygon, or \code{False}
1616  otherwise. If the first argument is a list of points, returns a Numeric
1617  array comprising the indices of the points in the list that lie inside the polygon.
1618  (If none of the points are inside, returns \code{zeros((0,), 'l')}.)
[2637]1619  Points on the edges of the polygon are regarded as inside if
[2560]1620  \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside.
1621  \end{funcdesc}
[2637]1622
[2560]1623  \begin{funcdesc}{outside_polygon}{points, polygon, closed = True, verbose = False}
[2637]1624  Module: \code{utilities.polygon}
1625
[2560]1626  Exactly like \code{inside_polygon}, but with the words `inside' and `outside' interchanged.
1627  \end{funcdesc}
[2637]1628
[2560]1629  \begin{funcdesc}{point_on_line}{x, y, x0, y0, x1, y1}
1630  Module: \code{utilities.polygon}
[2637]1631
[2560]1632  Returns \code{True} or \code{False}, depending on whether the point with coordinates
[2637]1633  \code{x, y} is on the line passing through the points with coordinates \code{x0, y0}
[2560]1634  and \code{x1, y1} (extended if necessary at either end).
1635  \end{funcdesc}
[2637]1636
[2560]1637  \begin{funcdesc}{separate_points_by_polygon}{points, polygon,
1638                               closed = True, verbose = False}\indexedcode{separate_points_by_polygon}
[2637]1639  Module: \code{utilities.polygon}
1640
[2560]1641  \end{funcdesc}
[2358]1642
[2560]1643
[2637]1644
[2358]1645\section{coordinate_transforms}
1646
1647\section{geo_spatial_data}
1648
[2637]1649This describes a class that represents arbitrary point data in UTM
[2497]1650coordinates along with named attribute values.
1651
1652TBA
1653
[2358]1654\section{pmesh GUI}
1655
1656\section{alpha_shape}
1657
1658
[2497]1659\section{utilities/numerical_tools} Do now.
[2358]1660
[2422]1661\begin{itemize}
[2412]1662  \item \indexedcode{ensure_numeric}
1663  \item \indexedcode{mean}
[2422]1664  \item
1665\end{itemize}
[2358]1666
[2463]1667%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1668
[2329]1669\chapter{Glossary}
[2274]1670
[2358]1671\begin{itemize}
[2463]1672    \item \indexedbold{\anuga} Name of software (joint development between ANU and GA)
[2274]1673
[2467]1674    \item \indexedbold{domain}
[2422]1675
[2358]1676    \item \indexedbold{Dirichlet boundary}
[2274]1677
[2467]1678    \item \indexedbold{elevation} - refers to bathymetry and topography
[2274]1679
[2467]1680    \item \indexedbold{bathymetry} offshore
[2274]1681
[2467]1682    \item \indexedbold{topography} onshore
[2328]1683
[2467]1684    \item \indexedbold{evolution} integration of the shallow water wave equations over time
[2328]1685
[2467]1686    \item \indexedbold{forcing term}
[2274]1687
[2463]1688    \item \indexedbold{IDLE} Development environment shipped with Python
[2274]1689
[2358]1690    \item \indexedbold{Manning friction coefficient}
[2422]1691
[2467]1692    \item \indexedbold{mesh}    Triangulation of domain
[2274]1693
[2467]1694    \item \indexedbold{meshfile}  [generic word for either .tsh or
[2457]1695    .msh file]
1696
[2467]1697    \item \indexedbold{points file}  [generic word for either .pts or
[2457]1698    .xya file]
1699
[2467]1700    \item \indexedbold{grid} evenly spaced
[2328]1701
[2358]1702    \item \indexedbold{NetCDF}
[2274]1703
[2358]1704    \item \indexedbold{pmesh} does this really need to be here? it's a class/module?
[2274]1705
[2358]1706    \item \indexedbold{pyvolution} does this really need to be here? it's a class/module?
[2274]1707
[2467]1708    \item \indexedbold{conserved quantity} conserved (state, x and y momentum)
[2274]1709
[2467]1710    \item \indexedbold{reflective boundary}
[2274]1711
[2467]1712    \item \indexedbold{smoothing} is this really needed?
[2274]1713
[2467]1714    \item \indexedbold{stage}
[2274]1715
[2467]1716%    \item \indexedbold{try this}
[2274]1717
[2467]1718    \item \indexedbold{swollen} visualisation tool
[2274]1719
[2467]1720    \item \indexedbold{time boundary} defined in the manual (flog from there)
[2274]1721
[2467]1722    \item \indexedbold{transmissive boundary} defined in the manual (flog from there)
1723
[2358]1724    \item \indexedbold{xmomentum} conserved quantity (note, two-dimensional SWW equations say only x and y and NOT z)
[2274]1725
[2422]1726    \item \indexedbold{ymomentum}  conserved quantity
1727
[2467]1728    \item \indexedbold{resolution}   The maximal area of a triangular cell in a mesh
[2422]1729
[2467]1730    \item \indexedbold{polygon} A sequence of points in the plane. (Arbitrary polygons can be created
[2463]1731    in this way.)
1732    \anuga represents a polygon in one of two ways. One way is to represent it as a
1733    list whose members are either Python tuples
1734    or Python lists of length 2. The unit square, for example, would be represented by the
1735    list
1736    [ [0,0], [1,0], [1,1], [0,1] ]. The alternative is to represent it as an
1737    $N \times 2$ Numeric array, where $N$ is the number of points.
[2422]1738
[2358]1739    NOTE: More can be read in the module utilities/polygon.py ....
[2274]1740
[2467]1741    \item \indexedbold{easting}
[2328]1742
[2467]1743    \item \indexedbold{northing}
[2328]1744
[2467]1745    \item \indexedbold{latitude}
[2328]1746
[2467]1747    \item \indexedbold{longitude}
[2328]1748
[2467]1749    \item \indexedbold{edge}
[2328]1750
[2467]1751    \item \indexedbold{vertex}
[2328]1752
[2467]1753    \item \indexedbold{finite volume}
[2328]1754
[2467]1755    \item \indexedbold{flux}
[2328]1756
[2422]1757    \item \indexedbold{Digital Elevation Model (DEM)}
[2328]1758
[2422]1759
[2358]1760\end{itemize}
[2274]1761
[2329]1762The \code{\e appendix} markup need not be repeated for additional
1763appendices.
1764
1765
1766%
1767%  The ugly "%begin{latexonly}" pseudo-environments are really just to
1768%  keep LaTeX2HTML quiet during the \renewcommand{} macros; they're
1769%  not really valuable.
1770%
1771%  If you don't want the Module Index, you can remove all of this up
1772%  until the second \input line.
1773%
1774
1775%begin{latexonly}
1776%\renewcommand{\indexname}{Module Index}
1777%end{latexonly}
[2422]1778%\input{mod\jobname.ind}        % Module Index
[2329]1779
1780%begin{latexonly}
1781\renewcommand{\indexname}{Index}
1782%end{latexonly}
[2422]1783\input{\jobname.ind}            % Index
[2329]1784
1785
1786
[2274]1787\end{document}
Note: See TracBrowser for help on using the repository browser.