Changeset 5353


Ignore:
Timestamp:
May 21, 2008, 3:49:12 PM (17 years ago)
Author:
ole
Message:

Rewrite of model section, references and other stuff

Location:
anuga_work/publications/anuga_2007
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/publications/anuga_2007/anuga-bibliography.bib

    r5341 r5353  
    4949
    5050
     51@ARTICLE{Stoker57,
     52  AUTHOR =       {J. J. Stoker},
     53  TITLE =        {Water waves: the mathematical theory with applications},
     54  JOURNAL =      {Interscience Publishers, New York},
     55  YEAR =         {1957},
     56  volume =       {},
     57  number =       {},
     58  pages =        {}
     59}
     60
     61
     62@ARTICLE{Peregrine67,
     63  AUTHOR =       {J. J. Stoker},
     64  TITLE =        {Long waves on a beach},
     65  JOURNAL =      {Journal of Fluid Mechanics},
     66  YEAR =         {1967},
     67  volume =       {27},
     68  number =       {4},
     69  pages =        {}
     70}
     71
     72
    5173
    5274
    … …  
    339361JOURNAL = {Coastal Engineering},
    340362}
     363
     364@ARTICLE{Zhang07,
     365  author =       {Y. J. Zhang and A. M. Baptista},
     366  title =        {An efficient and robust tsunami model on unstructured grids.
     367                  Part I: inundation benchmarks},
     368  journal =      {Pure and Applied Geophysics: Topical issue on Tsunamis},
     369  year =         {2007},
     370  month =        {},
     371  date =         {},
     372  volume =       {},
     373  pages =        {},
     374  note = {In press - wait for Phil Cummins to notify acceptance}
     375}
     376
     377

  • anuga_work/publications/anuga_2007/anuga_validation.tex

    r5341 r5353  
    5050\author[GA]{M.~J.~Sexton}
    5151\ead{Jane.Sexton@ga.gov.au}
     52\author[GA]{L.~Fountain}
     53\author[GA]{K.~VanPutten}
    5254\author[ANU]{S.~G.~Roberts}
    5355\ead{Stephen.Roberts@anu.edu.au}
    … …  
    158160present an exhaustive validation of the numerical model.
    159161Further to these tests, we will
    160 incorporate a test to verify friction values. The six tests are:
    161 
    162 (1) verification against the 1D analytical solution of Carrier and Greenspan;
    163 (2) testing against 1D (flume) data sets to verify wave height and velocity
    164 (3) determining friction values from 1D flume data sets;
    165 (4) validation against a genuinely 2D analytical solution of the model equations;
    166 (5) testing against the 2D Okushiri benchmark problem; and
    167 (6) testing against the 2D data sets modelling wave run-up around a circular island by Briggs et al.
    168 Throughout the paper, qualitative comparisons will be drawn against other models.
     162incorporate a test to verify friction values. The tests reported in
     163this paper are:
     164\begin{itemize}
     165  \item Verification against the 1D analytical solution of Carrier and
     166  Greenspan (Section \ref{sec:carrier})
     167  \item Testing against 1D (flume) data sets to verify wave height and velocity
     168  (Section \ref{sec:stage and velocity})
     169  \item Determining friction values from 1D flume data sets
     170  (Section \ref{sec:friction})
     171  \item Validation against a genuinely 2D analytical solution of the
     172  model equations
     173  (Section \ref{sec:XXX})
     174  \item Testing against the 2D Okushiri benchmark problem
     175  (Section \ref{sec:okushiri})   
     176  \item Testing against the 2D data sets modelling wave run-up around a circular island by Briggs et al.
     177  (Section \ref{sec:circular island})
     178\end{itemize}   
     179
     180
     181Throughout the paper, qualitative comparisons will be drawn against
     182other models.  Moreover, all source code necessary to reproduce the
     183results reported in this paper is available as part of the \ANUGA{}
     184distribution in the form of a test suite. It is thus possible for
     185anyone to readily verify that the implementation meets the
     186requirements set out by these benchmarks.
     187 
    169188
    170189%Hubbard and Dodd's model, OTT-2D, has some similarities to \ANUGA{}, and
    … …  
    178197conclusions outlined in section~\ref{sec:conclusions}.
    179198
    180 {\bf question - if the Okushiri result has already been presented in the
    181 MODSIM paper, how should it be presented in this paper? - simply refer to it I think}
    182 
    183 \section{Model}
     199
     200\section{Mathematical model, numerical scheme and implementation}
    184201\label{sec:model}
    185202
    186 The shallow water wave equations are a system of differential
    187 conservation equations which describe the flow of a thin layer of
    188 fluid over terrain. The form of the equations are:
     203The ANUGA model is based on the shallow water wave equations which are
     204widely regarded as suitable for modelling 2D flows subject to the
     205assumptions that horizontal scales (e.g. wave lengths) greatly exceed
     206the depth, vertical velocities are negligible and the fluid is treated
     207as inviscid and incompressible. See e.g. the classical texts
     208\citet{Stoker57} and \citet{Peregrine67} for the background or
     209\citet{Roberts1999} for more details on the mathematical model
     210used by \ANUGA{}.
     211
     212The conservation form of the shallow water wave
     213equations used in \ANUGA{} are:
    189214\[
    190215\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
    … …  
    196221$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
    197222entering the system are bed elevation $z$ and stage (absolute water
    198 level) $w$, where the relation $w = z + h$ holds true at all times.
     223level above a reference datum such as Mean Sea Level) $w$,
     224where the relation $w = z + h$ holds true at all times.
    199225The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
    200226by
    … …  
    227253in which $\eta$ is the Manning resistance coefficient.
    228254
    229 As demonstrated in our papers, \cite{modsim2005,Roberts1999} these
    230 equations provide an excellent model of flows associated with
    231 inundation such as dam breaks and tsunamis. Question - how do we
    232 know it is excellent?
     255%%As demonstrated in our papers, \cite{modsim2005,Roberts1999} these
     256%%equations provide an excellent model of flows associated with
     257%%inundation such as dam breaks and tsunamis. Question - how do we
     258%%know it is excellent?
    233259
    234260\ANUGA{} uses a finite-volume method as
    235 described in \citet{Roberts1999} where the study area is represented by an
    236 unstructured triangular mesh. The flexibility afforded by this approach
    237 is the ability for the user to refine the mesh in areas of interest.
    238 \ANUGA{} is mostly written in the object-oriented programming
    239 language Python with computationally intensive parts implemented
    240 as highly optimised shared objects written in C. The API is a
    241 Python script where the user sets up the scenario. This script
    242 defines the study area, mesh refinement as well as initial and boundary conditions.
    243 The user is free to update quantity values or boundary conditions through
    244 the simulation. Reference to user manual
    245 
    246 Could include here a brief overview of the numerical
    247 technique and reference the CTAC 2006 paper and has Steve written something
    248 that could be included?
    249 
    250 \section{Finite Volume Method}
    251 \label{sec:fvm}
    252 
    253 {\bf Jane: I don't think this section is needed here, but the
    254 content is referred to at the end of section 1}
    255 
    256 We use a finite-volume method for solving the shallow water wave
    257 equations \cite{Roberts1999}. The study area is represented by a mesh of
    258 triangular cells as in Figure~\ref{fig:mesh} in which the conserved
    259 quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
    260 in each volume are to be determined. The size of the triangles may
    261 be varied within the mesh to allow greater resolution in regions of
    262 particular interest.
    263 
    264 \begin{figure}
    265 \begin{center}
    266 %\includegraphics[width=5.0cm,keepaspectratio=true]{step-five}
    267 \caption{Triangular mesh used in our finite volume method. Conserved
    268 quantities $h$, $uh$ and $vh$ are associated with the centroid of
    269 each triangular cell.} \label{fig:mesh}
    270 \end{center}
    271 \end{figure}
    272 
    273 The equations constituting the finite-volume method are obtained by
    274 integrating the differential conservation equations over each
    275 triangular cell of the mesh. Introducing some notation we use $i$ to
    276 refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
    277 set of indices referring to the cells neighbouring the $i$th cell.
    278 Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
    279 the length of the edge between the $i$th and $j$th cells.
    280 
    281 By applying the divergence theorem we obtain for each volume an
    282 equation which describes the rate of change of the average of the
    283 conserved quantities within each cell, in terms of the fluxes across
    284 the edges of the cells and the effect of the source terms. In
    285 particular, rate equations associated with each cell have the form
    286 $$
    287  \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
    288 $$
    289 where
     261described in \citet{Roberts2007} where the study area is represented by an
     262unstructured triangular mesh in which the vector of conserved quantities
     263$\UU$ is maintained and updated over time. The flexibility afforded by
     264allowing unstructed meshes rather than fixed resolution grids
     265is the ability for the user to refine the mesh in areas of interest
     266while leaving other areas coarse and thereby conserving computational
     267resources.
     268
     269
     270The approach used in \ANUGA{} are distinguished from many
     271other implementations (e.g. \citet{Hubbard02} or \citet{Zhang07}) by the
     272following features:
    290273\begin{itemize}
    291 \item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
    292 \item $\SSS_i$ is the source term associated with the $i$th cell,
    293 and
    294 \item $\HH_{ij}$ is the outward normal flux of
    295 material across the \textit{ij}th edge.
    296 \end{itemize}
    297 
    298 
    299 %\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
    300 %cells
    301 %\item $m_{ij}$ is the midpoint of
    302 %the \textit{ij}th edge,
    303 %\item
    304 %$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
    305 %normal along the \textit{ij}th edge, and The
    306 
    307 The flux $\HH_{ij}$ is evaluated using a numerical flux function
    308 $\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
    309 water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
    310 $$
    311 H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
    312 $$
    313 
    314 Then
    315 $$
    316 \HH_{ij}  = \HH(\UU_i(m_{ij}),
    317 \UU_j(m_{ij}); \mathbf{n}_{ij})
    318 $$
    319 where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
    320 $\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
    321 \textit{ij}th edge. The function $\UU_i(x)$ for $x \in
    322 T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
    323 neighbouring  cells.
    324 
    325 We use a second order reconstruction to produce a piece-wise linear
    326 function construction of the conserved quantities for  all $x \in
    327 T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This
    328 function is allowed to be discontinuous across the edges of the
    329 cells, but the slope of this function is limited to avoid
    330 artificially introduced oscillations.
    331 
    332 Godunov's method (see \cite{Toro-92}) involves calculating the
    333 numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
    334 solving the corresponding one dimensional Riemann problem normal to
    335 the edge. We use the central-upwind scheme of \cite{KurNP2001} to
    336 calculate an approximation of the flux across each edge.
    337 
     274    \item The fluxes across each edge are computed using the semi-discrete
     275    central-upwind scheme for approximating the Riemann problem
     276    proposed by \citet{KurNP2001}. This scheme deals with different
     277    flow regimes such as shocks, rarefactions and sub to super
     278    critical flow transitions using one general approach. We have
     279    found this scheme to be pleasingly simple, robust and efficient.
     280    \item \ANUGA{} does not employ a shoreline detection algorithm as the
     281    central-upwind scheme is capable of resolving fluxes arising between
     282    wet and dry cells. \ANUGA{} does optionally bypass unnecessary
     283    computations for dry-dry cell boundaries for purely performance reasons.
     284    \item \ANUGA{} employs a second order spatial reconstruction of triangles
     285    to produce a piece-wise linear function construction of the conserved
     286    quantities. This function is allowed to be discontinuous across the
     287    edges of the cells, but the slope of this function is limited to avoid
     288    artificially introduced oscillations. This approach provides good
     289    approximation of steep gradients in the solution. However,
     290    where the depths are very small compared to the bed-slope a linear
     291    combination between second order and first order reconstructions is
     292    employed to guarantee numerical stability that may arise form very
     293    small depths.
     294\end{itemize}     
     295   
    338296\begin{figure}
    339297\begin{center}
    … …  
    347305
    348306In the computations presented in this paper we use an explicit Euler
    349 time stepping method with variable timestepping adapted to the
    350 observed CFL condition.
    351 
    352 
    353 \section{Software Implementation}
    354 \label{sec:software}
    355 
    356 {\bf Jane: I don't think this section is needed here, but the
    357 content is referred to at the end of section 1}
    358  
     307time stepping method with variable timestepping subject to the
     308CFL condition:
     309\[
     310  \delta t = \min_k \frac{r_k}{v_k} 
     311\]
     312where $r_k$ refers to the radius of the inscribed circle of triangle
     313$k$, $v_k$ refers to the maximal velocity calculated from fluxes
     314passing in or out of triangle $k$ and $\delta t$ is the resulting
     315'safe' timestep to be used for the next iteration.
     316
     317
     318\ANUGA{} utilises a general velocity limiter described in the
     319manual which guarantees a gradual compression of computed velocities
     320in the presence of very shallow depths:
     321\begin{equation}
     322  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h},
     323\end{equation}
     324where $h_0$ is a regularisation parameter that controls the minimal
     325magnitude of the denominator. The default value is $h_0 = 10^{-6}$.
     326
     327
    359328\ANUGA{} is mostly written in the object-oriented programming
    360329language Python with computationally intensive parts implemented
    … …  
    366335language syntax. In addition, Python's automatic memory management,
    367336dynamic typing, object model and vast number of libraries means that
    368 software can be produced quickly and can be readily adapted to
    369 changing requirements throughout its lifetime.
    370 
    371 The fundamental object in \ANUGA{} is the \code{Domain} which
    372 inherits functionality from a hierarchy of increasingly specialised
    373 classes starting with a basic structural Mesh to classes implementing
    374 the finite-volume scheme described in section \ref{sec:fvm}. Other classes
    375 are \code{Quantity} which represents values of one variable across the mesh
    376 along with their associated operations, \code{Geospatial_data} which
    377 represents georeferenced elevation data and a collection of \code{Boundary}
    378 classes which allows for a 'pluggable' way of driving the model.
    379 The conserved quantities updated automatically by the numerical scheme
    380 are stage (water level) $w$, $x$-momentum $uh$ and $y$-momentum
    381 $vh$. The quanitites elevation $z$ and friction $\eta$ are
    382 quantities that are not updated automatically but can be changed explicitly
    383 during run-time if the user wishes to do so.
    384 
    385 To set up a scenario the user specifies the study area along with any internal
    386 regions where increased mesh resolution is required. External edges may
    387 be labelled using symbolic tags which are subsequently used to bind
    388 boundary condition objects to tagged segments of the mesh boundary.
    389 The mesh is then generated using \ANUGA{}'s built-in mesh generator and
    390 converted into the \code{Domain} object which provides all methods used to
    391 setup and run the flow simulation. Figure \ref{fig:anuga mesh} shows an example of a mesh generated by \ANUGA{}.
    392 \begin{figure}
    393 \begin{center}
    394 \includegraphics[width=4in,keepaspectratio=true]{triangular-mesh}
    395 \caption{Triangular mesh used in our finite volume method. Conserved
    396 quantities $h$, $uh$ and $vh$ are associated with each triangular cell.}
    397 \label{fig:anuga mesh}
    398 \end{center}
    399 \end{figure}
    400 
    401 
    402 
    403 Next step is to setup initial conditions for each \code{Quantity} object. For
    404 the elevation $z$ this is typically obtained from bathymetric and
    405 topographic data sets. Setting initial values for quantities is done
    406 through the method \code{domain.set_quantity(name, X, location,
    407 region)} where name is the name of the quantity (e.g.\ 'stage',
    408 'xmomentum', 'ymomentum', 'elevation' or 'friction').  The variable X
    409 represents the source data for populating the quantity and may take
    410 one of the following forms:
    411 
    412 \begin{itemize}
    413 \item A constant value as in \code{domain.set_quantity('stage', 1)} which
    414 will set the initial water level to 1 m everywhere.
    415 \item Another quantity or a linear combination of quantities.  If \code{q1}
    416 and \code{q2} are two arbitrary quantities defined within the same domain,
    417 the expression \code{domain.set_quantity('stage', q1*(3*q2 + 5))} will set the stage
    418 quantity accordingly.  One common application of this would be to
    419 assign the stage as a constant depth above the bed elevation.
    420 \item An arbitrary function (or a callable object), \code{f(x, y)}, where \code{x}
    421 and \code{y} are assumed to be vectors. The quantity will be assigned values by
    422 evaluating \code{f} at each location within the mesh.
    423 \item An arbitrary set of points and associated values (wrapped into a
    424 Geospatial_data object). The points need not coincide with triangle
    425 vertices or centroids and a penalised least squares technique is
    426 employed to populate the quantity in a smooth and stable way.
    427 Since the least squares technique can be time consuming for large
    428 problems, \code{set_quantity} employs a caching technique which automatically
    429 decides whether to perform the computations or retrieve them from a
    430 cache.  This will typically speed up the build by several orders of
    431 magnitude after each computation has been performed once.
    432 \item A filename containing points and attributes.
    433 \item A Numerical Python array (or a list of numbers) ordered
    434 according to the internal data structure.
    435 \end{itemize}
    436 The parameter \code{location} determines whether the values should be
    437 assigned to triangle edge, midpoints or vertices and \code{region} allows the
    438 operation to be restricted to a region specified by a symbolic tag or
    439 a set of indices.
    440 
    441 Boundary conditions are bound to symbolic tags through the method
    442 \code{domain.set_boundary} which takes as input a lookup table (implemented
    443 as a Python dictionary) of the form \code{\{tag:~boundary_object\}}.
    444 The boundary objects are all assumed to be callable functions of vectors x
    445 and y.  Several predefined standard boundary objects are available and
    446 it is relatively straightforward to define problem-specific custom
    447 boundaries if needed.  The predefined boundary conditions include
    448 Dirichlet, Reflective, Transmissive, Temporal, and Spatio-Temporal
    449 boundaries.
    450 
    451 Forcing terms can be written according to a fixed protocol and added
    452 to the model using the idiom \code{domain.forcing_terms.append(F)} where \code{F} is
    453 assumed to be a user-defined callable object.
    454 
    455 When the simulation is running, the length of each time step is
    456 determined from the maximal speeds encountered and the sizes of
    457 triangles in order not to violate the CFL condition which specifies
    458 that no information should skip any triangles in one time step.  With
    459 large speeds and small triangles, time steps can become very small.
    460 In order to access the state of the simulation at regular time
    461 intervals, \ANUGA{} uses the method evolve:
    462 \begin{verbatim}
    463 For t in domain.evolve(yieldstep, duration):
    464     <model interrogation and modification>
    465 \end{verbatim}
    466 The parameter \code{duration} specifies the time period over which
    467 evolve operates, and control is passed to the body of the for-loop at
    468 each fixed time step called \code{yieldstep}.  The internal workings
    469 of the numerical scheme and its variable time stepping are thus
    470 decoupled from the fixed time stepping of the evolve loop.  This means
    471 that the user of the API may access the model at fixed timesteps to
    472 e.g.\ store model outputs, interrogate quantities or change the model
    473 itself at runtime. The evolve method has been implemented using a
    474 Python generator hence the reference to 'yield' in the parameter name.
    475 
    476 %Figure \ref{fig:beach runup} shows a simulation of water flowing onto a
    477 %hypothetical beach with obstacles.
    478 %A number of complex patterns are captured in this example including a shock where water reflected off the wall far (at the right hand side) meets the main flow. Other physical features are the standing waves and interference patterns.
    479 %See the \ANUGA{} User Manual at \url{http://sourceforge.net/projects/anuga} for more details and examples.
    480 %\begin{figure}
    481 %\begin{center}
    482 %\includegraphics[width=4in,keepaspectratio=true]{runup}
    483 %\caption{A hypothetical runup scenario.}
    484 %\label{fig:beach runup}
    485 %\end{center}
    486 %\end{figure}
     337\ANUGA scripts can be produced quickly and can be adapted fairly easily to
     338changing requirements.
    487339
    488340
    … …  
    503355
    504356\subsection{Stage and Velocity Validation in a Flume}
     357\label{sec:stage and velocity}
    505358This section will describe tilting flume tank experiments that were
    506359conducted at the Gordon McKay Hydraulics Laboratory at the University of
    … …  
    561414
    562415\subsection{1D flume tank to verify friction}
    563 
     416\label{sec:friction}
    564417The same tilting flume tank was used to validate stage and velocity
    565418was used to validate the ANUGA friction model. A ground slope of 1:20,
    … …  
    603456
    604457\subsection{Okushiri Wavetank Validation}
    605 
     458\label{sec:okushiri}
    606459As part of the Third International Workshop on Long-wave Runup
    607460Models in 2004 (\url{http://www.cee.cornell.edu/longwave}), four
    … …  
    618471This dataset has been used by to validate tsunami models by
    619472a number of tsunami scientists. Examples include Titov ... lit review
    620 here on who has used this example for verification
     473here on who has used this example for verification (Leharne?)
    621474
    622475\begin{figure}[htbp]
    … …  
    668521
    669522\subsection{Runup of solitary wave on circular island wavetank validation}
    670 
     523\label{sec:circular island}
    671524This section will describe the ANUGA results for the experiments
    672525conducted by Briggs et al (1995). Here, a 30x25m basin with a conical
    … …  
    716569
    717570\section{Conclusions}
    718 \label{sec:6}
     571\label{sec:conclusions}
    719572\ANUGA{} is a flexible and robust modelling system
    720573that simulates hydrodynamics by solving the shallow water wave
Note: See TracChangeset for help on using the changeset viewer.