Changeset 5353
 Timestamp:
 May 21, 2008, 3:49:12 PM (16 years ago)
 Location:
 anuga_work/publications/anuga_2007
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/publications/anuga_2007/anugabibliography.bib
r5341 r5353 49 49 50 50 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 51 73 52 74 … … 339 361 JOURNAL = {Coastal Engineering}, 340 362 } 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 50 50 \author[GA]{M.~J.~Sexton} 51 51 \ead{Jane.Sexton@ga.gov.au} 52 \author[GA]{L.~Fountain} 53 \author[GA]{K.~VanPutten} 52 54 \author[ANU]{S.~G.~Roberts} 53 55 \ead{Stephen.Roberts@anu.edu.au} … … 158 160 present an exhaustive validation of the numerical model. 159 161 Further 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 runup around a circular island by Briggs et al. 168 Throughout the paper, qualitative comparisons will be drawn against other models. 162 incorporate a test to verify friction values. The tests reported in 163 this 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 runup around a circular island by Briggs et al. 177 (Section \ref{sec:circular island}) 178 \end{itemize} 179 180 181 Throughout the paper, qualitative comparisons will be drawn against 182 other models. Moreover, all source code necessary to reproduce the 183 results reported in this paper is available as part of the \ANUGA{} 184 distribution in the form of a test suite. It is thus possible for 185 anyone to readily verify that the implementation meets the 186 requirements set out by these benchmarks. 187 169 188 170 189 %Hubbard and Dodd's model, OTT2D, has some similarities to \ANUGA{}, and … … 178 197 conclusions outlined in section~\ref{sec:conclusions}. 179 198 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} 184 201 \label{sec:model} 185 202 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: 203 The ANUGA model is based on the shallow water wave equations which are 204 widely regarded as suitable for modelling 2D flows subject to the 205 assumptions that horizontal scales (e.g. wave lengths) greatly exceed 206 the depth, vertical velocities are negligible and the fluid is treated 207 as 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 210 used by \ANUGA{}. 211 212 The conservation form of the shallow water wave 213 equations used in \ANUGA{} are: 189 214 \[ 190 215 \frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial … … 196 221 $h$, $x$momentum $uh$ and $y$momentum $vh$. Other quantities 197 222 entering the system are bed elevation $z$ and stage (absolute water 198 level) $w$, where the relation $w = z + h$ holds true at all times. 223 level above a reference datum such as Mean Sea Level) $w$, 224 where the relation $w = z + h$ holds true at all times. 199 225 The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given 200 226 by … … 227 253 in which $\eta$ is the Manning resistance coefficient. 228 254 229 As demonstrated in our papers, \cite{modsim2005,Roberts1999} these230 equations provide an excellent model of flows associated with231 inundation such as dam breaks and tsunamis. Question  how do we232 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? 233 259 234 260 \ANUGA{} uses a finitevolume 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 objectoriented 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 finitevolume 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]{stepfive} 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 finitevolume 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 261 described in \citet{Roberts2007} where the study area is represented by an 262 unstructured triangular mesh in which the vector of conserved quantities 263 $\UU$ is maintained and updated over time. The flexibility afforded by 264 allowing unstructed meshes rather than fixed resolution grids 265 is the ability for the user to refine the mesh in areas of interest 266 while leaving other areas coarse and thereby conserving computational 267 resources. 268 269 270 The approach used in \ANUGA{} are distinguished from many 271 other implementations (e.g. \citet{Hubbard02} or \citet{Zhang07}) by the 272 following features: 290 273 \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 piecewise 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{Toro92}) 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 centralupwind 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 semidiscrete 275 centralupwind 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 centralupwind scheme is capable of resolving fluxes arising between 282 wet and dry cells. \ANUGA{} does optionally bypass unnecessary 283 computations for drydry cell boundaries for purely performance reasons. 284 \item \ANUGA{} employs a second order spatial reconstruction of triangles 285 to produce a piecewise 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 bedslope 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 338 296 \begin{figure} 339 297 \begin{center} … … 347 305 348 306 In 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 307 time stepping method with variable timestepping subject to the 308 CFL condition: 309 \[ 310 \delta t = \min_k \frac{r_k}{v_k} 311 \] 312 where $r_k$ refers to the radius of the inscribed circle of triangle 313 $k$, $v_k$ refers to the maximal velocity calculated from fluxes 314 passing 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 319 manual which guarantees a gradual compression of computed velocities 320 in 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} 324 where $h_0$ is a regularisation parameter that controls the minimal 325 magnitude of the denominator. The default value is $h_0 = 10^{6}$. 326 327 359 328 \ANUGA{} is mostly written in the objectoriented programming 360 329 language Python with computationally intensive parts implemented … … 366 335 language syntax. In addition, Python's automatic memory management, 367 336 dynamic 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 finitevolume 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 runtime 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 builtin 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]{triangularmesh} 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 problemspecific custom 447 boundaries if needed. The predefined boundary conditions include 448 Dirichlet, Reflective, Transmissive, Temporal, and SpatioTemporal 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 userdefined 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 forloop 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 338 changing requirements. 487 339 488 340 … … 503 355 504 356 \subsection{Stage and Velocity Validation in a Flume} 357 \label{sec:stage and velocity} 505 358 This section will describe tilting flume tank experiments that were 506 359 conducted at the Gordon McKay Hydraulics Laboratory at the University of … … 561 414 562 415 \subsection{1D flume tank to verify friction} 563 416 \label{sec:friction} 564 417 The same tilting flume tank was used to validate stage and velocity 565 418 was used to validate the ANUGA friction model. A ground slope of 1:20, … … 603 456 604 457 \subsection{Okushiri Wavetank Validation} 605 458 \label{sec:okushiri} 606 459 As part of the Third International Workshop on Longwave Runup 607 460 Models in 2004 (\url{http://www.cee.cornell.edu/longwave}), four … … 618 471 This dataset has been used by to validate tsunami models by 619 472 a number of tsunami scientists. Examples include Titov ... lit review 620 here on who has used this example for verification 473 here on who has used this example for verification (Leharne?) 621 474 622 475 \begin{figure}[htbp] … … 668 521 669 522 \subsection{Runup of solitary wave on circular island wavetank validation} 670 523 \label{sec:circular island} 671 524 This section will describe the ANUGA results for the experiments 672 525 conducted by Briggs et al (1995). Here, a 30x25m basin with a conical … … 716 569 717 570 \section{Conclusions} 718 \label{sec: 6}571 \label{sec:conclusions} 719 572 \ANUGA{} is a flexible and robust modelling system 720 573 that simulates hydrodynamics by solving the shallow water wave
Note: See TracChangeset
for help on using the changeset viewer.