Adaptive Finite Volume Method for the Shallow Water Equations on Triangular Grids

This paper presents a numerical entropy production (NEP) scheme for two-dimensional shallow water equations on unstructured triangular grids. We implement NEP as the error indicator for adaptive mesh refinement or coarsening in solving the shallow water equations using a finite volumemethod. Numerical simulations show thatNEP is successful to be a refinement/coarsening indicator in the adaptive mesh finite volumemethod, as the method refines the mesh or grids around nonsmooth regions and coarsens them around smooth regions.


Introduction
The numerical entropy production (NEP) has been implemented successfully as smoothness indicator (error indicator) in adaptive finite volume methods used to solve compressible fluid problems [1][2][3][4].In the work of Mungkasi and Roberts [5,6], NEP was used for solving the onedimensional shallow water equations.We note that twodimensional problems are more complex and more applicable in real situations but challenging to solve.
The goal of this paper is to implement the NEP into an adaptive mesh finite volume method used to solve the two-dimensional shallow water equations on unstructured triangular grids.Recall that NEP diverges under grid refinement on nonsmooth regions, while it has the same order of magnitude of the local error on smooth ones [4].The difference between the degrees of accuracy between smooth and nonsmooth regions makes the NEP able to detect the smoothness of the numerical solution.Therefore, we are confident here to implement the NEP as smoothness indicator into the adaptive mesh finite volume method.
Our adaptive mesh finite volume method uses iFEM as the underlying data structure, an integrated finite element methods package in MATLAB.We note that iFEM and its data structure are not of our work, but they were developed by Chen [7].Even though iFEM was originally developed for finite element methods, we can still use it to achieve our goal for finite volume methods.This is because what we need is actually the data structure of iFEM.Therefore, we can adapt the original iFEM into our finite volume problems.
This paper is organised as follows.Section 2 summarises the iFEM data structure for triangulations to assist mesh refinement and coarsening techniques in the finite volume method.In Section 3, we recall the two-dimensional shallow water equations and present the scheme for the numerical entropy production.Some simulation results are reported in Section 4. Finally, we draw concluding remarks in Section 5.

Data Structure for Triangulations
We implement two-dimensional mesh triangulations of iFEM.In this section, we summarise the iFEM data structure for triangulations to assist mesh adaptation (refinement or coarsening) techniques in our finite volume method.More detailed explanations can be found in a technical report written by Chen [7] and a paper by Chen and Zhang [8].Note that in our work the spatial dimension is two.
Let  be the number of vertices and  be the number of elements, where the elements are triangles.A twodimensional triangulation is represented by the matrices node(1 : , 1 : 2) and elem(1 : , 1 : 3).The th row of 2   Two conditions on triangulations are imposed.The first condition is that a triangulation must be conforming.That is, the intersection of any two simplexes in the triangulation is either empty or a common lower dimensional simplex.The second condition is that all triangulations must be shaperegular.A triangulation T is called shape-regular if there exists a constant  0 such that max ∈ diam () 2  || ≤  0

Advances in Mathematical Physics
for all triangles  in T, where diam() is the diameter of  and || is the measure of .

Refinement: The Bisection Method.
A newest-vertex bisection method is implemented for grid refinement, but in the initial triangulation the longest edge is assigned as the refinement edge of each triangle.The function bisect refines a given triangulation by bisecting marked elements and minimal neighbouring elements to get a conforming and shape-regular triangulation.This means that a triangle is bisected based on the newest-vertex appearing on that triangle.In other words, the newest-vertex appearing on that triangle is used as the reference for the triangular bisection.A triangulation is labeled as follows.Note that global indices of three vertices of a triangle  are given by elem(, {1, 2, 3}).The newest-vertex of the triangle  is set as elem(, 1) and the refinement edge is elem(, {2, 3}).That is, if we label vertices of the parent  as  1 ,  2 , and  3 and the new node relating to a refinement of  is  4 , then the children of  have vertices  4 ,  1 , and  2 and  4 ,  3 , and  1 .The first child called the LEFT (L) child occupies the location used by .The second child called the RIGHT (R) child is appended as a new element to elem matrix.To be specific, the first (L) child is stored prior to the second (R) child in elem matrix.We call this storing technique the (, ) positioning.
Bisections over a collection of elements could result in a nonconforming triangulation.To enforce that we get a conforming triangulation, we do the following.Let cutEdge() be the refinement edge of  and () be the refinement neighbour of .The refinement neighbour means the neighbour sharing the refinement edge of .We note that a triangle  has at most three neighbours and we may have that cutEdge() ̸ = cutEdge(()).The rule which enforces conformity is that if cutEdge() is bisected, then cutEdge(()) must also be bisected.In iFEM, this rule is implemented in a loop and the loop terminates in a finite number of steps producing a conforming triangulation.
An example of the refinement procedure is illustrated in Figure 1.Let us consider an initial triangulation, as shown in Figure 1(a).Suppose that we want to refine triangle  1 .Triangle  1 is bisected with respect to the longest edge of  1 , as shown in Figure 1(b), leading to a nonconforming triangulation because a hanging node occurs.The neighbour of  1 having the common edge of  1 that is bisected is triangle  2 .Then triangle  2 is bisected with respect to the longest edge of  2 , as shown in Figure 1(c), leading to another nonconforming triangulation because the hanging node still occurs.Finally the hanging node is used as a reference for bisecting triangle  2.2 , which leads to a fine conforming triangulation.

Coarsening: A Nodal-Wise Algorithm.
A nodal-wise coarsening algorithm is used for adaptive grids obtained from newest-vertex bisections, with a note that the algorithm enforces not to coarsen the grids in the initial triangulation.A coarsening process is thought as the inverse of a conforming bisection and restricted to a conforming star.A conforming bisection means a bisection which results in a conforming triangulation, and a conforming star means a star-shaped domain consisting of some conforming triangles.
A  Finding good nodes can be described as follows.Let  and  0 be conforming triangulations, where  can be found by refinement of some triangles of  0 .A node of  is a good node if and only if (1) it is not a vertex of  0 , (2) it is the newest-vertex of all elements in the star of the node, (3) the number of elements in the star is four if the node is an interior node or two otherwise (a boundary node).
Therefore, if some triangles of  are marked to be coarsened, then we check the newest vertices, which satisfies the three conditions, of the marked elements.
After good nodes are identified, coarsening of conforming stars can be done by considering the (, ) type of positioning.If a good node is an interior node, then four conforming triangles are coarsened into two conforming ones.If a good node is a boundary node, then two conforming triangles are coarsened into one triangle.
An example of the coarsening procedure is illustrated in Figure 2. Let us consider the triangulation in Figure 1(d).Suppose that we want to coarsen the triangles  1.1 and  1.2 .Coarsening these two triangles results in a nonconforming triangulation, because a hanging node occurs, as shown in Figure 2(a).Then the two triangles  2.2.1 and  2.2.2 sharing the same hanging node are coarsened, which leads to a coarse conforming triangulation, as shown in Figure 2(b).

Conserved Quantities in Finite Volume
Methods.The data structure for conserved quantities follows from the structure for triangulations.In addition, if a parent-cell is bisected, then we do a quantitative splitting of the amount of the conserved quantity of the parent-cell into two parts for the resulting child-cells.If two conforming child-cells are coarsened then we do a quantitative averaging of the amount of the conserved quantity of the child-cells for the resulting parent-cell.The splitting and averaging are done in a proportional weight with respect to the areas of the base child-cells.
Remark 1.In the refinement and coarsening processes of iFEM, the longest edge is assigned as the refinement edge of each triangle in the initial triangulation T 0 .In iFEM, the subroutine label is needed only for the initial triangulation.Therefore, the subroutine label is called outside of the function bisect.This procedure leads to a shape-regular triangulation and so prevents the formation of triangles with very small angles (see Chen [7] for more details on this data structure of iFEM).

Numerical Entropy Production
Recall the two-dimensional shallow water (wave) equations where q = [ℎ ℎ ℎV]  is the vector of conserved quantities consisting of water depth ℎ, -momentum ℎ, and momentum ℎV.Here,  and V are velocities in the and direction; f(q) and g(q) are flux functions in the and direction given by the source term is where (, ) is the bed topography.The stage (absolute water level) is given by  fl  + ℎ.The constant  is the acceleration due to gravity.The free variables are the spaces  and  as well as time .More explanation about the shallow water equations can be found in the literature, such as Akyildiz [9], Peng [10], and Yue et al. [11].The entropy inequality is [12]   +   +   ≤ 0, where the entropy is and entropy fluxes are  (q (, , ) ,  (, )) = ( 1 2 ℎ ( 2 + V 2 ) + ℎ 2 + ℎ) ,  (q (, , ) ,  (, )) We solve the two-dimensional shallow water equations using a finite volume method described by Mungkasi and Roberts [13], but in the present paper, we improve the method adaptively.The finite volume method itself has been used to solve one-dimensional problems successfully [14][15][16].
For a smooth solution, (5) becomes an equality and is called entropy equation.To compute the NEP we need to evolve numerically the entropy based on the entropy residual, given by the discretisation of the left hand side of ( 5).This evolution is done together with the evolution of the conserved quantities.We may discretise the spatial domain into arbitrary polygonal cells, but in this paper we focus our work on triangular grids.Then for an arbitrary cell , we have a semidiscrete scheme for the entropy evolution The notation  is the edge (interface) between the  and  cells.
Here Θ  is an approximation of the averaged entropy over the  cell, Ψ  is the outward normal entropy flux across the edge , and   is the length of the edge .In addition, N() is the set of three neighbouring cells of the  cell and   is the area of the  cell.The flux Ψ  is evaluated using a consistent numerical flux function Ψ(⋅, ⋅) such that for all values of Θ Ψ (Θ, Θ) =  (Θ) . Furthermore, where   is the midpoint of the  edge.Note that in (8) we have used the sign "≲."With this sign, we mean that it is not guaranteed that inequality "<" in ( 8) holds.Note also that it has been proven that the entropy production, which is the left hand side of (8), is negative definite only for some cases.In general, the entropy production may have positive overshoots.More details can be found in the work of Puppo and Semplice [4].
Recall that, in evolving the conserved quantities [13], we employ rotations from global coordinates into local coordinates in order to compute the conserved quantity fluxes.The rotations are needed because the conserved quantities are vector-valued.However, unlike the conserved quantities, the entropy is scalar-valued.Therefore, in evolving the entropy numerically we do not need to do any rotation but directly compute the edge-flux at each interface between cells.
The NEP is formulated as In (11), Θ   is found from the numerical entropy scheme (8) with (Q −1 ) as the input; Q   is found from the finite volume scheme with Q −1 as the input; and Δ is the time step used in those two schemes.Following Puppo and Semplice [4], the NEP formula is expressed as Remark 2. The left hand side of ( 8) is not zero.It defines the entropy residual (in semidiscrete form), from which we then derive    in (11) and (12).The left hand side of (8) in fact is not exactly zero, even on smooth flows, because we have discretised the entropy relation in space.More discussions about entropy and related schemes are given comprehensively by a number of researchers, such as Tadmor et al. [17][18][19][20][21].

Numerical Simulations
In this section, we present two simulations to show the performance of the NEP as smoothness indicator in the adaptive mesh finite volume methods.The two simulations are planar and radial dam-break problems.A first-order method is used.Quantities are measured in SI units.The numerical flux due to Kurganov, Noelle, and Petrova (KNP) [22] is used.The acceleration due to gravity is taken as  = 9.81.( The bottom topography is (, ) = 0 for all  and .We take fixed time step Δ = 0.002 with 100 iterations.The maximum level of refinement/coarsening is 2 at each iteration.That is, a triangular cell can be refined or coarsened up to two levels finer or two levels coarser at each iteration.However, we restrict our refinement/coarsening only up to 8 levels of triangles from the smallest to the largest.The NEP tolerance is NEP tol = 0.25 max|NEP|.That is, if the NEP is larger than the tolerance, then the corresponding cell is refined in order to get more accurate results.After the refinement procedure if the NEP is lower than the tolerance, then the corresponding cell is coarsened in order to get faster computation.Representatives of simulation results are shown in Figures 3 and 4. We see that fine grids occur around shock and rarefaction regions in order to get accurate solutions at those regions.Coarser grids occur in smooth regions having flat water surface in this case.This means that the adaptive finite volume method works well as we expect and is consistent with the behaviour for the one-dimensional problems as in our Ph.D. thesis [23].At this time ( = 0.2) the number of cells is 1544, which is a smaller number than 2048, the initial number of cells.

Radial Dam-Break Problem. We consider a spatial domain defined by (𝑥, 𝑦
The bottom topography is (, ) = 0 for all  and .The domain is discretised into 2048 triangular cells initially.We take fixed time step Δ = 0.002 with 25 iterations.The maximum level of refinement/coarsening is 2 at each iteration, so that a triangular cell can be refined or coarsened up to two levels finer or two levels coarser at each iteration.Again, we restrict our refinement/coarsening only up to 8 levels of triangles from the smallest to the largest.The NEP tolerance is NEP tol = 0.25 max|NEP|.Cells having absolute NEP larger than NEP tol are candidates for refinement, while otherwise they are candidates for coarsening.
Representatives of simulation results are shown in Figures 5 and 6.Similar to the previous simulation, we see that the grids are small around the shock and rarefaction in order to get accurate solutions and large at flat water surface in order to reduce the computational cost.At this time ( = 0.05), the number of cells is 1968.Remark 3. In our simulations above, we have specified intuitively the criterion for choosing the NEP threshold (tolerance) used to trigger grid refinement and coarsening.
Obviously the tolerance should be between zero and the maximum value of the absolute NEP.If the tolerance is too large, then the grids will be too coarse so the wave will be too diffusive and inaccurate.If the tolerance is too small, then the grids will be too fine; that is, the computation will be very expensive even though the solution is accurate.This trade-off could be studied further but is out of the scope of this paper.

Conclusions
Numerical entropy production has been implemented into an adaptive mesh finite volume method used to solve the two-dimensional shallow water equations.The implementation is successful, as the method does refinement around nonsmooth regions and coarsening around smooth regions.Future work will focus on the error analysis and investigate the rate of convergence of the adaptive method.

Disclosure
Some parts of this work were conducted and written as a portion of the author's Ph.D. thesis at the Australian National University.Other parts were conducted and written at Sanata Dharma University after the author received the Ph.D. degree.

Figure 1 :
Figure 1: Illustration of the refinement procedure.

Figure 2 :
Figure 2: Illustration of the coarsening procedure.

Figure 3 :Figure 4 :
Figure 3: The water surface of the planar dam-break problem at  = 0.2 using the adaptive method.The base plane is the -plane.

Figure 5 :Figure 6 :
Figure5: The water surface of the radial dam-break problem at  = 0.05 using the adaptive method.The base plane is the -plane.
node  is called a good-for-coarsening node or a good node in short if there exists a conforming bisection   such that  is the middle point of , where   bisects (i) two conforming triangles into four conforming triangles if  is an interior edge or (ii) one triangle into two conforming triangles if  is a boundary edge.