[4123] | 1 | % Complete documentation on the extended LaTeX markup used for Python |
---|
| 2 | % documentation is available in ``Documenting Python'', which is part |
---|
| 3 | % of the standard documentation for Python. It may be found online |
---|
| 4 | % at: |
---|
| 5 | % |
---|
| 6 | % http://www.python.org/doc/current/doc/doc.html |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | %labels |
---|
| 10 | %Sections and subsections \label{sec: } |
---|
| 11 | %Chapters \label{ch: } |
---|
| 12 | %Equations \label{eq: } |
---|
| 13 | %Figures \label{fig: } |
---|
| 14 | |
---|
| 15 | % Is latex failing with; |
---|
| 16 | % `modanuga_user_manual.ind' not found? |
---|
| 17 | % try this command-line |
---|
| 18 | % makeindex modanuga_user_manual.idx |
---|
| 19 | % To produce the modanuga_user_manual.ind file. |
---|
| 20 | |
---|
| 21 | |
---|
[5730] | 22 | %%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%% |
---|
| 23 | % |
---|
| 24 | % ensure_geospatial |
---|
| 25 | % ensure_absolute |
---|
| 26 | % set_geo_reference |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | |
---|
[4123] | 31 | \documentclass{manual} |
---|
| 32 | |
---|
| 33 | \usepackage{graphicx} |
---|
[5744] | 34 | \usepackage[english]{babel} |
---|
[4123] | 35 | \usepackage{datetime} |
---|
| 36 | |
---|
| 37 | \input{definitions} |
---|
| 38 | |
---|
| 39 | \title{\anuga User Manual} |
---|
| 40 | \author{Geoscience Australia and the Australian National University} |
---|
| 41 | |
---|
| 42 | % Please at least include a long-lived email address; |
---|
| 43 | % the rest is at your discretion. |
---|
| 44 | \authoraddress{Geoscience Australia \\ |
---|
| 45 | Email: \email{ole.nielsen@ga.gov.au} |
---|
| 46 | } |
---|
| 47 | |
---|
| 48 | %Draft date |
---|
| 49 | |
---|
| 50 | % update before release! |
---|
| 51 | % Use an explicit date so that reformatting |
---|
| 52 | % doesn't cause a new date to be used. Setting |
---|
| 53 | % the date to \today can be used during draft |
---|
| 54 | % stages to make it easier to handle versions. |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | \longdate % Make date format long using datetime.sty |
---|
| 58 | %\settimeformat{xxivtime} % 24 hour Format |
---|
| 59 | \settimeformat{oclock} % Verbose |
---|
| 60 | \date{\today, \ \currenttime} |
---|
| 61 | %\hyphenation{set\_datadir} |
---|
| 62 | |
---|
| 63 | \ifhtml |
---|
| 64 | \date{\today} % latex2html does not know about datetime |
---|
| 65 | \fi |
---|
| 66 | |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | |
---|
[4785] | 70 | \input{version} % Get version info - this file may be modified by |
---|
[4953] | 71 | % update_anuga_user_manual.py - if not a dummy |
---|
[4785] | 72 | % will be used. |
---|
| 73 | |
---|
| 74 | %\release{1.0} % release version; this is used to define the |
---|
| 75 | % % \version macro |
---|
[4123] | 76 | |
---|
| 77 | \makeindex % tell \index to actually write the .idx file |
---|
| 78 | \makemodindex % If this contains a lot of module sections. |
---|
| 79 | |
---|
| 80 | \setcounter{tocdepth}{3} |
---|
| 81 | \setcounter{secnumdepth}{3} |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | \begin{document} |
---|
| 85 | \maketitle |
---|
| 86 | |
---|
| 87 | |
---|
[5744] | 88 | |
---|
[4123] | 89 | % This makes the contents more accessible from the front page of the HTML. |
---|
| 90 | \ifhtml |
---|
| 91 | \chapter*{Front Matter\label{front}} |
---|
| 92 | \fi |
---|
| 93 | |
---|
| 94 | %Subversion keywords: |
---|
| 95 | % |
---|
| 96 | %$LastChangedDate: 2008-09-18 01:03:45 +0000 (Thu, 18 Sep 2008) $ |
---|
| 97 | %$LastChangedRevision: 5765 $ |
---|
| 98 | %$LastChangedBy: rwilson $ |
---|
| 99 | |
---|
| 100 | \input{copyright} |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | \begin{abstract} |
---|
| 104 | \label{def:anuga} |
---|
| 105 | |
---|
| 106 | \noindent \anuga\index{\anuga} is a hydrodynamic modelling tool that |
---|
[5566] | 107 | allows users to model realistic flow problems in complex 2D geometries. |
---|
[4123] | 108 | Examples include dam breaks or the effects of natural hazards such |
---|
| 109 | as riverine flooding, storm surges and tsunami. |
---|
| 110 | |
---|
| 111 | The user must specify a study area represented by a mesh of triangular |
---|
| 112 | cells, the topography and bathymetry, frictional resistance, initial |
---|
| 113 | values for water level (called \emph{stage}\index{stage} within \anuga), |
---|
[5566] | 114 | boundary conditions and forces such as rainfall, stream flows, windstress or pressure gradients if applicable. |
---|
[4123] | 115 | |
---|
| 116 | \anuga tracks the evolution of water depth and horizontal momentum |
---|
| 117 | within each cell over time by solving the shallow water wave equation |
---|
| 118 | governing equation using a finite-volume method. |
---|
| 119 | |
---|
[4736] | 120 | \anuga also incorporates a mesh generator %, called \code{graphical |
---|
| 121 | %mesh generator}, |
---|
[4123] | 122 | that |
---|
| 123 | allows the user to set up the geometry of the problem interactively as |
---|
| 124 | well as tools for interpolation and surface fitting, and a number of |
---|
| 125 | auxiliary tools for visualising and interrogating the model output. |
---|
| 126 | |
---|
| 127 | Most \anuga components are written in the object-oriented programming |
---|
| 128 | language Python and most users will interact with \anuga by writing |
---|
| 129 | small Python programs based on the \anuga library |
---|
| 130 | functions. Computationally intensive components are written for |
---|
| 131 | efficiency in C routines working directly with the Numerical Python |
---|
| 132 | structures. |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | \end{abstract} |
---|
| 136 | |
---|
| 137 | \tableofcontents |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | \chapter{Introduction} |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | \section{Purpose} |
---|
| 144 | |
---|
| 145 | The purpose of this user manual is to introduce the new user to the |
---|
| 146 | inundation software, describe what it can do and give step-by-step |
---|
| 147 | instructions for setting up and running hydrodynamic simulations. |
---|
| 148 | |
---|
| 149 | \section{Scope} |
---|
| 150 | |
---|
| 151 | This manual covers only what is needed to operate the software after |
---|
| 152 | installation and configuration. It does not includes instructions |
---|
| 153 | for installing the software or detailed API documentation, both of |
---|
| 154 | which will be covered in separate publications and by documentation |
---|
| 155 | in the source code. |
---|
| 156 | |
---|
| 157 | \section{Audience} |
---|
| 158 | |
---|
[5744] | 159 | Readers are assumed to be familiar with the Python Programming language and |
---|
[5129] | 160 | its object oriented approach. |
---|
[5744] | 161 | Python tutorials include |
---|
[5129] | 162 | \url{http://docs.python.org/tut}, |
---|
| 163 | \url{http://www.sthurlow.com/python}, and |
---|
[5130] | 164 | %\url{http://datamining.anu.edu.au/\%7e ole/work/teaching/ctac2006/exercise1.pdf}. |
---|
| 165 | \url{http://datamining.anu.edu.au/\~{}ole/work/teaching/ctac2006/exercise1.pdf}. |
---|
[5129] | 166 | |
---|
| 167 | |
---|
| 168 | Readers also need to have a general understanding of scientific modelling, |
---|
| 169 | as well as |
---|
[4123] | 170 | enough programming experience to adapt the code to different |
---|
[5129] | 171 | requirements. |
---|
[4123] | 172 | |
---|
[5129] | 173 | |
---|
| 174 | |
---|
[4123] | 175 | \pagebreak |
---|
| 176 | \chapter{Background} |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | Modelling the effects on the built environment of natural hazards such |
---|
| 180 | as riverine flooding, storm surges and tsunami is critical for |
---|
| 181 | understanding their economic and social impact on our urban |
---|
| 182 | communities. Geoscience Australia and the Australian National |
---|
| 183 | University are developing a hydrodynamic inundation modelling tool |
---|
| 184 | called \anuga to help simulate the impact of these hazards. |
---|
| 185 | |
---|
| 186 | The core of \anuga is the fluid dynamics module, called \code{shallow\_water}, |
---|
| 187 | which is based on a finite-volume method for solving the Shallow Water |
---|
| 188 | Wave Equation. The study area is represented by a mesh of triangular |
---|
| 189 | cells. By solving the governing equation within each cell, water |
---|
| 190 | depth and horizontal momentum are tracked over time. |
---|
| 191 | |
---|
| 192 | A major capability of \anuga is that it can model the process of |
---|
| 193 | wetting and drying as water enters and leaves an area. This means |
---|
| 194 | that it is suitable for simulating water flow onto a beach or dry land |
---|
| 195 | and around structures such as buildings. \anuga is also capable |
---|
| 196 | of modelling hydraulic jumps due to the ability of the finite-volume |
---|
| 197 | method to accommodate discontinuities in the solution. |
---|
| 198 | |
---|
| 199 | To set up a particular scenario the user specifies the geometry |
---|
| 200 | (bathymetry and topography), the initial water level (stage), |
---|
| 201 | boundary conditions such as tide, and any forcing terms that may |
---|
[5506] | 202 | drive the system such as rain_fall, abstraction of water, wind stress or atmospheric pressure |
---|
[4123] | 203 | gradients. Gravity and frictional resistance from the different |
---|
| 204 | terrains in the model are represented by predefined forcing terms. |
---|
[5506] | 205 | See section \ref{sec:forcing terms} for details on forcing terms available in ANUGA. |
---|
[4123] | 206 | |
---|
[4673] | 207 | The built-in mesh generator, called \code{graphical\_mesh\_generator}, |
---|
[4123] | 208 | allows the user to set up the geometry |
---|
| 209 | of the problem interactively and to identify boundary segments and |
---|
| 210 | regions using symbolic tags. These tags may then be used to set the |
---|
| 211 | actual boundary conditions and attributes for different regions |
---|
| 212 | (e.g.\ the Manning friction coefficient) for each simulation. |
---|
| 213 | |
---|
| 214 | Most \anuga components are written in the object-oriented programming |
---|
| 215 | language Python. Software written in Python can be produced quickly |
---|
| 216 | and can be readily adapted to changing requirements throughout its |
---|
| 217 | lifetime. Computationally intensive components are written for |
---|
| 218 | efficiency in C routines working directly with the Numerical Python |
---|
| 219 | structures. The animation tool developed for \anuga is based on |
---|
| 220 | OpenSceneGraph, an Open Source Software (OSS) component allowing high |
---|
| 221 | level interaction with sophisticated graphics primitives. |
---|
| 222 | See \cite{nielsen2005} for more background on \anuga. |
---|
| 223 | |
---|
| 224 | \chapter{Restrictions and limitations on \anuga} |
---|
| 225 | \label{ch:limitations} |
---|
| 226 | |
---|
| 227 | Although a powerful and flexible tool for hydrodynamic modelling, \anuga has a |
---|
| 228 | number of limitations that any potential user need to be aware of. They are |
---|
| 229 | |
---|
| 230 | \begin{itemize} |
---|
[4209] | 231 | \item The mathematical model is the 2D shallow water wave equation. |
---|
| 232 | As such it cannot resolve vertical convection and consequently not breaking |
---|
[4123] | 233 | waves or 3D turbulence (e.g.\ vorticity). |
---|
[5566] | 234 | %\item The surface is assumed to be open, e.g. \anuga cannot model |
---|
| 235 | %flow under ceilings or in pipes |
---|
[4209] | 236 | \item All spatial coordinates are assumed to be UTM (meters). As such, |
---|
| 237 | ANUGA is unsuitable for modelling flows in areas larger than one UTM zone |
---|
| 238 | (6 degrees wide). |
---|
[5744] | 239 | \item Fluid is assumed to be inviscid - i.e.\ no kinematic viscosity included. |
---|
[4209] | 240 | \item The finite volume is a very robust and flexible numerical technique, |
---|
| 241 | but it is not the fastest method around. If the geometry is sufficiently |
---|
| 242 | simple and if there is no need for wetting or drying, a finite-difference |
---|
[4123] | 243 | method may be able to solve the problem faster than \anuga. |
---|
[4209] | 244 | %\item Mesh resolutions near coastlines with steep gradients need to be... |
---|
| 245 | \item Frictional resistance is implemented using Manning's formula, but |
---|
[4123] | 246 | \anuga has not yet been fully validated in regard to bottom roughness |
---|
[5566] | 247 | %\item ANUGA contains no tsunami-genic functionality relating to |
---|
| 248 | %earthquakes. |
---|
[4123] | 249 | \end{itemize} |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | \chapter{Getting Started} |
---|
| 254 | \label{ch:getstarted} |
---|
| 255 | |
---|
| 256 | This section is designed to assist the reader to get started with |
---|
| 257 | \anuga by working through some examples. Two examples are discussed; |
---|
[5566] | 258 | the first is a simple example to illustrate many of the concepts, and |
---|
[4123] | 259 | the second is a more realistic example. |
---|
| 260 | |
---|
| 261 | \section{A Simple Example} |
---|
| 262 | \label{sec:simpleexample} |
---|
| 263 | |
---|
| 264 | \subsection{Overview} |
---|
| 265 | |
---|
| 266 | What follows is a discussion of the structure and operation of a |
---|
| 267 | script called \file{runup.py}. |
---|
| 268 | |
---|
| 269 | This example carries out the solution of the shallow-water wave |
---|
| 270 | equation in the simple case of a configuration comprising a flat |
---|
| 271 | bed, sloping at a fixed angle in one direction and having a |
---|
| 272 | constant depth across each line in the perpendicular direction. |
---|
| 273 | |
---|
| 274 | The example demonstrates the basic ideas involved in setting up a |
---|
| 275 | complex scenario. In general the user specifies the geometry |
---|
| 276 | (bathymetry and topography), the initial water level, boundary |
---|
| 277 | conditions such as tide, and any forcing terms that may drive the |
---|
[5506] | 278 | system such as rain_fall, abstraction of water, wind stress or atmospheric pressure gradients. |
---|
[4123] | 279 | Frictional resistance from the different terrains in the model is |
---|
| 280 | represented by predefined forcing terms. In this example, the |
---|
| 281 | boundary is reflective on three sides and a time dependent wave on |
---|
| 282 | one side. |
---|
| 283 | |
---|
| 284 | The present example represents a simple scenario and does not |
---|
| 285 | include any forcing terms, nor is the data taken from a file as it |
---|
| 286 | would typically be. |
---|
| 287 | |
---|
| 288 | The conserved quantities involved in the |
---|
| 289 | problem are stage (absolute height of water surface), |
---|
| 290 | $x$-momentum and $y$-momentum. Other quantities |
---|
| 291 | involved in the computation are the friction and elevation. |
---|
| 292 | |
---|
| 293 | Water depth can be obtained through the equation |
---|
| 294 | |
---|
| 295 | \begin{tabular}{rcrcl} |
---|
| 296 | \code{depth} &=& \code{stage} &$-$& \code{elevation} |
---|
| 297 | \end{tabular} |
---|
| 298 | |
---|
| 299 | |
---|
| 300 | \subsection{Outline of the Program} |
---|
| 301 | |
---|
| 302 | In outline, \file{runup.py} performs the following steps: |
---|
| 303 | |
---|
| 304 | \begin{enumerate} |
---|
| 305 | |
---|
| 306 | \item Sets up a triangular mesh. |
---|
| 307 | |
---|
| 308 | \item Sets certain parameters governing the mode of |
---|
| 309 | operation of the model-specifying, for instance, where to store the model output. |
---|
| 310 | |
---|
| 311 | \item Inputs various quantities describing physical measurements, such |
---|
| 312 | as the elevation, to be specified at each mesh point (vertex). |
---|
| 313 | |
---|
| 314 | \item Sets up the boundary conditions. |
---|
| 315 | |
---|
| 316 | \item Carries out the evolution of the model through a series of time |
---|
| 317 | steps and outputs the results, providing a results file that can |
---|
| 318 | be visualised. |
---|
| 319 | |
---|
| 320 | \end{enumerate} |
---|
| 321 | |
---|
| 322 | \subsection{The Code} |
---|
| 323 | |
---|
| 324 | %FIXME: we are using the \code function here. |
---|
| 325 | %This should be used wherever possible |
---|
| 326 | For reference we include below the complete code listing for |
---|
| 327 | \file{runup.py}. Subsequent paragraphs provide a |
---|
| 328 | `commentary' that describes each step of the program and explains it |
---|
| 329 | significance. |
---|
| 330 | |
---|
| 331 | \verbatiminput{demos/runup.py} |
---|
| 332 | |
---|
| 333 | \subsection{Establishing the Mesh}\index{mesh, establishing} |
---|
| 334 | |
---|
| 335 | The first task is to set up the triangular mesh to be used for the |
---|
| 336 | scenario. This is carried out through the statement: |
---|
| 337 | |
---|
| 338 | {\small \begin{verbatim} |
---|
[4953] | 339 | points, vertices, boundary = rectangular_cross(10, 10) |
---|
[4123] | 340 | \end{verbatim}} |
---|
| 341 | |
---|
[4953] | 342 | The function \function{rectangular_cross} is imported from a module |
---|
[4123] | 343 | \module{mesh\_factory} defined elsewhere. (\anuga also contains |
---|
| 344 | several other schemes that can be used for setting up meshes, but we |
---|
| 345 | shall not discuss these.) The above assignment sets up a $10 \times |
---|
| 346 | 10$ rectangular mesh, triangulated in a regular way. The assignment |
---|
| 347 | |
---|
| 348 | {\small \begin{verbatim} |
---|
[4953] | 349 | points, vertices, boundary = rectangular_cross(m, n) |
---|
[4123] | 350 | \end{verbatim}} |
---|
| 351 | |
---|
| 352 | returns: |
---|
| 353 | |
---|
| 354 | \begin{itemize} |
---|
| 355 | |
---|
| 356 | \item a list \code{points} giving the coordinates of each mesh point, |
---|
| 357 | |
---|
| 358 | \item a list \code{vertices} specifying the three vertices of each triangle, and |
---|
| 359 | |
---|
| 360 | \item a dictionary \code{boundary} that stores the edges on |
---|
| 361 | the boundary and associates each with one of the symbolic tags \code{`left'}, \code{`right'}, |
---|
| 362 | \code{`top'} or \code{`bottom'}. |
---|
| 363 | |
---|
| 364 | \end{itemize} |
---|
| 365 | |
---|
| 366 | (For more details on symbolic tags, see page |
---|
| 367 | \pageref{ref:tagdescription}.) |
---|
| 368 | |
---|
| 369 | An example of a general unstructured mesh and the associated data |
---|
| 370 | structures \code{points}, \code{vertices} and \code{boundary} is |
---|
| 371 | given in Section \ref{sec:meshexample}. |
---|
| 372 | |
---|
| 373 | |
---|
| 374 | |
---|
| 375 | |
---|
| 376 | \subsection{Initialising the Domain} |
---|
| 377 | |
---|
| 378 | These variables are then used to set up a data structure |
---|
| 379 | \code{domain}, through the assignment: |
---|
| 380 | |
---|
| 381 | {\small \begin{verbatim} |
---|
| 382 | domain = Domain(points, vertices, boundary) |
---|
| 383 | \end{verbatim}} |
---|
| 384 | |
---|
| 385 | This creates an instance of the \class{Domain} class, which |
---|
| 386 | represents the domain of the simulation. Specific options are set at |
---|
| 387 | this point, including the basename for the output file and the |
---|
| 388 | directory to be used for data: |
---|
| 389 | |
---|
| 390 | {\small \begin{verbatim} |
---|
| 391 | domain.set_name('runup') |
---|
| 392 | \end{verbatim}} |
---|
| 393 | |
---|
| 394 | {\small \begin{verbatim} |
---|
| 395 | domain.set_datadir('.') |
---|
| 396 | \end{verbatim}} |
---|
| 397 | |
---|
| 398 | In addition, the following statement now specifies that the |
---|
| 399 | quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are |
---|
| 400 | to be stored: |
---|
| 401 | |
---|
| 402 | {\small \begin{verbatim} |
---|
| 403 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', |
---|
| 404 | 'ymomentum']) |
---|
| 405 | \end{verbatim}} |
---|
| 406 | |
---|
| 407 | |
---|
| 408 | \subsection{Initial Conditions} |
---|
| 409 | |
---|
| 410 | The next task is to specify a number of quantities that we wish to |
---|
| 411 | set for each mesh point. The class \class{Domain} has a method |
---|
| 412 | \method{set\_quantity}, used to specify these quantities. It is a |
---|
| 413 | flexible method that allows the user to set quantities in a variety |
---|
| 414 | of ways---using constants, functions, numeric arrays, expressions |
---|
| 415 | involving other quantities, or arbitrary data points with associated |
---|
| 416 | values, all of which can be passed as arguments. All quantities can |
---|
| 417 | be initialised using \method{set\_quantity}. For a conserved |
---|
| 418 | quantity (such as \code{stage, xmomentum, ymomentum}) this is called |
---|
| 419 | an \emph{initial condition}. However, other quantities that aren't |
---|
| 420 | updated by the equation are also assigned values using the same |
---|
| 421 | interface. The code in the present example demonstrates a number of |
---|
| 422 | forms in which we can invoke \method{set\_quantity}. |
---|
| 423 | |
---|
| 424 | |
---|
| 425 | \subsubsection{Elevation} |
---|
| 426 | |
---|
| 427 | The elevation, or height of the bed, is set using a function, |
---|
| 428 | defined through the statements below, which is specific to this |
---|
| 429 | example and specifies a particularly simple initial configuration |
---|
| 430 | for demonstration purposes: |
---|
| 431 | |
---|
| 432 | {\small \begin{verbatim} |
---|
| 433 | def f(x,y): |
---|
| 434 | return -x/2 |
---|
| 435 | \end{verbatim}} |
---|
| 436 | |
---|
| 437 | This simply associates an elevation with each point \code{(x, y)} of |
---|
| 438 | the plane. It specifies that the bed slopes linearly in the |
---|
| 439 | \code{x} direction, with slope $-\frac{1}{2}$, and is constant in |
---|
| 440 | the \code{y} direction. |
---|
| 441 | |
---|
| 442 | Once the function \function{f} is specified, the quantity |
---|
| 443 | \code{elevation} is assigned through the simple statement: |
---|
| 444 | |
---|
| 445 | {\small \begin{verbatim} |
---|
| 446 | domain.set_quantity('elevation', f) |
---|
| 447 | \end{verbatim}} |
---|
| 448 | |
---|
[4953] | 449 | NOTE: If using function to set \code{elevation} it must be vector |
---|
[4743] | 450 | compatible. For example square root will not work. |
---|
[4123] | 451 | |
---|
| 452 | \subsubsection{Friction} |
---|
| 453 | |
---|
| 454 | The assignment of the friction quantity (a forcing term) |
---|
| 455 | demonstrates another way we can use \method{set\_quantity} to set |
---|
| 456 | quantities---namely, assign them to a constant numerical value: |
---|
| 457 | |
---|
| 458 | {\small \begin{verbatim} |
---|
| 459 | domain.set_quantity('friction', 0.1) |
---|
| 460 | \end{verbatim}} |
---|
| 461 | |
---|
| 462 | This specifies that the Manning friction coefficient is set to 0.1 |
---|
| 463 | at every mesh point. |
---|
| 464 | |
---|
| 465 | \subsubsection{Stage} |
---|
| 466 | |
---|
| 467 | The stage (the height of the water surface) is related to the |
---|
| 468 | elevation and the depth at any time by the equation |
---|
| 469 | |
---|
| 470 | {\small \begin{verbatim} |
---|
| 471 | stage = elevation + depth |
---|
| 472 | \end{verbatim}} |
---|
| 473 | |
---|
| 474 | |
---|
| 475 | For this example, we simply assign a constant value to \code{stage}, |
---|
| 476 | using the statement |
---|
| 477 | |
---|
| 478 | {\small \begin{verbatim} |
---|
| 479 | domain.set_quantity('stage', -.4) |
---|
| 480 | \end{verbatim}} |
---|
| 481 | |
---|
| 482 | which specifies that the surface level is set to a height of $-0.4$, |
---|
| 483 | i.e. 0.4 units (m) below the zero level. |
---|
| 484 | |
---|
| 485 | Although it is not necessary for this example, it may be useful to |
---|
| 486 | digress here and mention a variant to this requirement, which allows |
---|
| 487 | us to illustrate another way to use \method{set\_quantity}---namely, |
---|
| 488 | incorporating an expression involving other quantities. Suppose, |
---|
| 489 | instead of setting a constant value for the stage, we wished to |
---|
| 490 | specify a constant value for the \emph{depth}. For such a case we |
---|
| 491 | need to specify that \code{stage} is everywhere obtained by adding |
---|
| 492 | that value to the value already specified for \code{elevation}. We |
---|
| 493 | would do this by means of the statements: |
---|
| 494 | |
---|
| 495 | {\small \begin{verbatim} |
---|
| 496 | h = 0.05 # Constant depth |
---|
| 497 | domain.set_quantity('stage', expression = 'elevation + %f' %h) |
---|
| 498 | \end{verbatim}} |
---|
| 499 | |
---|
| 500 | That is, the value of \code{stage} is set to $\code{h} = 0.05$ plus |
---|
| 501 | the value of \code{elevation} already defined. |
---|
| 502 | |
---|
| 503 | The reader will probably appreciate that this capability to |
---|
| 504 | incorporate expressions into statements using \method{set\_quantity} |
---|
[5508] | 505 | greatly expands its power.) See Section \ref{sec:initial conditions} for more |
---|
[4123] | 506 | details. |
---|
| 507 | |
---|
| 508 | \subsection{Boundary Conditions}\index{boundary conditions} |
---|
| 509 | |
---|
| 510 | The boundary conditions are specified as follows: |
---|
| 511 | |
---|
| 512 | {\small \begin{verbatim} |
---|
| 513 | Br = Reflective_boundary(domain) |
---|
| 514 | |
---|
| 515 | Bt = Transmissive_boundary(domain) |
---|
| 516 | |
---|
| 517 | Bd = Dirichlet_boundary([0.2,0.,0.]) |
---|
| 518 | |
---|
| 519 | Bw = Time_boundary(domain=domain, |
---|
| 520 | f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0]) |
---|
| 521 | \end{verbatim}} |
---|
| 522 | |
---|
| 523 | The effect of these statements is to set up a selection of different |
---|
| 524 | alternative boundary conditions and store them in variables that can be |
---|
| 525 | assigned as needed. Each boundary condition specifies the |
---|
| 526 | behaviour at a boundary in terms of the behaviour in neighbouring |
---|
| 527 | elements. The boundary conditions introduced here may be briefly described as |
---|
| 528 | follows: |
---|
| 529 | |
---|
| 530 | \begin{itemize} |
---|
| 531 | \item \textbf{Reflective boundary}\label{def:reflective boundary} Returns same \code{stage} as |
---|
| 532 | as present in its neighbour volume but momentum vector |
---|
| 533 | reversed 180 degrees (reflected). |
---|
| 534 | Specific to the shallow water equation as it works with the |
---|
| 535 | momentum quantities assumed to be the second and third conserved |
---|
| 536 | quantities. A reflective boundary condition models a solid wall. |
---|
[5765] | 537 | \item \textbf{Transmissive boundary}\label{def:transmissive boundary} |
---|
| 538 | Returns same conserved quantities as |
---|
[4123] | 539 | those present in its neighbour volume. This is one way of modelling |
---|
| 540 | outflow from a domain, but it should be used with caution if flow is |
---|
| 541 | not steady state as replication of momentum at the boundary |
---|
[5765] | 542 | may cause numerical instabilities propagating into the domain and |
---|
| 543 | eventually causing ANUGA to crash. If this occurs, |
---|
| 544 | consider using e.g. a Dirichlet boundary condition with a stage value |
---|
| 545 | less than the elevation at the boundary. |
---|
[4123] | 546 | \item \textbf{Dirichlet boundary}\label{def:dirichlet boundary} Specifies |
---|
| 547 | constant values for stage, $x$-momentum and $y$-momentum at the boundary. |
---|
| 548 | \item \textbf{Time boundary}\label{def:time boundary} Like a Dirichlet |
---|
| 549 | boundary but with behaviour varying with time. |
---|
| 550 | \end{itemize} |
---|
| 551 | |
---|
| 552 | \label{ref:tagdescription}Before describing how these boundary |
---|
| 553 | conditions are assigned, we recall that a mesh is specified using |
---|
| 554 | three variables \code{points}, \code{vertices} and \code{boundary}. |
---|
| 555 | In the code we are discussing, these three variables are returned by |
---|
| 556 | the function \code{rectangular}; however, the example given in |
---|
| 557 | Section \ref{sec:realdataexample} illustrates another way of |
---|
| 558 | assigning the values, by means of the function |
---|
| 559 | \code{create\_mesh\_from\_regions}. |
---|
| 560 | |
---|
| 561 | These variables store the data determining the mesh as follows. (You |
---|
| 562 | may find that the example given in Section \ref{sec:meshexample} |
---|
| 563 | helps to clarify the following discussion, even though that example |
---|
| 564 | is a \emph{non-rectangular} mesh.) |
---|
| 565 | |
---|
| 566 | \begin{itemize} |
---|
| 567 | \item The variable \code{points} stores a list of 2-tuples giving the |
---|
| 568 | coordinates of the mesh points. |
---|
| 569 | |
---|
| 570 | \item The variable \code{vertices} stores a list of 3-tuples of |
---|
| 571 | numbers, representing vertices of triangles in the mesh. In this |
---|
| 572 | list, the triangle whose vertices are \code{points[i]}, |
---|
| 573 | \code{points[j]}, \code{points[k]} is represented by the 3-tuple |
---|
| 574 | \code{(i, j, k)}. |
---|
| 575 | |
---|
| 576 | \item The variable \code{boundary} is a Python dictionary that |
---|
| 577 | not only stores the edges that make up the boundary but also assigns |
---|
| 578 | symbolic tags to these edges to distinguish different parts of the |
---|
| 579 | boundary. An edge with endpoints \code{points[i]} and |
---|
| 580 | \code{points[j]} is represented by the 2-tuple \code{(i, j)}. The |
---|
| 581 | keys for the dictionary are the 2-tuples \code{(i, j)} corresponding |
---|
| 582 | to boundary edges in the mesh, and the values are the tags are used |
---|
| 583 | to label them. In the present example, the value \code{boundary[(i, |
---|
| 584 | j)]} assigned to \code{(i, j)]} is one of the four tags |
---|
| 585 | \code{`left'}, \code{`right'}, \code{`top'} or \code{`bottom'}, |
---|
| 586 | depending on whether the boundary edge represented by \code{(i, j)} |
---|
| 587 | occurs at the left, right, top or bottom of the rectangle bounding |
---|
| 588 | the mesh. The function \code{rectangular} automatically assigns |
---|
| 589 | these tags to the boundary edges when it generates the mesh. |
---|
| 590 | \end{itemize} |
---|
| 591 | |
---|
| 592 | The tags provide the means to assign different boundary conditions |
---|
| 593 | to an edge depending on which part of the boundary it belongs to. |
---|
| 594 | (In Section \ref{sec:realdataexample} we describe an example that |
---|
[4673] | 595 | uses different boundary tags --- in general, the possible tags are entirely selectable by the user when generating the mesh and not |
---|
| 596 | limited to `left', `right', `top' and `bottom' as in this example.) |
---|
| 597 | All segments in bounding polygon must be tagged. If a tag is not supplied, the default tag name 'exterior' will be assigned by ANUGA. |
---|
[4123] | 598 | |
---|
[4673] | 599 | |
---|
[4123] | 600 | Using the boundary objects described above, we assign a boundary |
---|
| 601 | condition to each part of the boundary by means of a statement like |
---|
| 602 | |
---|
| 603 | {\small \begin{verbatim} |
---|
| 604 | domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br}) |
---|
| 605 | \end{verbatim}} |
---|
| 606 | |
---|
[4673] | 607 | It is critical that all tags are assoctiated with a boundary conditing in this statement. If not the program will halt with a statement like |
---|
| 608 | |
---|
| 609 | \begin{verbatim} |
---|
| 610 | |
---|
| 611 | Traceback (most recent call last): |
---|
| 612 | File "mesh_test.py", line 114, in ? |
---|
| 613 | domain.set_boundary({'west': Bi, 'east': Bo, 'north': Br, 'south': Br}) |
---|
| 614 | File "X:\inundation\sandpits\onielsen\anuga_core\source\anuga\abstract_2d_finite_volumes\domain.py", line 505, in set_boundary |
---|
| 615 | raise msg |
---|
| 616 | ERROR (domain.py): Tag "exterior" has not been bound to a boundary object. |
---|
| 617 | All boundary tags defined in domain must appear in the supplied dictionary. |
---|
| 618 | The tags are: ['ocean', 'east', 'north', 'exterior', 'south'] |
---|
[4953] | 619 | \end{verbatim} |
---|
[4673] | 620 | |
---|
| 621 | |
---|
| 622 | The command \code{set\_boundary} stipulates that, in the current example, the right |
---|
[4123] | 623 | boundary varies with time, as defined by the lambda function, while the other |
---|
| 624 | boundaries are all reflective. |
---|
| 625 | |
---|
| 626 | The reader may wish to experiment by varying the choice of boundary |
---|
| 627 | types for one or more of the boundaries. (In the case of \code{Bd} |
---|
| 628 | and \code{Bw}, the three arguments in each case represent the |
---|
| 629 | \code{stage}, $x$-momentum and $y$-momentum, respectively.) |
---|
| 630 | |
---|
| 631 | {\small \begin{verbatim} |
---|
| 632 | Bw = Time_boundary(domain=domain, |
---|
| 633 | f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0]) |
---|
| 634 | \end{verbatim}} |
---|
| 635 | |
---|
| 636 | |
---|
| 637 | |
---|
| 638 | \subsection{Evolution}\index{evolution} |
---|
| 639 | |
---|
| 640 | The final statement \nopagebreak[3] |
---|
| 641 | {\small \begin{verbatim} |
---|
| 642 | for t in domain.evolve(yieldstep = 0.1, duration = 4.0): |
---|
| 643 | print domain.timestepping_statistics() |
---|
| 644 | \end{verbatim}} |
---|
| 645 | |
---|
| 646 | causes the configuration of the domain to `evolve', over a series of |
---|
| 647 | steps indicated by the values of \code{yieldstep} and |
---|
| 648 | \code{duration}, which can be altered as required. The value of |
---|
| 649 | \code{yieldstep} controls the time interval between successive model |
---|
| 650 | outputs. Behind the scenes more time steps are generally taken. |
---|
| 651 | |
---|
| 652 | |
---|
| 653 | \subsection{Output} |
---|
| 654 | |
---|
| 655 | The output is a NetCDF file with the extension \code{.sww}. It |
---|
[4209] | 656 | contains stage and momentum information and can be used with the |
---|
| 657 | ANUGA viewer \code{animate} (see Section \ref{sec:animate}) |
---|
[4123] | 658 | visualisation package |
---|
| 659 | to generate a visual display. See Section \ref{sec:file formats} |
---|
| 660 | (page \pageref{sec:file formats}) for more on NetCDF and other file |
---|
| 661 | formats. |
---|
| 662 | |
---|
| 663 | The following is a listing of the screen output seen by the user |
---|
| 664 | when this example is run: |
---|
| 665 | |
---|
| 666 | \verbatiminput{examples/runupoutput.txt} |
---|
| 667 | |
---|
| 668 | |
---|
| 669 | \section{How to Run the Code} |
---|
| 670 | |
---|
| 671 | The code can be run in various ways: |
---|
| 672 | |
---|
| 673 | \begin{itemize} |
---|
| 674 | \item{from a Windows or Unix command line} as in\ \code{python runup.py} |
---|
| 675 | \item{within the Python IDLE environment} |
---|
| 676 | \item{within emacs} |
---|
| 677 | \item{within Windows, by double-clicking the \code{runup.py} |
---|
| 678 | file.} |
---|
| 679 | \end{itemize} |
---|
| 680 | |
---|
| 681 | |
---|
| 682 | \section{Exploring the Model Output} |
---|
| 683 | |
---|
| 684 | The following figures are screenshots from the \anuga visualisation |
---|
| 685 | tool \code{animate}. Figure \ref{fig:runupstart} shows the domain |
---|
| 686 | with water surface as specified by the initial condition, $t=0$. |
---|
| 687 | Figure \ref{fig:runup2} shows later snapshots for $t=2.3$ and |
---|
| 688 | $t=4$ where the system has been evolved and the wave is encroaching |
---|
| 689 | on the previously dry bed. All figures are screenshots from an |
---|
| 690 | interactive animation tool called animate which is part of \anuga and |
---|
| 691 | distributed as in the package anuga\_viewer. |
---|
| 692 | Animate is described in more detail is Section \ref{sec:animate}. |
---|
| 693 | |
---|
| 694 | \begin{figure}[hbt] |
---|
| 695 | |
---|
| 696 | \centerline{\includegraphics[width=75mm, height=75mm] |
---|
| 697 | {graphics/bedslopestart.jpg}} |
---|
| 698 | |
---|
| 699 | \caption{Runup example viewed with the ANUGA viewer} |
---|
| 700 | \label{fig:runupstart} |
---|
| 701 | \end{figure} |
---|
| 702 | |
---|
| 703 | |
---|
| 704 | \begin{figure}[hbt] |
---|
| 705 | |
---|
| 706 | \centerline{ |
---|
| 707 | \includegraphics[width=75mm, height=75mm]{graphics/bedslopeduring.jpg} |
---|
| 708 | \includegraphics[width=75mm, height=75mm]{graphics/bedslopeend.jpg} |
---|
| 709 | } |
---|
| 710 | |
---|
| 711 | \caption{Runup example viewed with ANGUA viewer} |
---|
| 712 | \label{fig:runup2} |
---|
| 713 | \end{figure} |
---|
| 714 | |
---|
| 715 | |
---|
| 716 | |
---|
| 717 | \clearpage |
---|
| 718 | |
---|
| 719 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 720 | |
---|
| 721 | \section{A slightly more complex example} |
---|
| 722 | \label{sec:channelexample} |
---|
| 723 | |
---|
| 724 | \subsection{Overview} |
---|
| 725 | |
---|
[4209] | 726 | The next example is about waterflow in a channel with varying boundary conditions and |
---|
[4123] | 727 | more complex topograhies. These examples build on the |
---|
| 728 | concepts introduced through the \file{runup.py} in Section \ref{sec:simpleexample}. |
---|
| 729 | The example will be built up through three progressively more complex scripts. |
---|
| 730 | |
---|
| 731 | \subsection{Overview} |
---|
| 732 | As in the case of \file{runup.py}, the actions carried |
---|
| 733 | out by the program can be organised according to this outline: |
---|
| 734 | |
---|
| 735 | \begin{enumerate} |
---|
| 736 | |
---|
| 737 | \item Set up a triangular mesh. |
---|
| 738 | |
---|
| 739 | \item Set certain parameters governing the mode of |
---|
| 740 | operation of the model---specifying, for instance, where to store the |
---|
| 741 | model output. |
---|
| 742 | |
---|
| 743 | \item Set up initial conditions for various quantities such as the elevation, to be specified at each mesh point (vertex). |
---|
| 744 | |
---|
| 745 | \item Set up the boundary conditions. |
---|
| 746 | |
---|
| 747 | \item Carry out the evolution of the model through a series of time |
---|
| 748 | steps and output the results, providing a results file that can be |
---|
| 749 | visualised. |
---|
| 750 | |
---|
| 751 | \end{enumerate} |
---|
| 752 | |
---|
| 753 | |
---|
| 754 | \subsection{The Code} |
---|
| 755 | |
---|
| 756 | Here is the code for the first version of the channel flow \file{channel1.py}: |
---|
| 757 | |
---|
| 758 | \verbatiminput{demos/channel1.py} |
---|
| 759 | |
---|
| 760 | In discussing the details of this example, we follow the outline |
---|
| 761 | given above, discussing each major step of the code in turn. |
---|
| 762 | |
---|
| 763 | \subsection{Establishing the Mesh}\index{mesh, establishing} |
---|
| 764 | |
---|
[4209] | 765 | In this example we use a similar simple structured triangular mesh as in \code{runup.py} |
---|
| 766 | for simplicity, but this time we will use a symmetric one and also |
---|
[4123] | 767 | change the physical extent of the domain. The assignment |
---|
| 768 | |
---|
| 769 | {\small \begin{verbatim} |
---|
[4209] | 770 | points, vertices, boundary = rectangular_cross(m, n, |
---|
[4123] | 771 | len1=length, len2=width) |
---|
| 772 | \end{verbatim}} |
---|
[4209] | 773 | returns a m x n mesh similar to the one used in the previous example, except that now the |
---|
| 774 | extent in the x and y directions are given by the value of \code{length} and \code{width} |
---|
[4123] | 775 | respectively. |
---|
| 776 | |
---|
[4209] | 777 | Defining m and n in terms of the extent as in this example provides a convenient way of |
---|
[4123] | 778 | controlling the resolution: By defining dx and dy to be the desired size of each hypothenuse in the mesh we can write the mesh generation as follows: |
---|
| 779 | |
---|
| 780 | {\small \begin{verbatim} |
---|
| 781 | length = 10. |
---|
| 782 | width = 5. |
---|
| 783 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 784 | |
---|
| 785 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
---|
| 786 | len1=length, len2=width) |
---|
| 787 | \end{verbatim}} |
---|
| 788 | which yields a mesh of length=10m, width=5m with 1m spacings. To increase the resolution, as we will later in this example, one merely decrease the values of dx and dy. |
---|
| 789 | |
---|
| 790 | The rest of this script is as in the previous example. |
---|
| 791 | % except for an application of the 'expression' form of \code{set\_quantity} where we use the value of \code{elevation} to define the (dry) initial condition for \code{stage}: |
---|
| 792 | %{\small \begin{verbatim} |
---|
| 793 | % domain.set_quantity('stage', expression='elevation') |
---|
| 794 | %\end{verbatim}} |
---|
| 795 | |
---|
| 796 | \section{Model Output} |
---|
| 797 | |
---|
| 798 | The following figure is a screenshot from the \anuga visualisation |
---|
| 799 | tool \code{animate} of output from this example. |
---|
| 800 | \begin{figure}[hbt] |
---|
| 801 | \centerline{\includegraphics[height=75mm] |
---|
| 802 | {graphics/channel1.png}}% |
---|
| 803 | |
---|
| 804 | \caption{Simple channel example viewed with the ANUGA viewer.} |
---|
| 805 | \label{fig:channel1} |
---|
| 806 | \end{figure} |
---|
| 807 | |
---|
| 808 | |
---|
| 809 | \subsection{Changing boundary conditions on the fly} |
---|
[4205] | 810 | \label{sec:change boundary} |
---|
[4123] | 811 | |
---|
| 812 | Here is the code for the second version of the channel flow \file{channel2.py}: |
---|
| 813 | \verbatiminput{demos/channel2.py} |
---|
[4209] | 814 | This example differs from the first version in that a constant outflow boundary condition has |
---|
| 815 | been defined |
---|
[4123] | 816 | {\small \begin{verbatim} |
---|
| 817 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
| 818 | \end{verbatim}} |
---|
| 819 | and that it is applied to the right hand side boundary when the water level there exceeds 0m. |
---|
| 820 | {\small \begin{verbatim} |
---|
| 821 | for t in domain.evolve(yieldstep = 0.2, finaltime = 40.0): |
---|
| 822 | domain.write_time() |
---|
| 823 | |
---|
[4209] | 824 | if domain.get_quantity('stage').get_values(interpolation_points=[[10, 2.5]]) > 0: |
---|
[4123] | 825 | print 'Stage > 0: Changing to outflow boundary' |
---|
| 826 | domain.set_boundary({'right': Bo}) |
---|
| 827 | \end{verbatim}} |
---|
[4206] | 828 | \label{sec:change boundary code} |
---|
[4123] | 829 | |
---|
| 830 | The if statement in the timestepping loop (evolve) gets the quantity |
---|
| 831 | \code{stage} and obtain the interpolated value at the point (10m, |
---|
| 832 | 2.5m) which is on the right boundary. If the stage exceeds 0m a |
---|
| 833 | message is printed and the old boundary condition at tag 'right' is |
---|
[4209] | 834 | replaced by the outflow boundary using the method |
---|
| 835 | {\small \begin{verbatim} |
---|
[4123] | 836 | domain.set_boundary({'right': Bo}) |
---|
| 837 | \end{verbatim}} |
---|
[4209] | 838 | This type of dynamically varying boundary could for example be |
---|
| 839 | used to model the |
---|
| 840 | breakdown of a sluice door when water exceeds a certain level. |
---|
[4123] | 841 | |
---|
| 842 | \subsection{Output} |
---|
| 843 | |
---|
| 844 | The text output from this example looks like this |
---|
[4209] | 845 | {\small \begin{verbatim} |
---|
[4123] | 846 | ... |
---|
| 847 | Time = 15.4000, delta t in [0.03789902, 0.03789916], steps=6 (6) |
---|
| 848 | Time = 15.6000, delta t in [0.03789896, 0.03789908], steps=6 (6) |
---|
| 849 | Time = 15.8000, delta t in [0.03789891, 0.03789903], steps=6 (6) |
---|
| 850 | Stage > 0: Changing to outflow boundary |
---|
| 851 | Time = 16.0000, delta t in [0.02709050, 0.03789898], steps=6 (6) |
---|
| 852 | Time = 16.2000, delta t in [0.03789892, 0.03789904], steps=6 (6) |
---|
| 853 | ... |
---|
| 854 | \end{verbatim}} |
---|
| 855 | |
---|
| 856 | |
---|
| 857 | \subsection{Flow through more complex topograhies} |
---|
| 858 | |
---|
| 859 | Here is the code for the third version of the channel flow \file{channel3.py}: |
---|
| 860 | \verbatiminput{demos/channel3.py} |
---|
| 861 | |
---|
[4209] | 862 | This example differs from the first two versions in that the topography |
---|
[4123] | 863 | contains obstacles. |
---|
| 864 | |
---|
| 865 | This is accomplished here by defining the function \code{topography} as follows |
---|
| 866 | {\small \begin{verbatim} |
---|
| 867 | def topography(x,y): |
---|
| 868 | """Complex topography defined by a function of vectors x and y |
---|
| 869 | """ |
---|
| 870 | |
---|
[4209] | 871 | z = -x/10 |
---|
| 872 | |
---|
[4123] | 873 | N = len(x) |
---|
| 874 | for i in range(N): |
---|
| 875 | |
---|
| 876 | # Step |
---|
| 877 | if 10 < x[i] < 12: |
---|
[4209] | 878 | z[i] += 0.4 - 0.05*y[i] |
---|
| 879 | |
---|
[4123] | 880 | # Constriction |
---|
| 881 | if 27 < x[i] < 29 and y[i] > 3: |
---|
[4209] | 882 | z[i] += 2 |
---|
| 883 | |
---|
[4123] | 884 | # Pole |
---|
| 885 | if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2: |
---|
| 886 | z[i] += 2 |
---|
| 887 | |
---|
| 888 | return z |
---|
| 889 | \end{verbatim}} |
---|
| 890 | |
---|
| 891 | In addition, changing the resolution to dx=dy=0.1 creates a finer mesh resolving the new featurse better. |
---|
| 892 | |
---|
| 893 | A screenshot of this model at time == 15s is |
---|
| 894 | \begin{figure}[hbt] |
---|
| 895 | |
---|
| 896 | \centerline{\includegraphics[height=75mm] |
---|
| 897 | {graphics/channel3.png}} |
---|
| 898 | |
---|
| 899 | \caption{More complex flow in a channel} |
---|
| 900 | \label{fig:channel3} |
---|
| 901 | \end{figure} |
---|
| 902 | |
---|
| 903 | |
---|
| 904 | |
---|
| 905 | |
---|
| 906 | \clearpage |
---|
| 907 | |
---|
| 908 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 909 | |
---|
| 910 | \section{An Example with Real Data} |
---|
| 911 | \label{sec:realdataexample} The following discussion builds on the |
---|
| 912 | concepts introduced through the \file{runup.py} example and |
---|
| 913 | introduces a second example, \file{runcairns.py}. This refers to |
---|
[4953] | 914 | a {\bf hypothetical} scenario using real-life data, |
---|
[4875] | 915 | in which the domain of interest surrounds the |
---|
[4123] | 916 | Cairns region. Two scenarios are given; firstly, a |
---|
| 917 | hypothetical tsunami wave is generated by a submarine mass failure |
---|
| 918 | situated on the edge of the continental shelf, and secondly, a fixed wave |
---|
| 919 | of given amplitude and period is introduced through the boundary. |
---|
| 920 | |
---|
[4875] | 921 | {\bf |
---|
| 922 | Each scenario has been designed to generate a tsunami which will |
---|
[4953] | 923 | inundate the Cairns region. To achieve this, suitably large |
---|
[4875] | 924 | parameters were chosen and were not based on any known tsunami sources |
---|
| 925 | or realistic amplitudes. |
---|
| 926 | } |
---|
| 927 | |
---|
[4123] | 928 | \subsection{Overview} |
---|
| 929 | As in the case of \file{runup.py}, the actions carried |
---|
| 930 | out by the program can be organised according to this outline: |
---|
| 931 | |
---|
| 932 | \begin{enumerate} |
---|
| 933 | |
---|
| 934 | \item Set up a triangular mesh. |
---|
| 935 | |
---|
| 936 | \item Set certain parameters governing the mode of |
---|
| 937 | operation of the model---specifying, for instance, where to store the |
---|
| 938 | model output. |
---|
| 939 | |
---|
| 940 | \item Input various quantities describing physical measurements, such |
---|
| 941 | as the elevation, to be specified at each mesh point (vertex). |
---|
| 942 | |
---|
| 943 | \item Set up the boundary conditions. |
---|
| 944 | |
---|
| 945 | \item Carry out the evolution of the model through a series of time |
---|
| 946 | steps and output the results, providing a results file that can be |
---|
| 947 | visualised. |
---|
| 948 | |
---|
| 949 | \end{enumerate} |
---|
| 950 | |
---|
| 951 | |
---|
| 952 | |
---|
| 953 | \subsection{The Code} |
---|
| 954 | |
---|
| 955 | Here is the code for \file{runcairns.py}: |
---|
| 956 | |
---|
| 957 | \verbatiminput{demos/cairns/runcairns.py} |
---|
| 958 | |
---|
| 959 | In discussing the details of this example, we follow the outline |
---|
| 960 | given above, discussing each major step of the code in turn. |
---|
| 961 | |
---|
| 962 | \subsection{Establishing the Mesh}\index{mesh, establishing} |
---|
| 963 | |
---|
| 964 | One obvious way that the present example differs from |
---|
| 965 | \file{runup.py} is in the use of a more complex method to |
---|
| 966 | create the mesh. Instead of imposing a mesh structure on a |
---|
| 967 | rectangular grid, the technique used for this example involves |
---|
| 968 | building mesh structures inside polygons specified by the user, |
---|
[4736] | 969 | using a mesh-generator. |
---|
[4123] | 970 | |
---|
[4736] | 971 | In its simplest form, the mesh-generator creates the mesh within a single |
---|
[4123] | 972 | polygon whose vertices are at geographical locations specified by |
---|
| 973 | the user. The user specifies the \emph{resolution}---that is, the |
---|
| 974 | maximal area of a triangle used for triangulation---and a triangular |
---|
| 975 | mesh is created inside the polygon using a mesh generation engine. |
---|
[4209] | 976 | On any given platform, the same mesh will be returned. |
---|
[4123] | 977 | %Figure |
---|
| 978 | %\ref{fig:pentagon} shows a simple example of this, in which the |
---|
| 979 | %triangulation is carried out within a pentagon. |
---|
| 980 | |
---|
| 981 | |
---|
| 982 | %\begin{figure}[hbt] |
---|
| 983 | |
---|
| 984 | % \caption{Mesh points are created inside the polygon} |
---|
| 985 | %\label{fig:pentagon} |
---|
| 986 | %\end{figure} |
---|
| 987 | |
---|
| 988 | Boundary tags are not restricted to \code{`left'}, \code{`bottom'}, |
---|
| 989 | \code{`right'} and \code{`top'}, as in the case of |
---|
| 990 | \file{runup.py}. Instead the user specifies a list of |
---|
| 991 | tags appropriate to the configuration being modelled. |
---|
| 992 | |
---|
[4736] | 993 | In addition, the mesh-generator provides a way to adapt to geographic or |
---|
[4123] | 994 | other features in the landscape, whose presence may require an |
---|
| 995 | increase in resolution. This is done by allowing the user to specify |
---|
| 996 | a number of \emph{interior polygons}, each with a specified |
---|
| 997 | resolution. It is also |
---|
| 998 | possible to specify one or more `holes'---that is, areas bounded by |
---|
| 999 | polygons in which no triangulation is required. |
---|
| 1000 | |
---|
| 1001 | %\begin{figure}[hbt] |
---|
| 1002 | % \caption{Interior meshes with individual resolution} |
---|
| 1003 | % \label{fig:interior meshes} |
---|
| 1004 | %\end{figure} |
---|
| 1005 | |
---|
[4736] | 1006 | In its general form, the mesh-generator takes for its input a bounding |
---|
[4123] | 1007 | polygon and (optionally) a list of interior polygons. The user |
---|
| 1008 | specifies resolutions, both for the bounding polygon and for each of |
---|
[4736] | 1009 | the interior polygons. Given this data, the mesh-generator first creates a |
---|
[4123] | 1010 | triangular mesh with varying resolution. |
---|
| 1011 | |
---|
| 1012 | The function used to implement this process is |
---|
| 1013 | \function{create\_mesh\_from\_regions}. Its arguments include the |
---|
| 1014 | bounding polygon and its resolution, a list of boundary tags, and a |
---|
| 1015 | list of pairs \code{[polygon, resolution]}, specifying the interior |
---|
| 1016 | polygons and their resolutions. |
---|
| 1017 | |
---|
| 1018 | The resulting mesh is output to a \emph{mesh file}\index{mesh |
---|
| 1019 | file}\label{def:mesh file}. This term is used to describe a file of |
---|
| 1020 | a specific format used to store the data specifying a mesh. (There |
---|
| 1021 | are in fact two possible formats for such a file: it can either be a |
---|
| 1022 | binary file, with extension \code{.msh}, or an ASCII file, with |
---|
| 1023 | extension \code{.tsh}. In the present case, the binary file format |
---|
| 1024 | \code{.msh} is used. See Section \ref{sec:file formats} (page |
---|
| 1025 | \pageref{sec:file formats}) for more on file formats.) |
---|
| 1026 | |
---|
| 1027 | In practice, the details of the polygons used are read from a |
---|
| 1028 | separate file \file{project.py}. Here is a complete listing of |
---|
| 1029 | \file{project.py}: |
---|
| 1030 | |
---|
| 1031 | \verbatiminput{demos/cairns/project.py} |
---|
| 1032 | |
---|
| 1033 | Figure \ref{fig:cairns3d} illustrates the landscape of the region |
---|
| 1034 | for the Cairns example. Understanding the landscape is important in |
---|
[4209] | 1035 | determining the location and resolution of interior polygons. The |
---|
[4123] | 1036 | supporting data is found in the ASCII grid, \code{cairns.asc}, which |
---|
| 1037 | has been sourced from the publically available Australian Bathymetry |
---|
| 1038 | and Topography Grid 2005, \cite{grid250}. The required resolution |
---|
[4209] | 1039 | for inundation modelling will depend on the underlying topography and |
---|
[4123] | 1040 | bathymetry; as the terrain becomes more complex, the desired resolution |
---|
| 1041 | would decrease to the order of tens of metres. |
---|
| 1042 | |
---|
| 1043 | \begin{figure}[hbt] |
---|
| 1044 | \centerline{\includegraphics[scale=0.5]{graphics/cairns3.jpg}} |
---|
| 1045 | \caption{Landscape of the Cairns scenario.} |
---|
| 1046 | \label{fig:cairns3d} |
---|
| 1047 | |
---|
| 1048 | \end{figure} |
---|
| 1049 | The following statements are used to read in the specific polygons |
---|
| 1050 | from \code{project.cairns} and assign a defined resolution to |
---|
| 1051 | each polygon. |
---|
| 1052 | |
---|
| 1053 | {\small \begin{verbatim} |
---|
| 1054 | islands_res = 100000 |
---|
| 1055 | cairns_res = 100000 |
---|
| 1056 | shallow_res = 500000 |
---|
| 1057 | interior_regions = [[project.poly_cairns, cairns_res], |
---|
| 1058 | [project.poly_island0, islands_res], |
---|
| 1059 | [project.poly_island1, islands_res], |
---|
| 1060 | [project.poly_island2, islands_res], |
---|
| 1061 | [project.poly_island3, islands_res], |
---|
| 1062 | [project.poly_shallow, shallow_res]] |
---|
| 1063 | \end{verbatim}} |
---|
| 1064 | |
---|
| 1065 | Figure \ref{fig:cairnspolys} |
---|
[4209] | 1066 | illustrates the polygons used for the Cairns scenario. |
---|
[4123] | 1067 | |
---|
| 1068 | \begin{figure}[hbt] |
---|
| 1069 | |
---|
| 1070 | \centerline{\includegraphics[scale=0.5] |
---|
| 1071 | {graphics/cairnsmodel.jpg}} |
---|
| 1072 | \caption{Interior and bounding polygons for the Cairns example.} |
---|
| 1073 | \label{fig:cairnspolys} |
---|
| 1074 | \end{figure} |
---|
| 1075 | |
---|
| 1076 | The statement |
---|
| 1077 | |
---|
| 1078 | |
---|
| 1079 | {\small \begin{verbatim} |
---|
[4209] | 1080 | remainder_res = 10000000 |
---|
[4123] | 1081 | create_mesh_from_regions(project.bounding_polygon, |
---|
| 1082 | boundary_tags={'top': [0], |
---|
| 1083 | 'ocean_east': [1], |
---|
| 1084 | 'bottom': [2], |
---|
| 1085 | 'onshore': [3]}, |
---|
| 1086 | maximum_triangle_area=remainder_res, |
---|
| 1087 | filename=meshname, |
---|
| 1088 | interior_regions=interior_regions, |
---|
| 1089 | use_cache=True, |
---|
| 1090 | verbose=True) |
---|
| 1091 | \end{verbatim}} |
---|
| 1092 | is then used to create the mesh, taking the bounding polygon to be |
---|
| 1093 | the polygon \code{bounding\_polygon} specified in \file{project.py}. |
---|
| 1094 | The argument \code{boundary\_tags} assigns a dictionary, whose keys |
---|
| 1095 | are the names of the boundary tags used for the bounding |
---|
| 1096 | polygon---\code{`top'}, \code{`ocean\_east'}, \code{`bottom'}, and |
---|
| 1097 | \code{`onshore'}--- and whose values identify the indices of the |
---|
[4953] | 1098 | segments associated with each of these tags. |
---|
| 1099 | The polygon may be arranged either clock-wise or counter clock-wise and the |
---|
[4818] | 1100 | indices refer to edges in the order they appear: Edge 0 connects vertex 0 and vertex 1, edge 1 connects vertex 1 and 2; and so forth. |
---|
| 1101 | (Here, the values associated with each boundary tag are one-element lists, but they can have as many indices as there are edges) |
---|
[4673] | 1102 | If polygons intersect, or edges coincide the resolution may be undefined in some regions. |
---|
[4953] | 1103 | Use the underlying mesh interface for such cases. See Section |
---|
[4674] | 1104 | \ref{sec:mesh interface}. |
---|
[4123] | 1105 | |
---|
[4953] | 1106 | Note that every point on each polygon defining the mesh will be used as vertices in triangles. |
---|
| 1107 | Consequently, polygons with points very close together will cause triangles with very small |
---|
[4691] | 1108 | areas to be generated irrespective of the requested resolution. |
---|
[4953] | 1109 | Make sure points on polygons are spaced to be no closer than the smallest resolution requested. |
---|
[4123] | 1110 | |
---|
| 1111 | |
---|
| 1112 | \subsection{Initialising the Domain} |
---|
| 1113 | |
---|
| 1114 | As with \file{runup.py}, once we have created the mesh, the next |
---|
| 1115 | step is to create the data structure \code{domain}. We did this for |
---|
| 1116 | \file{runup.py} by inputting lists of points and triangles and |
---|
| 1117 | specifying the boundary tags directly. However, in the present case, |
---|
| 1118 | we use a method that works directly with the mesh file |
---|
| 1119 | \code{meshname}, as follows: |
---|
| 1120 | |
---|
| 1121 | |
---|
| 1122 | {\small \begin{verbatim} |
---|
| 1123 | domain = Domain(meshname, use_cache=True, verbose=True) |
---|
| 1124 | \end{verbatim}} |
---|
| 1125 | |
---|
| 1126 | Providing a filename instead of the lists used in \file{runup.py} |
---|
| 1127 | above causes \code{Domain} to convert a mesh file \code{meshname} |
---|
| 1128 | into an instance of \code{Domain}, allowing us to use methods like |
---|
| 1129 | \method{set\_quantity} to set quantities and to apply other |
---|
| 1130 | operations. |
---|
| 1131 | |
---|
| 1132 | %(In principle, the |
---|
| 1133 | %second argument of \function{pmesh\_to\_domain\_instance} can be any |
---|
| 1134 | %subclass of \class{Domain}, but for applications involving the |
---|
| 1135 | %shallow-water wave equation, the second argument of |
---|
| 1136 | %\function{pmesh\_to\_domain\_instance} can always be set simply to |
---|
| 1137 | %\class{Domain}.) |
---|
| 1138 | |
---|
| 1139 | The following statements specify a basename and data directory, and |
---|
| 1140 | identify quantities to be stored. For the first two, values are |
---|
| 1141 | taken from \file{project.py}. |
---|
| 1142 | |
---|
| 1143 | {\small \begin{verbatim} |
---|
| 1144 | domain.set_name(project.basename) |
---|
| 1145 | domain.set_datadir(project.outputdir) |
---|
| 1146 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', |
---|
| 1147 | 'ymomentum']) |
---|
| 1148 | \end{verbatim}} |
---|
| 1149 | |
---|
| 1150 | |
---|
| 1151 | \subsection{Initial Conditions} |
---|
| 1152 | Quantities for \file{runcairns.py} are set |
---|
| 1153 | using similar methods to those in \file{runup.py}. However, |
---|
| 1154 | in this case, many of the values are read from the auxiliary file |
---|
| 1155 | \file{project.py} or, in the case of \code{elevation}, from an |
---|
| 1156 | ancillary points file. |
---|
| 1157 | |
---|
| 1158 | |
---|
| 1159 | |
---|
| 1160 | \subsubsection{Stage} |
---|
| 1161 | |
---|
| 1162 | For the scenario we are modelling in this case, we use a callable |
---|
| 1163 | object \code{tsunami\_source}, assigned by means of a function |
---|
| 1164 | \function{slide\_tsunami}. This is similar to how we set elevation in |
---|
| 1165 | \file{runup.py} using a function---however, in this case the |
---|
| 1166 | function is both more complex and more interesting. |
---|
| 1167 | |
---|
| 1168 | The function returns the water displacement for all \code{x} and |
---|
| 1169 | \code{y} in the domain. The water displacement is a double Gaussian |
---|
| 1170 | function that depends on the characteristics of the slide (length, |
---|
| 1171 | width, thickness, slope, etc), its location (origin) and the depth at that |
---|
[4209] | 1172 | location. For this example, we choose to apply the slide function |
---|
[4875] | 1173 | at a specified time into the simulation. {\bf Note, the parameters used |
---|
[4953] | 1174 | in this example have been deliberately chosen to generate a suitably |
---|
[4875] | 1175 | large amplitude tsunami which would inundate the Cairns region.} |
---|
[4123] | 1176 | |
---|
| 1177 | \subsubsection{Friction} |
---|
| 1178 | |
---|
| 1179 | We assign the friction exactly as we did for \file{runup.py}: |
---|
| 1180 | |
---|
| 1181 | {\small \begin{verbatim} |
---|
| 1182 | domain.set_quantity('friction', 0.0) |
---|
| 1183 | \end{verbatim}} |
---|
| 1184 | |
---|
| 1185 | |
---|
| 1186 | \subsubsection{Elevation} |
---|
| 1187 | |
---|
| 1188 | The elevation is specified by reading data from a file: |
---|
| 1189 | |
---|
| 1190 | {\small \begin{verbatim} |
---|
| 1191 | domain.set_quantity('elevation', |
---|
| 1192 | filename = project.dem_name + '.pts', |
---|
| 1193 | use_cache = True, |
---|
| 1194 | verbose = True) |
---|
| 1195 | \end{verbatim}} |
---|
| 1196 | |
---|
| 1197 | %However, before this step can be executed, some preliminary steps |
---|
| 1198 | %are needed to prepare the file from which the data is taken. Two |
---|
| 1199 | %source files are used for this data---their names are specified in |
---|
| 1200 | %the file \file{project.py}, in the variables \code{coarsedemname} |
---|
| 1201 | %and \code{finedemname}. They contain `coarse' and `fine' data, |
---|
| 1202 | %respectively---that is, data sampled at widely spaced points over a |
---|
| 1203 | %large region and data sampled at closely spaced points over a |
---|
| 1204 | %smaller subregion. The data in these files is combined through the |
---|
| 1205 | %statement |
---|
| 1206 | |
---|
| 1207 | %{\small \begin{verbatim} |
---|
| 1208 | %combine_rectangular_points_files(project.finedemname + '.pts', |
---|
| 1209 | % project.coarsedemname + '.pts', |
---|
| 1210 | % project.combineddemname + '.pts') |
---|
| 1211 | %\end{verbatim}} |
---|
| 1212 | %The effect of this is simply to combine the datasets by eliminating |
---|
| 1213 | %any coarse data associated with points inside the smaller region |
---|
| 1214 | %common to both datasets. The name to be assigned to the resulting |
---|
| 1215 | %dataset is also derived from the name stored in the variable |
---|
| 1216 | %\code{combinedname} in the file \file{project.py}. |
---|
| 1217 | |
---|
| 1218 | \subsection{Boundary Conditions}\index{boundary conditions} |
---|
| 1219 | |
---|
| 1220 | Setting boundaries follows a similar pattern to the one used for |
---|
| 1221 | \file{runup.py}, except that in this case we need to associate a |
---|
| 1222 | boundary type with each of the |
---|
| 1223 | boundary tag names introduced when we established the mesh. In place of the four |
---|
| 1224 | boundary types introduced for \file{runup.py}, we use the reflective |
---|
| 1225 | boundary for each of the |
---|
| 1226 | eight tagged segments defined by \code{create_mesh_from_regions}: |
---|
| 1227 | |
---|
| 1228 | {\small \begin{verbatim} |
---|
| 1229 | Bd = Dirichlet_boundary([0.0,0.0,0.0]) |
---|
| 1230 | domain.set_boundary( {'ocean_east': Bd, 'bottom': Bd, 'onshore': Bd, |
---|
| 1231 | 'top': Bd} ) |
---|
| 1232 | \end{verbatim}} |
---|
| 1233 | |
---|
| 1234 | \subsection{Evolution} |
---|
| 1235 | |
---|
| 1236 | With the basics established, the running of the `evolve' step is |
---|
| 1237 | very similar to the corresponding step in \file{runup.py}. For the slide |
---|
[4209] | 1238 | scenario, |
---|
[4123] | 1239 | the simulation is run for 5000 seconds with the output stored every ten seconds. |
---|
| 1240 | For this example, we choose to apply the slide at 60 seconds into the simulation. |
---|
| 1241 | |
---|
| 1242 | {\small \begin{verbatim} |
---|
| 1243 | import time t0 = time.time() |
---|
| 1244 | |
---|
[4209] | 1245 | |
---|
| 1246 | for t in domain.evolve(yieldstep = 10, finaltime = 60): |
---|
[4123] | 1247 | domain.write_time() |
---|
[4209] | 1248 | domain.write_boundary_statistics(tags = 'ocean_east') |
---|
| 1249 | |
---|
[4123] | 1250 | # add slide |
---|
| 1251 | thisstagestep = domain.get_quantity('stage') |
---|
| 1252 | if allclose(t, 60): |
---|
| 1253 | slide = Quantity(domain) |
---|
| 1254 | slide.set_values(tsunami_source) |
---|
| 1255 | domain.set_quantity('stage', slide + thisstagestep) |
---|
[4209] | 1256 | |
---|
| 1257 | for t in domain.evolve(yieldstep = 10, finaltime = 5000, |
---|
[4123] | 1258 | skip_initial_step = True): |
---|
| 1259 | domain.write_time() |
---|
| 1260 | domain.write_boundary_statistics(tags = 'ocean_east') |
---|
| 1261 | \end{verbatim}} |
---|
| 1262 | |
---|
[4209] | 1263 | For the fixed wave scenario, the simulation is run to 10000 seconds, |
---|
[4123] | 1264 | with the first half of the simulation stored at two minute intervals, |
---|
| 1265 | and the second half of the simulation stored at ten second intervals. |
---|
| 1266 | This functionality is especially convenient as it allows the detailed |
---|
| 1267 | parts of the simulation to be viewed at higher time resolution. |
---|
| 1268 | |
---|
| 1269 | |
---|
| 1270 | {\small \begin{verbatim} |
---|
| 1271 | |
---|
| 1272 | # save every two mins leading up to wave approaching land |
---|
[4209] | 1273 | for t in domain.evolve(yieldstep = 120, finaltime = 5000): |
---|
[4123] | 1274 | domain.write_time() |
---|
[4209] | 1275 | domain.write_boundary_statistics(tags = 'ocean_east') |
---|
[4123] | 1276 | |
---|
| 1277 | # save every 30 secs as wave starts inundating ashore |
---|
[4209] | 1278 | for t in domain.evolve(yieldstep = 10, finaltime = 10000, |
---|
[4123] | 1279 | skip_initial_step = True): |
---|
| 1280 | domain.write_time() |
---|
| 1281 | domain.write_boundary_statistics(tags = 'ocean_east') |
---|
[4209] | 1282 | |
---|
[4123] | 1283 | \end{verbatim}} |
---|
| 1284 | |
---|
| 1285 | \section{Exploring the Model Output} |
---|
| 1286 | |
---|
| 1287 | Now that the scenario has been run, the user can view the output in a number of ways. |
---|
| 1288 | As described earlier, the user may run animate to view a three-dimensional representation |
---|
| 1289 | of the simulation. |
---|
| 1290 | |
---|
| 1291 | The user may also be interested in a maximum inundation map. This simply shows the |
---|
[4209] | 1292 | maximum water depth over the domain and is achieved with the function sww2dem (described in |
---|
[4207] | 1293 | Section \ref{sec:basicfileconversions}). |
---|
[4123] | 1294 | \file{ExportResults.py} demonstrates how this function can be used: |
---|
| 1295 | |
---|
| 1296 | \verbatiminput{demos/cairns/ExportResults.py} |
---|
| 1297 | |
---|
| 1298 | The script generates an maximum water depth ASCII grid at a defined |
---|
| 1299 | resolution (here 100 m$^2$) which can then be viewed in a GIS environment, for |
---|
| 1300 | example. The parameters used in the function are defined in \file{project.py}. |
---|
| 1301 | Figures \ref{fig:maxdepthcairnsslide} and \ref{fig:maxdepthcairnsfixedwave} show |
---|
[4209] | 1302 | the maximum water depth within the defined region for the slide and fixed wave scenario |
---|
[4875] | 1303 | respectively. {\bf Note, these inundation maps have been based on purely hypothetical |
---|
| 1304 | scenarios and were designed explicitly for demonstration purposes only.} |
---|
[4123] | 1305 | The user could develop a maximum absolute momentum or other expressions which can be |
---|
[4209] | 1306 | derived from the quantities. |
---|
[4869] | 1307 | It must be noted here that depth is more meaningful when the elevation is positive |
---|
[4870] | 1308 | (\code{depth} = \code{stage} $-$ \code{elevation}) as it describes the water height |
---|
[4869] | 1309 | above the available elevation. When the elevation is negative, depth is meauring the |
---|
[4953] | 1310 | water height from the sea floor. With this in mind, maximum inundation maps are |
---|
| 1311 | typically "clipped" to the coastline. However, the data input here did not contain a |
---|
[4869] | 1312 | coastline. |
---|
[4123] | 1313 | |
---|
| 1314 | \begin{figure}[hbt] |
---|
| 1315 | \centerline{\includegraphics[scale=0.5]{graphics/slidedepth.jpg}} |
---|
[4875] | 1316 | \caption{Maximum inundation map for the Cairns slide scenario. \bf Note, this |
---|
| 1317 | inundaiton map has been based on a purely hypothetical scenario which was |
---|
| 1318 | designed explictiy for demonstration purposes only.} |
---|
[4123] | 1319 | \label{fig:maxdepthcairnsslide} |
---|
| 1320 | \end{figure} |
---|
| 1321 | |
---|
| 1322 | \begin{figure}[hbt] |
---|
| 1323 | \centerline{\includegraphics[scale=0.5]{graphics/fixedwavedepth.jpg}} |
---|
[4953] | 1324 | \caption{Maximum inundation map for the Cairns fixed wave scenario. |
---|
[4875] | 1325 | \bf Note, this |
---|
| 1326 | inundaiton map has been based on a purely hypothetical scenario which was |
---|
| 1327 | designed explictiy for demonstration purposes only.} |
---|
[4123] | 1328 | \label{fig:maxdepthcairnsfixedwave} |
---|
| 1329 | \end{figure} |
---|
| 1330 | |
---|
| 1331 | The user may also be interested in interrogating the solution at a particular spatial |
---|
| 1332 | location to understand the behaviour of the system through time. To do this, the user |
---|
| 1333 | must first define the locations of interest. A number of locations have been |
---|
[4209] | 1334 | identified for the Cairns scenario, as shown in Figure \ref{fig:cairnsgauges}. |
---|
[4123] | 1335 | |
---|
| 1336 | \begin{figure}[hbt] |
---|
| 1337 | \centerline{\includegraphics[scale=0.5]{graphics/cairnsgauges.jpg}} |
---|
| 1338 | \caption{Point locations to show time series information for the Cairns scenario.} |
---|
| 1339 | \label{fig:cairnsgauges} |
---|
| 1340 | \end{figure} |
---|
| 1341 | |
---|
| 1342 | These locations |
---|
| 1343 | must be stored in either a .csv or .txt file. The corresponding .csv file for |
---|
| 1344 | the gauges shown in Figure \ref{fig:cairnsgauges} is \file{gauges.csv} |
---|
| 1345 | |
---|
[4871] | 1346 | \verbatiminput{demos/cairns/gauges.csv} |
---|
[4123] | 1347 | |
---|
| 1348 | Header information has been included to identify the location in terms of eastings and |
---|
| 1349 | northings, and each gauge is given a name. The elevation column can be zero here. |
---|
[4945] | 1350 | This information is then passed to the function \code{sww2csv_gauges} (shown in |
---|
[4953] | 1351 | \file{GetTimeseries.py} which generates the csv files for each point location. The csv files |
---|
| 1352 | can then be used in \code{csv2timeseries_graphs} to create the timeseries plot for each desired |
---|
[4949] | 1353 | quantity. \code{csv2timeseries_graphs} relies on \code{pylab} to be installed which is not part |
---|
[4945] | 1354 | of the standard \code{anuga} release, however it can be downloaded and installed from \code{http://matplotlib.sourceforge.net/} |
---|
[4123] | 1355 | |
---|
| 1356 | \verbatiminput{demos/cairns/GetTimeseries.py} |
---|
| 1357 | |
---|
[4945] | 1358 | Here, the time series for the quantities stage, depth and speed will be generated for |
---|
| 1359 | each gauge defined in the gauge file. As described earlier, depth is more meaningful |
---|
[4953] | 1360 | for onshore gauges, and stage is more appropriate for offshore gauges. |
---|
[4123] | 1361 | |
---|
[4869] | 1362 | As an example output, |
---|
[4953] | 1363 | Figure \ref{fig:reef} shows the time series for the quantity stage for the |
---|
| 1364 | Elford Reef location for each scenario (the elevation at this location is negative, |
---|
[4949] | 1365 | therefore stage is the more appropriate quantity to plot). Note the large negative stage value when the slide was |
---|
[4875] | 1366 | introduced. This is due to the double gaussian form of the initial surface |
---|
| 1367 | displacement of the slide. By contrast, the time series for depth is shown for the onshore location of the Cairns |
---|
[4945] | 1368 | Airport in Figure \ref{fig:airportboth}. |
---|
[4869] | 1369 | |
---|
[4123] | 1370 | \begin{figure}[hbt] |
---|
[4948] | 1371 | \centerline{\includegraphics[scale=0.5]{graphics/gaugeElfordReefstage.png}} |
---|
[4953] | 1372 | \caption{Time series information of the quantity stage for the Elford Reef location for the |
---|
[4869] | 1373 | fixed wave and slide scenario.} |
---|
[4123] | 1374 | \label{fig:reef} |
---|
| 1375 | \end{figure} |
---|
| 1376 | |
---|
| 1377 | \begin{figure}[hbt] |
---|
[4948] | 1378 | \centerline{\includegraphics[scale=0.5]{graphics/gaugeCairnsAirportdepth.png}} |
---|
[4953] | 1379 | \caption{Time series information of the quantity depth for the Cairns Airport |
---|
[4869] | 1380 | location for the slide and fixed wave scenario.} |
---|
| 1381 | \label{fig:airportboth} |
---|
[4123] | 1382 | \end{figure} |
---|
| 1383 | |
---|
| 1384 | |
---|
| 1385 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 1386 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 1387 | |
---|
| 1388 | \chapter{\anuga Public Interface} |
---|
| 1389 | \label{ch:interface} |
---|
| 1390 | |
---|
| 1391 | This chapter gives an overview of the features of \anuga available |
---|
| 1392 | to the user at the public interface. These are grouped under the |
---|
| 1393 | following headings, which correspond to the outline of the examples |
---|
| 1394 | described in Chapter \ref{ch:getstarted}: |
---|
| 1395 | |
---|
| 1396 | \begin{itemize} |
---|
[5508] | 1397 | \item Establishing the Mesh: Section \ref{sec:establishing the mesh} |
---|
| 1398 | \item Initialising the Domain: Section \ref{sec:initialising the domain} |
---|
| 1399 | \item Specifying the Quantities: Section \ref{sec:quantitis} |
---|
| 1400 | \item Initial Conditions: Section \ref{sec:initial conditions} |
---|
| 1401 | \item Boundary Conditions: Section \ref{sec:boundary conditions} |
---|
| 1402 | \item Forcing Terms: Section \ref{sec:forcing terms} |
---|
| 1403 | \item Evolution: Section \ref{sec:evolution} |
---|
[4123] | 1404 | \end{itemize} |
---|
| 1405 | |
---|
| 1406 | The listings are intended merely to give the reader an idea of what |
---|
| 1407 | each feature is, where to find it and how it can be used---they do |
---|
| 1408 | not give full specifications; for these the reader |
---|
| 1409 | may consult the code. The code for every function or class contains |
---|
| 1410 | a documentation string, or `docstring', that specifies the precise |
---|
| 1411 | syntax for its use. This appears immediately after the line |
---|
| 1412 | introducing the code, between two sets of triple quotes. |
---|
| 1413 | |
---|
| 1414 | Each listing also describes the location of the module in which |
---|
| 1415 | the code for the feature being described can be found. All modules |
---|
| 1416 | are in the folder \file{inundation} or one of its subfolders, and the |
---|
| 1417 | location of each module is described relative to \file{inundation}. Rather |
---|
| 1418 | than using pathnames, whose syntax depends on the operating system, |
---|
| 1419 | we use the format adopted for importing the function or class for |
---|
| 1420 | use in Python code. For example, suppose we wish to specify that the |
---|
| 1421 | function \function{create\_mesh\_from\_regions} is in a module called |
---|
| 1422 | \module{mesh\_interface} in a subfolder of \module{inundation} called |
---|
| 1423 | \code{pmesh}. In Linux or Unix syntax, the pathname of the file |
---|
| 1424 | containing the function, relative to \file{inundation}, would be |
---|
| 1425 | |
---|
| 1426 | \begin{center} |
---|
| 1427 | % \code{pmesh/mesh\_interface.py} |
---|
| 1428 | \code{pmesh}$\slash$\code{mesh\_interface.py} |
---|
| 1429 | \end{center} |
---|
[4674] | 1430 | \label{sec:mesh interface} |
---|
[4123] | 1431 | |
---|
| 1432 | while in Windows syntax it would be |
---|
| 1433 | |
---|
| 1434 | \begin{center} |
---|
| 1435 | \code{pmesh}$\backslash$\code{mesh\_interface.py} |
---|
| 1436 | \end{center} |
---|
| 1437 | |
---|
| 1438 | Rather than using either of these forms, in this chapter we specify |
---|
| 1439 | the location simply as \code{pmesh.mesh\_interface}, in keeping with |
---|
| 1440 | the usage in the Python statement for importing the function, |
---|
| 1441 | namely: |
---|
| 1442 | \begin{center} |
---|
| 1443 | \code{from pmesh.mesh\_interface import create\_mesh\_from\_regions} |
---|
| 1444 | \end{center} |
---|
| 1445 | |
---|
| 1446 | Each listing details the full set of parameters for the class or |
---|
| 1447 | function; however, the description is generally limited to the most |
---|
| 1448 | important parameters and the reader is again referred to the code |
---|
| 1449 | for more details. |
---|
| 1450 | |
---|
| 1451 | The following parameters are common to many functions and classes |
---|
| 1452 | and are omitted from the descriptions given below: |
---|
| 1453 | |
---|
| 1454 | %\begin{center} |
---|
| 1455 | \begin{tabular}{ll} %\hline |
---|
| 1456 | %\textbf{Name } & \textbf{Description}\\ |
---|
| 1457 | %\hline |
---|
| 1458 | \emph{use\_cache} & Specifies whether caching is to be used for improved performance. See Section \ref{sec:caching} for details on the underlying caching functionality\\ |
---|
| 1459 | \emph{verbose} & If \code{True}, provides detailed terminal output |
---|
| 1460 | to the user\\ % \hline |
---|
| 1461 | \end{tabular} |
---|
| 1462 | %\end{center} |
---|
| 1463 | |
---|
[5508] | 1464 | \section{Mesh Generation}\index{Mesh!generation} |
---|
| 1465 | \label{sec:establishing the mesh} |
---|
[4123] | 1466 | Before discussing the part of the interface relating to mesh |
---|
| 1467 | generation, we begin with a description of a simple example of a |
---|
| 1468 | mesh and use it to describe how mesh data is stored. |
---|
| 1469 | |
---|
| 1470 | \label{sec:meshexample} Figure \ref{fig:simplemesh} represents a |
---|
| 1471 | very simple mesh comprising just 11 points and 10 triangles. |
---|
| 1472 | |
---|
| 1473 | |
---|
| 1474 | \begin{figure}[h] |
---|
| 1475 | \begin{center} |
---|
| 1476 | \includegraphics[width=90mm, height=90mm]{triangularmesh.jpg} |
---|
| 1477 | \end{center} |
---|
| 1478 | |
---|
| 1479 | \caption{A simple mesh} |
---|
| 1480 | \label{fig:simplemesh} |
---|
| 1481 | \end{figure} |
---|
| 1482 | |
---|
| 1483 | |
---|
[5673] | 1484 | The variables \code{points}, \code{triangles} and \code{boundary} |
---|
[4123] | 1485 | represent the data displayed in Figure \ref{fig:simplemesh} as |
---|
| 1486 | follows. The list \code{points} stores the coordinates of the |
---|
| 1487 | points, and may be displayed schematically as in Table |
---|
| 1488 | \ref{tab:points}. |
---|
| 1489 | |
---|
| 1490 | |
---|
| 1491 | \begin{table} |
---|
| 1492 | \begin{center} |
---|
| 1493 | \begin{tabular}[t]{|c|cc|} \hline |
---|
| 1494 | index & \code{x} & \code{y}\\ \hline |
---|
| 1495 | 0 & 1 & 1\\ |
---|
| 1496 | 1 & 4 & 2\\ |
---|
| 1497 | 2 & 8 & 1\\ |
---|
| 1498 | 3 & 1 & 3\\ |
---|
| 1499 | 4 & 5 & 5\\ |
---|
| 1500 | 5 & 8 & 6\\ |
---|
| 1501 | 6 & 11 & 5\\ |
---|
| 1502 | 7 & 3 & 6\\ |
---|
| 1503 | 8 & 1 & 8\\ |
---|
| 1504 | 9 & 4 & 9\\ |
---|
| 1505 | 10 & 10 & 7\\ \hline |
---|
| 1506 | \end{tabular} |
---|
| 1507 | \end{center} |
---|
| 1508 | |
---|
| 1509 | \caption{Point coordinates for mesh in |
---|
| 1510 | Figure \protect \ref{fig:simplemesh}} |
---|
| 1511 | \label{tab:points} |
---|
| 1512 | \end{table} |
---|
| 1513 | |
---|
[5673] | 1514 | The list \code{triangles} specifies the triangles that make up the |
---|
[4123] | 1515 | mesh. It does this by specifying, for each triangle, the indices |
---|
| 1516 | (the numbers shown in the first column above) that correspond to the |
---|
[5673] | 1517 | three points at the triangles vertices, taken in an anti-clockwise order |
---|
[4123] | 1518 | around the triangle. Thus, in the example shown in Figure |
---|
[5673] | 1519 | \ref{fig:simplemesh}, the variable \code{triangles} contains the |
---|
| 1520 | entries shown in Table \ref{tab:triangles}. The starting point is |
---|
[4123] | 1521 | arbitrary so triangle $(0,1,3)$ is considered the same as $(1,3,0)$ |
---|
| 1522 | and $(3,0,1)$. |
---|
| 1523 | |
---|
| 1524 | |
---|
| 1525 | \begin{table} |
---|
| 1526 | \begin{center} |
---|
| 1527 | \begin{tabular}{|c|ccc|} \hline |
---|
[5673] | 1528 | index & \multicolumn{3}{c|}{\code{points}}\\ \hline |
---|
[4123] | 1529 | 0 & 0 & 1 & 3\\ |
---|
| 1530 | 1 & 1 & 2 & 4\\ |
---|
| 1531 | 2 & 2 & 5 & 4\\ |
---|
| 1532 | 3 & 2 & 6 & 5\\ |
---|
| 1533 | 4 & 4 & 5 & 9\\ |
---|
| 1534 | 5 & 4 & 9 & 7\\ |
---|
| 1535 | 6 & 3 & 4 & 7\\ |
---|
| 1536 | 7 & 7 & 9 & 8\\ |
---|
| 1537 | 8 & 1 & 4 & 3\\ |
---|
| 1538 | 9 & 5 & 10 & 9\\ \hline |
---|
| 1539 | \end{tabular} |
---|
| 1540 | \end{center} |
---|
| 1541 | |
---|
[5673] | 1542 | \caption{Triangles for mesh in Figure \protect \ref{fig:simplemesh}} |
---|
| 1543 | \label{tab:triangles} |
---|
[4123] | 1544 | \end{table} |
---|
| 1545 | |
---|
| 1546 | Finally, the variable \code{boundary} identifies the boundary |
---|
| 1547 | triangles and associates a tag with each. |
---|
| 1548 | |
---|
| 1549 | \refmodindex[pmesh.meshinterface]{pmesh.mesh\_interface}\label{sec:meshgeneration} |
---|
| 1550 | |
---|
| 1551 | \begin{funcdesc} {create\_mesh\_from\_regions}{bounding_polygon, |
---|
| 1552 | boundary_tags, |
---|
| 1553 | maximum_triangle_area, |
---|
| 1554 | filename=None, |
---|
| 1555 | interior_regions=None, |
---|
| 1556 | poly_geo_reference=None, |
---|
| 1557 | mesh_geo_reference=None, |
---|
| 1558 | minimum_triangle_angle=28.0} |
---|
| 1559 | Module: \module{pmesh.mesh\_interface} |
---|
| 1560 | |
---|
| 1561 | This function allows a user to initiate the automatic creation of a |
---|
| 1562 | mesh inside a specified polygon (input \code{bounding_polygon}). |
---|
| 1563 | Among the parameters that can be set are the \emph{resolution} |
---|
| 1564 | (maximal area for any triangle in the mesh) and the minimal angle |
---|
| 1565 | allowable in any triangle. The user can specify a number of internal |
---|
[4736] | 1566 | polygons within each of which the resolution of the mesh can be |
---|
| 1567 | specified. \code{interior_regions} is a paired list containing the |
---|
| 1568 | interior polygon and its resolution. Additionally, the user specifies |
---|
| 1569 | a list of boundary tags, one for each edge of the bounding polygon. |
---|
[4123] | 1570 | |
---|
| 1571 | \textbf{WARNING}. Note that the dictionary structure used for the |
---|
| 1572 | parameter \code{boundary\_tags} is different from that used for the |
---|
| 1573 | variable \code{boundary} that occurs in the specification of a mesh. |
---|
| 1574 | In the case of \code{boundary}, the tags are the \emph{values} of |
---|
| 1575 | the dictionary, whereas in the case of \code{boundary_tags}, the |
---|
| 1576 | tags are the \emph{keys} and the \emph{value} corresponding to a |
---|
| 1577 | particular tag is a list of numbers identifying boundary edges |
---|
| 1578 | labelled with that tag. Because of this, it is theoretically |
---|
| 1579 | possible to assign the same edge to more than one tag. However, an |
---|
| 1580 | attempt to do this will cause an error. |
---|
[4738] | 1581 | |
---|
| 1582 | \textbf{WARNING}. Do not have polygon lines cross or be on-top of each |
---|
| 1583 | other. This can result in regions of unspecified resolutions. Do |
---|
| 1584 | not have polygon close to each other. This can result in the area |
---|
| 1585 | between the polygons having small triangles. For more control |
---|
| 1586 | over the mesh outline use the methods described below. |
---|
[4953] | 1587 | |
---|
[4123] | 1588 | \end{funcdesc} |
---|
| 1589 | |
---|
| 1590 | |
---|
| 1591 | |
---|
| 1592 | \subsection{Advanced mesh generation} |
---|
| 1593 | |
---|
| 1594 | For more control over the creation of the mesh outline, use the |
---|
[4209] | 1595 | methods of the class \class{Mesh}. |
---|
[4123] | 1596 | |
---|
| 1597 | |
---|
| 1598 | \begin{classdesc} {Mesh}{userSegments=None, |
---|
| 1599 | userVertices=None, |
---|
| 1600 | holes=None, |
---|
| 1601 | regions=None} |
---|
| 1602 | Module: \module{pmesh.mesh} |
---|
| 1603 | |
---|
| 1604 | A class used to build a mesh outline and generate a two-dimensional |
---|
| 1605 | triangular mesh. The mesh outline is used to describe features on the |
---|
| 1606 | mesh, such as the mesh boundary. Many of this classes methods are used |
---|
| 1607 | to build a mesh outline, such as \code{add\_vertices} and |
---|
| 1608 | \code{add\_region\_from\_polygon}. |
---|
| 1609 | |
---|
| 1610 | \end{classdesc} |
---|
| 1611 | |
---|
| 1612 | |
---|
| 1613 | \subsubsection{Key Methods of Class Mesh} |
---|
| 1614 | |
---|
| 1615 | |
---|
| 1616 | \begin{methoddesc} {add\_hole}{x,y} |
---|
| 1617 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1618 | |
---|
| 1619 | This method is used to build the mesh outline. It defines a hole, |
---|
| 1620 | when the boundary of the hole has already been defined, by selecting a |
---|
[4209] | 1621 | point within the boundary. |
---|
[4123] | 1622 | |
---|
| 1623 | \end{methoddesc} |
---|
| 1624 | |
---|
| 1625 | |
---|
| 1626 | \begin{methoddesc} {add\_hole\_from\_polygon}{self, polygon, tags=None} |
---|
| 1627 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1628 | |
---|
| 1629 | This method is used to add a `hole' within a region ---that is, to |
---|
| 1630 | define a interior region where the triangular mesh will not be |
---|
| 1631 | generated---to a \class{Mesh} instance. The region boundary is described by |
---|
| 1632 | the polygon passed in. Additionally, the user specifies a list of |
---|
| 1633 | boundary tags, one for each edge of the bounding polygon. |
---|
| 1634 | \end{methoddesc} |
---|
| 1635 | |
---|
| 1636 | |
---|
| 1637 | \begin{methoddesc} {add\_points_and_segments}{self, points, segments, |
---|
| 1638 | segment\_tags=None} |
---|
| 1639 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1640 | |
---|
| 1641 | This method is used to build the mesh outline. It adds points and |
---|
[5673] | 1642 | segments connecting the points. Points is a list of points. Segments |
---|
| 1643 | is a list of segments. Each segment is defined by the start and end |
---|
| 1644 | of the line by it's point index, e.g. use \code{segments = |
---|
| 1645 | [[0,1],[1,2]]} to make a polyline between points 0, 1 and 2. A tag for |
---|
| 1646 | each segment can optionally be added. |
---|
[4123] | 1647 | |
---|
| 1648 | \end{methoddesc} |
---|
| 1649 | |
---|
| 1650 | \begin{methoddesc} {add\_region}{x,y} |
---|
| 1651 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1652 | |
---|
| 1653 | This method is used to build the mesh outline. It defines a region, |
---|
| 1654 | when the boundary of the region has already been defined, by selecting |
---|
| 1655 | a point within the boundary. A region instance is returned. This can |
---|
| 1656 | be used to set the resolution. |
---|
| 1657 | |
---|
| 1658 | \end{methoddesc} |
---|
| 1659 | |
---|
[4953] | 1660 | \begin{methoddesc} {add\_region\_from\_polygon}{self, polygon, |
---|
[4888] | 1661 | segment_tags=None, region_tag=None |
---|
[4123] | 1662 | max_triangle_area=None} |
---|
| 1663 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1664 | |
---|
| 1665 | This method is used to build the mesh outline. It adds a region to a |
---|
| 1666 | \class{Mesh} instance. Regions are commonly used to describe an area |
---|
| 1667 | with an increased density of triangles, by setting |
---|
| 1668 | \code{max_triangle_area}. The |
---|
| 1669 | region boundary is described by the input \code{polygon}. Additionally, the |
---|
| 1670 | user specifies a list of segment tags, one for each edge of the |
---|
[4888] | 1671 | bounding polygon. The regional tag is set using \code{region}. |
---|
[4123] | 1672 | |
---|
| 1673 | \end{methoddesc} |
---|
| 1674 | |
---|
| 1675 | |
---|
| 1676 | |
---|
| 1677 | |
---|
| 1678 | |
---|
| 1679 | \begin{methoddesc} {add\_vertices}{point_data} |
---|
| 1680 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1681 | |
---|
| 1682 | Add user vertices. The point_data can be a list of (x,y) values, a numeric |
---|
[4209] | 1683 | array or a geospatial_data instance. |
---|
[4123] | 1684 | \end{methoddesc} |
---|
| 1685 | |
---|
| 1686 | \begin{methoddesc} {auto\_segment}{raw_boundary=raw_boundary, |
---|
| 1687 | remove_holes=remove_holes, |
---|
| 1688 | smooth_indents=smooth_indents, |
---|
| 1689 | expand_pinch=expand_pinch} |
---|
| 1690 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1691 | |
---|
| 1692 | Add segments between some of the user vertices to give the vertices an |
---|
| 1693 | outline. The outline is an alpha shape. This method is |
---|
| 1694 | useful since a set of user vertices need to be outlined by segments |
---|
| 1695 | before generate_mesh is called. |
---|
[4209] | 1696 | |
---|
[4123] | 1697 | \end{methoddesc} |
---|
| 1698 | |
---|
| 1699 | \begin{methoddesc} {export\_mesh_file}{self,ofile} |
---|
| 1700 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1701 | |
---|
| 1702 | This method is used to save the mesh to a file. \code{ofile} is the |
---|
| 1703 | name of the mesh file to be written, including the extension. Use |
---|
| 1704 | the extension \code{.msh} for the file to be in NetCDF format and |
---|
| 1705 | \code{.tsh} for the file to be ASCII format. |
---|
| 1706 | \end{methoddesc} |
---|
| 1707 | |
---|
| 1708 | \begin{methoddesc} {generate\_mesh}{self, |
---|
| 1709 | maximum_triangle_area=None, |
---|
| 1710 | minimum_triangle_angle=28.0, |
---|
| 1711 | verbose=False} |
---|
| 1712 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1713 | |
---|
| 1714 | This method is used to generate the triangular mesh. The maximal |
---|
| 1715 | area of any triangle in the mesh can be specified, which is used to |
---|
| 1716 | control the triangle density, along with the |
---|
| 1717 | minimum angle in any triangle. |
---|
| 1718 | \end{methoddesc} |
---|
| 1719 | |
---|
| 1720 | |
---|
| 1721 | |
---|
[5744] | 1722 | \begin{methoddesc} {import_ungenerate_file}{self,ofile, tag=None, |
---|
[5207] | 1723 | region_tag=None} |
---|
[4123] | 1724 | Module: \module{pmesh.mesh}, Class: \class{Mesh} |
---|
| 1725 | |
---|
[5207] | 1726 | This method is used to import a polygon file in the ungenerate format, |
---|
| 1727 | which is used by arcGIS. The polygons from the file are converted to |
---|
[4123] | 1728 | vertices and segments. \code{ofile} is the name of the polygon file. |
---|
| 1729 | \code{tag} is the tag given to all the polygon's segments. |
---|
[5207] | 1730 | \code{region_tag} is the tag given to all the polygon's segments. If |
---|
| 1731 | it is a string the one value will be assigned to all regions. If it |
---|
| 1732 | is a list the first value in the list will be applied to the first |
---|
[5650] | 1733 | polygon etc. If \code{tag} is not given a value it defaults to None, |
---|
| 1734 | which means the segement will not effect the water flow, it will only |
---|
| 1735 | effect the mesh generation. |
---|
[4123] | 1736 | |
---|
| 1737 | This function can be used to import building footprints. |
---|
| 1738 | \end{methoddesc} |
---|
| 1739 | |
---|
| 1740 | %%%%%% |
---|
[5508] | 1741 | \section{Initialising the Domain}\index{Initialising the Domain} |
---|
| 1742 | \label{sec:initialising the domain} |
---|
[4123] | 1743 | |
---|
| 1744 | %Include description of the class Domain and the module domain. |
---|
| 1745 | |
---|
| 1746 | %FIXME (Ole): This is also defined in a later chapter |
---|
| 1747 | %\declaremodule{standard}{...domain} |
---|
| 1748 | |
---|
| 1749 | \begin{classdesc} {Domain} {source=None, |
---|
| 1750 | triangles=None, |
---|
| 1751 | boundary=None, |
---|
| 1752 | conserved_quantities=None, |
---|
| 1753 | other_quantities=None, |
---|
| 1754 | tagged_elements=None, |
---|
| 1755 | use_inscribed_circle=False, |
---|
| 1756 | mesh_filename=None, |
---|
| 1757 | use_cache=False, |
---|
| 1758 | verbose=False, |
---|
| 1759 | full_send_dict=None, |
---|
| 1760 | ghost_recv_dict=None, |
---|
| 1761 | processor=0, |
---|
| 1762 | numproc=1} |
---|
| 1763 | Module: \refmodule{abstract_2d_finite_volumes.domain} |
---|
| 1764 | |
---|
| 1765 | This class is used to create an instance of a data structure used to |
---|
| 1766 | store and manipulate data associated with a mesh. The mesh is |
---|
| 1767 | specified either by assigning the name of a mesh file to |
---|
| 1768 | \code{source} or by specifying the points, triangle and boundary of the |
---|
| 1769 | mesh. |
---|
| 1770 | \end{classdesc} |
---|
| 1771 | |
---|
| 1772 | \subsection{Key Methods of Domain} |
---|
| 1773 | |
---|
| 1774 | \begin{methoddesc} {set\_name}{name} |
---|
[4209] | 1775 | Module: \refmodule{abstract\_2d\_finite\_volumes.domain}, |
---|
| 1776 | page \pageref{mod:domain} |
---|
[4123] | 1777 | |
---|
| 1778 | Assigns the name \code{name} to the domain. |
---|
| 1779 | \end{methoddesc} |
---|
| 1780 | |
---|
| 1781 | \begin{methoddesc} {get\_name}{} |
---|
| 1782 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 1783 | |
---|
| 1784 | Returns the name assigned to the domain by \code{set\_name}. If no name has been |
---|
| 1785 | assigned, returns \code{`domain'}. |
---|
| 1786 | \end{methoddesc} |
---|
| 1787 | |
---|
| 1788 | \begin{methoddesc} {set\_datadir}{name} |
---|
| 1789 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 1790 | |
---|
[4209] | 1791 | Specifies the directory used for SWW files, assigning it to the |
---|
[4123] | 1792 | pathname \code{name}. The default value, before |
---|
| 1793 | \code{set\_datadir} has been run, is the value \code{default\_datadir} |
---|
| 1794 | specified in \code{config.py}. |
---|
| 1795 | |
---|
| 1796 | Since different operating systems use different formats for specifying pathnames, |
---|
| 1797 | it is necessary to specify path separators using the Python code \code{os.sep}, rather than |
---|
| 1798 | the operating-specific ones such as `$\slash$' or `$\backslash$'. |
---|
| 1799 | For this to work you will need to include the statement \code{import os} |
---|
| 1800 | in your code, before the first appearance of \code{set\_datadir}. |
---|
| 1801 | |
---|
| 1802 | For example, to set the data directory to a subdirectory |
---|
| 1803 | \code{data} of the directory \code{project}, you could use |
---|
| 1804 | the statements: |
---|
| 1805 | |
---|
| 1806 | {\small \begin{verbatim} |
---|
| 1807 | import os |
---|
| 1808 | domain.set_datadir{'project' + os.sep + 'data'} |
---|
| 1809 | \end{verbatim}} |
---|
| 1810 | \end{methoddesc} |
---|
| 1811 | |
---|
| 1812 | \begin{methoddesc} {get\_datadir}{} |
---|
| 1813 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 1814 | |
---|
| 1815 | Returns the data directory set by \code{set\_datadir} or, |
---|
| 1816 | if \code{set\_datadir} has not |
---|
| 1817 | been run, returns the value \code{default\_datadir} specified in |
---|
| 1818 | \code{config.py}. |
---|
| 1819 | \end{methoddesc} |
---|
| 1820 | |
---|
[4258] | 1821 | |
---|
| 1822 | \begin{methoddesc} {set\_minimum_allowed_height}{} |
---|
| 1823 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 1824 | |
---|
[4377] | 1825 | Set the minimum depth (in meters) that will be recognised in |
---|
[4258] | 1826 | the numerical scheme (including limiters and flux computations) |
---|
[4377] | 1827 | |
---|
| 1828 | Default value is $10^{-3}$ m, but by setting this to a greater value, |
---|
| 1829 | e.g.\ for large scale simulations, the computation time can be |
---|
| 1830 | significantly reduced. |
---|
[4258] | 1831 | \end{methoddesc} |
---|
| 1832 | |
---|
| 1833 | |
---|
[4123] | 1834 | \begin{methoddesc} {set\_minimum_storable_height}{} |
---|
| 1835 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 1836 | |
---|
| 1837 | Sets the minimum depth that will be recognised when writing |
---|
| 1838 | to an sww file. This is useful for removing thin water layers |
---|
| 1839 | that seems to be caused by friction creep. |
---|
| 1840 | \end{methoddesc} |
---|
| 1841 | |
---|
| 1842 | |
---|
| 1843 | \begin{methoddesc} {set\_maximum_allowed_speed}{} |
---|
| 1844 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 1845 | |
---|
| 1846 | Set the maximum particle speed that is allowed in water |
---|
| 1847 | shallower than minimum_allowed_height. This is useful for |
---|
| 1848 | controlling speeds in very thin layers of water and at the same time |
---|
| 1849 | allow some movement avoiding pooling of water. |
---|
| 1850 | \end{methoddesc} |
---|
| 1851 | |
---|
| 1852 | |
---|
| 1853 | \begin{methoddesc} {set\_time}{time=0.0} |
---|
| 1854 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 1855 | |
---|
| 1856 | Sets the initial time, in seconds, for the simulation. The |
---|
| 1857 | default is 0.0. |
---|
| 1858 | \end{methoddesc} |
---|
| 1859 | |
---|
| 1860 | \begin{methoddesc} {set\_default\_order}{n} |
---|
| 1861 | Sets the default (spatial) order to the value specified by |
---|
| 1862 | \code{n}, which must be either 1 or 2. (Assigning any other value |
---|
| 1863 | to \code{n} will cause an error.) |
---|
| 1864 | \end{methoddesc} |
---|
| 1865 | |
---|
| 1866 | |
---|
[4471] | 1867 | \begin{methoddesc} {set\_store\_vertices\_uniquely}{flag} |
---|
[4123] | 1868 | Decide whether vertex values should be stored uniquely as |
---|
| 1869 | computed in the model or whether they should be reduced to one |
---|
[4471] | 1870 | value per vertex using averaging. |
---|
[4782] | 1871 | |
---|
[4953] | 1872 | Triangles stored in the sww file can be discontinuous reflecting |
---|
| 1873 | the internal representation of the finite-volume scheme |
---|
| 1874 | (this is a feature allowing for arbitrary steepness). |
---|
| 1875 | However, for visual purposes and also for use with \code{Field\_boundary} |
---|
| 1876 | (and \code{File\_boundary}) it is often desirable to store triangles |
---|
| 1877 | with values at each vertex point as the average of the potentially |
---|
| 1878 | discontinuous numbers found at vertices of different triangles sharing the |
---|
| 1879 | same vertex location. |
---|
| 1880 | |
---|
| 1881 | Storing one way or the other is controlled in ANUGA through the method |
---|
[4782] | 1882 | \code{domain.store\_vertices\_uniquely}. Options are |
---|
[4953] | 1883 | \begin{itemize} |
---|
| 1884 | \item \code{domain.store\_vertices\_uniquely(True)}: Allow discontinuities in the sww file |
---|
| 1885 | \item \code{domain.store\_vertices\_uniquely(False)}: (Default). |
---|
| 1886 | Average values |
---|
| 1887 | to ensure continuity in sww file. The latter also makes for smaller |
---|
[4782] | 1888 | sww files. |
---|
[4953] | 1889 | \end{itemize} |
---|
[4782] | 1890 | |
---|
[4123] | 1891 | \end{methoddesc} |
---|
| 1892 | |
---|
| 1893 | |
---|
| 1894 | % Structural methods |
---|
| 1895 | \begin{methoddesc}{get\_nodes}{absolute=False} |
---|
| 1896 | Return x,y coordinates of all nodes in mesh. |
---|
| 1897 | |
---|
| 1898 | The nodes are ordered in an Nx2 array where N is the number of nodes. |
---|
| 1899 | This is the same format they were provided in the constructor |
---|
| 1900 | i.e. without any duplication. |
---|
| 1901 | |
---|
| 1902 | Boolean keyword argument absolute determines whether coordinates |
---|
| 1903 | are to be made absolute by taking georeference into account |
---|
| 1904 | Default is False as many parts of ANUGA expects relative coordinates. |
---|
| 1905 | \end{methoddesc} |
---|
| 1906 | |
---|
| 1907 | |
---|
| 1908 | \begin{methoddesc}{get\_vertex_coordinates}{absolute=False} |
---|
[4209] | 1909 | |
---|
| 1910 | Return vertex coordinates for all triangles. |
---|
| 1911 | |
---|
[4123] | 1912 | Return all vertex coordinates for all triangles as a 3*M x 2 array |
---|
| 1913 | where the jth vertex of the ith triangle is located in row 3*i+j and |
---|
| 1914 | M the number of triangles in the mesh. |
---|
| 1915 | |
---|
| 1916 | Boolean keyword argument absolute determines whether coordinates |
---|
| 1917 | are to be made absolute by taking georeference into account |
---|
| 1918 | Default is False as many parts of ANUGA expects relative coordinates. |
---|
| 1919 | \end{methoddesc} |
---|
[4209] | 1920 | |
---|
| 1921 | |
---|
[4123] | 1922 | \begin{methoddesc}{get\_triangles}{indices=None} |
---|
| 1923 | |
---|
| 1924 | Return Mx3 integer array where M is the number of triangles. |
---|
| 1925 | Each row corresponds to one triangle and the three entries are |
---|
| 1926 | indices into the mesh nodes which can be obtained using the method |
---|
| 1927 | get\_nodes() |
---|
| 1928 | |
---|
| 1929 | Optional argument, indices is the set of triangle ids of interest. |
---|
| 1930 | \end{methoddesc} |
---|
[4209] | 1931 | |
---|
[4123] | 1932 | \begin{methoddesc}{get\_disconnected\_triangles}{} |
---|
| 1933 | |
---|
| 1934 | Get mesh based on nodes obtained from get_vertex_coordinates. |
---|
| 1935 | |
---|
| 1936 | Return array Mx3 array of integers where each row corresponds to |
---|
| 1937 | a triangle. A triangle is a triplet of indices into |
---|
| 1938 | point coordinates obtained from get_vertex_coordinates and each |
---|
| 1939 | index appears only once.\\ |
---|
| 1940 | |
---|
| 1941 | This provides a mesh where no triangles share nodes |
---|
| 1942 | (hence the name disconnected triangles) and different |
---|
| 1943 | nodes may have the same coordinates.\\ |
---|
| 1944 | |
---|
| 1945 | This version of the mesh is useful for storing meshes with |
---|
| 1946 | discontinuities at each node and is e.g. used for storing |
---|
| 1947 | data in sww files.\\ |
---|
| 1948 | |
---|
| 1949 | The triangles created will have the format |
---|
| 1950 | |
---|
[4209] | 1951 | {\small \begin{verbatim} |
---|
[4123] | 1952 | [[0,1,2], |
---|
| 1953 | [3,4,5], |
---|
| 1954 | [6,7,8], |
---|
| 1955 | ... |
---|
[4209] | 1956 | [3*M-3 3*M-2 3*M-1]] |
---|
| 1957 | \end{verbatim}} |
---|
[4123] | 1958 | \end{methoddesc} |
---|
| 1959 | |
---|
| 1960 | |
---|
| 1961 | |
---|
| 1962 | %%%%%% |
---|
[5508] | 1963 | \section{Initial Conditions}\index{Initial Conditions} |
---|
| 1964 | \label{sec:initial conditions} |
---|
[4123] | 1965 | In standard usage of partial differential equations, initial conditions |
---|
| 1966 | refers to the values associated to the system variables (the conserved |
---|
| 1967 | quantities here) for \code{time = 0}. In setting up a scenario script |
---|
| 1968 | as described in Sections \ref{sec:simpleexample} and \ref{sec:realdataexample}, |
---|
| 1969 | \code{set_quantity} is used to define the initial conditions of variables |
---|
| 1970 | other than the conserved quantities, such as friction. Here, we use the terminology |
---|
| 1971 | of initial conditions to refer to initial values for variables which need |
---|
| 1972 | prescription to solve the shallow water wave equation. Further, it must be noted |
---|
| 1973 | that \code{set_quantity} does not necessarily have to be used in the initial |
---|
| 1974 | condition setting; it can be used at any time throughout the simulation. |
---|
| 1975 | |
---|
| 1976 | \begin{methoddesc}{set\_quantity}{name, |
---|
| 1977 | numeric = None, |
---|
| 1978 | quantity = None, |
---|
| 1979 | function = None, |
---|
| 1980 | geospatial_data = None, |
---|
| 1981 | filename = None, |
---|
| 1982 | attribute_name = None, |
---|
| 1983 | alpha = None, |
---|
| 1984 | location = 'vertices', |
---|
| 1985 | indices = None, |
---|
| 1986 | verbose = False, |
---|
| 1987 | use_cache = False} |
---|
| 1988 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 1989 | (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values}) |
---|
| 1990 | |
---|
| 1991 | This function is used to assign values to individual quantities for a |
---|
| 1992 | domain. It is very flexible and can be used with many data types: a |
---|
| 1993 | statement of the form \code{domain.set\_quantity(name, x)} can be used |
---|
| 1994 | to define a quantity having the name \code{name}, where the other |
---|
| 1995 | argument \code{x} can be any of the following: |
---|
| 1996 | |
---|
| 1997 | \begin{itemize} |
---|
| 1998 | \item a number, in which case all vertices in the mesh gets that for |
---|
| 1999 | the quantity in question. |
---|
| 2000 | \item a list of numbers or a Numeric array ordered the same way as the mesh vertices. |
---|
| 2001 | \item a function (e.g.\ see the samples introduced in Chapter 2) |
---|
| 2002 | \item an expression composed of other quantities and numbers, arrays, lists (for |
---|
| 2003 | example, a linear combination of quantities, such as |
---|
| 2004 | \code{domain.set\_quantity('stage','elevation'+x))} |
---|
| 2005 | \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. |
---|
[4209] | 2006 | \item a geospatial dataset (See Section \ref{sec:geospatial}). |
---|
[4123] | 2007 | Optional argument attribute\_name applies here as with files. |
---|
| 2008 | \end{itemize} |
---|
| 2009 | |
---|
| 2010 | |
---|
| 2011 | Exactly one of the arguments |
---|
| 2012 | numeric, quantity, function, points, filename |
---|
| 2013 | must be present. |
---|
| 2014 | |
---|
| 2015 | |
---|
| 2016 | Set quantity will look at the type of the second argument (\code{numeric}) and |
---|
| 2017 | determine what action to take. |
---|
| 2018 | |
---|
| 2019 | Values can also be set using the appropriate keyword arguments. |
---|
| 2020 | If 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)} |
---|
| 2021 | are all equivalent. |
---|
| 2022 | |
---|
| 2023 | |
---|
| 2024 | Other optional arguments are |
---|
| 2025 | \begin{itemize} |
---|
| 2026 | \item \code{indices} which is a list of ids of triangles to which set\_quantity should apply its assignment of values. |
---|
| 2027 | \item \code{location} determines which part of the triangles to assign |
---|
| 2028 | to. Options are 'vertices' (default), 'edges', 'unique vertices', and 'centroids'. |
---|
| 2029 | \end{itemize} |
---|
| 2030 | |
---|
| 2031 | %%% |
---|
| 2032 | \anuga provides a number of predefined initial conditions to be used |
---|
[4209] | 2033 | with \code{set\_quantity}. See for example callable object |
---|
[4123] | 2034 | \code{slump\_tsunami} below. |
---|
| 2035 | |
---|
| 2036 | \end{methoddesc} |
---|
| 2037 | |
---|
| 2038 | |
---|
| 2039 | |
---|
| 2040 | |
---|
| 2041 | \begin{funcdesc}{set_region}{tag, quantity, X, location='vertices'} |
---|
| 2042 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
[4209] | 2043 | |
---|
[4123] | 2044 | (see also \module{abstract\_2d\_finite\_volumes.quantity.set\_values}) |
---|
[4209] | 2045 | |
---|
[4123] | 2046 | This function is used to assign values to individual quantities given |
---|
[4209] | 2047 | a regional tag. It is similar to \code{set\_quantity}. |
---|
[4736] | 2048 | For example, if in the mesh-generator a regional tag of 'ditch' was |
---|
[4123] | 2049 | used, set\_region can be used to set elevation of this region to |
---|
| 2050 | -10m. X is the constant or function to be applied to the quantity, |
---|
| 2051 | over the tagged region. Location describes how the values will be |
---|
| 2052 | applied. Options are 'vertices' (default), 'edges', 'unique |
---|
| 2053 | vertices', and 'centroids'. |
---|
| 2054 | |
---|
| 2055 | This method can also be called with a list of region objects. This is |
---|
| 2056 | useful for adding quantities in regions, and having one quantity |
---|
| 2057 | value based on another quantity. See \module{abstract\_2d\_finite\_volumes.region} for |
---|
| 2058 | more details. |
---|
| 2059 | \end{funcdesc} |
---|
| 2060 | |
---|
| 2061 | |
---|
| 2062 | |
---|
| 2063 | |
---|
| 2064 | \begin{funcdesc}{slump_tsunami}{length, depth, slope, width=None, thickness=None, |
---|
| 2065 | x0=0.0, y0=0.0, alpha=0.0, |
---|
| 2066 | gravity=9.8, gamma=1.85, |
---|
| 2067 | massco=1, dragco=1, frictionco=0, psi=0, |
---|
| 2068 | dx=None, kappa=3.0, kappad=0.8, zsmall=0.01, |
---|
| 2069 | domain=None, |
---|
| 2070 | verbose=False} |
---|
| 2071 | Module: \module{shallow\_water.smf} |
---|
| 2072 | |
---|
| 2073 | This function returns a callable object representing an initial water |
---|
| 2074 | displacement generated by a submarine sediment failure. These failures can take the form of |
---|
| 2075 | a submarine slump or slide. In the case of a slide, use \code{slide_tsunami} instead. |
---|
| 2076 | |
---|
| 2077 | The arguments include as a minimum, the slump or slide length, the water depth to the centre of sediment |
---|
| 2078 | mass, and the bathymetric slope. Other slump or slide parameters can be included if they are known. |
---|
| 2079 | \end{funcdesc} |
---|
| 2080 | |
---|
| 2081 | |
---|
| 2082 | %%% |
---|
| 2083 | \begin{funcdesc}{file\_function}{filename, |
---|
| 2084 | domain = None, |
---|
| 2085 | quantities = None, |
---|
| 2086 | interpolation_points = None, |
---|
| 2087 | verbose = False, |
---|
| 2088 | use_cache = False} |
---|
| 2089 | Module: \module{abstract\_2d\_finite\_volumes.util} |
---|
| 2090 | |
---|
| 2091 | Reads the time history of spatial data for |
---|
| 2092 | specified interpolation points from a NetCDF file (\code{filename}) |
---|
| 2093 | and returns |
---|
[5555] | 2094 | a callable object. \code{filename} could be a \code{sww} or \code{sts} file. |
---|
[4123] | 2095 | Returns interpolated values based on the input |
---|
| 2096 | file using the underlying \code{interpolation\_function}. |
---|
| 2097 | |
---|
| 2098 | \code{quantities} is either the name of a single quantity to be |
---|
| 2099 | interpolated or a list of such quantity names. In the second case, the resulting |
---|
| 2100 | function will return a tuple of values---one for each quantity. |
---|
| 2101 | |
---|
| 2102 | \code{interpolation\_points} is a list of absolute coordinates or a |
---|
| 2103 | geospatial object |
---|
| 2104 | for points at which values are sought. |
---|
| 2105 | |
---|
[5555] | 2106 | \code{boundary_polygon} is a list of coordinates specifying the vertices of the boundary. This must be the same polygon as used when calling \code{create_mesh_from_regions}. This argument can only be used when reading boundary data from the STS format. |
---|
| 2107 | |
---|
[4123] | 2108 | The model time stored within the file function can be accessed using |
---|
| 2109 | the method \code{f.get\_time()} |
---|
| 2110 | |
---|
| 2111 | |
---|
| 2112 | The underlying algorithm used is as follows:\\ |
---|
| 2113 | Given a time series (i.e.\ a series of values associated with |
---|
[5555] | 2114 | different times), whose values are either just numbers, a set of |
---|
[4123] | 2115 | numbers defined at the vertices of a triangular mesh (such as those |
---|
[5555] | 2116 | stored in SWW files) or a set of |
---|
| 2117 | numbers defined at a number of points on the boundary (such as those |
---|
| 2118 | stored in STS files), \code{Interpolation\_function} is used to |
---|
[4123] | 2119 | create a callable object that interpolates a value for an arbitrary |
---|
| 2120 | time \code{t} within the model limits and possibly a point \code{(x, |
---|
| 2121 | y)} within a mesh region. |
---|
| 2122 | |
---|
| 2123 | The actual time series at which data is available is specified by |
---|
| 2124 | means of an array \code{time} of monotonically increasing times. The |
---|
| 2125 | quantities containing the values to be interpolated are specified in |
---|
| 2126 | an array---or dictionary of arrays (used in conjunction with the |
---|
| 2127 | optional argument \code{quantity\_names}) --- called |
---|
| 2128 | \code{quantities}. The optional arguments \code{vertex\_coordinates} |
---|
| 2129 | and \code{triangles} represent the spatial mesh associated with the |
---|
[5555] | 2130 | quantity arrays. If omitted the function must be created using an STS file |
---|
| 2131 | or a TMS file. |
---|
[4123] | 2132 | |
---|
| 2133 | Since, in practice, values need to be computed at specified points, |
---|
| 2134 | the syntax allows the user to specify, once and for all, a list |
---|
| 2135 | \code{interpolation\_points} of points at which values are required. |
---|
| 2136 | In this case, the function may be called using the form \code{f(t, |
---|
| 2137 | id)}, where \code{id} is an index for the list |
---|
| 2138 | \code{interpolation\_points}. |
---|
| 2139 | |
---|
| 2140 | |
---|
| 2141 | \end{funcdesc} |
---|
| 2142 | |
---|
[5719] | 2143 | % FIXME (OLE): Why has this been commented out? |
---|
[4123] | 2144 | %%% |
---|
| 2145 | %% \begin{classdesc}{Interpolation\_function}{self, |
---|
| 2146 | %% time, |
---|
| 2147 | %% quantities, |
---|
| 2148 | %% quantity_names = None, |
---|
| 2149 | %% vertex_coordinates = None, |
---|
| 2150 | %% triangles = None, |
---|
| 2151 | %% interpolation_points = None, |
---|
| 2152 | %% verbose = False} |
---|
| 2153 | %% Module: \module{abstract\_2d\_finite\_volumes.least\_squares} |
---|
| 2154 | |
---|
| 2155 | %% Given a time series (i.e.\ a series of values associated with |
---|
| 2156 | %% different times), whose values are either just numbers or a set of |
---|
| 2157 | %% numbers defined at the vertices of a triangular mesh (such as those |
---|
| 2158 | %% stored in SWW files), \code{Interpolation\_function} is used to |
---|
| 2159 | %% create a callable object that interpolates a value for an arbitrary |
---|
| 2160 | %% time \code{t} within the model limits and possibly a point \code{(x, |
---|
| 2161 | %% y)} within a mesh region. |
---|
| 2162 | |
---|
| 2163 | %% The actual time series at which data is available is specified by |
---|
| 2164 | %% means of an array \code{time} of monotonically increasing times. The |
---|
| 2165 | %% quantities containing the values to be interpolated are specified in |
---|
| 2166 | %% an array---or dictionary of arrays (used in conjunction with the |
---|
| 2167 | %% optional argument \code{quantity\_names}) --- called |
---|
| 2168 | %% \code{quantities}. The optional arguments \code{vertex\_coordinates} |
---|
| 2169 | %% and \code{triangles} represent the spatial mesh associated with the |
---|
| 2170 | %% quantity arrays. If omitted the function created by |
---|
| 2171 | %% \code{Interpolation\_function} will be a function of \code{t} only. |
---|
| 2172 | |
---|
| 2173 | %% Since, in practice, values need to be computed at specified points, |
---|
| 2174 | %% the syntax allows the user to specify, once and for all, a list |
---|
| 2175 | %% \code{interpolation\_points} of points at which values are required. |
---|
| 2176 | %% In this case, the function may be called using the form \code{f(t, |
---|
| 2177 | %% id)}, where \code{id} is an index for the list |
---|
| 2178 | %% \code{interpolation\_points}. |
---|
| 2179 | |
---|
| 2180 | %% \end{classdesc} |
---|
| 2181 | |
---|
| 2182 | %%% |
---|
| 2183 | %\begin{funcdesc}{set\_region}{functions} |
---|
| 2184 | %[Low priority. Will be merged into set\_quantity] |
---|
| 2185 | |
---|
| 2186 | %Module:\module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2187 | %\end{funcdesc} |
---|
| 2188 | |
---|
| 2189 | |
---|
| 2190 | |
---|
| 2191 | %%%%%% |
---|
| 2192 | \section{Boundary Conditions}\index{boundary conditions} |
---|
[5508] | 2193 | \label{sec:boundary conditions} |
---|
[4123] | 2194 | |
---|
| 2195 | \anuga provides a large number of predefined boundary conditions, |
---|
| 2196 | represented by objects such as \code{Reflective\_boundary(domain)} and |
---|
| 2197 | \code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, described in the examples |
---|
| 2198 | in Chapter 2. Alternatively, you may prefer to ``roll your own'', |
---|
| 2199 | following the method explained in Section \ref{sec:roll your own}. |
---|
| 2200 | |
---|
| 2201 | These boundary objects may be used with the function \code{set\_boundary} described below |
---|
| 2202 | to assign boundary conditions according to the tags used to label boundary segments. |
---|
| 2203 | |
---|
| 2204 | \begin{methoddesc}{set\_boundary}{boundary_map} |
---|
| 2205 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2206 | |
---|
| 2207 | This function allows you to assign a boundary object (corresponding to a |
---|
| 2208 | pre-defined or user-specified boundary condition) to every boundary segment that |
---|
| 2209 | has been assigned a particular tag. |
---|
| 2210 | |
---|
| 2211 | This is done by specifying a dictionary \code{boundary\_map}, whose values are the boundary objects |
---|
| 2212 | and whose keys are the symbolic tags. |
---|
| 2213 | |
---|
| 2214 | \end{methoddesc} |
---|
| 2215 | |
---|
| 2216 | \begin{methoddesc} {get\_boundary\_tags}{} |
---|
| 2217 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2218 | |
---|
| 2219 | Returns a list of the available boundary tags. |
---|
| 2220 | \end{methoddesc} |
---|
| 2221 | |
---|
| 2222 | %%% |
---|
| 2223 | \subsection{Predefined boundary conditions} |
---|
| 2224 | |
---|
| 2225 | \begin{classdesc}{Reflective\_boundary}{Boundary} |
---|
| 2226 | Module: \module{shallow\_water} |
---|
| 2227 | |
---|
| 2228 | Reflective boundary returns same conserved quantities as those present in |
---|
| 2229 | the neighbouring volume but reflected. |
---|
| 2230 | |
---|
| 2231 | This class is specific to the shallow water equation as it works with the |
---|
| 2232 | momentum quantities assumed to be the second and third conserved quantities. |
---|
| 2233 | \end{classdesc} |
---|
| 2234 | |
---|
| 2235 | %%% |
---|
| 2236 | \begin{classdesc}{Transmissive\_boundary}{domain = None} |
---|
| 2237 | Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions} |
---|
| 2238 | |
---|
| 2239 | A transmissive boundary returns the same conserved quantities as |
---|
| 2240 | those present in the neighbouring volume. |
---|
| 2241 | |
---|
| 2242 | The underlying domain must be specified when the boundary is instantiated. |
---|
| 2243 | \end{classdesc} |
---|
| 2244 | |
---|
| 2245 | %%% |
---|
| 2246 | \begin{classdesc}{Dirichlet\_boundary}{conserved_quantities=None} |
---|
| 2247 | Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions} |
---|
| 2248 | |
---|
| 2249 | A Dirichlet boundary returns constant values for each of conserved |
---|
| 2250 | quantities. In the example of \code{Dirichlet\_boundary([0.2, 0.0, 0.0])}, |
---|
| 2251 | the \code{stage} value at the boundary is 0.2 and the \code{xmomentum} and |
---|
| 2252 | \code{ymomentum} at the boundary are set to 0.0. The list must contain |
---|
| 2253 | a value for each conserved quantity. |
---|
| 2254 | \end{classdesc} |
---|
| 2255 | |
---|
| 2256 | %%% |
---|
| 2257 | \begin{classdesc}{Time\_boundary}{domain = None, f = None} |
---|
| 2258 | Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions} |
---|
| 2259 | |
---|
| 2260 | A time-dependent boundary returns values for the conserved |
---|
| 2261 | quantities as a function \code{f(t)} of time. The user must specify |
---|
| 2262 | the domain to get access to the model time. |
---|
| 2263 | \end{classdesc} |
---|
| 2264 | |
---|
| 2265 | %%% |
---|
| 2266 | \begin{classdesc}{File\_boundary}{Boundary} |
---|
| 2267 | Module: \module{abstract\_2d\_finite\_volumes.generic\_boundary\_conditions} |
---|
| 2268 | |
---|
[5555] | 2269 | This method may be used if the user wishes to apply a SWW file, STS file or |
---|
| 2270 | a time series file (TMS) to a boundary segment or segments. |
---|
[4123] | 2271 | The boundary values are obtained from a file and interpolated to the |
---|
| 2272 | appropriate segments for each conserved quantity. |
---|
[5657] | 2273 | |
---|
[5672] | 2274 | Optional argument \code{default\_boundary} can be used to specify another boundary object to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}. |
---|
[5744] | 2275 | The \code{default\_boundary} could be a simple Dirichlet condition or |
---|
| 2276 | even another \code{File\_boundary} |
---|
| 2277 | typically using data pertaining to another time interval. |
---|
[4123] | 2278 | \end{classdesc} |
---|
| 2279 | |
---|
[5672] | 2280 | \begin{classdesc}{Field\_boundary}{Boundary} |
---|
| 2281 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
[4123] | 2282 | |
---|
[5744] | 2283 | This method works in the same way as \code{File\_boundary} except that it |
---|
| 2284 | allows for the value of stage to be offset by a constant specified in the |
---|
[5672] | 2285 | keyword argument \code{mean\_stage}. |
---|
[4123] | 2286 | |
---|
[5672] | 2287 | This functionality allows for models to be run for a range of tides using |
---|
| 2288 | the same boundary information (from .sts, .sww or .tms files). The tidal value |
---|
[5744] | 2289 | for each run would then be specified in the keyword argument |
---|
| 2290 | \code{mean\_stage}. |
---|
| 2291 | If \code{mean\_stage} = 0.0, \code{Field\_boundary} and \code{File\_boundary} |
---|
| 2292 | behave identically. |
---|
[5672] | 2293 | |
---|
| 2294 | |
---|
| 2295 | Note that if the optional argument \code{default\_boundary} is specified |
---|
| 2296 | it's stage value will be adjusted by \code{mean\_stage} just like the values |
---|
[5744] | 2297 | obtained from the file. |
---|
[5672] | 2298 | |
---|
| 2299 | See \code{File\_boundary} for further details. |
---|
| 2300 | \end{classdesc} |
---|
| 2301 | |
---|
[4123] | 2302 | %%% |
---|
| 2303 | \begin{classdesc}{Transmissive\_Momentum\_Set\_Stage\_boundary}{Boundary} |
---|
| 2304 | Module: \module{shallow\_water} |
---|
| 2305 | |
---|
| 2306 | This boundary returns same momentum conserved quantities as |
---|
| 2307 | those present in its neighbour volume but sets stage as in a Time\_boundary. |
---|
| 2308 | The underlying domain must be specified when boundary is instantiated |
---|
| 2309 | |
---|
| 2310 | This type of boundary is useful when stage is known at the boundary as a |
---|
| 2311 | function of time, but momenta (or speeds) aren't. |
---|
| 2312 | |
---|
| 2313 | This class is specific to the shallow water equation as it works with the |
---|
| 2314 | momentum quantities assumed to be the second and third conserved quantities. |
---|
| 2315 | \end{classdesc} |
---|
| 2316 | |
---|
| 2317 | |
---|
| 2318 | \begin{classdesc}{Dirichlet\_Discharge\_boundary}{Boundary} |
---|
| 2319 | Module: \module{shallow\_water} |
---|
| 2320 | |
---|
| 2321 | Sets stage (stage0) |
---|
| 2322 | Sets momentum (wh0) in the inward normal direction. |
---|
| 2323 | \end{classdesc} |
---|
| 2324 | |
---|
| 2325 | |
---|
| 2326 | |
---|
| 2327 | \subsection{User-defined boundary conditions} |
---|
| 2328 | \label{sec:roll your own} |
---|
| 2329 | |
---|
| 2330 | All boundary classes must inherit from the generic boundary class |
---|
[4209] | 2331 | \code{Boundary} and have a method called \code{evaluate} which must |
---|
[4123] | 2332 | take as inputs \code{self, vol\_id, edge\_id} where self refers to the |
---|
| 2333 | object itself and vol\_id and edge\_id are integers referring to |
---|
[4209] | 2334 | particular edges. The method must return a list of three floating point |
---|
| 2335 | numbers representing values for \code{stage}, |
---|
[4123] | 2336 | \code{xmomentum} and \code{ymomentum}, respectively. |
---|
| 2337 | |
---|
[4209] | 2338 | The constructor of a particular boundary class may be used to specify |
---|
[4123] | 2339 | particular values or flags to be used by the \code{evaluate} method. |
---|
[4209] | 2340 | Please refer to the source code for the existing boundary conditions |
---|
[4123] | 2341 | for examples of how to implement boundary conditions. |
---|
| 2342 | |
---|
| 2343 | |
---|
| 2344 | |
---|
[5508] | 2345 | \section{Forcing Terms}\index{Forcing terms} |
---|
| 2346 | \label{sec:forcing terms} |
---|
[4123] | 2347 | |
---|
[5506] | 2348 | \anuga provides a number of predefined forcing functions to be used with simulations. |
---|
[5744] | 2349 | Gravity and friction are always calculated using the elevation and friction quantities, but the user may additionally add forcing terms to the list |
---|
| 2350 | \code{domain.forcing\_terms} and have them affect the model. |
---|
| 2351 | |
---|
[5508] | 2352 | Currently, predefined forcing terms are |
---|
[4123] | 2353 | |
---|
[5507] | 2354 | \begin{funcdesc}{General\_forcing}{} |
---|
[5506] | 2355 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
[4123] | 2356 | |
---|
[5506] | 2357 | This is a general class to modify any quantity according to a given rate of change. |
---|
| 2358 | Other specific forcing terms are based on this class but it can be used by itself as well (e.g.\ to modify momentum). |
---|
[5744] | 2359 | |
---|
[5506] | 2360 | The General\_forcing class takes as input: |
---|
[5744] | 2361 | \begin{itemize} |
---|
[5506] | 2362 | \item \code{domain}: a reference to the domain being evolved |
---|
| 2363 | \item \code{quantity\_name}: The name of the quantity that will be affected by this forcing term |
---|
| 2364 | \item \code{rate}: The rate at which the quantity should change. The parameter \code{rate} can be eithe a constant or a |
---|
[5744] | 2365 | function of time. Positive values indicate increases, |
---|
[5506] | 2366 | negative values indicate decreases. |
---|
| 2367 | The parametr \code{rate} can be \code{None} at initialisation but must be specified |
---|
| 2368 | before forcing term is applied (i.e. simulation has started). |
---|
| 2369 | The default value is 0.0 - i.e.\ no forcing. |
---|
| 2370 | \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius. |
---|
| 2371 | \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon. |
---|
| 2372 | \end{itemize} |
---|
| 2373 | Note specifying both center, radius and polygon will cause an exception to be thrown. |
---|
[5566] | 2374 | Moreover, if the specified polygon or circle does not lie fully within the mesh boundary an Exception will be thrown. |
---|
[5507] | 2375 | |
---|
[5744] | 2376 | \bigskip |
---|
[5506] | 2377 | Example: |
---|
[5744] | 2378 | {\scriptsize \begin{verbatim} |
---|
[5508] | 2379 | P = [[x0, y0], [x1, y0], [x1, y1], [x0, y1]] # Square polygon |
---|
[5744] | 2380 | |
---|
[5508] | 2381 | xmom = General_forcing(domain, 'xmomentum', polygon=P) |
---|
| 2382 | ymom = General_forcing(domain, 'ymomentum', polygon=P) |
---|
[4123] | 2383 | |
---|
[5508] | 2384 | xmom.rate = f |
---|
| 2385 | ymom.rate = g |
---|
[5744] | 2386 | |
---|
[5508] | 2387 | domain.forcing_terms.append(xmom) |
---|
| 2388 | domain.forcing_terms.append(ymom) |
---|
[5507] | 2389 | \end{verbatim}} |
---|
[5506] | 2390 | Here, \code{f}, \code{g} are assumed to be defined as functions of time providing a time dependent rate of change for xmomentum and ymomentum respectively. |
---|
[5744] | 2391 | P is assumed to be polygon, specified as a list of points. |
---|
[5506] | 2392 | |
---|
[5744] | 2393 | \end{funcdesc} |
---|
[5506] | 2394 | |
---|
[5744] | 2395 | |
---|
[5506] | 2396 | \begin{funcdesc}{Inflow}{} |
---|
| 2397 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 2398 | |
---|
[5566] | 2399 | This is a general class for inflow and abstraction of water according to a given rate of change. |
---|
[5506] | 2400 | This class will always modify the \code{stage} quantity. |
---|
[5744] | 2401 | |
---|
[5506] | 2402 | Inflow is based on the General_forcing class so the functionality is similar. |
---|
[5744] | 2403 | |
---|
[5506] | 2404 | The Inflow class takes as input: |
---|
[5744] | 2405 | \begin{itemize} |
---|
[5506] | 2406 | \item \code{domain}: a reference to the domain being evolved |
---|
| 2407 | \item \code{rate}: The flow rate in $m^3/s$ at which the stage should change. The parameter \code{rate} can be eithe a constant or a |
---|
[5744] | 2408 | function of time. Positive values indicate inflow, |
---|
[5506] | 2409 | negative values indicate outflow. |
---|
| 2410 | |
---|
| 2411 | Note: The specified flow will be divided by the area of |
---|
| 2412 | the inflow region and then applied to update the |
---|
[5744] | 2413 | stage quantity. |
---|
[5506] | 2414 | \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius. |
---|
| 2415 | \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon. |
---|
| 2416 | \end{itemize} |
---|
[5507] | 2417 | |
---|
[5744] | 2418 | \bigskip |
---|
[5506] | 2419 | Example: |
---|
[5744] | 2420 | {\scriptsize \begin{verbatim} |
---|
[5506] | 2421 | hydrograph = Inflow(center=(320, 300), radius=10, |
---|
| 2422 | rate=file_function('QPMF_Rot_Sub13.tms')) |
---|
| 2423 | |
---|
| 2424 | domain.forcing_terms.append(hydrograph) |
---|
[5507] | 2425 | \end{verbatim}} |
---|
| 2426 | Here, \code{'QPMF_Rot_Sub13.tms'} is assumed to be a NetCDF file in the format \code{tms} defining a timeseries for a hydrograph. |
---|
[5744] | 2427 | \end{funcdesc} |
---|
[5506] | 2428 | |
---|
| 2429 | |
---|
| 2430 | \begin{funcdesc}{Rainfall}{} |
---|
| 2431 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 2432 | |
---|
| 2433 | This is a general class for implementing rainfall over the domain, possibly restricted to a given circle or polygon. |
---|
| 2434 | This class will always modify the \code{stage} quantity. |
---|
[5744] | 2435 | |
---|
[5506] | 2436 | Rainfall is based on the General_forcing class so the functionality is similar. |
---|
[5744] | 2437 | |
---|
[5506] | 2438 | The Rainfall class takes as input: |
---|
[5744] | 2439 | \begin{itemize} |
---|
[5506] | 2440 | \item \code{domain}: a reference to the domain being evolved |
---|
[5744] | 2441 | \item \code{rate}: Total rain rate over the specified domain. |
---|
[5506] | 2442 | Note: Raingauge Data needs to reflect the time step. |
---|
| 2443 | For example: if rain gauge is mm read every \code{dt} seconds, then the input |
---|
| 2444 | here is as \code{mm/dt} so 10 mm in 5 minutes becomes |
---|
| 2445 | 10/(5x60) = 0.0333mm/s. |
---|
| 2446 | |
---|
| 2447 | This parameter can be either a constant or a |
---|
[5744] | 2448 | function of time. Positive values indicate rain being added (or be used for general infiltration), |
---|
[5506] | 2449 | negative values indicate outflow at the specified rate (presumably this could model evaporation or abstraction). |
---|
| 2450 | \item \code{center, radius}: Optionally restrict forcing to a circle with given center and radius. |
---|
| 2451 | \item \code{polygon}: Optionally restrict forcing to an area enclosed by given polygon. |
---|
| 2452 | \end{itemize} |
---|
[5744] | 2453 | |
---|
| 2454 | \bigskip |
---|
[5506] | 2455 | Example: |
---|
[5744] | 2456 | {\scriptsize \begin{verbatim} |
---|
| 2457 | |
---|
| 2458 | catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) |
---|
[5506] | 2459 | domain.forcing_terms.append(catchmentrainfall) |
---|
| 2460 | |
---|
[5507] | 2461 | \end{verbatim}} |
---|
| 2462 | Here, \code{'Q100_2hr_Rain.tms'} is assumed to be a NetCDF file in the format \code{tms} defining a timeseries for the rainfall. |
---|
[5744] | 2463 | \end{funcdesc} |
---|
[5506] | 2464 | |
---|
| 2465 | |
---|
| 2466 | |
---|
| 2467 | \begin{funcdesc}{Culvert\_flow}{} |
---|
| 2468 | Module: \module{culver\_flows.culvert\_class} |
---|
| 2469 | |
---|
| 2470 | This is a general class for implementing flow through a culvert. |
---|
| 2471 | This class modifies the quantities \code{stage, xmomentum, ymomentum} in areas at both ends of the culvert. |
---|
[5744] | 2472 | |
---|
[5507] | 2473 | The Culvert\_flow forcing term uses \code{Inflow} and {General\_forcing} to update the quantities. The flow drection is determined on-the-fly so |
---|
[5506] | 2474 | openings are referenced simple as opening0 and opening1 with either being able to take the role as Inflow and Outflow. |
---|
[5744] | 2475 | |
---|
[5506] | 2476 | The Culvert\_flow class takes as input: |
---|
[5744] | 2477 | \begin{itemize} |
---|
[5506] | 2478 | \item \code{domain}: a reference to the domain being evolved |
---|
| 2479 | \item \code{label}: Short text naming the culvert |
---|
| 2480 | \item \code{description}: Text describing it |
---|
[5744] | 2481 | \item \code{end_point0}: Coordinates of one opening |
---|
[5506] | 2482 | \item \code{end_point1}: Coordinates of other opening |
---|
[5744] | 2483 | \item \code{width}: |
---|
[5506] | 2484 | \item \code{height}: |
---|
| 2485 | \item \code{diameter}: |
---|
| 2486 | \item \code{manning}: Mannings Roughness for Culvert |
---|
| 2487 | \item \code{invert_level0}: Invert level if not the same as the Elevation on the Domain |
---|
| 2488 | \item \code{invert_level1}: Invert level if not the same as the Elevation on the Domain |
---|
| 2489 | \item \code{culvert_routine}: Function specifying the calculation of flow based on energy difference between the two openings (see below) |
---|
| 2490 | \end{itemize} |
---|
| 2491 | |
---|
[5566] | 2492 | The user can specify different culvert routines. Hower ANUGA currently provides only one, namely the \code{boyd\_generalised\_culvert\_model} as used in the example below. |
---|
[5744] | 2493 | |
---|
| 2494 | \bigskip |
---|
[5506] | 2495 | Example: |
---|
[5744] | 2496 | {\scriptsize \begin{verbatim} |
---|
[5506] | 2497 | from anuga.culvert_flows.culvert_class import Culvert_flow |
---|
[5744] | 2498 | from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
---|
[5506] | 2499 | |
---|
| 2500 | culvert1 = Culvert_flow(domain, |
---|
| 2501 | label='Culvert No. 1', |
---|
[5744] | 2502 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
| 2503 | end_point0=[9.0, 2.5], |
---|
[5506] | 2504 | end_point1=[13.0, 2.5], |
---|
| 2505 | width=1.20,height=0.75, |
---|
[5744] | 2506 | culvert_routine=boyd_generalised_culvert_model, |
---|
[5506] | 2507 | verbose=True) |
---|
| 2508 | |
---|
| 2509 | culvert2 = Culvert_flow(domain, |
---|
| 2510 | label='Culvert No. 2', |
---|
[5744] | 2511 | description='This culvert is a circular test with d=1.2m', |
---|
| 2512 | end_point0=[9.0, 1.5], |
---|
[5506] | 2513 | end_point1=[30.0, 3.5], |
---|
| 2514 | diameter=1.20, |
---|
| 2515 | invert_level0=7, |
---|
[5744] | 2516 | culvert_routine=boyd_generalised_culvert_model, |
---|
[5506] | 2517 | verbose=True) |
---|
[5744] | 2518 | |
---|
[5506] | 2519 | domain.forcing_terms.append(culvert1) |
---|
| 2520 | domain.forcing_terms.append(culvert2) |
---|
| 2521 | |
---|
[5744] | 2522 | |
---|
[5507] | 2523 | \end{verbatim}} |
---|
[5744] | 2524 | \end{funcdesc} |
---|
[5506] | 2525 | |
---|
| 2526 | |
---|
| 2527 | |
---|
| 2528 | |
---|
| 2529 | |
---|
| 2530 | |
---|
[4123] | 2531 | \section{Evolution}\index{evolution} |
---|
[5508] | 2532 | \label{sec:evolution} |
---|
[4123] | 2533 | |
---|
| 2534 | \begin{methoddesc}{evolve}{yieldstep = None, finaltime = None, duration = None, skip_initial_step = False} |
---|
| 2535 | |
---|
| 2536 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2537 | |
---|
| 2538 | This function (a method of \class{domain}) is invoked once all the |
---|
| 2539 | preliminaries have been completed, and causes the model to progress |
---|
| 2540 | through successive steps in its evolution, storing results and |
---|
| 2541 | outputting statistics whenever a user-specified period |
---|
| 2542 | \code{yieldstep} is completed (generally during this period the |
---|
| 2543 | model will evolve through several steps internally |
---|
| 2544 | as the method forces the water speed to be calculated |
---|
| 2545 | on successive new cells). The user |
---|
| 2546 | specifies the total time period over which the evolution is to take |
---|
| 2547 | place, by specifying values (in seconds) for either \code{duration} |
---|
| 2548 | or \code{finaltime}, as well as the interval in seconds after which |
---|
| 2549 | results are to be stored and statistics output. |
---|
| 2550 | |
---|
| 2551 | You can include \method{evolve} in a statement of the type: |
---|
| 2552 | |
---|
| 2553 | {\small \begin{verbatim} |
---|
| 2554 | for t in domain.evolve(yieldstep, finaltime): |
---|
| 2555 | <Do something with domain and t> |
---|
| 2556 | \end{verbatim}} |
---|
| 2557 | |
---|
| 2558 | \end{methoddesc} |
---|
| 2559 | |
---|
| 2560 | |
---|
| 2561 | |
---|
| 2562 | \subsection{Diagnostics} |
---|
[4554] | 2563 | \label{sec:diagnostics} |
---|
[4123] | 2564 | |
---|
| 2565 | |
---|
| 2566 | \begin{funcdesc}{statistics}{} |
---|
| 2567 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2568 | |
---|
| 2569 | \end{funcdesc} |
---|
| 2570 | |
---|
| 2571 | \begin{funcdesc}{timestepping\_statistics}{} |
---|
| 2572 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2573 | |
---|
| 2574 | Returns a string of the following type for each |
---|
| 2575 | timestep: |
---|
| 2576 | |
---|
| 2577 | \code{Time = 0.9000, delta t in [0.00598964, 0.01177388], steps=12 |
---|
| 2578 | (12)} |
---|
| 2579 | |
---|
| 2580 | Here the numbers in \code{steps=12 (12)} indicate the number of steps taken and |
---|
[4369] | 2581 | the number of first-order steps, respectively.\\ |
---|
[4377] | 2582 | |
---|
| 2583 | The optional keyword argument \code{track_speeds=True} will |
---|
| 2584 | generate a histogram of speeds generated by each triangle. The |
---|
| 2585 | speeds relate to the size of the timesteps used by ANUGA and |
---|
| 2586 | this diagnostics may help pinpoint problem areas where excessive speeds |
---|
| 2587 | are generated. |
---|
| 2588 | |
---|
[4123] | 2589 | \end{funcdesc} |
---|
| 2590 | |
---|
| 2591 | |
---|
| 2592 | \begin{funcdesc}{boundary\_statistics}{quantities = None, tags = None} |
---|
| 2593 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2594 | |
---|
| 2595 | Returns a string of the following type when \code{quantities = 'stage'} and \code{tags = ['top', 'bottom']}: |
---|
| 2596 | |
---|
| 2597 | {\small \begin{verbatim} |
---|
| 2598 | Boundary values at time 0.5000: |
---|
| 2599 | top: |
---|
| 2600 | stage in [ -0.25821218, -0.02499998] |
---|
| 2601 | bottom: |
---|
| 2602 | stage in [ -0.27098821, -0.02499974] |
---|
| 2603 | \end{verbatim}} |
---|
| 2604 | |
---|
| 2605 | \end{funcdesc} |
---|
| 2606 | |
---|
| 2607 | |
---|
| 2608 | \begin{funcdesc}{get\_quantity}{name, location='vertices', indices = None} |
---|
| 2609 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
[4209] | 2610 | |
---|
[4123] | 2611 | Allow access to individual quantities and their methods |
---|
| 2612 | |
---|
| 2613 | \end{funcdesc} |
---|
| 2614 | |
---|
[4953] | 2615 | |
---|
[4705] | 2616 | \begin{funcdesc}{set\_quantities\_to\_be\_monitored}{} |
---|
| 2617 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
[4123] | 2618 | |
---|
[4953] | 2619 | Selects quantities and derived quantities for which extrema attained at internal timesteps |
---|
[4705] | 2620 | will be collected. |
---|
[4953] | 2621 | |
---|
| 2622 | Information can be tracked in the evolve loop by printing \code{quantity\_statistics} and |
---|
[4705] | 2623 | collected data will be stored in the sww file. |
---|
| 2624 | |
---|
[4953] | 2625 | Optional parameters \code{polygon} and \code{time\_interval} may be specified to restrict the |
---|
[4705] | 2626 | extremum computation. |
---|
[4953] | 2627 | \end{funcdesc} |
---|
| 2628 | |
---|
[4705] | 2629 | \begin{funcdesc}{quantity\_statistics}{} |
---|
| 2630 | Module: \module{abstract\_2d\_finite\_volumes.domain} |
---|
| 2631 | |
---|
| 2632 | Reports on extrema attained by selected quantities. |
---|
[4953] | 2633 | |
---|
[4705] | 2634 | Returns a string of the following type for each |
---|
| 2635 | timestep: |
---|
| 2636 | |
---|
[4953] | 2637 | \begin{verbatim} |
---|
[4705] | 2638 | Monitored quantities at time 1.0000: |
---|
| 2639 | stage-elevation: |
---|
| 2640 | values since time = 0.00 in [0.00000000, 0.30000000] |
---|
| 2641 | minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333) |
---|
| 2642 | maximum attained at time = 0.00000000, location = (0.83333333, 0.16666667) |
---|
| 2643 | ymomentum: |
---|
| 2644 | values since time = 0.00 in [0.00000000, 0.06241221] |
---|
| 2645 | minimum attained at time = 0.00000000, location = (0.33333333, 0.16666667) |
---|
| 2646 | maximum attained at time = 0.22472667, location = (0.83333333, 0.66666667) |
---|
| 2647 | xmomentum: |
---|
| 2648 | values since time = 0.00 in [-0.06062178, 0.47886313] |
---|
| 2649 | minimum attained at time = 0.00000000, location = (0.16666667, 0.33333333) |
---|
| 2650 | maximum attained at time = 0.35103646, location = (0.83333333, 0.16666667) |
---|
[4953] | 2651 | \end{verbatim} |
---|
[4705] | 2652 | |
---|
[4953] | 2653 | The quantities (and derived quantities) listed here must be selected at model |
---|
[4705] | 2654 | initialisation using the method \code{domain.set_quantities_to_be_monitored}.\\ |
---|
[4953] | 2655 | |
---|
[4705] | 2656 | The optional keyword argument \code{precision='\%.4f'} will |
---|
| 2657 | determine the precision used for floating point values in the output. |
---|
[4953] | 2658 | This diagnostics helps track extrema attained by the selected quantities |
---|
[4705] | 2659 | at every internal timestep. |
---|
| 2660 | |
---|
| 2661 | These values are also stored in the sww file for post processing. |
---|
| 2662 | |
---|
| 2663 | \end{funcdesc} |
---|
| 2664 | |
---|
[4953] | 2665 | |
---|
| 2666 | |
---|
[4123] | 2667 | \begin{funcdesc}{get\_values}{location='vertices', indices = None} |
---|
| 2668 | Module: \module{abstract\_2d\_finite\_volumes.quantity} |
---|
| 2669 | |
---|
| 2670 | Extract values for quantity as an array |
---|
| 2671 | |
---|
| 2672 | \end{funcdesc} |
---|
| 2673 | |
---|
[4209] | 2674 | |
---|
[5744] | 2675 | |
---|
[5719] | 2676 | \begin{funcdesc}{set\_values}{location='vertices', indices = None} |
---|
| 2677 | Module: \module{abstract\_2d\_finite\_volumes.quantity} |
---|
| 2678 | |
---|
| 2679 | Assign values to a quantity object. |
---|
| 2680 | This method works the same way as \code{set\_quantity} except that it doesn't take |
---|
| 2681 | a quantity name as the first argument. The reason to use \code{set\_values} is for |
---|
[5744] | 2682 | example to assign values to a new quantity that has been created but which is |
---|
[5719] | 2683 | not part of the domain's predefined quantities. |
---|
| 2684 | |
---|
[5744] | 2685 | The method \code{set\_values} is always called by \code{set\_quantity} |
---|
| 2686 | behind the scenes. |
---|
| 2687 | |
---|
[5719] | 2688 | \end{funcdesc} |
---|
| 2689 | |
---|
[5744] | 2690 | |
---|
| 2691 | |
---|
[4123] | 2692 | \begin{funcdesc}{get\_integral}{} |
---|
| 2693 | Module: \module{abstract\_2d\_finite\_volumes.quantity} |
---|
| 2694 | |
---|
| 2695 | Return computed integral over entire domain for this quantity |
---|
| 2696 | |
---|
| 2697 | \end{funcdesc} |
---|
| 2698 | |
---|
| 2699 | |
---|
[4209] | 2700 | |
---|
| 2701 | |
---|
[4123] | 2702 | \begin{funcdesc}{get\_maximum\_value}{indices = None} |
---|
| 2703 | Module: \module{abstract\_2d\_finite\_volumes.quantity} |
---|
| 2704 | |
---|
| 2705 | Return maximum value of quantity (on centroids) |
---|
| 2706 | |
---|
[4209] | 2707 | Optional argument indices is the set of element ids that |
---|
[4123] | 2708 | the operation applies to. If omitted all elements are considered. |
---|
| 2709 | |
---|
| 2710 | We do not seek the maximum at vertices as each vertex can |
---|
[4209] | 2711 | have multiple values - one for each triangle sharing it. |
---|
[4123] | 2712 | \end{funcdesc} |
---|
| 2713 | |
---|
| 2714 | |
---|
[4209] | 2715 | |
---|
[4123] | 2716 | \begin{funcdesc}{get\_maximum\_location}{indices = None} |
---|
| 2717 | Module: \module{abstract\_2d\_finite\_volumes.quantity} |
---|
| 2718 | |
---|
| 2719 | Return location of maximum value of quantity (on centroids) |
---|
| 2720 | |
---|
[4209] | 2721 | Optional argument indices is the set of element ids that |
---|
[4123] | 2722 | the operation applies to. |
---|
| 2723 | |
---|
| 2724 | We do not seek the maximum at vertices as each vertex can |
---|
| 2725 | have multiple values - one for each triangle sharing it. |
---|
| 2726 | |
---|
| 2727 | If there are multiple cells with same maximum value, the |
---|
[4209] | 2728 | first cell encountered in the triangle array is returned. |
---|
[4123] | 2729 | \end{funcdesc} |
---|
| 2730 | |
---|
| 2731 | |
---|
[4209] | 2732 | |
---|
[4123] | 2733 | \begin{funcdesc}{get\_wet\_elements}{indices=None} |
---|
[4209] | 2734 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
[4123] | 2735 | |
---|
| 2736 | Return indices for elements where h $>$ minimum_allowed_height |
---|
| 2737 | Optional argument indices is the set of element ids that the operation applies to. |
---|
| 2738 | \end{funcdesc} |
---|
| 2739 | |
---|
| 2740 | |
---|
| 2741 | \begin{funcdesc}{get\_maximum\_inundation\_elevation}{indices=None} |
---|
[4209] | 2742 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
[4123] | 2743 | |
---|
| 2744 | Return highest elevation where h $>$ 0.\\ |
---|
| 2745 | Optional argument indices is the set of element ids that the operation applies to.\\ |
---|
[4209] | 2746 | |
---|
| 2747 | Example to find maximum runup elevation:\\ |
---|
| 2748 | z = domain.get_maximum_inundation_elevation() |
---|
[4123] | 2749 | \end{funcdesc} |
---|
| 2750 | |
---|
| 2751 | |
---|
| 2752 | \begin{funcdesc}{get\_maximum\_inundation\_location}{indices=None} |
---|
[4209] | 2753 | Module: \module{shallow\_water.shallow\_water\_domain} |
---|
| 2754 | |
---|
[4123] | 2755 | Return location (x,y) of highest elevation where h $>$ 0.\\ |
---|
| 2756 | Optional argument indices is the set of element ids that the operation applies to.\\ |
---|
| 2757 | |
---|
| 2758 | Example to find maximum runup location:\\ |
---|
[4209] | 2759 | x, y = domain.get_maximum_inundation_location() |
---|
[4123] | 2760 | \end{funcdesc} |
---|
| 2761 | |
---|
[4209] | 2762 | |
---|
[4953] | 2763 | \section{Queries of SWW model output files} |
---|
| 2764 | After a model has been run, it is often useful to extract various information from the sww |
---|
[4554] | 2765 | output file (see Section \ref{sec:sww format}). This is typically more convenient than using the |
---|
[4953] | 2766 | diagnostics described in Section \ref{sec:diagnostics} which rely on the model running - something |
---|
| 2767 | that can be very time consuming. The sww files are easy and quick to read and offer much information |
---|
| 2768 | about the model results such as runup heights, time histories of selected quantities, |
---|
[4554] | 2769 | flow through cross sections and much more. |
---|
[4209] | 2770 | |
---|
[4953] | 2771 | \begin{funcdesc}{get\_maximum\_inundation\_elevation}{filename, polygon=None, |
---|
[4554] | 2772 | time_interval=None, verbose=False} |
---|
| 2773 | Module: \module{shallow\_water.data\_manager} |
---|
| 2774 | |
---|
[4556] | 2775 | Return highest elevation where depth is positive ($h > 0$) |
---|
[4554] | 2776 | |
---|
[4953] | 2777 | Example to find maximum runup elevation:\\ |
---|
[4554] | 2778 | max_runup = get_maximum_inundation_elevation(filename, |
---|
| 2779 | polygon=None, |
---|
| 2780 | time_interval=None, |
---|
| 2781 | verbose=False) |
---|
| 2782 | |
---|
[4953] | 2783 | |
---|
| 2784 | filename is a NetCDF sww file containing ANUGA model output. |
---|
[4554] | 2785 | Optional arguments polygon and time_interval restricts the maximum runup calculation |
---|
| 2786 | to a points that lie within the specified polygon and time interval. |
---|
| 2787 | |
---|
| 2788 | If no inundation is found within polygon and time_interval the return value |
---|
| 2789 | is None signifying "No Runup" or "Everything is dry". |
---|
| 2790 | |
---|
| 2791 | See doc string for general function get_maximum_inundation_data for details. |
---|
| 2792 | \end{funcdesc} |
---|
| 2793 | |
---|
| 2794 | |
---|
[4953] | 2795 | \begin{funcdesc}{get\_maximum\_inundation\_location}{filename, polygon=None, |
---|
[4554] | 2796 | time_interval=None, verbose=False} |
---|
| 2797 | Module: \module{shallow\_water.data\_manager} |
---|
| 2798 | |
---|
[4556] | 2799 | Return location (x,y) of highest elevation where depth is positive ($h > 0$) |
---|
[4554] | 2800 | |
---|
[4953] | 2801 | Example to find maximum runup location:\\ |
---|
[4554] | 2802 | max_runup_location = get_maximum_inundation_location(filename, |
---|
| 2803 | polygon=None, |
---|
| 2804 | time_interval=None, |
---|
| 2805 | verbose=False) |
---|
| 2806 | |
---|
[4953] | 2807 | |
---|
| 2808 | filename is a NetCDF sww file containing ANUGA model output. |
---|
[4554] | 2809 | Optional arguments polygon and time_interval restricts the maximum runup calculation |
---|
| 2810 | to a points that lie within the specified polygon and time interval. |
---|
| 2811 | |
---|
| 2812 | If no inundation is found within polygon and time_interval the return value |
---|
| 2813 | is None signifying "No Runup" or "Everything is dry". |
---|
| 2814 | |
---|
| 2815 | See doc string for general function get_maximum_inundation_data for details. |
---|
| 2816 | \end{funcdesc} |
---|
| 2817 | |
---|
| 2818 | |
---|
[4953] | 2819 | \begin{funcdesc}{sww2timeseries}{swwfiles, gauge_filename, production_dirs, report = None, reportname = None, |
---|
| 2820 | plot_quantity = None, generate_fig = False, surface = None, time_min = None, time_max = None, time_thinning = 1, |
---|
[4746] | 2821 | time_unit = None, title_on = None, use_cache = False, verbose = False} |
---|
| 2822 | |
---|
| 2823 | Module: \module{anuga.abstract\_2d\_finite\_volumes.util} |
---|
[4953] | 2824 | |
---|
[4746] | 2825 | Return csv files for the location in the \code{gauge_filename} and can also return plots of them |
---|
[4953] | 2826 | |
---|
[4746] | 2827 | See doc string for general function sww2timeseries for details. |
---|
| 2828 | |
---|
[4554] | 2829 | \end{funcdesc} |
---|
| 2830 | |
---|
[4953] | 2831 | |
---|
[5288] | 2832 | \begin{funcdesc}{get\_flow\_through\_cross\_section}{filename, cross\_section, verbose=False} |
---|
| 2833 | Module: \module{shallow\_water.data\_manager} |
---|
[4953] | 2834 | |
---|
[5566] | 2835 | Obtain flow $[m^3/s]$ perpendicular to specified cross section. |
---|
[5288] | 2836 | |
---|
| 2837 | Inputs: |
---|
[5744] | 2838 | \begin{itemize} |
---|
[5288] | 2839 | \item filename: Name of sww file containing ANUGA model output. |
---|
| 2840 | \item polyline: Representation of desired cross section - it may contain multiple |
---|
| 2841 | sections allowing for complex shapes. Assume absolute UTM coordinates. |
---|
[5744] | 2842 | \end{itemize} |
---|
[5288] | 2843 | |
---|
| 2844 | Output: |
---|
| 2845 | \begin{itemize} |
---|
| 2846 | \item time: All stored times in sww file |
---|
| 2847 | \item Q: Hydrograph of total flow across given segments for all stored times. |
---|
[5744] | 2848 | \end{itemize} |
---|
| 2849 | |
---|
[5288] | 2850 | The normal flow is computed for each triangle intersected by the polyline and |
---|
| 2851 | added up. Multiple segments at different angles are specified the normal flows |
---|
| 2852 | may partially cancel each other. |
---|
[5744] | 2853 | |
---|
[5288] | 2854 | Example to find flow through cross section: |
---|
[5744] | 2855 | |
---|
| 2856 | \begin{verbatim} |
---|
[5288] | 2857 | cross_section = [[x, 0], [x, width]] |
---|
| 2858 | time, Q = get_flow_through_cross_section(filename, |
---|
| 2859 | cross_section, |
---|
| 2860 | verbose=False) |
---|
[5744] | 2861 | \end{verbatim} |
---|
[5288] | 2862 | |
---|
| 2863 | |
---|
| 2864 | See doc string for general function get_maximum_inundation_data for details. |
---|
[5744] | 2865 | |
---|
[5288] | 2866 | \end{funcdesc} |
---|
| 2867 | |
---|
| 2868 | |
---|
| 2869 | |
---|
[4123] | 2870 | \section{Other} |
---|
| 2871 | |
---|
| 2872 | \begin{funcdesc}{domain.create\_quantity\_from\_expression}{string} |
---|
| 2873 | |
---|
| 2874 | Handy for creating derived quantities on-the-fly as for example |
---|
[4209] | 2875 | \begin{verbatim} |
---|
[4123] | 2876 | Depth = domain.create_quantity_from_expression('stage-elevation') |
---|
| 2877 | |
---|
[4209] | 2878 | exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5') |
---|
[4123] | 2879 | Absolute_momentum = domain.create_quantity_from_expression(exp) |
---|
[4209] | 2880 | \end{verbatim} |
---|
| 2881 | |
---|
[4123] | 2882 | %See also \file{Analytical\_solution\_circular\_hydraulic\_jump.py} for an example of use. |
---|
| 2883 | \end{funcdesc} |
---|
| 2884 | |
---|
| 2885 | |
---|
| 2886 | |
---|
| 2887 | |
---|
| 2888 | |
---|
| 2889 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 2890 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 2891 | |
---|
| 2892 | \chapter{\anuga System Architecture} |
---|
| 2893 | |
---|
| 2894 | |
---|
| 2895 | \section{File Formats} |
---|
| 2896 | \label{sec:file formats} |
---|
| 2897 | |
---|
| 2898 | \anuga makes use of a number of different file formats. The |
---|
| 2899 | following table lists all these formats, which are described in more |
---|
| 2900 | detail in the paragraphs below. |
---|
| 2901 | |
---|
| 2902 | \bigskip |
---|
| 2903 | |
---|
| 2904 | \begin{center} |
---|
| 2905 | |
---|
| 2906 | \begin{tabular}{|ll|} \hline |
---|
| 2907 | |
---|
| 2908 | \textbf{Extension} & \textbf{Description} \\ |
---|
| 2909 | \hline\hline |
---|
| 2910 | |
---|
[5555] | 2911 | \code{.sww} & NetCDF format for storing model output with mesh information |
---|
[4123] | 2912 | \code{f(t,x,y)}\\ |
---|
| 2913 | |
---|
[5555] | 2914 | \code{.sts} & NetCDF format for storing model ouput \code{f(t,x,y)} without any mesh information\\ |
---|
| 2915 | |
---|
[4123] | 2916 | \code{.tms} & NetCDF format for storing time series \code{f(t)}\\ |
---|
| 2917 | |
---|
[4662] | 2918 | \code{.csv/.txt} & ASCII format called points csv for storing |
---|
| 2919 | arbitrary points and associated attributes\\ |
---|
[4123] | 2920 | |
---|
| 2921 | \code{.pts} & NetCDF format for storing arbitrary points and |
---|
| 2922 | associated attributes\\ |
---|
| 2923 | |
---|
| 2924 | \code{.asc} & ASCII format of regular DEMs as output from ArcView\\ |
---|
| 2925 | |
---|
| 2926 | \code{.prj} & Associated ArcView file giving more metadata for |
---|
| 2927 | \code{.asc} format\\ |
---|
| 2928 | |
---|
| 2929 | \code{.ers} & ERMapper header format of regular DEMs for ArcView\\ |
---|
| 2930 | |
---|
| 2931 | \code{.dem} & NetCDF representation of regular DEM data\\ |
---|
| 2932 | |
---|
| 2933 | \code{.tsh} & ASCII format for storing meshes and associated |
---|
| 2934 | boundary and region info\\ |
---|
| 2935 | |
---|
| 2936 | \code{.msh} & NetCDF format for storing meshes and associated |
---|
| 2937 | boundary and region info\\ |
---|
| 2938 | |
---|
| 2939 | \code{.nc} & Native ferret NetCDF format\\ |
---|
| 2940 | |
---|
| 2941 | \code{.geo} & Houdinis ASCII geometry format (?) \\ \par \hline |
---|
| 2942 | %\caption{File formats used by \anuga} |
---|
| 2943 | \end{tabular} |
---|
| 2944 | |
---|
| 2945 | |
---|
| 2946 | \end{center} |
---|
| 2947 | |
---|
| 2948 | The above table shows the file extensions used to identify the |
---|
| 2949 | formats of files. However, typically, in referring to a format we |
---|
| 2950 | capitalise the extension and omit the initial full stop---thus, we |
---|
| 2951 | refer, for example, to `SWW files' or `PRJ files'. |
---|
| 2952 | |
---|
| 2953 | \bigskip |
---|
| 2954 | |
---|
| 2955 | A typical dataflow can be described as follows: |
---|
| 2956 | |
---|
| 2957 | \subsection{Manually Created Files} |
---|
| 2958 | |
---|
| 2959 | \begin{tabular}{ll} |
---|
| 2960 | ASC, PRJ & Digital elevation models (gridded)\\ |
---|
| 2961 | NC & Model outputs for use as boundary conditions (e.g. from MOST) |
---|
| 2962 | \end{tabular} |
---|
| 2963 | |
---|
| 2964 | \subsection{Automatically Created Files} |
---|
| 2965 | |
---|
| 2966 | \begin{tabular}{ll} |
---|
| 2967 | ASC, PRJ $\rightarrow$ DEM $\rightarrow$ PTS & Convert |
---|
| 2968 | DEMs to native \code{.pts} file\\ |
---|
| 2969 | |
---|
| 2970 | NC $\rightarrow$ SWW & Convert MOST boundary files to |
---|
| 2971 | boundary \code{.sww}\\ |
---|
| 2972 | |
---|
| 2973 | PTS + TSH $\rightarrow$ TSH with elevation & Least squares fit\\ |
---|
| 2974 | |
---|
| 2975 | TSH $\rightarrow$ SWW & Convert TSH to \code{.sww}-viewable using |
---|
| 2976 | \code{animate}\\ |
---|
| 2977 | |
---|
| 2978 | TSH + Boundary SWW $\rightarrow$ SWW & Simulation using |
---|
| 2979 | \code{\anuga}\\ |
---|
| 2980 | |
---|
| 2981 | Polygonal mesh outline $\rightarrow$ & TSH or MSH |
---|
| 2982 | \end{tabular} |
---|
| 2983 | |
---|
| 2984 | |
---|
| 2985 | |
---|
| 2986 | |
---|
| 2987 | \bigskip |
---|
| 2988 | |
---|
[5555] | 2989 | \subsection{SWW, STS and TMS Formats} |
---|
[4554] | 2990 | \label{sec:sww format} |
---|
[4123] | 2991 | |
---|
[5555] | 2992 | The SWW, STS and TMS formats are all NetCDF formats, and are of key |
---|
[4123] | 2993 | importance for \anuga. |
---|
| 2994 | |
---|
| 2995 | An SWW file is used for storing \anuga output and therefore pertains |
---|
| 2996 | to a set of points and a set of times at which a model is evaluated. |
---|
| 2997 | It contains, in addition to dimension information, the following |
---|
| 2998 | variables: |
---|
| 2999 | |
---|
| 3000 | \begin{itemize} |
---|
| 3001 | \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays |
---|
| 3002 | \item \code{elevation}, a Numeric array storing bed-elevations |
---|
| 3003 | \item \code{volumes}, a list specifying the points at the vertices of each of the |
---|
| 3004 | triangles |
---|
| 3005 | % Refer here to the example to be provided in describing the simple example |
---|
| 3006 | \item \code{time}, a Numeric array containing times for model |
---|
| 3007 | evaluation |
---|
| 3008 | \end{itemize} |
---|
| 3009 | |
---|
| 3010 | |
---|
[4209] | 3011 | The contents of an SWW file may be viewed using the anuga viewer |
---|
[4123] | 3012 | \code{animate}, which creates an on-screen geometric |
---|
| 3013 | representation. See section \ref{sec:animate} (page |
---|
| 3014 | \pageref{sec:animate}) in Appendix \ref{ch:supportingtools} for more |
---|
| 3015 | on \code{animate}. |
---|
| 3016 | |
---|
| 3017 | Alternatively, there are tools, such as \code{ncdump}, that allow |
---|
| 3018 | you to convert an NetCDF file into a readable format such as the |
---|
| 3019 | Class Definition Language (CDL). The following is an excerpt from a |
---|
| 3020 | CDL representation of the output file \file{runup.sww} generated |
---|
| 3021 | from running the simple example \file{runup.py} of |
---|
| 3022 | Chapter \ref{ch:getstarted}: |
---|
| 3023 | |
---|
| 3024 | \verbatiminput{examples/bedslopeexcerpt.cdl} |
---|
| 3025 | |
---|
| 3026 | The SWW format is used not only for output but also serves as input |
---|
| 3027 | for functions such as \function{file\_boundary} and |
---|
| 3028 | \function{file\_function}, described in Chapter \ref{ch:interface}. |
---|
| 3029 | |
---|
[5555] | 3030 | An STS file is used for storing a set of points and and associated set of times. |
---|
| 3031 | It contains, in addition to dimension information, the following |
---|
| 3032 | variables: |
---|
| 3033 | \begin{itemize} |
---|
| 3034 | \item \code{x} and \code{y}: coordinates of the points, represented as Numeric arrays |
---|
[5745] | 3035 | \item \code{permutation}: Original indices of the points as specified by |
---|
| 3036 | the optional \code{ordering\_file} |
---|
| 3037 | (see the function \code{urs2sts} in Section \ref{sec:basicfileconversions}). |
---|
| 3038 | |
---|
[5555] | 3039 | \item \code{elevation}, a Numeric array storing bed-elevations |
---|
| 3040 | % Refer here to the example to be provided in describing the simple example |
---|
| 3041 | \item \code{time}, a Numeric array containing times for model |
---|
| 3042 | evaluation |
---|
| 3043 | \end{itemize} |
---|
| 3044 | The only difference between the STS format and the SWW format is the former does not contain a list specifying the points at the vertices of each of the triangles (\code{volumes}). Consequenlty information (arrays) stored within an STS file such as \code{elevation} can be accessed in exactly the same way as it would be extraced from an SWW file. |
---|
| 3045 | |
---|
[4123] | 3046 | A TMS file is used to store time series data that is independent of |
---|
| 3047 | position. |
---|
| 3048 | |
---|
| 3049 | |
---|
| 3050 | \subsection{Mesh File Formats} |
---|
| 3051 | |
---|
| 3052 | A mesh file is a file that has a specific format suited to |
---|
| 3053 | triangular meshes and their outlines. A mesh file can have one of |
---|
| 3054 | two formats: it can be either a TSH file, which is an ASCII file, or |
---|
| 3055 | an MSH file, which is a NetCDF file. A mesh file can be generated |
---|
| 3056 | from the function \function{create\_mesh\_from\_regions} (see |
---|
| 3057 | Section \ref{sec:meshgeneration}) and used to initialise a domain. |
---|
| 3058 | |
---|
| 3059 | A mesh file can define the outline of the mesh---the vertices and |
---|
| 3060 | line segments that enclose the region in which the mesh is |
---|
| 3061 | created---and the triangular mesh itself, which is specified by |
---|
| 3062 | listing the triangles and their vertices, and the segments, which |
---|
| 3063 | are those sides of the triangles that are associated with boundary |
---|
| 3064 | conditions. |
---|
| 3065 | |
---|
| 3066 | In addition, a mesh file may contain `holes' and/or `regions'. A |
---|
| 3067 | hole represents an area where no mesh is to be created, while a |
---|
| 3068 | region is a labelled area used for defining properties of a mesh, |
---|
| 3069 | such as friction values. A hole or region is specified by a point |
---|
| 3070 | and bounded by a number of segments that enclose that point. |
---|
| 3071 | |
---|
| 3072 | A mesh file can also contain a georeference, which describes an |
---|
| 3073 | offset to be applied to $x$ and $y$ values---eg to the vertices. |
---|
| 3074 | |
---|
| 3075 | |
---|
| 3076 | \subsection{Formats for Storing Arbitrary Points and Attributes} |
---|
| 3077 | |
---|
| 3078 | |
---|
[4472] | 3079 | A CSV/TXT file is used to store data representing |
---|
| 3080 | arbitrary numerical attributes associated with a set of points. |
---|
| 3081 | |
---|
| 3082 | The format for an CSV/TXT file is:\\ |
---|
| 3083 | %\begin{verbatim} |
---|
| 3084 | |
---|
| 3085 | first line: \code{[column names]}\\ |
---|
| 3086 | other lines: \code{[x value], [y value], [attributes]}\\ |
---|
| 3087 | |
---|
| 3088 | for example:\\ |
---|
| 3089 | \code{x, y, elevation, friction}\\ |
---|
| 3090 | \code{0.6, 0.7, 4.9, 0.3}\\ |
---|
| 3091 | \code{1.9, 2.8, 5, 0.3}\\ |
---|
| 3092 | \code{2.7, 2.4, 5.2, 0.3} |
---|
| 3093 | |
---|
| 3094 | The delimiter is a comma. The first two columns are assumed to |
---|
[4953] | 3095 | be x, y coordinates. |
---|
[4472] | 3096 | |
---|
| 3097 | |
---|
[4662] | 3098 | A PTS file is a NetCDF representation of the data held in an points CSV |
---|
[4123] | 3099 | file. If the data is associated with a set of $N$ points, then the |
---|
| 3100 | data is stored using an $N \times 2$ Numeric array of float |
---|
| 3101 | variables for the points and an $N \times 1$ Numeric array for each |
---|
| 3102 | attribute. |
---|
| 3103 | |
---|
| 3104 | %\end{verbatim} |
---|
| 3105 | |
---|
| 3106 | \subsection{ArcView Formats} |
---|
| 3107 | |
---|
| 3108 | Files of the three formats ASC, PRJ and ERS are all associated with |
---|
| 3109 | data from ArcView. |
---|
| 3110 | |
---|
| 3111 | An ASC file is an ASCII representation of DEM output from ArcView. |
---|
| 3112 | It contains a header with the following format: |
---|
| 3113 | |
---|
| 3114 | \begin{tabular}{l l} |
---|
| 3115 | \code{ncols} & \code{753}\\ |
---|
| 3116 | \code{nrows} & \code{766}\\ |
---|
| 3117 | \code{xllcorner} & \code{314036.58727982}\\ |
---|
| 3118 | \code{yllcorner} & \code{6224951.2960092}\\ |
---|
| 3119 | \code{cellsize} & \code{100}\\ |
---|
| 3120 | \code{NODATA_value} & \code{-9999} |
---|
| 3121 | \end{tabular} |
---|
| 3122 | |
---|
| 3123 | The remainder of the file contains the elevation data for each grid point |
---|
| 3124 | in the grid defined by the above information. |
---|
| 3125 | |
---|
| 3126 | A PRJ file is an ArcView file used in conjunction with an ASC file |
---|
| 3127 | to represent metadata for a DEM. |
---|
| 3128 | |
---|
| 3129 | |
---|
| 3130 | \subsection{DEM Format} |
---|
| 3131 | |
---|
[5619] | 3132 | A DEM file in \anuga is a NetCDF representation of regular DEM data. |
---|
[4123] | 3133 | |
---|
| 3134 | |
---|
| 3135 | \subsection{Other Formats} |
---|
| 3136 | |
---|
| 3137 | |
---|
| 3138 | |
---|
| 3139 | |
---|
| 3140 | \subsection{Basic File Conversions} |
---|
| 3141 | \label{sec:basicfileconversions} |
---|
| 3142 | |
---|
| 3143 | \begin{funcdesc}{sww2dem}{basename_in, basename_out = None, |
---|
| 3144 | quantity = None, |
---|
| 3145 | timestep = None, |
---|
| 3146 | reduction = None, |
---|
| 3147 | cellsize = 10, |
---|
[5632] | 3148 | number_of_decimal_places = None, |
---|
[4123] | 3149 | NODATA_value = -9999, |
---|
| 3150 | easting_min = None, |
---|
| 3151 | easting_max = None, |
---|
| 3152 | northing_min = None, |
---|
| 3153 | northing_max = None, |
---|
| 3154 | expand_search = False, |
---|
| 3155 | verbose = False, |
---|
| 3156 | origin = None, |
---|
| 3157 | datum = 'WGS84', |
---|
| 3158 | format = 'ers'} |
---|
| 3159 | Module: \module{shallow\_water.data\_manager} |
---|
| 3160 | |
---|
| 3161 | Takes data from an SWW file \code{basename_in} and converts it to DEM format (ASC or |
---|
[5744] | 3162 | ERS) of a desired grid size \code{cellsize} in metres. The user can select how |
---|
| 3163 | many the decimal places the output data can be written to using \code{number_of_decimal_places}, |
---|
[5632] | 3164 | with the default being 3. |
---|
| 3165 | The easting and northing values are used if the user wished to determine the output |
---|
| 3166 | within a specified rectangular area. The \code{reduction} input refers to a function |
---|
[4123] | 3167 | to reduce the quantities over all time step of the SWW file, example, maximum. |
---|
| 3168 | \end{funcdesc} |
---|
| 3169 | |
---|
| 3170 | |
---|
| 3171 | \begin{funcdesc}{dem2pts}{basename_in, basename_out=None, |
---|
| 3172 | easting_min=None, easting_max=None, |
---|
| 3173 | northing_min=None, northing_max=None, |
---|
| 3174 | use_cache=False, verbose=False} |
---|
| 3175 | Module: \module{shallow\_water.data\_manager} |
---|
| 3176 | |
---|
| 3177 | Takes DEM data (a NetCDF file representation of data from a regular Digital |
---|
| 3178 | Elevation Model) and converts it to PTS format. |
---|
| 3179 | \end{funcdesc} |
---|
| 3180 | |
---|
[5744] | 3181 | \begin{funcdesc}{urs2sts}{basename_in, basename_out=None, |
---|
| 3182 | weights=None, verbose=False, |
---|
| 3183 | origin=None,mean_stage=0.0, |
---|
[5555] | 3184 | zscale=1.0, ordering_filename=None} |
---|
| 3185 | Module: \module{shallow\_water.data\_manager} |
---|
[4123] | 3186 | |
---|
[5745] | 3187 | Takes urs data in (timeseries data in mux2 format) and converts it to STS format. The optional filename \code{ordering\_filename} specifies the permutation of indices of points to select along with their longitudes and latitudes. This permutation will also be |
---|
| 3188 | stored in the STS file. If absent, all points are taken and the permutation will be trivial, i.e. $0, 1, \ldots, N-1$, where $N$ is the total number of points stored. |
---|
[5555] | 3189 | \end{funcdesc} |
---|
[4377] | 3190 | |
---|
[5555] | 3191 | |
---|
| 3192 | |
---|
[4265] | 3193 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3194 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
[4123] | 3195 | |
---|
[4265] | 3196 | \chapter{\anuga mathematical background} |
---|
| 3197 | \label{cd:mathematical background} |
---|
| 3198 | |
---|
| 3199 | \section{Introduction} |
---|
| 3200 | |
---|
| 3201 | This chapter outlines the mathematics underpinning \anuga. |
---|
| 3202 | |
---|
[4377] | 3203 | |
---|
| 3204 | |
---|
[4265] | 3205 | \section{Model} |
---|
| 3206 | \label{sec:model} |
---|
| 3207 | |
---|
| 3208 | The shallow water wave equations are a system of differential |
---|
| 3209 | conservation equations which describe the flow of a thin layer of |
---|
| 3210 | fluid over terrain. The form of the equations are: |
---|
| 3211 | \[ |
---|
| 3212 | \frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial |
---|
| 3213 | x}+\frac{\partial \GG}{\partial y}=\SSS |
---|
| 3214 | \] |
---|
| 3215 | where $\UU=\left[ {{\begin{array}{*{20}c} |
---|
| 3216 | h & {uh} & {vh} \\ |
---|
| 3217 | \end{array} }} \right]^T$ is the vector of conserved quantities; water depth |
---|
| 3218 | $h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities |
---|
| 3219 | entering the system are bed elevation $z$ and stage (absolute water |
---|
| 3220 | level) $w$, where the relation $w = z + h$ holds true at all times. |
---|
| 3221 | The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given |
---|
| 3222 | by |
---|
| 3223 | \[ |
---|
| 3224 | \EE=\left[ {{\begin{array}{*{20}c} |
---|
| 3225 | {uh} \hfill \\ |
---|
| 3226 | {u^2h+gh^2/2} \hfill \\ |
---|
| 3227 | {uvh} \hfill \\ |
---|
| 3228 | \end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c} |
---|
| 3229 | {vh} \hfill \\ |
---|
| 3230 | {vuh} \hfill \\ |
---|
| 3231 | {v^2h+gh^2/2} \hfill \\ |
---|
| 3232 | \end{array} }} \right] |
---|
| 3233 | \] |
---|
| 3234 | and the source term (which includes gravity and friction) is given |
---|
| 3235 | by |
---|
| 3236 | \[ |
---|
| 3237 | \SSS=\left[ {{\begin{array}{*{20}c} |
---|
| 3238 | 0 \hfill \\ |
---|
| 3239 | -{gh(z_{x} + S_{fx} )} \hfill \\ |
---|
| 3240 | -{gh(z_{y} + S_{fy} )} \hfill \\ |
---|
| 3241 | \end{array} }} \right] |
---|
| 3242 | \] |
---|
| 3243 | where $S_f$ is the bed friction. The friction term is modelled using |
---|
| 3244 | Manning's resistance law |
---|
| 3245 | \[ |
---|
| 3246 | S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy} |
---|
| 3247 | =\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}} |
---|
| 3248 | \] |
---|
| 3249 | in which $\eta$ is the Manning resistance coefficient. |
---|
[5566] | 3250 | The model does not currently include consideration of kinematic viscosity. |
---|
[4265] | 3251 | |
---|
[4377] | 3252 | As demonstrated in our papers, \cite{ZR1999,nielsen2005} these |
---|
[5744] | 3253 | equations and their implementation in \anuga provide a reliable |
---|
| 3254 | model of general flows associated with inundation such as dam breaks |
---|
[5566] | 3255 | and tsunamis. |
---|
[4265] | 3256 | |
---|
| 3257 | \section{Finite Volume Method} |
---|
| 3258 | \label{sec:fvm} |
---|
| 3259 | |
---|
| 3260 | We use a finite-volume method for solving the shallow water wave |
---|
[4377] | 3261 | equations \cite{ZR1999}. The study area is represented by a mesh of |
---|
[4265] | 3262 | triangular cells as in Figure~\ref{fig:mesh} in which the conserved |
---|
| 3263 | quantities of water depth $h$, and horizontal momentum $(uh, vh)$, |
---|
| 3264 | in each volume are to be determined. The size of the triangles may |
---|
| 3265 | be varied within the mesh to allow greater resolution in regions of |
---|
| 3266 | particular interest. |
---|
| 3267 | |
---|
| 3268 | \begin{figure} |
---|
| 3269 | \begin{center} |
---|
[4377] | 3270 | \includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-five} |
---|
[4265] | 3271 | \caption{Triangular mesh used in our finite volume method. Conserved |
---|
| 3272 | quantities $h$, $uh$ and $vh$ are associated with the centroid of |
---|
| 3273 | each triangular cell.} \label{fig:mesh} |
---|
| 3274 | \end{center} |
---|
| 3275 | \end{figure} |
---|
| 3276 | |
---|
| 3277 | The equations constituting the finite-volume method are obtained by |
---|
| 3278 | integrating the differential conservation equations over each |
---|
| 3279 | triangular cell of the mesh. Introducing some notation we use $i$ to |
---|
| 3280 | refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the |
---|
| 3281 | set of indices referring to the cells neighbouring the $i$th cell. |
---|
| 3282 | Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is |
---|
| 3283 | the length of the edge between the $i$th and $j$th cells. |
---|
| 3284 | |
---|
| 3285 | By applying the divergence theorem we obtain for each volume an |
---|
| 3286 | equation which describes the rate of change of the average of the |
---|
| 3287 | conserved quantities within each cell, in terms of the fluxes across |
---|
| 3288 | the edges of the cells and the effect of the source terms. In |
---|
| 3289 | particular, rate equations associated with each cell have the form |
---|
| 3290 | $$ |
---|
| 3291 | \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i |
---|
| 3292 | $$ |
---|
| 3293 | where |
---|
| 3294 | \begin{itemize} |
---|
| 3295 | \item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell, |
---|
| 3296 | \item $\SSS_i$ is the source term associated with the $i$th cell, |
---|
| 3297 | and |
---|
| 3298 | \item $\HH_{ij}$ is the outward normal flux of |
---|
| 3299 | material across the \textit{ij}th edge. |
---|
| 3300 | \end{itemize} |
---|
| 3301 | |
---|
| 3302 | |
---|
| 3303 | %\item $l_{ij}$ is the length of the edge between the $i$th and $j$th |
---|
| 3304 | %cells |
---|
| 3305 | %\item $m_{ij}$ is the midpoint of |
---|
| 3306 | %the \textit{ij}th edge, |
---|
| 3307 | %\item |
---|
| 3308 | %$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing |
---|
| 3309 | %normal along the \textit{ij}th edge, and The |
---|
| 3310 | |
---|
| 3311 | The flux $\HH_{ij}$ is evaluated using a numerical flux function |
---|
| 3312 | $\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow |
---|
| 3313 | water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$ |
---|
| 3314 | $$ |
---|
| 3315 | H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 . |
---|
| 3316 | $$ |
---|
| 3317 | |
---|
| 3318 | Then |
---|
| 3319 | $$ |
---|
| 3320 | \HH_{ij} = \HH(\UU_i(m_{ij}), |
---|
| 3321 | \UU_j(m_{ij}); \mathbf{n}_{ij}) |
---|
| 3322 | $$ |
---|
| 3323 | where $m_{ij}$ is the midpoint of the \textit{ij}th edge and |
---|
| 3324 | $\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the |
---|
| 3325 | \textit{ij}th edge. The function $\UU_i(x)$ for $x \in |
---|
| 3326 | T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and |
---|
| 3327 | neighbouring cells. |
---|
| 3328 | |
---|
| 3329 | We use a second order reconstruction to produce a piece-wise linear |
---|
| 3330 | function construction of the conserved quantities for all $x \in |
---|
| 3331 | T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This |
---|
| 3332 | function is allowed to be discontinuous across the edges of the |
---|
| 3333 | cells, but the slope of this function is limited to avoid |
---|
| 3334 | artificially introduced oscillations. |
---|
| 3335 | |
---|
[4377] | 3336 | Godunov's method (see \cite{Toro1992}) involves calculating the |
---|
[4265] | 3337 | numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly |
---|
| 3338 | solving the corresponding one dimensional Riemann problem normal to |
---|
| 3339 | the edge. We use the central-upwind scheme of \cite{KurNP2001} to |
---|
| 3340 | calculate an approximation of the flux across each edge. |
---|
| 3341 | |
---|
| 3342 | \begin{figure} |
---|
| 3343 | \begin{center} |
---|
[4377] | 3344 | \includegraphics[width=8.0cm,keepaspectratio=true]{graphics/step-reconstruct} |
---|
[4265] | 3345 | \caption{From the values of the conserved quantities at the centroid |
---|
| 3346 | of the cell and its neighbouring cells, a discontinuous piecewise |
---|
| 3347 | linear reconstruction of the conserved quantities is obtained.} |
---|
| 3348 | \label{fig:mesh:reconstruct} |
---|
| 3349 | \end{center} |
---|
| 3350 | \end{figure} |
---|
| 3351 | |
---|
| 3352 | In the computations presented in this paper we use an explicit Euler |
---|
| 3353 | time stepping method with variable timestepping adapted to the |
---|
[5618] | 3354 | observed CFL condition: |
---|
[4265] | 3355 | |
---|
[5744] | 3356 | \begin{equation} |
---|
[5621] | 3357 | \Delta t = \min_{k,i=[0,1,2]} \min \left( \frac{r_k}{v_{k,i}}, \frac{r_{n_{k,i}}}{v_{k,i}} \right ) |
---|
[5618] | 3358 | \label{eq:CFL condition} |
---|
[5744] | 3359 | \end{equation} |
---|
| 3360 | where $r_k$ is the radius of the $k$'th triangle and $v_{k,i}$ is the maximal velocity across |
---|
| 3361 | edge joining triangle $k$ and it's $i$'th neighbour, triangle $n_{k,i}$, as calculated by the |
---|
| 3362 | numerical flux function |
---|
| 3363 | using the central upwind scheme of \cite{KurNP2001}. The symbol $r_{n_{k,i}}$ denotes the radius |
---|
| 3364 | of the $i$'th neighbour of triangle $k$. The radii are calculated as radii of the inscribed circles |
---|
[5620] | 3365 | of each triangle. |
---|
[4265] | 3366 | |
---|
| 3367 | \section{Flux limiting} |
---|
| 3368 | |
---|
[4377] | 3369 | The shallow water equations are solved numerically using a |
---|
[4265] | 3370 | finite volume method on unstructured triangular grid. |
---|
[4377] | 3371 | The upwind central scheme due to Kurganov and Petrova is used as an |
---|
[4265] | 3372 | approximate Riemann solver for the computation of inviscid flux functions. |
---|
[4377] | 3373 | This makes it possible to handle discontinuous solutions. |
---|
[4265] | 3374 | |
---|
[4377] | 3375 | To alleviate the problems associated with numerical instabilities due to |
---|
[4265] | 3376 | small water depths near a wet/dry boundary we employ a new flux limiter that |
---|
| 3377 | ensures that unphysical fluxes are never encounted. |
---|
| 3378 | |
---|
| 3379 | |
---|
[4377] | 3380 | Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction, |
---|
[4265] | 3381 | $w$ the absolute water level (stage) and |
---|
[4377] | 3382 | $z$ the bed elevation. The latter are assumed to be relative to the |
---|
| 3383 | same height datum. |
---|
| 3384 | The conserved quantities tracked by ANUGA are momentum in the |
---|
| 3385 | $x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$) |
---|
[4265] | 3386 | and depth ($h = w-z$). |
---|
| 3387 | |
---|
[4377] | 3388 | The flux calculation requires access to the velocity vector $(u, v)$ |
---|
[4265] | 3389 | where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively. |
---|
[4377] | 3390 | In the presence of very small water depths, these calculations become |
---|
[4265] | 3391 | numerically unreliable and will typically cause unphysical speeds. |
---|
| 3392 | |
---|
[4377] | 3393 | We have employed a flux limiter which replaces the calculations above with |
---|
[4265] | 3394 | the limited approximations. |
---|
| 3395 | \begin{equation} |
---|
[4377] | 3396 | \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h}, |
---|
[4265] | 3397 | \end{equation} |
---|
[4377] | 3398 | where $h_0$ is a regularisation parameter that controls the minimal |
---|
[4265] | 3399 | magnitude of the denominator. Taking the limits we have for $\hat{u}$ |
---|
| 3400 | \[ |
---|
[4377] | 3401 | \lim_{h \rightarrow 0} \hat{u} = |
---|
[4265] | 3402 | \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0 |
---|
| 3403 | \] |
---|
[4377] | 3404 | and |
---|
[4265] | 3405 | \[ |
---|
[4377] | 3406 | \lim_{h \rightarrow \infty} \hat{u} = |
---|
[4265] | 3407 | \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u |
---|
| 3408 | \] |
---|
| 3409 | with similar results for $\hat{v}$. |
---|
| 3410 | |
---|
| 3411 | The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator) |
---|
| 3412 | \[ |
---|
| 3413 | 1 - h_0/h^2 = 0 |
---|
| 3414 | \] |
---|
| 3415 | or |
---|
| 3416 | \[ |
---|
| 3417 | h_0 = h^2 |
---|
| 3418 | \] |
---|
| 3419 | |
---|
| 3420 | |
---|
[4377] | 3421 | ANUGA has a global parameter $H_0$ that controls the minimal depth which |
---|
[4265] | 3422 | is considered in the various equations. This parameter is typically set to |
---|
| 3423 | $10^{-3}$. Setting |
---|
| 3424 | \[ |
---|
| 3425 | h_0 = H_0^2 |
---|
| 3426 | \] |
---|
| 3427 | provides a reasonable balance between accurracy and stability. In fact, |
---|
| 3428 | setting $h=H_0$ will scale the predicted speed by a factor of $0.5$: |
---|
| 3429 | \[ |
---|
| 3430 | \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0} |
---|
| 3431 | \] |
---|
[4377] | 3432 | In general, for multiples of the minimal depth $N H_0$ one obtains |
---|
[4265] | 3433 | \[ |
---|
[4377] | 3434 | \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} = |
---|
[4265] | 3435 | \frac{\mu}{H_0 (1 + 1/N^2)} |
---|
| 3436 | \] |
---|
[4377] | 3437 | which converges quadratically to the true value with the multiple N. |
---|
[4265] | 3438 | |
---|
| 3439 | |
---|
| 3440 | %The developed numerical model has been applied to several test cases as well as to real flows. Numerical tests prove the robustness and accuracy of the model. |
---|
| 3441 | |
---|
| 3442 | |
---|
| 3443 | |
---|
| 3444 | |
---|
| 3445 | |
---|
| 3446 | \section{Slope limiting} |
---|
| 3447 | A multidimensional slope-limiting technique is employed to achieve second-order spatial accuracy and to prevent spurious oscillations. This is using the MinMod limiter and is documented elsewhere. |
---|
| 3448 | |
---|
| 3449 | However close to the bed, the limiter must ensure that no negative depths occur. On the other hand, in deep water, the bed topography should be ignored for the purpose of the limiter. |
---|
| 3450 | |
---|
| 3451 | |
---|
| 3452 | Let $w, z, h$ be the stage, bed elevation and depth at the centroid and |
---|
| 3453 | let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$. |
---|
| 3454 | Define the minimal depth across all vertices as $\hmin$ as |
---|
| 3455 | \[ |
---|
| 3456 | \hmin = \min_i h_i |
---|
| 3457 | \] |
---|
| 3458 | |
---|
[4377] | 3459 | Let $\tilde{w_i}$ be the stage obtained from a gradient limiter |
---|
[4265] | 3460 | limiting on stage only. The corresponding depth is then defined as |
---|
| 3461 | \[ |
---|
| 3462 | \tilde{h_i} = \tilde{w_i} - z_i |
---|
| 3463 | \] |
---|
[4377] | 3464 | We would use this limiter in deep water which we will define (somewhat boldly) |
---|
[4265] | 3465 | as |
---|
| 3466 | \[ |
---|
| 3467 | \hmin \ge \epsilon |
---|
| 3468 | \] |
---|
| 3469 | |
---|
| 3470 | |
---|
[4377] | 3471 | Similarly, let $\bar{w_i}$ be the stage obtained from a gradient |
---|
[4265] | 3472 | limiter limiting on depth respecting the bed slope. |
---|
| 3473 | The corresponding depth is defined as |
---|
| 3474 | \[ |
---|
| 3475 | \bar{h_i} = \bar{w_i} - z_i |
---|
| 3476 | \] |
---|
| 3477 | |
---|
| 3478 | |
---|
| 3479 | We introduce the concept of a balanced stage $w_i$ which is obtained as |
---|
| 3480 | the linear combination |
---|
| 3481 | |
---|
| 3482 | \[ |
---|
| 3483 | w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i} |
---|
| 3484 | \] |
---|
| 3485 | or |
---|
| 3486 | \[ |
---|
| 3487 | w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} |
---|
| 3488 | \] |
---|
| 3489 | where $\alpha \in [0, 1]$. |
---|
| 3490 | |
---|
[4377] | 3491 | Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope |
---|
| 3492 | is ignored we have immediately that |
---|
[4265] | 3493 | \[ |
---|
| 3494 | \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0 |
---|
| 3495 | \] |
---|
| 3496 | %where the maximal bed elevation range $dz$ is defined as |
---|
| 3497 | %\[ |
---|
| 3498 | % dz = \max_i |z_i - z| |
---|
| 3499 | %\] |
---|
| 3500 | |
---|
| 3501 | If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that |
---|
| 3502 | no negative depths occur. Formally, we will require that |
---|
| 3503 | \[ |
---|
| 3504 | \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i |
---|
| 3505 | \] |
---|
[4377] | 3506 | or |
---|
| 3507 | \begin{equation} |
---|
[4265] | 3508 | \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i |
---|
| 3509 | \label{eq:limiter bound} |
---|
[4377] | 3510 | \end{equation} |
---|
[4265] | 3511 | |
---|
| 3512 | There are two cases: |
---|
[4377] | 3513 | \begin{enumerate} |
---|
| 3514 | \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage) |
---|
| 3515 | vertex is at least as far away from the bed than the shallow water |
---|
[4265] | 3516 | (limited using depth). In this case we won't need any contribution from |
---|
[4404] | 3517 | $\bar{h_i}$ and can accept any $\alpha$. |
---|
[4377] | 3518 | |
---|
[4265] | 3519 | E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to |
---|
| 3520 | \[ |
---|
[4377] | 3521 | \tilde{h_i} > \epsilon |
---|
[4265] | 3522 | \] |
---|
| 3523 | whereas $\alpha=0$ yields |
---|
| 3524 | \[ |
---|
[4377] | 3525 | \bar{h_i} > \epsilon |
---|
[4265] | 3526 | \] |
---|
| 3527 | all well and good. |
---|
[4377] | 3528 | \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is |
---|
| 3529 | closer to the bed than the shallow water vertex or even below the bed. |
---|
[4404] | 3530 | In this case we need to find an $\alpha$ that will ensure a positive depth. |
---|
[4377] | 3531 | Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one |
---|
[4265] | 3532 | obtains the bound |
---|
| 3533 | \[ |
---|
| 3534 | \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i |
---|
[4377] | 3535 | \] |
---|
| 3536 | \end{enumerate} |
---|
[4265] | 3537 | |
---|
[4377] | 3538 | Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one |
---|
[4265] | 3539 | arrives at the definition |
---|
| 3540 | \[ |
---|
| 3541 | \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}} |
---|
| 3542 | \] |
---|
| 3543 | which will guarantee that no vertex 'cuts' through the bed. Finally, should |
---|
[4377] | 3544 | $\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting |
---|
[4404] | 3545 | $\alpha=0$ and similarly capping $\alpha$ at 1 just in case. |
---|
[4265] | 3546 | |
---|
| 3547 | %Furthermore, |
---|
[4377] | 3548 | %dropping the $\epsilon$ ensures that alpha is always positive and also |
---|
[4265] | 3549 | %provides a numerical safety {??) |
---|
| 3550 | |
---|
| 3551 | |
---|
| 3552 | |
---|
[4377] | 3553 | |
---|
| 3554 | |
---|
[4123] | 3555 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3556 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3557 | |
---|
| 3558 | \chapter{Basic \anuga Assumptions} |
---|
| 3559 | |
---|
| 3560 | |
---|
| 3561 | Physical model time cannot be earlier than 1 Jan 1970 00:00:00. |
---|
| 3562 | If one wished to recreate scenarios prior to that date it must be done |
---|
| 3563 | using some relative time (e.g. 0). |
---|
| 3564 | |
---|
| 3565 | |
---|
| 3566 | All spatial data relates to the WGS84 datum (or GDA94) and has been |
---|
| 3567 | projected into UTM with false easting of 500000 and false northing of |
---|
| 3568 | 1000000 on the southern hemisphere (0 on the northern). |
---|
| 3569 | |
---|
[4953] | 3570 | It is assumed that all computations take place within one UTM zone and |
---|
| 3571 | all locations must consequently be specified in Cartesian coordinates |
---|
[4543] | 3572 | (eastings, northings) or (x,y) where the unit is metres. |
---|
[4123] | 3573 | |
---|
| 3574 | DEMs, meshes and boundary conditions can have different origins within |
---|
| 3575 | one UTM zone. However, the computation will use that of the mesh for |
---|
| 3576 | numerical stability. |
---|
| 3577 | |
---|
| 3578 | When generating a mesh it is assumed that polygons do not cross. |
---|
| 3579 | Having polygons tht cross can cause the mesh generation to fail or bad |
---|
| 3580 | meshes being produced. |
---|
| 3581 | |
---|
| 3582 | |
---|
| 3583 | %OLD |
---|
| 3584 | %The dataflow is: (See data_manager.py and from scenarios) |
---|
| 3585 | % |
---|
| 3586 | % |
---|
| 3587 | %Simulation scenarios |
---|
| 3588 | %--------------------% |
---|
| 3589 | %% |
---|
| 3590 | % |
---|
| 3591 | %Sub directories contain scrips and derived files for each simulation. |
---|
| 3592 | %The directory ../source_data contains large source files such as |
---|
| 3593 | %DEMs provided externally as well as MOST tsunami simulations to be used |
---|
| 3594 | %as boundary conditions. |
---|
| 3595 | % |
---|
| 3596 | %Manual steps are: |
---|
| 3597 | % Creation of DEMs from argcview (.asc + .prj) |
---|
| 3598 | % Creation of mesh from pmesh (.tsh) |
---|
| 3599 | % Creation of tsunami simulations from MOST (.nc) |
---|
| 3600 | %% |
---|
| 3601 | % |
---|
| 3602 | %Typical scripted steps are% |
---|
| 3603 | % |
---|
| 3604 | % prepare_dem.py: Convert asc and prj files supplied by arcview to |
---|
| 3605 | % native dem and pts formats% |
---|
| 3606 | % |
---|
| 3607 | % prepare_pts.py: Convert netcdf output from MOST to an sww file suitable |
---|
| 3608 | % as boundary condition% |
---|
| 3609 | % |
---|
| 3610 | % prepare_mesh.py: Merge DEM (pts) and mesh (tsh) using least squares |
---|
| 3611 | % smoothing. The outputs are tsh files with elevation data.% |
---|
| 3612 | % |
---|
| 3613 | % run_simulation.py: Use the above together with various parameters to |
---|
| 3614 | % run inundation simulation. |
---|
| 3615 | |
---|
| 3616 | |
---|
| 3617 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3618 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 3619 | |
---|
| 3620 | \appendix |
---|
| 3621 | |
---|
| 3622 | \chapter{Supporting Tools} |
---|
| 3623 | \label{ch:supportingtools} |
---|
| 3624 | |
---|
| 3625 | This section describes a number of supporting tools, supplied with \anuga, that offer a |
---|
| 3626 | variety of types of functionality and enhance the basic capabilities of \anuga. |
---|
| 3627 | |
---|
| 3628 | \section{caching} |
---|
| 3629 | \label{sec:caching} |
---|
| 3630 | |
---|
[4209] | 3631 | The \code{cache} function is used to provide supervised caching of function |
---|
[4123] | 3632 | results. A Python function call of the form |
---|
| 3633 | |
---|
| 3634 | {\small \begin{verbatim} |
---|
| 3635 | result = func(arg1,...,argn) |
---|
| 3636 | \end{verbatim}} |
---|
| 3637 | |
---|
| 3638 | can be replaced by |
---|
| 3639 | |
---|
| 3640 | {\small \begin{verbatim} |
---|
| 3641 | from caching import cache |
---|
| 3642 | result = cache(func,(arg1,...,argn)) |
---|
| 3643 | \end{verbatim}} |
---|
| 3644 | |
---|
| 3645 | which returns the same output but reuses cached |
---|
| 3646 | results if the function has been computed previously in the same context. |
---|
| 3647 | \code{result} and the arguments can be simple types, tuples, list, dictionaries or |
---|
| 3648 | objects, but not unhashable types such as functions or open file objects. |
---|
| 3649 | The function \code{func} may be a member function of an object or a module. |
---|
| 3650 | |
---|
| 3651 | This type of caching is particularly useful for computationally intensive |
---|
| 3652 | functions with few frequently used combinations of input arguments. Note that |
---|
| 3653 | if the inputs or output are very large caching may not save time because |
---|
| 3654 | disc access may dominate the execution time. |
---|
| 3655 | |
---|
| 3656 | If the function definition changes after a result has been cached, this will be |
---|
| 3657 | detected by examining the functions \code{bytecode (co\_code, co\_consts, |
---|
| 3658 | func\_defaults, co\_argcount)} and the function will be recomputed. |
---|
| 3659 | However, caching will not detect changes in modules used by \code{func}. |
---|
| 3660 | In this case cache must be cleared manually. |
---|
| 3661 | |
---|
| 3662 | Options are set by means of the function \code{set\_option(key, value)}, |
---|
| 3663 | where \code{key} is a key associated with a |
---|
| 3664 | Python dictionary \code{options}. This dictionary stores settings such as the name of |
---|
| 3665 | the directory used, the maximum |
---|
| 3666 | number of cached files allowed, and so on. |
---|
| 3667 | |
---|
| 3668 | The \code{cache} function allows the user also to specify a list of dependent files. If any of these |
---|
| 3669 | have been changed, the function is recomputed and the results stored again. |
---|
| 3670 | |
---|
| 3671 | %Other features include support for compression and a capability to \ldots |
---|
| 3672 | |
---|
| 3673 | |
---|
| 3674 | \textbf{USAGE:} \nopagebreak |
---|
| 3675 | |
---|
| 3676 | {\small \begin{verbatim} |
---|
| 3677 | result = cache(func, args, kwargs, dependencies, cachedir, verbose, |
---|
| 3678 | compression, evaluate, test, return_filename) |
---|
| 3679 | \end{verbatim}} |
---|
| 3680 | |
---|
| 3681 | |
---|
| 3682 | \section{ANUGA viewer - animate} |
---|
| 3683 | \label{sec:animate} |
---|
| 3684 | The output generated by \anuga may be viewed by |
---|
| 3685 | means of the visualisation tool \code{animate}, which takes the |
---|
| 3686 | \code{SWW} file output by \anuga and creates a visual representation |
---|
| 3687 | of the data. Examples may be seen in Figures \ref{fig:runupstart} |
---|
| 3688 | and \ref{fig:runup2}. To view an \code{SWW} file with |
---|
| 3689 | \code{animate} in the Windows environment, you can simply drag the |
---|
| 3690 | icon representing the file over an icon on the desktop for the |
---|
| 3691 | \code{animate} executable file (or a shortcut to it), or set up a |
---|
| 3692 | file association to make files with the extension \code{.sww} open |
---|
| 3693 | with \code{animate}. Alternatively, you can operate \code{animate} |
---|
| 3694 | from the command line, in both Windows and Linux environments. |
---|
| 3695 | |
---|
| 3696 | On successful operation, you will see an interactive moving-picture |
---|
| 3697 | display. You can use keys and the mouse to slow down, speed up or |
---|
| 3698 | stop the display, change the viewing position or carry out a number |
---|
| 3699 | of other simple operations. Help is also displayed when you press |
---|
| 3700 | the \code{h} key. |
---|
| 3701 | |
---|
| 3702 | The main keys operating the interactive screen are:\\ |
---|
| 3703 | |
---|
| 3704 | \begin{center} |
---|
| 3705 | \begin{tabular}{|ll|} \hline |
---|
| 3706 | |
---|
| 3707 | \code{w} & toggle wireframe \\ |
---|
| 3708 | |
---|
| 3709 | space bar & start/stop\\ |
---|
| 3710 | |
---|
| 3711 | up/down arrows & increase/decrease speed\\ |
---|
| 3712 | |
---|
| 3713 | left/right arrows & direction in time \emph{(when running)}\\ |
---|
| 3714 | & step through simulation \emph{(when stopped)}\\ |
---|
| 3715 | |
---|
| 3716 | left mouse button & rotate\\ |
---|
| 3717 | |
---|
| 3718 | middle mouse button & pan\\ |
---|
| 3719 | |
---|
| 3720 | right mouse button & zoom\\ \hline |
---|
| 3721 | |
---|
| 3722 | \end{tabular} |
---|
| 3723 | \end{center} |
---|
| 3724 | |
---|
| 3725 | \vfill |
---|
| 3726 | |
---|
| 3727 | The following table describes how to operate animate from the command line: |
---|
| 3728 | |
---|
| 3729 | Usage: \code{animate [options] swwfile \ldots}\\ \nopagebreak |
---|
| 3730 | Options:\\ \nopagebreak |
---|
| 3731 | \begin{tabular}{ll} |
---|
| 3732 | \code{--display <type>} & \code{MONITOR | POWERWALL | REALITY\_CENTER |}\\ |
---|
| 3733 | & \code{HEAD\_MOUNTED\_DISPLAY}\\ |
---|
| 3734 | \code{--rgba} & Request a RGBA colour buffer visual\\ |
---|
| 3735 | \code{--stencil} & Request a stencil buffer visual\\ |
---|
| 3736 | \code{--stereo} & Use default stereo mode which is \code{ANAGLYPHIC} if not \\ |
---|
| 3737 | & overridden by environmental variable\\ |
---|
| 3738 | \code{--stereo <mode>} & \code{ANAGLYPHIC | QUAD\_BUFFER | HORIZONTAL\_SPLIT |}\\ |
---|
| 3739 | & \code{VERTICAL\_SPLIT | LEFT\_EYE | RIGHT\_EYE |}\\ |
---|
| 3740 | & \code{ON | OFF} \\ |
---|
| 3741 | \code{-alphamax <float 0-1>} & Maximum transparency clamp value\\ |
---|
| 3742 | \code{-alphamin <float 0-1>} & Transparency value at \code{hmin}\\ |
---|
| 3743 | \end{tabular} |
---|
| 3744 | |
---|
| 3745 | \begin{tabular}{ll} |
---|
| 3746 | \code{-cullangle <float angle 0-90>} & Cull triangles steeper than this value\\ |
---|
| 3747 | \code{-help} & Display this information\\ |
---|
| 3748 | \code{-hmax <float>} & Height above which transparency is set to |
---|
| 3749 | \code{alphamax}\\ |
---|
| 3750 | \end{tabular} |
---|
| 3751 | |
---|
| 3752 | \begin{tabular}{ll} |
---|
| 3753 | |
---|
| 3754 | \code{-hmin <float>} & Height below which transparency is set to |
---|
| 3755 | zero\\ |
---|
| 3756 | \end{tabular} |
---|
| 3757 | |
---|
| 3758 | \begin{tabular}{ll} |
---|
| 3759 | \code{-lightpos <float>,<float>,<float>} & $x,y,z$ of bedslope directional light ($z$ is |
---|
| 3760 | up, default is overhead)\\ |
---|
| 3761 | \end{tabular} |
---|
| 3762 | |
---|
| 3763 | \begin{tabular}{ll} |
---|
| 3764 | \code{-loop} & Repeated (looped) playback of \code{.swm} files\\ |
---|
| 3765 | |
---|
| 3766 | \end{tabular} |
---|
| 3767 | |
---|
| 3768 | \begin{tabular}{ll} |
---|
| 3769 | \code{-movie <dirname>} & Save numbered images to named directory and |
---|
| 3770 | quit\\ |
---|
| 3771 | |
---|
| 3772 | \code{-nosky} & Omit background sky\\ |
---|
| 3773 | |
---|
| 3774 | |
---|
| 3775 | \code{-scale <float>} & Vertical scale factor\\ |
---|
| 3776 | \code{-texture <file>} & Image to use for bedslope topography\\ |
---|
| 3777 | \code{-tps <rate>} & Timesteps per second\\ |
---|
| 3778 | \code{-version} & Revision number and creation (not compile) |
---|
| 3779 | date\\ |
---|
| 3780 | \end{tabular} |
---|
| 3781 | |
---|
| 3782 | \section{utilities/polygons} |
---|
| 3783 | |
---|
| 3784 | \declaremodule{standard}{utilities.polygon} |
---|
| 3785 | \refmodindex{utilities.polygon} |
---|
| 3786 | |
---|
[5088] | 3787 | \begin{classdesc}{Polygon\_function}{regions, default=0.0, geo_reference=None} |
---|
[4123] | 3788 | Module: \code{utilities.polygon} |
---|
| 3789 | |
---|
| 3790 | Creates a callable object that returns one of a specified list of values when |
---|
| 3791 | evaluated at a point \code{x, y}, depending on which polygon, from a specified list of polygons, the |
---|
| 3792 | point belongs to. The parameter \code{regions} is a list of pairs |
---|
| 3793 | \code{(P, v)}, where each \code{P} is a polygon and each \code{v} |
---|
| 3794 | is either a constant value or a function of coordinates \code{x} |
---|
| 3795 | and \code{y}, specifying the return value for a point inside \code{P}. The |
---|
[5744] | 3796 | optional parameter \code{default} may be used to specify a value |
---|
[5088] | 3797 | (or a function) |
---|
[4123] | 3798 | for a point not lying inside any of the specified polygons. When a |
---|
| 3799 | point lies in more than one polygon, the return value is taken to |
---|
| 3800 | be the value for whichever of these polygon appears later in the |
---|
| 3801 | list. |
---|
| 3802 | %FIXME (Howard): CAN x, y BE VECTORS? |
---|
[5090] | 3803 | The optional parameter geo\_reference refers to the status of points |
---|
[5744] | 3804 | that are passed into the function. Typically they will be relative to |
---|
[5090] | 3805 | some origin. In ANUGA, a typical call will take the form: |
---|
| 3806 | {\small \begin{verbatim} |
---|
[5744] | 3807 | set_quantity('elevation', |
---|
[5090] | 3808 | Polygon_function([(P1, v1), (P2, v2)], |
---|
[5091] | 3809 | default=v3, |
---|
| 3810 | geo_reference=domain.geo_reference)) |
---|
[5090] | 3811 | \end{verbatim}} |
---|
[4123] | 3812 | |
---|
[5744] | 3813 | |
---|
[4123] | 3814 | \end{classdesc} |
---|
| 3815 | |
---|
| 3816 | \begin{funcdesc}{read\_polygon}{filename} |
---|
| 3817 | Module: \code{utilities.polygon} |
---|
| 3818 | |
---|
| 3819 | Reads the specified file and returns a polygon. Each |
---|
| 3820 | line of the file must contain exactly two numbers, separated by a comma, which are interpreted |
---|
| 3821 | as coordinates of one vertex of the polygon. |
---|
| 3822 | \end{funcdesc} |
---|
| 3823 | |
---|
| 3824 | \begin{funcdesc}{populate\_polygon}{polygon, number_of_points, seed = None, exclude = None} |
---|
| 3825 | Module: \code{utilities.polygon} |
---|
| 3826 | |
---|
| 3827 | Populates the interior of the specified polygon with the specified number of points, |
---|
| 3828 | selected by means of a uniform distribution function. |
---|
| 3829 | \end{funcdesc} |
---|
| 3830 | |
---|
| 3831 | \begin{funcdesc}{point\_in\_polygon}{polygon, delta=1e-8} |
---|
| 3832 | Module: \code{utilities.polygon} |
---|
| 3833 | |
---|
| 3834 | Returns a point inside the specified polygon and close to the edge. The distance between |
---|
| 3835 | the returned point and the nearest point of the polygon is less than $\sqrt{2}$ times the |
---|
| 3836 | second argument \code{delta}, which is taken as $10^{-8}$ by default. |
---|
| 3837 | \end{funcdesc} |
---|
| 3838 | |
---|
| 3839 | \begin{funcdesc}{inside\_polygon}{points, polygon, closed = True, verbose = False} |
---|
| 3840 | Module: \code{utilities.polygon} |
---|
| 3841 | |
---|
| 3842 | Used to test whether the members of a list of points |
---|
| 3843 | are inside the specified polygon. Returns a Numeric |
---|
| 3844 | array comprising the indices of the points in the list that lie inside the polygon. |
---|
| 3845 | (If none of the points are inside, returns \code{zeros((0,), 'l')}.) |
---|
| 3846 | Points on the edges of the polygon are regarded as inside if |
---|
| 3847 | \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside. |
---|
| 3848 | \end{funcdesc} |
---|
| 3849 | |
---|
| 3850 | \begin{funcdesc}{outside\_polygon}{points, polygon, closed = True, verbose = False} |
---|
| 3851 | Module: \code{utilities.polygon} |
---|
| 3852 | |
---|
| 3853 | Exactly like \code{inside\_polygon}, but with the words `inside' and `outside' interchanged. |
---|
| 3854 | \end{funcdesc} |
---|
| 3855 | |
---|
| 3856 | \begin{funcdesc}{is\_inside\_polygon}{point, polygon, closed=True, verbose=False} |
---|
| 3857 | Module: \code{utilities.polygon} |
---|
| 3858 | |
---|
| 3859 | Returns \code{True} if \code{point} is inside \code{polygon} or |
---|
| 3860 | \code{False} otherwise. Points on the edges of the polygon are regarded as inside if |
---|
| 3861 | \code{closed} is set to \code{True} or omitted; otherwise they are regarded as outside. |
---|
| 3862 | \end{funcdesc} |
---|
| 3863 | |
---|
| 3864 | \begin{funcdesc}{is\_outside\_polygon}{point, polygon, closed=True, verbose=False} |
---|
| 3865 | Module: \code{utilities.polygon} |
---|
| 3866 | |
---|
| 3867 | Exactly like \code{is\_outside\_polygon}, but with the words `inside' and `outside' interchanged. |
---|
| 3868 | \end{funcdesc} |
---|
| 3869 | |
---|
| 3870 | \begin{funcdesc}{point\_on\_line}{x, y, x0, y0, x1, y1} |
---|
| 3871 | Module: \code{utilities.polygon} |
---|
| 3872 | |
---|
| 3873 | Returns \code{True} or \code{False}, depending on whether the point with coordinates |
---|
| 3874 | \code{x, y} is on the line passing through the points with coordinates \code{x0, y0} |
---|
| 3875 | and \code{x1, y1} (extended if necessary at either end). |
---|
| 3876 | \end{funcdesc} |
---|
| 3877 | |
---|
| 3878 | \begin{funcdesc}{separate\_points\_by\_polygon}{points, polygon, closed = True, verbose = False} |
---|
| 3879 | \indexedcode{separate\_points\_by\_polygon} |
---|
| 3880 | Module: \code{utilities.polygon} |
---|
| 3881 | |
---|
| 3882 | \end{funcdesc} |
---|
| 3883 | |
---|
| 3884 | \begin{funcdesc}{polygon\_area}{polygon} |
---|
| 3885 | Module: \code{utilities.polygon} |
---|
| 3886 | |
---|
| 3887 | Returns area of arbitrary polygon (reference http://mathworld.wolfram.com/PolygonArea.html) |
---|
| 3888 | \end{funcdesc} |
---|
| 3889 | |
---|
[5484] | 3890 | \begin{funcdesc}{plot\_polygons}{polygons, style, figname, verbose = False} |
---|
| 3891 | Module: \code{utilities.polygon} |
---|
[5744] | 3892 | |
---|
[5484] | 3893 | Plots each polygon contained in input polygon list, e.g. |
---|
| 3894 | \code{polygons = [poly1, poly2, poly3]} where \code{poly1 = [[x11,y11],[x12,y12],[x13,y13]]} |
---|
| 3895 | etc. Each polygon can be closed for plotting purposes by assigning the style type to each |
---|
| 3896 | polygon in the list, e.g. \code{style = ['line','line','line']}. The default will be a line |
---|
| 3897 | type when \code{style = None}. |
---|
| 3898 | The subsequent plot will be saved to \code{figname} or defaulted to \code{test_image.png}. |
---|
| 3899 | The function returns a list containing the minimum and maximum of \code{x} and \code{y}, |
---|
| 3900 | i.e. \code{[x_{min}, x_{max}, y_{min}, y_{max}]}. |
---|
[4123] | 3901 | \end{funcdesc} |
---|
| 3902 | |
---|
| 3903 | \section{coordinate\_transforms} |
---|
| 3904 | |
---|
| 3905 | \section{geospatial\_data} |
---|
| 3906 | \label{sec:geospatial} |
---|
| 3907 | |
---|
| 3908 | This describes a class that represents arbitrary point data in UTM |
---|
| 3909 | coordinates along with named attribute values. |
---|
| 3910 | |
---|
| 3911 | %FIXME (Ole): This gives a LaTeX error |
---|
| 3912 | %\declaremodule{standard}{geospatial_data} |
---|
| 3913 | %\refmodindex{geospatial_data} |
---|
| 3914 | |
---|
| 3915 | |
---|
| 3916 | |
---|
| 3917 | \begin{classdesc}{Geospatial\_data} |
---|
| 3918 | {data_points = None, |
---|
| 3919 | attributes = None, |
---|
| 3920 | geo_reference = None, |
---|
| 3921 | default_attribute_name = None, |
---|
| 3922 | file_name = None} |
---|
| 3923 | Module: \code{geospatial\_data} |
---|
| 3924 | |
---|
| 3925 | This class is used to store a set of data points and associated |
---|
| 3926 | attributes, allowing these to be manipulated by methods defined for |
---|
| 3927 | the class. |
---|
| 3928 | |
---|
| 3929 | The data points are specified either by reading them from a NetCDF |
---|
[4662] | 3930 | or CSV file, identified through the parameter \code{file\_name}, or |
---|
[4123] | 3931 | by providing their \code{x}- and \code{y}-coordinates in metres, |
---|
| 3932 | either as a sequence of 2-tuples of floats or as an $M \times 2$ |
---|
| 3933 | Numeric array of floats, where $M$ is the number of points. |
---|
| 3934 | Coordinates are interpreted relative to the origin specified by the |
---|
| 3935 | object \code{geo\_reference}, which contains data indicating the UTM |
---|
| 3936 | zone, easting and northing. If \code{geo\_reference} is not |
---|
| 3937 | specified, a default is used. |
---|
| 3938 | |
---|
| 3939 | Attributes are specified through the parameter \code{attributes}, |
---|
| 3940 | set either to a list or array of length $M$ or to a dictionary whose |
---|
| 3941 | keys are the attribute names and whose values are lists or arrays of |
---|
| 3942 | length $M$. One of the attributes may be specified as the default |
---|
| 3943 | attribute, by assigning its name to \code{default\_attribute\_name}. |
---|
| 3944 | If no value is specified, the default attribute is taken to be the |
---|
| 3945 | first one. |
---|
[5566] | 3946 | |
---|
[5744] | 3947 | Note that the Geospatial\_data object currently reads entire datasets |
---|
| 3948 | into memory i.e.\ no memomry blocking takes place. |
---|
[5619] | 3949 | For this we refer to the set\_quantity method which will read .pts and .csv files into \anuga using memory blocking allowing large files to be used. |
---|
[4123] | 3950 | \end{classdesc} |
---|
| 3951 | |
---|
| 3952 | |
---|
| 3953 | \begin{methoddesc}{import\_points\_file}{delimiter = None, verbose = False} |
---|
| 3954 | |
---|
| 3955 | \end{methoddesc} |
---|
| 3956 | |
---|
|
---|