Anuga Software for Numerical Simulations of Shallow Water Flows

Shallow water flows are governed by the shallow water wave equations, also known as the Saint-Venant system. This paper presents a finite volume method used to solve the two-dimensional shallow water wave equations and how the finite volume method is implemented in ANUGA software. This finite volume method is the numerical method underlying the software. ANUGA is open source software developed by Australian National University (ANU) and Geoscience Australia (GA). This software uses the finite volume method with triangular domain discretisation for the computation. Four test cases are considered in order to evaluate the performance of the software. Overall, ANUGA is a robust software to simulate two-dimensional shallow water flows. Arus air dangkal diatur dalam persamaan gelombang air dangkal, dikenal sebagai sistem Saint-Venant. Penelitian ini menyajikan metode finite volumeyang digunakan untuk menyelesaikan persamaan gelombang air dangkal dua dimensi dan bagaimana metode finite volumediimplementasikan dalam perangkat lunak ANUGA. Metode finite volumeadalah metode numerik yang mendasari perangkat lunakANUGA. ANUGA sendiri adalah perangkat lunak open source yang dikembangkan oleh Australian National University(ANU) dan Geoscience Australia (GA). Perangkat lunak ini menggunakan metode finite volumedengan diskritisasi domain segitiga dalam proseskomputasi. Empat uji kasus digunakan untuk mengevaluasi kinerja perangkat lunak. Secara keseluruhan, ANUGA adalah perangkat lunak yang robust untuk mensimulasikan dua dimensi aliran arus air dangkal.


Introduction
Since the nineteenth century the shallow water wave equations have emerged as important mathematical models to describe water flows. Some of its developments and applications are given by Roberts et al. [1][2][3][4][5], and Stoker [6].
In 1892, Ritter [7] studied dam break modelling analytically. This old work has had a big impact as the dam break problem is now becoming a standard benchmark to test the performance of numerical methods. Even though this problem is originally modelled for onedimensional shallow water wave equations, it can This paper is the extended version from paper titled "A finite volume method for shallow water flows on triangular computational grids" that has been published in Proceeding of ICACSIS 2012. also be treated as a two-dimensional problem, which we present in this paper as a planar dam break problem. In addition, it can also be extended to a circular two-dimensional problem.
Oscillation of water on a paraboloid canal developed by Thacker [8] is another benchmark problem, which is interesting, and usually considered by shallow water researchers to test their numerical methods. This problem involves a wetting and drying process, a process that is wellknown to be very difficult to resolve [9][10][11][12][13]. Yoonand Cho [14] used this type of problem to test their numerical method.
To complete testing a numerical method, a steady flow problem must also be performed in addition to unsteady flow test problems. In this paper, we do perform a steady flow involving a shock over a parabolic obstruction. This test case is usually used in testing a one-dimensional shallow water solver.
One open-source free software developed to simulate two-dimensional shallow water flows is ANUGA, named after Australian National University (ANU) and Geoscience Australia (GA). This software uses a finite volume method as the underlying mathematical background. We will present this finite volume method and test the performance of ANUGA using the planar dam break problem, circular dam break problem, oscillation on a paraboloid channel, and a steady flow over an obstruction. ANUGA uses triangular domain discretisation for the computation.
The remainder of this paper is organised as follows. In Section 2, we present the shallow water wave equations governing water flows, the finite volume method and the ANUGA software. Section 3 contains four numerical simulations involving unsteady and steady flow problems. Finally, Section 4 concludes the paper with some remarks.

Methodology
Shallow water wave equations,the onedimensional shallow water wave equations can be written as The notations used here are described as follows: M L >DQD? Í is the vector of conserved quantities consisting of water depth D and momentum QD, B:M; L >QDQ 6 D E C å D 6 t? Í is the flux function, O L >r F C å D:V ë E ð;? Í is the source vector with ð is some defined friction terms, T represents the distance variable along the flow, P represents the time variable, V:T; is the fixed water bed, D:Tá P; is the height of the water at point Tand at time P, Q:Tá P; denotes the velocity of the water flow at point T and at time P, and C å is a constant denoting the acceleration due to gravity. In addition, a new variable called the stage can be defined as S:Tá P; L V:T; E D:Tá P; which is the absolute water level.
More general equations governing shallow water flows are the two-dimensional shallow water wave equations [15±20] and The source term including gravity and friction is where V:Tá U; is the bed topography, and 5 Ù L §5 Ùë 6 E 5 Ùì 6 is the bed friction modelled using Manning's resistance law and in which ß is the Manning resistance coefficient. It should be noted that the stage (absolute water level) S is given by S V E D. In this paper, we focus on the two-dimensional shallow water wave equations.
Integrating equation 2 over an arbitrary closed and connected spatial domain 3 having boundary and applying the Gauss divergence theorem to the flux terms, we get the integral form ò òP ± M @: The rotational invariance property of the shallow water wave equations implies that (ä J L 6 ?5 B:6M; (9) where 6 is the transformation matrix Therefore, equation 8 can be rewritten as ! !ç ì M @: Finite volume method,a finite volume method is based on subdividing the spatial domain into cells and keeping track of an approximation to the integral of the quantity over each of these cells. In this paper, we limit our presentation for triangular computational grids (cells).
After the spatial domain is discretized, we have the equation constituting the finite volume method over each triangular cell of the grids [16] whereM Ü is the vector of conserved quantities averaged over the Eth cell, O Ü is the source term associated with the Eth cell, * ÜÝ is the outward normal flux of material across the EFth edge, and H ÜÝ is the length of the EFth edge. Here, the EFth edge is the interface between the Eth and Fth cells. The flux * ÜÝ is evaluated using a numerical flux function *:ä á ä â ä ; such that for all conservation vectors M and normal vectors J Furthermore, where I ÜÝ is the midpoint of the EFth edge and J ÜÝ is the outward normal vector, with respect to the Here, we recall the algorithm to solve the two-dimensional shallow water wave equations given by Guinot [15]. In the semi-discrete framework, four steps are considered as follows: For each interface :Eá F;, transform the quantity M Ü and M Ý in the global coordinate system :Tá U; into the quantity M Ü Ü and M Ü Ý in the local coordinate system system :T Üá U Ü;. The water depth D is unchanged as it is a scalar variable, while the velocities Q and R are transformed into Q Ü and R Ü. Therefore, the new quantities in the local coordinate system are j.
Compute the flux B at the interface :Eá F;. In the local coordinate system, the problem is merely a one-dimensional Riemann problem, as the flux vector is parallel with the normal vector of the interface. Therefore, the equations to be solved are with initial condition given by M Ü Ü on one side of the interface :Eá F; and M Ü Ý on the other side of the interface :Eá F;. It should be stressed that at this step, we do not need to integrate equation19 for the quantity M Ü, but all we need is the value of the flux B . Hence, applying either the exact or approximate Riemann solvers to compute the flux B at the interface :Eá F; is enough. Note that here the source term does not need to be transformed in the local coordinate system, since it involves scalar variables only.
Transform the flux B back to the global coordinate system :Tá U;,so the flux at the midpoint of the interface :Eá F;in the global coordinate system is Here, B is the flux computed with a Riemann solver for the one-dimensional problem stated in Step 2. Finally, solve equation12 where è:E; L <rását=, as triangular grids are considered, for M Ü . We multiply * ÜÝ with H ÜÝ in order to get the flux over the interface :Eá F;. This is because the flux at the midpoint of the interface :Eá F; is * ÜÝ , and the lengthof the interface H ÜÝ .
In this part, we briefly describe an opensource free software developed using mathematical models described in the previous sections called ANUGA.Uppercased word 'ANUGA' refers to the name of the software, whereas lowercased word 'anuga' refers to the name of ANUGA module in programming. See ANUGA User Manual for the details of anuga [14].
ANUGA is an open source (free software) developed by Australian National University (ANU) and Geoscience Australia (GA) to be used for simulating shallow water flows, such as floods and tsunamis. The mathematical background underlying the software is the two-dimensional shallow water wave equations and finite volume method presented in Sections 2 and 3 above. The interface of this software is in Python language, but the mostly expensive computation parts, such as flux computation, are written in C language. A combination of these two languages provides two different advantages: Python has the flexibility in terms of software engineering, while C gives a very fast computation.
As has been described in the ANUGA User Manual [16], a simple simulation using ANUGA generally has five steps in the interface code. These five steps are importing necessary modules, setting-up the computational domain, setting-up initial conditions, setting-up boundary conditions, and evolving the system through time. The necessary modules to be imported are anuga ii and some standard libraries such as numpy, scipy, pylab, and time. The simplest creation of the triangular computational domain is in the rectangular-cross framework, which returns points, vertices, and boundary for the computation. The initial condition includes the definition of the topography, friction, and stage. The boundary conditions can be chosen in several ways depending on the needs: reflective boundary is for solid reflective wall; transmissive boundary is for continuing all values on the boundary; Dirichlet boundary is for constant boundary values; and time boundary is for time dependent boundary.A thorough description of this software is available at http://anuga.anu.edu.au.

Results and Analysis
This section presents some computational experiments using the finite volume method with ANUGA as the software used to perform the simulation. For simplicity of the experiments, frictionless topography is assumed. Four different test cases are considered: a planar dam break problem, circular dam break problem, water oscillation on a paraboloid channel, and steady flow over a parabolic obstruction.
The numerical settings are as follows. Even though ANUGA is capable to do the computations with unstructured triangular grids, we limit our presentation only to the structured triangular grids for brevity with rectangular-cross as the basis. The spatial reconstruction and temporal discretisation are both second order with edge limiter. A reflective boundary is imposed. The numerical flux is due to Kurganov et al. [21]. SI units are used in the measurements of the quantities.   Consider a flat topography without friction on a spatial domain given by <:Tá U; :Tá U; Ð >Fwráwr? H >Fsrásr?=. Suppose that we are given a planar dam break problem with wet/dry interface bewtween the dam, where we have sr meters water depth on the left of the dam, but dry area on the right. The dam is at <:rá U; U Ð >Fsrásr?=. The simulation is then done to predict the flow of water after the dam is removed. Note that this dam break problem is actually a one-dimensional dam-break problem, but here we treat this as a two-dimensional problem.    The domain setting is as follows. The spatial domain is discretised into srr by tr rectangularcrosses, where each rectangular cross has v uniform triangles. This means that we have z " sr 7 structured triangles as the discretisation of the spatial domain.
The water surface and topography are shown in figure1, 2, 3, and 4. In particular, figure 1 shows the initial condition, that is, the water profile at timeP L rär s. After the dam is broken, the profiles at time P L räwá särá säw s are illustrated in figure 2, figure 3, and figure 4 respectively. At säw seconds after dam break, the water flows to the right for a distance almost ur m. This agrees with the analytical solution (ur m) for the one-dimensional dam break problem. Now suppose that we have a circular dam break problem with wet/wet interface between the dam. Consider a circular dam, where the point of origin is the centre, having tr m as its radius. The water depth inside the circular dam is sr m, while outside is s m. The spatial domain of interest is<:Tá U; :Tá U; Ð >Fwráwr? H >Fwráwr?=. Similar to the planar dam break problem, the simulation of this test case is then done to predict the flow of water, after the dam is removed.
The domain setting is as follows. The spatial domain is discretised into srr bysrr rectangularcrosses, where each rectangular cross has v uniform triangles. This means that we havev " sr 8 structured triangles as the discretisation of the spatial domain.
The water surface and topography are shown in figure 5, 6, 7, and 8. In particular, figure 5 shows the initial condition, that is, the water profile at time P L rär s. After the dam is broken, the profiles at time P L räwá särá säw s are illustrated in figure 6, figure 7, and figure 8 respectively. At säw seconds after dam break, the water floods to the surrounding area approximately tw m outwards.
This test case is adapted from the work of Thacker [8] and Yoon and Chou [14]. Suppose that we have a spatial domain given by<:Tá U; :Tá U; Ð >Fvrrrávrrr? H >Fvrrrávrrr?=.
We define some parameters & 4 L srrrá . L twrrá 4 4 L trrrsuch that are the amplitude of oscillation at the origin, the angular frequency, and the period of oscillation. Furthermore, the topography is given by and the initial water level (the initial wet region) is The domain setting is as follows. The spatial domain is discretised into wr by wr rectangularcrosses, where each rectangular cross has v uniform triangles. This means that we have sr 8 structured triangles as the discretisation of the spatial domain.
The simulation is then done to perdict the flow of the water. According to Thacker [8], the flow will be a periodic oscillation. The simulation, using the numerical method and software discussed in this paper, indeed shows this periodic motion. The water profiles at time P L rá 6 tá 6 are shown in figure 9, figure 10, and figure 11 respectively. In addition, the periodic motion of the water surface corresponding to the origin is depicted in figure 12. From the numerical results, we see a minor damping of the oscillation. However, a perfect method or software should not produce this damping. Although it may be impossible to have a perfect method, we take this for future research for the ANUGA software to minimise or eliminate this kind of damping. We also see from the numerical results that a negligible amount of water is left over in the supposed-dry area.
Steady flow over a parabolic obstruction,this test case is for checking if ANUGA can produce a steady shock. Suppose that we have a spatial domain given by<:Tá U; :Tá U; Ð >rásw? H >rásw?=.
We consider the topography as horizontal bed except for a defined parabolic obstruction V:Tá U; L rät F rärw:T F sr; 6 (27) We use Dirichlet boundaries for the left and right boundaries. The left boundary is set up as water height rävt, Tmomentum 0.18, and U-momentum 0.0. The right boundary is set up as water height räuu, Tmomentum 0.18, and U-momentum 0.0. The initial condition is just a quiet lake at rest with the stage level as räuu.    The domain setting is as follows. The spatial domain is discretised into str by str rectangular-crosses, where each rectangular cross has v uniform triangles. This means that we have wyáxrrstructured triangles as the discretisation of the spatial domain.
The summary of the simualtion results is as follows. Figures 13 and 14 show the water surface at initial condition and at time P L twä After water is flowing for 25 seconds from left boundary to the right, a steady shock occurs on the parabolic obstruction. This agrees with the simulation using purely one-dimensional model using equation1 in our previous work [22]. (It is unfortunate that due to the small water height, the shock and the bottom topography are not clearly shownin figure 14. However, these can be clearly seen using ANUGA Viewer as parts of the picture can be magnified.).

Conclusion
We have presented some mathematical (both analytical and numerical) background underlying the ANUGA software. ANUGA was tested using four numerical simulations. We infer from the simulation results that ANUGA can solve nicely the wetting problem. For the drying problem, some very small amount of water, which is negligible, may be left over on the supposed-dry area. Regardless of what we have investigated, we need to note that different parameters or numerical settings generally lead to different results. We also note that ANUGA can simulate both steady and unsteady state problems.