Abstract
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 volume method. Numerical simulations show that NEP is successful to be a refinement/coarsening indicator in the adaptive mesh finite volume method, as the method refines the mesh or grids around nonsmooth regions and coarsens them around smooth regions.
1. 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–4]. In the work of Mungkasi and Roberts [5, 6], NEP was used for solving the one-dimensional shallow water equations. We note that two-dimensional 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.
2. 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 two-dimensional triangulation is represented by the matrices node and elem. The th row of the matrix node represents the coordinates of the th node. The th row of the matrix elem stores the three nodes (global indices of vertices) of the th element. The three nodes are stored counterclockwise.
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 shape-regular. A triangulation is called shape-regular if there exists a constant such thatfor all triangles in , where is the diameter of and is the measure of .
2.1. 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. The newest-vertex of the triangle is set as elem and the refinement edge is elem. That is, if we label vertices of the parent as , , and and the new node relating to a refinement of is , then the children of have vertices , , and and , , and . 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 cutEdgecutEdge. 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 . Triangle is bisected with respect to the longest edge of , as shown in Figure 1(b), leading to a nonconforming triangulation because a hanging node occurs. The neighbour of having the common edge of that is bisected is triangle . Then triangle is bisected with respect to the longest edge of , 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 , which leads to a fine conforming triangulation.
(a) Initial triangles
(b) After a triangle is bisected
(c) Hanging node still occurs
(d) Conforming triangles
2.2. 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 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.
Finding good nodes can be described as follows. Let and be conforming triangulations, where can be found by refinement of some triangles of . A node of is a good node if and only if(1)it is not a vertex of ,(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 and . Coarsening these two triangles results in a nonconforming triangulation, because a hanging node occurs, as shown in Figure 2(a). Then the two triangles and sharing the same hanging node are coarsened, which leads to a coarse conforming triangulation, as shown in Figure 2(b).
(a) After two triangles are coarsened
(b) Conforming triangles
2.3. 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 . 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).
3. Numerical Entropy Production
Recall the two-dimensional shallow water (wave) equationswhere is the vector of conserved quantities consisting of water depth , -momentum , and -momentum . Here, and are velocities in the - and -direction; and are flux functions in the - and -direction given bythe source term iswhere is the bed topography. The stage (absolute water level) is given by . 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]where the entropy isand entropy fluxes are
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–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 evolutionThe 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, 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 asIn (11), is found from the numerical entropy scheme (8) with as the input; is found from the finite volume scheme with 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–21].
4. 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 .
4.1. Planar Dam-Break Problem
We consider a spatial domain defined by . The domain is discretised into triangular cells initially. The initial condition isThe bottom topography is for all and . We take fixed time step 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 . 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 () the number of cells is , which is a smaller number than 2048, the initial number of cells.
4.2. Radial Dam-Break Problem
We consider a spatial domain defined by . The initial condition isThe bottom topography is for all and .
The domain is discretised into triangular cells initially. We take fixed time step 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 . Cells having absolute NEP larger than 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 (), the number of cells is .
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.
5. 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.
Competing Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was financially supported by the Australian National University (ANU), Canberra, Australia, and Sanata Dharma University, Yogyakarta, Indonesia. The author thanks Professor Stephen Gwyn Roberts at the Mathematical Sciences Institute of the ANU for the advice and discussions.