Adaptive Central-Upwind Scheme on Triangular Grids for the Shallow Water Model with variable density

Thuong Nguyen¹¹1Department of Mathematics, The University of Utah, Salt Lake City, UT 84112, USA; [email protected]

Abstract

In this paper, we construct a robust adaptive central-upwind scheme on unstructured triangular grids for two-dimensional shallow water equations with variable density. The method is well-balanced, positivity-preserving, and oscillation free at the curve where two types of fluid merge. The proposed approach is an extension of the adaptive well-balanced, positivity-preserving scheme developed in Epshteyn and Nguyen (arXiv preprint arXiv:2011.06143, 2020). In particular, to preserve “lake-at-rest” steady states, we utilize the Riemann Solver with appropriately rotated coordinates to obtain the point values in neighborhood of the fluid interface. In addition, to improve the efficiency of an adaptive method in the multifluid flow, the curve of density discontinuity is reconstructed by using the level set method and volume fraction method. To demonstrate the accuracy, high-resolution, and efficiency of the new adaptive central-upwind scheme, several challenging tests for Shallow water models with variable density are performed.

keywords:

Shallow water equations with variable density, central-upwind scheme, well-balanced and positivity-preserving scheme, adaptive algorithm, interface tracking, Riemann solver, weak local residual error estimator, unstructured triangular grid

AMS subject classification: 76M12, 65M08, 35L65, 86-08, 86A05

1 Introduction

The main goal of this paper is to develop an adaptive well-balanced positivity-preserving central-upwind scheme on triangular grids for shallow water equations with variable density (SWEDs). The two-dimensional (2-D) system of SWEDs can be written as,


	$\displaystyle w_{t}+(hu)_{x}+(hv)_{y}=0,$		(1.1a)
	$\displaystyle(hu)_{t}+\Big{(}hu^{2}+\frac{g}{2\rho_{0}}h^{2}\rho\Big{)}_{x}+(huv)_{y}=-\frac{g}{\rho_{0}}h\rho B_{x},$		(1.1b)
	$\displaystyle(hv)_{t}+(huv)_{x}+\Big{(}hv^{2}+\frac{g}{2\rho_{0}}h^{2}\rho\Big{)}_{y}=-\frac{g}{\rho_{0}}h\rho B_{y},$		(1.1c)
	$\displaystyle(h\rho)_{t}+(hu\rho)_{x}+(hv\rho)_{y}=0,$		(1.1d)

where $t$ is the time, $x$ and $y$ are spatial coordinates ( $(x,y)\in\Omega$ ), $h(x,y,t)$ is the water height, $\rho(x,y,t)$ is the density, $u(x,y,t)$ and $v(x,y,t)$ are the $x$ - and $y$ -components of the flow velocity, $B(x,y)$ is the bottom topography, $g$ is the constant gravitational acceleration, and $\rho_{0}$ is the reference density. The system (1.1a)–(1.1d) was proposed in [10, 35, 36, 15] as a variation of the Saint-Venant equations to model multi-phase flows in estuaries or deep ocean currents. The derivation of the system is based on hydrostatic approximation which eliminates the variability in the z-direction. The design of robust and accurate numerical algorithms for computing the solutions of SWEDs system is an important and challenging problem that has been extensively studied in the recent years.

A number of numerical schemes for balance laws have been introduced in recent years, [28, 24, 20, 5, 27, 7, 6, 8, 23, 22, 2, 26, 3, 4, 39, 29]. Most of them utilize a Riemann problem solver for the upwind evolution of the calculated solution. However, as discussed in [9], the eigensystem of the system (1.1a)–(1.1d) may be incomplete due to the resonance phenomenon. Hence, it may be very difficult to design a reliable upwind scheme for the SWEDs. In our paper, we therefore use central-upwind schemes which are Riemann-problem-solver free methods, [21, 25, 29]. Central-upwind schemes have been referred to “black-box” solvers for general multidimensional systems of hyperbolic systems of conservation laws. In our prior work [13], we have derived a successful adaptive central-upwind method for Saint-Venant system on triangular grids.

Similar to the Saint-Venant system, a good method for SWEDs system should preserve the non-negativity of $h$ and $\rho$ , which is called the positivity-preserving property. In addition, the scheme must ensure a well-balanced property obtained when the numerical method preserve “lake-at-rest” steady-state solutions. Otherwise, the numerical method may lead to significant oscillations. Note that, the system (1.1a)–(1.1d) admits the following two “lake-at-rest” steady-state solutions, [9]:

w=\max\big{\{}C,B(x,y)\big{\}},\quad C=\rm{Const},\quad\rho=P\equiv Const,\quad u\equiv v\equiv 0,

(1.2)

and

B\equiv\rm{Const},\quad h^{2}\rho\equiv Const,\quad u\equiv v\equiv 0.

(1.3)

Preserving the solution (1.3) is a big challenge for numerically solving the system (1.1a)–(1.1d) since using the conventional central-upwind methods may not ensure the variable $h^{2}\rho$ , so-called variable pressure, to be constant at the contact waves. Therefore, for our adaptive scheme of SWEDs, we will utilize the central-upwind scheme which is derived from [9] for shallow water model with horizontal temperature variable (the system in [9] has similar properties with the SWEDs). In [9], the proposed second-order semi-discrete central-upwind scheme is capable of preserving the “lake at rest” steady state (1.2) and (1.3) as well as the positivity of the water depth and the temperature (the temperature variable is equivalent to the variable density in our work). In particular, to preserve the second type of “lake at rest” steady state (1.3) and suppress the pressure oscillations across the interface, an efficient interface tracking method is performed. The main idea of the interface approach in [9] is to completely avoid to use the information from the cells where two types of fluids are numerically mixed, so-called “mixed” cells when evolving the solution in the neighborhood of the interface. The data in the “mixed” cells is replaced by the interpolated values that are calculated using the reliable information from the nearby “single fluid” cells. Namely, the point values in “mixed” cells are obtained by using the approximated solution of the 1-D Riemann problems between the reliable “single fluid” cell averages. However, the central-upwind method and the interface tracking in [9] are designed on structured rectangular grids. In practice, one needs to deal with complicated geometries, where the use of triangular grids could be advantageous or even unavoidable. Hence, in this study, we extend the interface tracking method from [9] to unstructured triangular grids by utilizing the idea of the rotated coordinates proposed in [1]. The 1-D Riemann solver [9] is then performed in the direction of normal vectors on each edge of “mixed” cells. This extended central-upwind scheme maintains the well-balanced and positivity-preserving properties.

In addition to achieving the accuracy of the solution, minimizing computational cost is also one of the major challenge in the modeling and numerical analysis of hydrodynamics. The traditional numerical schemes for the system (1.1a)–(1.1d) are based on the use of very fine fixed meshes to reconstruct delicate features of the solution. However, this can lead to high computational cost, as well as poor resolution of all small scale features of the problem. In many engineering and scientific applications, it is beneficial to use adaptive meshes for improving the accuracy of the approximation at a much lower cost. Therefore, another goal of this work is to design an adaptive numerical algorithm for shallow water equations with variable density. There is some very recent effort on the design of adaptive well-balanced and positivity-preserving central-upwind schemes on quad-tree grids for shallow water models [37, 15, 14], but no research has been done for the development of such adaptive schemes for the shallow water with variable density on unstructured triangular grids.

As a part of the mesh reconstruction, we will need to project the data from the old mesh onto a new adaptive mesh. Since we avoid using the cell averages of mixed cells, the data of a new cell which is reconstructed from an old mixed cell is calculated based on the location of the density jumps and the data from the nearby “reliable” cells. Therefore, the interface separating different fluid phases must be accurately tracked or captured. In each “mixed” cell, we will reconstruct an approximate interface, which is called “interface reconstruction”. Many methods have been developed for this purpose in the last two decades, [31, 38, 17, 34, 42, 41]. In our work, we consider the interface reconstruction method derived in [45] due to the great advantages of this approach such as high-resolution, efficiency and simplicity.

This paper is organized as follows. In Section 2, we present a well-balanced positivity-preserving central-upwind scheme on unstructured triangular grids for SWEDs which serves as the underlying discretization for the developed adaptive algorithm. Next, in Section 3.1, we present the procedure to detect mixed cells. The interface approximation will be briefly reviewed in Section 3.2. In Section 3.3, we discuss the cell averages correction, which is used to prevent the density diffusion around the density jumps. We summarize the adaptive central-upwind method in Section 4.1. We discuss a strategy of adaptive mesh refinement in Section 4.2. In Section 4.3, we present an adaptive second-order strong stability preserving Runge-Kutta method, employed as a part of time evolution for the adaptive central-upwind scheme. We develop a local posteriori error estimator in Section 4.4 which is used as a robust indicator for the adaptive mesh refinement in our work. Finally, in Section 5, we illustrate the high accuracy and efficiency of the developed adaptive central-upwind scheme on a number of challenging tests for multi-fluid shallow water models.

2 The Central Upwind Scheme for SWEDs in Triangular Mesh

In this section, we focus on developing the central-upwind method for the SWEDs system (1.1a)-(1.1d), [15]. In particular, we will extend the adaptive scheme from the Saint-Venant system in [13] to the SWEDs system. The developed scheme will eliminate the oscillation appearing at the interface and ensure the well-balanced and positivity-preserving properties.

We first rewrite the SWEDs system (1.1a)-(1.1d) into the vector form as,

\mathbf{U}_{t}+\mathbf{F}(\mathbf{U},B)_{x}+\mathbf{G}(\mathbf{U},B)_{y}=\mathbf{S}(\mathbf{q},B),

(2.1)

where

U=(w,hu,hv,h\rho),

and the fluxes and source term are:

$\displaystyle\mathbf{F}(\mathbf{U},B)$	$\displaystyle=(hu,\dfrac{(hu)^{2}}{w-B}+\dfrac{g}{2\rho_{0}}\rho(w-B)^{2},\dfrac{(hu)(hv)}{w-B},hu\rho)^{T},$	(2.2)
$\displaystyle\mathbf{G}(\mathbf{q},B)$	$\displaystyle=(hv,\dfrac{(hu)(hv)}{w-B},\dfrac{(hv)^{2}}{w-B}+\dfrac{g}{2\rho_{0}}\rho(w-B)^{2},hv\rho)^{T},$
$\displaystyle\mathbf{S}(\mathbf{U},B)$	$\displaystyle=(0,-g\frac{\rho}{\rho_{0}}(w-B)B_{x},-g\frac{\rho}{\rho_{0}}(w-B)B_{y},0)^{T}.$

In our research, we consider the triangular mesh illustrated in Fig. 2.1 with the following notations.

Refer to caption — Figure 2.1: A typical triangular cell with three neighbors.

${\mathcal{T}}:=\bigcup_{j}T_{j}$ is an unstructured triangulation of the computational domain $\Omega$ ;

$T_{j}\in\mathcal{T}$ is a triangular cell of size $|T_{j}|$ with the barycenter $(x_{j},y_{j})$ ;

$V_{j\kappa}=(\widetilde{x}_{j\kappa},\widetilde{y}_{j\kappa}),~{}\kappa=12,23,31$ are the three vertices of $T_{j}$ ;

$T_{jk},~{}k=1,2,3$ are the neighboring triangles that share a common side with $T_{j}$ ;

$\ell_{jk}$ is the length of the common side of $T_{j}$ and $T_{jk}$ , and $M_{jk}$ is its midpoint;

$\bm{n}_{jk}:=(\cos(\theta_{jk}),\sin(\theta_{jk}))^{\top}$ is the outer unit normal to the $k$ th side of $T_{j}$ .

Next, the bottom topography $B$ is replaced with its continuous piecewise linear approximation $\widetilde{B}$ given by

\left|\begin{array}[]{ccc}x-\widetilde{x}_{j12}&y-\widetilde{y}_{j12}&{\widetilde{B}}(x,y)-\widehat{B}_{j12}\\[2.15277pt] \widetilde{x}_{j23}-\widetilde{x}_{j12}&\widetilde{y}_{j23}-\widetilde{y}_{j12}&\widehat{B}_{j23}-\widehat{B}_{j12}\\[2.15277pt] \widetilde{x}_{j13}-\widetilde{x}_{j12}&\widetilde{y}_{j13}-\widetilde{y}_{j12}&\widehat{B}_{j13}-\widehat{B}_{j12}\end{array}\right|=0,\quad(x,y)\in T_{j},

where, in the case of continuous bottom topography, $\widehat{B}_{j\kappa}:=B(V_{j\kappa})~{}\kappa=12,23,31$ . Then, denote:

B_{jk}:=\widetilde{B}(M_{jk}),\quad B_{j}:=\widetilde{B}(x_{j},y_{j})=\frac{1}{3}(\widehat{B}_{j12}+\widehat{B}_{j23}+\widehat{B}_{j13}).

At time $t$ , define by $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}(t)$ the approximation of the cell averages of the solution,

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}(t)\approx\frac{1}{|T_{j}|}\iint\limits_{T_{j}}\bm{U}(x,y,t)\,dxdy.

From now on, in order to shorten the formula, we suppress the time-dependence in the algorithm. All indexed quantities used in the following formula will be computed at time $t$ . Then, as shown in [29, 4], the semi-discrete second-order central-upwind scheme for the Saint-Venant system (1.1a)-(1.1d) on triangular grids is written as the following system of ODEs,

\displaystyle\frac{d\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}}{dt}=

\displaystyle-\frac{1}{|T_{j}|}\big{[}\bm{H}_{j1}+\bm{H}_{j2}+\bm{H}_{j3}\big{]}+\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}_{j},

(2.3)

where the numerical fluxes through the edges of the triangular cell $T_{j}$ are

$\displaystyle\bm{H}_{jk}=$	$\displaystyle\frac{\ell_{jk}\cos(\theta_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}{a_{jk}^{\rm in}\bm{F}(\bm{U}_{jk}(M_{jk}),B_{jk})+a_{jk}^{\rm out}\bm{F}(\bm{U}_{j}(M_{jk}),B_{jk})}\Big{]}$	(2.4)
$\displaystyle+$	$\displaystyle\frac{\ell_{jk}\sin(\theta_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}{a_{jk}^{\rm in}\bm{G}(\bm{U}_{jk}(M_{jk}),B_{jk})+a_{jk}^{\rm out}\bm{G}(\bm{U}_{j}(M_{jk}),B_{jk})}\Big{]}$	(2.5)
$\displaystyle-$	$\displaystyle\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\big{[}\bm{U}_{jk}(M_{jk})-\bm{U}_{j}(M_{jk})\big{]},\quad k=1,2,3.$	(2.6)

Here, $\bm{U}_{j}(M_{jk})$ and $\bm{U}_{jk}(M_{jk})$ are the reconstructed point values of $\bm{U}$ at the middle points of the edges $M_{jk}$ . To obtain these values, first, a piecewise linear reconstruction of the variables $\bm{\Upsilon}:=(w,u,v,\rho)^{\top}$ is computed as,

\widetilde{\bm{\Upsilon}}(x,y)=\sum\limits_{j}\bm{\Upsilon}_{j}(x,y){\Large{\mbox{$\chi$}}}_{T_{j}},\quad\bm{\Upsilon}_{j}(x,y):=\bm{\Upsilon}_{j}+(\widehat{\bm{\Upsilon}}_{x})_{j}(x-x_{j})+(\widehat{\bm{\Upsilon}}_{y})_{j}(y-y_{j}),

(2.7)

where ${\Large{\mbox{$\chi$}}}_{T_{j}}$ is the characteristic function of the cell $T_{j}$ , $\bm{\Upsilon}_{j}$ are the point values of $\bm{\Upsilon}$ at the cell centers, and $(\widehat{\bm{\Upsilon}}_{x})_{j}$ and $(\widehat{\bm{\Upsilon}}_{y})_{j}$ are the limited partial derivatives (see Section 3 in [29] for more details of the computation). In order to compute the cell center point values of the density, $\rho_{j}\approx\rho(x_{j},y_{j},t)$ in cell $T_{j}$ and velocities $u_{j}\approx u(x_{j},y_{j},t)$ and $v_{j}\approx v(x_{j},y_{j},t)$ , we use the desingularization procedure presented in [29]. After that, the second and third components of the point values $\bm{U}_{j}(M_{jk})$ and $\bm{U}_{jk}(M_{jk})$ are obtained from, $\bm{\Upsilon}_{j}(M_{jk})$ and $\bm{\Upsilon}_{jk}(M_{jk})$ ,

		$\displaystyle(hu)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,u_{j}(M_{jk}),$		$\displaystyle(hu)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,u_{jk}(M_{jk}),$
		$\displaystyle(hv)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,v_{j}(M_{jk}),$		$\displaystyle(hv)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,v_{jk}(M_{jk}).$
		$\displaystyle(h\rho)_{j}(M_{jk})=(w_{j}(M_{jk})-B_{jk})\,\rho_{j}(M_{jk}),$		$\displaystyle(h\rho)_{jk}(M_{jk})=(w_{jk}(M_{jk})-B_{jk})\,\rho_{jk}(M_{jk}).$

However, in some numerical examples, the oscillation may appear when we use the piecewise linear approximation (2.7) to obtain the point values in some cells. It may occur due to the appearance of local extrema values at middle points or as a consequence of the density discontinuity. To prevent this oscillation, we will use some techniques that will be discussed later in Sections 2.1 and 2.2.

In (2.6), $a_{jk}^{\rm in}$ and $a_{jk}^{\rm out}$ are the one-sided local speeds of propagation in the directions $\pm\bm{n}_{jk}$ . These speeds are related to the largest and smallest eigenvalues of the Jacobian matrix $J_{jk}=\cos(\theta_{jk})\,\frac{\partial{\bm{F}}}{\partial{\bm{U}}}+\sin(\theta_{jk})\,\frac{\partial{\bm{G}}}{\partial{\bm{U}}}$ , denoted by $\lambda_{+}[J_{jk}]$ and $\lambda_{-}[J_{jk}]$ , respectively, and are defined by

		$\displaystyle a^{\rm in}_{jk}=-\min\{\lambda_{-}[J_{jk}(\bm{U}_{j}(M_{jk}))],\,\lambda_{-}[J_{jk}(\bm{U}_{jk}(M_{jk})],\,0\},$		(2.8)
		$\displaystyle a^{\rm out}_{jk}=\max\{\lambda_{+}[J_{jk}(\bm{U}_{j}(M_{jk}))],\,\lambda_{+}[J_{jk}(\bm{U}_{jk}(M_{jk})],\,0\},$		(2.8)

where

		$\displaystyle\lambda_{\pm}[J_{jk}(\bm{U}_{j}(M_{jk}))]=\cos(\theta_{jk})u_{j}(M_{jk})+\sin(\theta_{jk})v_{j}(M_{jk})\pm\sqrt{\frac{g}{\rho_{0}}h_{j}(M_{jk})\rho_{j}(M_{jk})},$
		$\displaystyle\lambda_{\pm}[J_{jk}(\bm{U}_{jk}(M_{jk}))]=\cos(\theta_{jk})u_{jk}(M_{jk})+\sin(\theta_{jk})v_{jk}(M_{jk})\pm\sqrt{\frac{g}{\rho_{0}}h_{jk}(M_{jk})\rho_{jk}(M_{jk})}.$

Remark 2.1

In order to avoid division by 0 (or by a very small positive number), the numerical flux (2.6) is replaced with

	$\displaystyle\bm{H}_{jk}=$	$\displaystyle\frac{\ell_{jk}\cos(\theta_{jk})}{2}\left[\bm{F}(\bm{U}_{jk}(M_{jk}),B_{jk})+\bm{F}(\bm{U}_{j}(M_{jk}),B_{jk})\right]$
	$\displaystyle+$	$\displaystyle\frac{\ell_{jk}\sin(\theta_{jk})}{2}\left[\bm{G}(\bm{U}_{jk}(M_{jk}),B_{jk})+\bm{G}(\bm{U}_{j}(M_{jk}),B_{jk})\right]$

wherever $a_{jk}^{\rm in}+a_{jk}^{\rm out}<\sigma$ . In all of the reported numerical examples in Section 5, we have taken $\sigma=10^{-6}$ .

Finally, the cell average of the source term $\bm{S}_{j}$ in (2.3),

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}_{j}(t)\approx\frac{1}{|T_{j}|}\iint\limits_{T_{j}}\bm{S}\big{(}\bm{U}(x,y,t),B(x,y)\big{)}\,dxdy,

has to be discretized in a well-balanced manner presented later in Section 2.4

2.1 The Limiter for Piecewise Linear Approximation

As discussed above, we may observe the oscillation when using piecewise linear approximation (2.7) due to the appearance of local extrema at midpoints $M_{jk}$ , see [29]. Hence, in order to maintain the numerical stability of the scheme, the following monotonicity condition must be satisfied.

\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}\leq\bm{\Upsilon}^{(i)}_{j}(M_{jk})\leq\max\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)},\quad k=1,2,3,

(2.9)

where $\bm{\Upsilon}^{(i)}_{j}(M_{jk})$ is the $i$ -th component of the point values at midpoint $M_{jk}$ which is obtained by the linear reconstruction (2.7). Here, $\bm{\Upsilon}^{(i)}_{j}$ and $\bm{\Upsilon}^{(i)}_{jk}$ are the center point value in cell $T_{j}$ and the neighboring cell $T_{jk}$ , see [19]. To ensure that the reconstructed point value $\bm{\Upsilon}^{(i)}_{j}(M_{jk})$ is between two cell averages $\bm{\Upsilon}^{(i)}_{j}$ and $\bm{\Upsilon}^{(i)}_{jk}$ , in [29, 3], the gradients $\widehat{\bm{\Upsilon}}_{x}$ and $\widehat{\bm{\Upsilon}}_{y}$ are set to zero for cells violating at least one of the inequalities (2.9). However, using zero gradients may reduce the convergence rate in numerical simulation. Hence, in our research, instead of zero gradients, we consider the idea of the positivity-preserving limiter developed in [32, 43] to correct the piecewise linear reconstruction. The details for this correction are presented as follows.

Consider an arbitrary cell $T_{j}$ with the polynomial $\bm{\Upsilon}^{(i)}_{j}(x,y)$ that does not satisfy at least one of local extrema conditions (2.9). In that case, we will replace $\bm{\Upsilon}^{(i)}_{j}(x,y)$ by

(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)=\theta(\bm{\Upsilon}^{(i)}_{j}(x,y)-\bm{\Upsilon}^{(i)}_{j})+\bm{\Upsilon}^{(i)}_{j},

(2.10)

where $\theta\in[0,1]$ such that $(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)$ satisfies the inequalities (2.9) for all $k=1,2,3$ . We then need to solve for the constant $\theta$ . Plugging (2.10) into the inequalities (2.9), we have

\theta\in\left[\min(\alpha_{jk},\beta_{jk}),\max(\alpha_{jk},\beta_{jk})\right],\quad k=1,2,3

(2.11)

\theta\in[\max_{k}(\min(\alpha_{jk},\beta_{jk})),\min_{k}(\max(\alpha_{jk},\beta_{jk}))]

(2.12)

where $\alpha_{jk}=\dfrac{\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}-\bm{\Upsilon}^{(i)}_{j}}{\bm{\Upsilon}^{(i)}_{j}(M_{jk})-\bm{\Upsilon}^{(i)}_{j}},\quad\beta_{jk}=\dfrac{\max\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}-\bm{\Upsilon}^{(i)}_{j}}{\bm{\Upsilon}^{(i)}_{j}(M_{jk})-\bm{\Upsilon}^{(i)}_{j}}.$

Hence we can take

\theta=\min(\min_{k}(\max(\alpha_{jk},\beta_{jk})),1),

(2.13)

Note that $\theta\geq 0$ since $\theta=0$ , which corresponds to the constant approximation for $(\bm{\Upsilon}^{(i)}_{j})^{*}(x,y)$ , is one solution of the inequalities (2.9). Hence, with this optimized limiter, the point values stay within the given range. In addition, the approximated values of water depth and density at middle points are also positive as $\bm{\Upsilon}^{(i)}_{j}(M_{jk})\geq\min\big{(}\bm{\Upsilon}^{(i)}_{j},\bm{\Upsilon}^{(i)}_{jk}\big{)}>0$ . Therefore, the stability is ensured without using zero-gradients.

2.2 Riemann Solver for Mixed Cells

As mentioned in Section 1 and [9], a good scheme should not develop spurious pressure oscillations in the neighborhood of density jumps. The oscillations appear when we numerically solve the compressible multifluid systems by using conventional Godunov-type finite-volume methods. We define the cells which contains the curves of density discontinuity as “mixed” cells. Note that the quantities in mixed cells are calculated as an artificial numerical mixture of two different fluids. Hence, the cell averages of mixed cells may have no or very little physical sense and become ’unreliable’ [9]. Consequently, using the cell averages in mixed cells to obtain the linear approximation (2.7) in the neighborhood of the interface will lead to unexpected pressure oscillations. In [9], to eliminate the oscillation, instead of the linear approximation (2.7), the point values in the “mixed” cells are calculated by using the approximated solution of 1-D Riemann problems. Namely, according to this approach, the 1-D Riemann Solver is applied in the $x$ -direction to approximate the point values on the vertical boundaries of each rectangular mixed cell. Similarly, for the point values on horizontal sides of the rectangular, 1-D Riemann problem is considered in $y$ -direction. This leads to the idea of applying Riemann Solver approach in the direction of the normal vector, so called “normal direction”, to get the middle point values on each side of the mixed cells in triangular meshes. We can understand this idea as rotating the reference coordinate. Recently, a method using rotated coordinate frame to solve multi-dimensional Riemann problems has been proposed in [1]. In such approach, the 2-D Riemann Problem is converted to 1-D problem in a particular direction in order to easily obtain the corresponding solution in 2-D. Therefore, in our research, we will calculate the point values in triangles based on the intermediate states of 1-D Riemann problem in the normal direction as follows.

Assume that triangle $T_{j}$ is a mixed cell which has the outward unit normal vector $\mathbf{n}_{jk}=(\cos(\theta_{jk}),\sin(\theta_{jk}))$ on side $k$ . We will consider a new reference frame $(x^{\prime},y^{\prime})$ such that the new horizontal axis is in the direction of the outward normal vector $\mathbf{n}_{jk}$ and the origin $(0,0)$ is at the middle point $M_{jk}$ of side $k$ as shown in Fig. 2.2.

Suppose $T_{L}$ and $T_{R}$ to be the two closest single-fluid cells to $T_{j}$ such that $T_{L}$ and $T_{R}$ stay on different sides of the $k$ -th edge of cell $T_{j}$ . We assume that $T_{R}$ stays on the side where the outward normal vector $\mathbf{n}_{jk}$ is pointing. Let $\bm{\Upsilon}_{L}=(w_{L},u_{L},v_{L},\rho_{L})$ and $\bm{\Upsilon}_{R}=(w_{R},u_{R},v_{R},\rho_{R})$ be the center point values of single-fluid cells $T_{L}$ and $T_{R}$ . Note that $(u_{L},v_{L})$ and $(u_{R},v_{R})$ are the values of velocities in the Cartesian coordinate system $(x,y)$ . The projections of velocities $(u_{L},v_{L})$ and $(u_{R},v_{R})$ onto the rotated frame $(x^{\prime},y^{\prime})$ are

u^{\prime}_{L}=u_{L}\cos(\theta_{jk})+v_{L}\sin(\theta_{jk}),\quad v^{\prime}_{L}=-u_{L}\sin(\theta_{jk})+v_{L}\cos(\theta_{jk}),

u^{\prime}_{R}=u_{R}\cos(\theta_{jk})+v_{R}\sin(\theta_{jk}),\quad v^{\prime}_{R}=-u_{R}\sin(\theta_{jk})+v_{R}\cos(\theta_{jk}).

Similar to [1, 9], to compute the point value at midpoint $M_{jk}$ on side $k$ of triangle $T_{j}$ , we have to solve the following 1-D Riemann problems between states $\bm{\Upsilon}_{L}$ an $\bm{\Upsilon}_{R}$ in direction $\mathbf{n}_{jk}$ .

\begin{cases}w_{t}+(hu)_{x^{\prime}}=0\\ (hu)_{t}+\Big{(}\dfrac{(hu)^{2}}{w-B}+\frac{g}{2}\frac{\rho}{\rho_{0}}(w-B)\Big{)}_{x^{\prime}}=-g\dfrac{\rho}{\rho_{0}}(w-B)B_{x^{\prime}},\\ (hv)_{t}+\left(\dfrac{(hu)(hv)}{w-B}\right)_{x^{\prime}}=0,\\ (h\frac{\rho}{\rho_{0}})_{t}+\left(uh\dfrac{\rho}{\rho_{0}}\right)_{x^{\prime}}=0,\end{cases}

(2.14)

subject to the following initial condition

(w,u,v,\rho,B)(x^{\prime},0)=\begin{cases}(\bm{\Upsilon}_{L})^{\prime}:=(w_{L},u^{\prime}_{L},v^{\prime}_{L},\rho_{L},B_{L})\quad\mbox{if}\quad x^{\prime}<0,\\ (\bm{\Upsilon}_{R})^{\prime}:=(w_{R},u^{\prime}_{R},v^{\prime}_{R},\rho_{R},B_{R})\quad\mbox{if}\quad x^{\prime}>0.\end{cases}

(2.15)

The details of the approximate Riemann solver for the above Riemann problem can be seen in [9]. Suppose $(\bm{\Upsilon}^{*}_{L})^{\prime}:=(w^{*}_{L},(u^{*}_{L})^{\prime},(v^{*}_{L})^{\prime},\rho^{*}_{L})$ and $(\bm{\Upsilon}^{*}_{R})^{\prime}:=(w^{*}_{R},(u^{*}_{R})^{\prime},(v^{*}_{R})^{\prime},\rho^{*}_{R})$ are the intermediate states calculated by Riemann solver for (2.14)-(2.15) in $(x^{\prime},y^{\prime})$ coordinates. We recall that in [9], the point values at midpoints on left and right boundaries of a rectangular mixed cell are obtained respectively based on the left and right Riemann intermediate states, $(\bm{\Upsilon}^{*}_{L})^{\prime}$ and $(\bm{\Upsilon}^{*}_{R})^{\prime}$ . Here, in triangular grids, we can choose either $(\bm{\Upsilon}^{*}_{L})^{\prime}$ or $(\bm{\Upsilon}^{*}_{R})^{\prime}$ to calculate the point values at the middle point $M_{jk}$ of cell $T_{j}$ . Note that $T_{L}$ and $T_{R}$ are both single-fluid cells nearby mixed cell $T_{j}$ . However, if the neighboring cell $T_{jk}$ shown in Fig. 2.2 is a single-fluid cell, then $T_{R}\equiv T_{jk}$ and $M_{jk}\in T_{R}$ , which means $T_{R}$ is closer to $M_{jk}$ compared to $T_{L}$ . Hence, in our work, we select the right intermediate values $(\bm{\Upsilon}^{*}_{R})^{\prime}$ and the information from $T_{R}$ to compute the point values $\bm{\Upsilon}_{j}(M_{jk})$ . From [9], we first calculate the sound of speed $c^{*}_{R}=h^{2}\dfrac{g\rho}{2\rho_{0}}$ . If $h^{*}_{R}>0,\rho^{*}_{R}>0$ and $(u^{*})^{\prime}_{R}-c^{*}_{R}<0$ , then we replace the piecwise linear approximation in (2.7) by $\bm{\Upsilon}_{j}(M_{jk})=(\bm{\Upsilon}^{*}_{R})^{\prime}$ , otherwise $\bm{\Upsilon}_{j}(M_{jk})=(\bm{\Upsilon}_{R})^{\prime}$ (more details of the point values correction can be seen in [9]). Finally, we convert the computed solution back to original Cartesian coordinate $(x,y)$ as

\bm{\Upsilon}_{j}(M_{jk})=\begin{cases}\bm{\Upsilon}^{*}_{R},\quad&\mbox{if}\quad h^{*}_{R}>0,\rho^{*}_{R}>0,(u^{*})^{\prime}_{R}-c^{*}_{R}<0,\\ \bm{\Upsilon}_{R},\quad&\mbox{otherwise},\end{cases}

(2.16)

where $\bm{\Upsilon}^{*}_{R}:=(w^{*}_{R},(u^{*}_{R})^{\prime}\cos(\theta_{jk})-(v^{*}_{R})^{\prime}\sin(\theta_{jk}),(u^{*}_{R})^{\prime}\sin(\theta_{jk})+(v^{*}_{R})^{\prime}\cos(\theta_{jk}),\rho^{*}_{R})$ .

From [9], in “lake at rest” solutions (1.2) and (1.3), we have $\bm{\Upsilon}^{*}_{R}=\bm{\Upsilon}_{R}$ and the steady state solutions are then preserved by using the discretization of the source term presented in Section 2.4.

2.3 Positivity-preserving condition

Next, we will discuss the restriction of time step to preserve the positivity of water depth $h$ and the variable $h\rho$ . Consider the Forward Euler (FE) equation

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt${U}$\kern-0.50003pt}}}_{j}^{n+1}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt${U}$\kern-0.50003pt}}}_{j}^{n}-\frac{1}{|T_{j}|}\sum_{k=1}^{3}\Delta t_{jk}{H}_{jk}+\Delta t\,\bar{{S}}_{j}.

From [3], the time step condition to guarantee the positivity of water height $h$ is

\Delta t<\dfrac{1}{6a}\min_{jk}(r_{jk}).

(2.17)

One can clearly see that the time step size (2.17) will also achieve nonnegative water height $h$ in the density shallow water model (1.1a)-(1.1d), see [3, 29]. We now use the idea from [3] to find the CFL-type condition for the positivity-preserving property of the last variable $h\rho$ . To this end, we apply the forward Euler discretization to the last component of the scheme.

$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n+1}_{j}=(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n}_{j}-$	$\displaystyle\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\cos(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}(h\rho u)_{jk}(M_{jk})+a_{jk}^{\rm out}(h\rho u)_{j}(M_{jk})\Big{]}$	(2.18)
$\displaystyle-$	$\displaystyle\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\sin(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}(h\rho v)_{jk}(M_{jk})+a_{jk}^{\rm out}(h\rho v)_{j}(M_{jk})\Big{]}$
$\displaystyle+$	$\displaystyle\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\big{[}(h\rho)_{jk}(M_{jk})-(h\rho)_{j}(M_{jk})\big{]}.$

Note that

(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n}_{j}=h_{j}\rho_{j}=\left(\dfrac{1}{3}\sum_{m=1}^{3}h_{j}(M_{jm})\right)\left(\dfrac{1}{3}\sum_{s=1}^{3}\rho_{j}(M_{js})\right)=\dfrac{1}{9}\sum_{m,s}h_{j}(M_{jm})\rho_{j}(M_{js}).

(2.19)

Plug (2.19) to (2.18), we have

	$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n+1}_{j}=$	$\displaystyle\dfrac{1}{9}\sum_{m,s}h_{j}(M_{jm})\rho_{j}(M_{js})$
		$\displaystyle-\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\cos(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})u_{jk}(M_{jk})+a_{jk}^{\rm out}h_{j}(M_{jk})\rho_{j}(M_{jk})u_{j}(M_{jk})\Big{]}$
		$\displaystyle-\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\sin(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})v_{jk}(M_{jk})+a_{jk}^{\rm out}h_{j}(M_{jk})\rho_{j}(M_{jk})v_{j}(M_{jk})\Big{]}$
		$\displaystyle+\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\big{[}h_{jk}(M_{jk})\rho_{jk}(M_{jk})-h_{j}(M_{jk})\rho_{j}(M_{jk})\big{]}$
	$\displaystyle=$	$\displaystyle\dfrac{1}{9}\sum_{m\neq s}h_{j}(M_{jm})\rho_{j}(M_{js})$
		$\displaystyle+\sum_{k=1}^{3}h_{j}(M_{jk})\rho_{j}(M_{jk})\Big{[}\dfrac{1}{9}-\dfrac{\Delta t}{\|T_{j}\|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm in}+u^{\perp}_{j}(M_{jk})\right)\Big{]}$
		$\displaystyle+\sum_{k=1}^{3}h_{jk}(M_{jk})\rho_{jk}(M_{jk})\dfrac{\Delta t}{\|T_{j}\|}\frac{\ell_{jk}a_{jk}^{\rm in}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm out}-u^{\perp}_{jk}(M_{jk})\right).$

From the definitions of the local speeds (2.8) we obtain that $u^{\perp}_{j}(M_{jk})=\cos(\theta_{jk})u_{j}(M_{jk})+\sin(\theta_{jk})v_{j}(M_{jk})\leq a_{jk}^{\rm out}$ , $0\leq a_{jk}^{\rm in}$ , and $0\leq a_{jk}^{\rm out}$ . Besides, the corrected reconstruction in Section 2.1 guarantees that $h_{jk}(M_{jk})\geq 0$ and $\rho_{jk}(M_{jk})\geq 0$ for all $j$ and $k=1,2,3$ . Therefore, all terms in the first and third sum on the RHS of (2.3) are nonnegative. To enforce the second sum on the RHS of (2.3) to be nonegative, we then need:

\dfrac{\Delta t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm in}+u^{\perp}_{j}(M_{jk})\right)\leq\dfrac{\Delta t}{|T_{j}|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm in}+a_{jk}^{\rm out}\right)=\dfrac{\Delta t}{|T_{j}|}l_{jk}a_{jk}^{\rm out}<\dfrac{1}{9},\quad\forall j\mbox{ and }k=1,2,3.

Hence we conclude that all terms in the second sum on the RHS of (2.3) are also nonnegative under the following CFL-type condition.

\Delta t<\dfrac{1}{18a}\min_{jk}(r_{jk})=\dfrac{\min_{jk}\left(\dfrac{2|T_{j}|}{l_{jk}}\right)}{18\max_{jk}(a_{jk}^{\rm in},a_{jk}^{\rm out})}<\dfrac{1}{9}\min_{jk}\left(\dfrac{|T_{j}|}{l_{jk}a^{\rm out}_{jk}}\right),

(2.20)

where $a=\max_{jk}(a_{jk}^{\rm in},a_{jk}^{\rm out})$ and $r_{jk}=\dfrac{2|T_{j}|}{l_{jk}}$ is the $k$ -th altitude of triangle $T_{j}$ . This completes the proof. This proof is still valid if one uses a higher-order SSP ODE solver (either the Runge-Kutta or the multistep one), because such solvers can be written as a convex combination of several forward Euler steps.

2.4 Well-balanced discretization of source term.

In this section, we first develop a well-balanced discretization of the source terms, which guarantees the first type of “lake at rest” steady-state solutions,

w=\max\big{\{}C,B(x,y)\big{\}},\quad C=\rm{Const},\quad\rho=P\equiv Const,\quad u\equiv v\equiv 0,

(2.21)

are exactly preserved by the resulting central-upwind scheme. This means that the source discretization $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\mathbf{S}$\kern-0.50003pt}}}_{j}$ should exactly balance the numerical fluxes so that the right-hand side (RHS) of (1.1a)-(1.1d) vanishes at “lake at rest” steady states.

To this end, we substitute the “lake at rest” state (2.21) into FE scheme of the system (it is still correct with the adaptive Runge-Kutta scheme) and conclude that a well-balanced quadrature for $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\mathbf{S}$\kern-0.50003pt}}}_{j}$ should satisfy the following two conditions:

-\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\cos(\theta_{jk})\cdot\frac{\Delta t^{(\rho)}_{jk}}{\Delta t}\cdot\dfrac{P}{\rho_{0}}\big{[}C-B(M_{jk})\big{]}^{2}+\bar{S}_{j}^{\,(2)}=0,

(2.22)

-\frac{g}{2|T_{j}|}\sum_{k=1}^{3}\ell_{jk}\sin(\theta_{jk})\cdot\frac{\Delta t^{(\rho)}_{jk}}{\Delta t}\cdot\dfrac{P}{\rho_{0}}\big{[}C-B(M_{jk})\big{]}^{2}+\bar{S}_{j}^{\,(3)}=0,

(2.23)

where

\bar{S}_{j}^{\,(2)}\approx-\frac{g}{|T_{j}|\rho_{0}}\iint\limits_{T_{j}}P(C-B(x,y))B_{x}(x,y)\,dxdy,\quad\bar{S}_{j}^{\,(3)}\approx-\frac{g}{|T_{j}|\rho_{0}}\iint\limits_{T_{j}}P(C-B(x,y))B_{y}(x,y)\,dxdy.

In order to derive a well-balanced quadrature, similar to [3, 29], we first apply Green’s formula, $\iint_{T_{j}}{\rm div}\,\mathbf{{\cal G}}\,dxdy=\int_{\partial T_{j}}\mathbf{{\cal G}}\cdot\mathbf{n}\,ds$ , to the vector field $\mathbf{{\cal G}}=(\dfrac{1}{2}\rho(x,y)(w(x,y)-B(x,y))^{2},0)$ and obtain

$\displaystyle-\iint\limits_{T_{j}}\rho(x,y)(w(x,y)-B(x,y))B_{x}(x,y)\,dxdy=$	$\displaystyle\sum_{k=1}^{3}\int\limits_{(\partial{T_{j}})_{k}}\rho(x,y)\frac{(w(x,y)-B(x,y))^{2}}{2}\cos(\theta_{jk})\,ds$	(2.24)
	$\displaystyle-\iint\limits_{T_{j}}\rho(x,y)(w(x,y)-B(x,y))w_{x}(x,y)\,dxdy$
	$\displaystyle-\iint\limits_{T_{j}}\rho_{x}(x,y)\frac{(w(x,y)-B(x,y))^{2}}{2}\,dxdy,$

where $(\partial{T_{j}})_{k}$ is the $k$ -th side of the triangle $T_{j}$ , $k=1,2,3$ . The double integrals are approximated using the trapezoidal rule. This results in the following quadrature for $\bar{S}^{\,(2)}_{j}$ :

$\displaystyle\bar{S}^{\,(2)}_{j}$	$\displaystyle=\frac{g}{2\|T_{j}\|}\sum_{k=1}^{3}\ell_{jk}\cos(\theta_{jk})\cdot\frac{\Delta t^{(\rho)}_{jk}}{\Delta t}\cdot\dfrac{\rho(M_{jk})}{\rho_{0}}\big{[}w(M_{jk})-B(M_{jk})\big{]}^{2}$	(2.25)
	$\displaystyle-\frac{g}{3\rho_{0}}\Big{[}\rho_{j12}(w_{j12}-\widehat{B}_{j12})\,w_{x}(V_{j12})+\rho_{j23}(w_{j23}-\widehat{B}_{j23})\,w_{x}(V_{j23})+\rho_{j13}(w_{j13}-\widehat{B}_{j13})\,w_{x}(V_{j13})\Big{]}.$
	$\displaystyle-\frac{g}{6\rho_{0}}\Big{[}\rho_{x}(V_{j12})(w_{j12}-\widehat{B}_{j12})^{2}\,+\rho_{x}(V_{j23})(w_{j23}-\widehat{B}_{j23})^{2}\,+\rho_{x}(V_{j13})(w_{j13}-\widehat{B}_{j13})^{2}\,\Big{]}.$

A similar quadrature for $\bar{S}^{\,(3)}_{j}$ is

$\displaystyle\bar{S}^{\,(3)}_{j}$	$\displaystyle=\frac{g}{2\|T_{j}\|}\sum_{k=1}^{3}\ell_{jk}\sin(\theta_{jk})\cdot\frac{\Delta t^{(\rho)}_{jk}}{\Delta t}\cdot\dfrac{\rho(M_{jk})}{\rho_{0}}\big{[}w(M_{jk})-B(M_{jk})\big{]}^{2}$	(2.26)
	$\displaystyle-\frac{g}{3\rho_{0}}\Big{[}\rho_{j12}(w_{j12}-\widehat{B}_{j12})\,w_{y}(V_{j12})+\rho_{j23}(w_{j23}-\widehat{B}_{j23})\,w_{y}(V_{j23})+\rho_{j13}(w_{j13}-\widehat{B}_{j13})\,w_{y}(V_{j13})\Big{]}.$
	$\displaystyle-\frac{g}{6\rho_{0}}\Big{[}\rho_{y}(V_{j12})(w_{j12}-\widehat{B}_{j12})^{2}\,+\rho_{y}(V_{j23})(w_{j23}-\widehat{B}_{j23})^{2}\,+\rho_{y}(V_{j13})(w_{j13}-\widehat{B}_{j13})^{2}\,\Big{]}.$

Notice that the piecewise linear reconstruction procedure ensures that at the steady state (2.21), $\nabla\rho(V_{j\kappa})=0,u_{x}=v_{y}=0$ and $\nabla w(V_{j\kappa})=0$ throughout the entire computational domain. This implies that $(w_{j\kappa}-\widehat{B}_{j\kappa})\,w_{x}(V_{j\kappa})\equiv(w_{j\kappa}-\widehat{B}_{j\kappa})\,w_{y}(V_{j\kappa})\equiv 0$ and $\rho_{x}(V_{j\kappa})(w_{j\kappa}-\widehat{B}_{j\kappa})^{2}\equiv\rho_{y}(V_{j\kappa})(w_{j\kappa}-\widehat{B}_{j\kappa})^{2}\equiv 0$ . Therefore, the quadratures (2.25) and (2.26) satisfy the desired well-balanced requirements (2.22) and (2.23).

Next, we consider the “lake at rest” situation

B\equiv\rm{Const},\quad h^{2}\rho=Q\equiv Const,\quad u\equiv v\equiv 0.

(2.27)

Note that in the region occupied by only one fluid, the density $\rho$ is a constant and water surface $w$ is then also a constant based on (2.27). Hence, in single-fluid cells, the steady state (2.27) is equivalent to the solution (2.21) which is maintained by the discretization of source term (2.25)-(2.26). In addition, in mixed cells, the point values are obtained by the data from nearby single-fluid cells thanks to the Riemann Solver approximation presented in Section 2.2, see also [9]. Therefore, the discretization of source term (2.25)-(2.26), is also capable of preserving the steady state solution (2.27).

3 Interface Tracking

To achieve an efficient adaptive scheme for multifluid flows, it is essential to accurately capture the curve where two types of fluid join simultaneously with the flow field evolution. The interface tracking is important in preventing excessive numerical diffusion of variable density, see Section 3.3. We also need the location of density jumps to exactly update the cell averages in the adaptive mesh reconstruction, Section 4. Therefore, we desire a simple and effective interface reconstruction for the system (1.1a)-(1.1d). Over the last two decades, many methods have been proposed for this purpose such as the level set method [31, 38], the volume of fluid method [17, 34], and the front tracking method [42, 41]. Among various versions of interface tracking, we have considered the approach described in [45, 44] for our work due to its simplicity, robustness, and high efficiency. In particular, the interface is obtained by using both the level set and the volume fraction functions. The details of the interface reconstruction is described as follows.

3.1 Mixed Cell Detecting by Level Set Function

In the first part of this section, we will discuss the method of detecting the mixed cells in the triangular grids. For this end, we consider the level set function $\phi$ which is defined such that it is positive in one fluid, is negative in the other fluid, and has zero value at the interface. Very often, the level set value of each grid point is initialized by the signed distance from that point to the curve of density discontinuities, see [29, 31, 38, 45]. As discussed in [29, 31, 38, 45], level set function is evolved by the velocity $(u,v)$ of the flow field as

\phi_{t}+u\phi_{x}+v\phi_{y}=0.

(3.1)

The equation (3.1) can be rewritten in conservation form as follows.

\phi_{t}+(u\phi)_{x}+(v\phi)_{y}=(u_{x}+v_{y})\phi.

(3.2)

We have used the central-upwind scheme presented in Section 2 to solve for $\phi$ , where the point values $\phi(M_{jk}),k=1,2,3$ in triangle $T_{j}$ are approximated by the piecewise linear reconstruction (2.7). The integral of the source term in (3.2) can be approximated by midpoint rule. Note that using the central-upwind scheme to solve (3.1) only give us the cell averages of $\phi$ in each cell and does not point out which cells contain the interface $\phi=0$ . We have to provide a numerical method for detecting mixed cells. However, a triangle is a single-fluid cell, also called as “reliable” cell, if the level set values at its vertices are either all positive or all negative. Hence, we will locate each node in the grid based on its point value of level set as presented below.

Without the loss of generality, we assume that $\phi(x,y)$ is positive in fluid 1 and $\phi(x,y)$ is negative in fluid 2. The point value of level set function $\phi^{*}_{j\kappa}$ at each vertex $V_{j\kappa}$ of cell $T_{j}$ is approximated by extrapolation [9]

\phi^{*}_{j\kappa}=\dfrac{\displaystyle\sum_{i=0}^{m^{*}_{\kappa}}c_{i}\phi^{i}_{j\kappa}}{\displaystyle\sum_{i=0}^{m^{*}_{\kappa}}c_{i}},

(3.3)

where $\phi^{i}_{j\kappa},$ is the cell averages of level set function in cells $T^{i}_{j\kappa}$ which have common vertex at $V_{j\kappa}$ , and the weight $c_{i}$ is inversely proportional to the distance between the center of cells $T^{i}_{j\kappa},i=0,1,...,m^{*}_{\kappa}$ and vertex $V_{j\kappa}$ . If $\phi^{*}_{j\kappa}>0$ , we have the vertex $V_{j\kappa}$ staying in fluid 1. Otherwise, if $\phi^{*}_{j\kappa}<0$ , the vertex $V_{j\kappa}$ is in fluid 2. Finally, based on the physical meaning of the level set function, a triangle $T_{j}$ is naturally marked as a mixed cell if it has two vertices staying in different fluids.

3.2 Interface Reconstruction

In most works of interface tracking, see [31, 38, 17, 34, 42, 41], the curve of density discontinuity in a mixed cell is approximated as a linear line segment of the form $\mathbf{n}\cdot x=\alpha$ , where $\mathbf{n}$ is the normal vector of the interface, $x$ is the location of a point, and $\alpha$ is the line constant. A variety of methods for computing normal vector $\mathbf{n}$ and parameter $\alpha$ have been briefly introduced and compared in [45, 33]. Based on the advantages and disadvantages of conventional interface tracking methods, an innovative approach has been proposed in [45]. The approach employs both the level set and volume of fluid methods to denote the density segment in each flagged cell by its two endpoints. Namely, the interface normal vector is calculated from the level set function while the exact location of two endpoints of the segment is determined by enforcing mass conservation based on the volume fraction. In our work, we will apply this interface tracking method in the adaptive triangular mesh due to its simplicity, accuracy, and versatility. The details of this method can be seen in [45]. In the following, we will briefly present the main steps of the interface reconstruction.

Once the mixed cells are detected, we will first compute the interface normal vector in each flagged cell based on the level set function and a least square problem. Namely, for each mixed cell $T_{j}$ with center at $(x_{j},y_{j})$ , we consider a stencil that consists of all centers $(x_{i},y_{i}),i=1,2,...,N$ of $N$ cells that share at least a vertex with $T_{j}$ (to simplify the computation, we do not place the vertices of cell $T_{j}$ in the stencil as in [45]). A quadratic function for function $\phi$ of a point $(x,y)$ is then given in the generic form

\phi=ax^{\prime 2}+bx^{\prime}y^{\prime}+cy^{\prime 2}+dx^{\prime}+ey^{\prime}+g,

where $(x^{\prime}:=x-x_{j},y^{\prime}:=x-x_{j})$ is the location of point $(x,y)$ in the local coordinate frame $x^{\prime}$ and $y^{\prime}$ obtained by shifting the original coordinate such that the new origin is at the center of $T_{j}$ . The coefficients $a,b,c,d,e,$ and $g$ are determined by using the least squares method with the linear system

Qs=r,

where

Q=\begin{bmatrix}x^{\prime 2}_{1}&x^{\prime}_{1}y^{\prime}_{1}&y^{\prime 2}_{1}&x^{\prime}_{1}&y^{\prime}_{1}&1\\ x^{\prime 2}_{2}&x^{\prime}_{2}y^{\prime}_{2}&y^{\prime 2}_{2}&x^{\prime}_{2}&y^{\prime}_{2}&1\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ x^{\prime 2}_{N}&x^{\prime}_{N}y^{\prime}_{N}&y^{\prime 2}_{N}&x^{\prime}_{N}&y^{\prime}_{N}&1\\ \end{bmatrix},\quad s=\begin{bmatrix}a\\ b\\ c\\ d\\ e\\ g\end{bmatrix},\quad\phi=\begin{bmatrix}\phi_{1}\\ \phi_{2}\\ \vdots\\ \phi_{n}\end{bmatrix},

(3.4)

$(x^{\prime}_{i}=x_{i}-x_{j},y^{\prime}_{i}=y_{i}-y_{j}),i=1,...,N,$ are the local coordinates of the $i-$ th node in the stencil and $\phi_{i},i=1,...,N,$ is the cell average of level set function in cell $T_{i}$ corresponding to node $i$ -th. Once the coefficients are known, we compute the unit normal vector by $\mathbf{n}=\left(\dfrac{d}{\sqrt{d^{2}+e^{2}}},\dfrac{e}{\sqrt{d^{2}+e^{2}}}\right).$ The system (3.4) is solvable as explained in [45]. In addition, since the level set function is continuous, the normal calculation is second-order accurate.

Next, we continue to determine the exact coordinates of two endpoints of the interface by using the volume of fluid approach (VOF). The VOF method has been widely used and developed for interface tracking in [17, 34, 44, 45]. In the VOF method, the volume fraction function, denoted by $f$ , is defined as the ratio of the volume of one fluid, called fluid 1, in each cell to the total volume of the cell. Hence, $f$ is unity if the cell is a single-fluid cell in fluid 1, and is zero if the cell only contains fluid 2. For mixed cells, we have $0<f<1$ . The function $f$ is advected by

f_{t}+uf_{x}+vf_{y}=0.

(3.5)

The main idea of the VOF method proposed in [44, 45] is to reconstruct the interface segment by determining the coordinates of its endpoints. In particular, the two endpoints must form a segment which has the normal vector $\mathbf{n}$ as computed above and the interface truncates the cell with the given volume fraction. Fig. 3.1 is an illustration for two cases of a mixed cell $T_{j}$ where the interface $I_{1}I_{2}$ splits $T_{j}$ into two parts, $T^{(1)}_{j}$ and $T^{(2)}_{j}$ , respectively occupied by fluid 1 and fluid 2.

We then have

I_{1}I_{2}\cdot\mathbf{n}=0,\quad\dfrac{T^{(1)}_{j}}{T_{j}}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$f$\kern-0.50003pt}}}_{j},

where $\mathbf{n}$ is the normal vector of the interface computed by (3.4) and $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$f$\kern-0.50003pt}}}_{j}$ is the cell average of volume fraction in mixed cell $T_{j}$ obtained by using the central-upwind method, see Section 2, to solve (3.5). This interface reconstruction is explicit, accurate, and capable of conserving the volume of the fluid. The details of endpoint calculation can be seen in [45]. Finally, the reconstructed interface will be used in the adaptive mesh reconstruction as discussed in the following sections.

3.3 The Cell Averages Correction

The idea of correcting the solution in the neighborhood of the interface has been used in numerical methods for multifluid flows [9, 45] to improve the computed solution. This technique prevents the diffusion of density emerging when we use the central-upwind method to solve compressible systems. In our work, we will also perform the correction procedure such that the local conservation is ensured based on the location of density jumps. Namely, suppose we have two types of fluid, fluid 1 and fluid 2, in the flow. After determining the single-fluid cells and mixed cells by using the point values of the level set function, see Section 3.1, the cell average correction will proceed as follows.

•

If a cell $T_{j}$ is a single fluid cell in fluid $i=1,2$ , we correct the cell averages of $h\rho$ by $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{j}=h_{j}\rho^{(i)}_{j}$ , where $\rho^{(i)}_{j}$ is the value of density in fluid $i=1,2$ . $h_{j}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}_{j}-B_{j}$ is the cell average of water depth in $T_{j}$ . To preserve the mass of $h\rho$ in the domain, we equally split the change of $h\rho$ in $T_{j}$ into the nearby mixed cells. Namely, if the neighboring cell $T_{jk},k=1,2,3$ of $T_{j}$ is a mixed cell, we obtain the new cell average of $h\rho$ in $T_{jk}$ by

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{jk}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{jk}+\dfrac{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{new}_{j}}{n_{mix}},

(3.6)

where $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{j}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}_{jk}$ are the old cell averages computed by using the central-upwind scheme in triangles $T_{j}$ and $T_{jk}$ , and $n_{mix}$ is the number of mixed cells surrounding the cell $T_{j}$ .

•

If $T_{j}$ is a mixed cell, except the update (3.6) from its nearby single-fluid cells, no further correction is needed.

Due to the fact that the interface normally moves from one cell to its adjacent cells, only cells surrounding the interface are updated and the loss of mass from this correction is negligible. In addition, for cells in the “lake at rest” area, this correction will not change the existing solution therefore maintaining the well-balanced properties.

4 Adaptive Central-Upwind Scheme

The traditional numerical methods for system (1.1a)–(1.1d) consider very fine fixed meshes to reconstruct delicate features of the solution. However, this can lead to high computational cost, as well as to a poor accuracy of small scale characteristics of the problem. Therefore, we will use adaptive meshes to improve the accuracy of the approximation at a much lower cost. In this section, we will review an efficient and accurate adaptive central-upwind algorithm from [13] which we adapt to the system (1.1a)-(1.1d).

4.1 Adaptive Central-Upwind Algorithm

From [13], the adaptive central-upwind algorithm for the system (1.1a)- (1.1d) is described briefly by the following steps.

Step 0. At time $t=t^{0}$ , generate the initial uniform grid $\mathcal{T}^{0,0}$ .

Step 1. On mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ , evolve the cell averages $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n}$ of the solution from time $t^{n}$ to $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n+1}$ at the next time level $t^{n+1}$ using adaptive central-upwind scheme (4.5a)-(4.5b), see Section 4:

•

At time $t^{n}$ , determine the level $l=0,1,...,L$ of each triangle $T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ , (4.1), Section 4.3.
•

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$ , perform the piecewise polynomial reconstruction (2.7) for single-fluid cells and apply the Riemann Solver (2.14) on the mixed cells to calculate the point values, Section 2, Section 3.3.
•

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$ , compute the one-sided local speeds of propagation using (2.8), Section 2.
•

At time $t^{n}$ , calculate the reference time step $\Delta t$ using (4.2), Section 4.3.
•

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$ , compute the local time step for each cell level, (4.4), Section 4.3.
•

At each time level $t^{n,p}_{l},p=0,1,...,\mathcal{P}_{l}-1$ , compute numerical fluxes and source term in the adaptive central-upwind scheme (4.5a)-(4.5b), (2.6), (2.25)-(2.26), Section 2, Section 4.3.

Step 2. On mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ , compute WLR error using (4.9) in Section 4.4 and determine the refinement/de-refinement status for each cell, Section 4.4.

Step 3. Generate the new adaptive mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ at $t^{n+1}$ , Section 4.2. This step includes coarsening of some cells, refinement of some cells, and the appropriate projection of the cell averages from the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ at $t^{n}$ onto a new adaptive mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ at time $t^{n+1}$ , Section 4.2.

Step 4. Repeat Step 1 - Step 3 until final time.

4.2 Adaptive Mesh Refinement/Coarsening

The purpose of mesh reconstruction is to obtain adaptive grids which delivers higher accuracy with a lower computational cost. An efficient adaptive mesh should have small cells in the regions with large errors and large cells in the other regions. In particular, at time $t^{n}$ , we start with the given mesh, denoted as $\mathcal{T}^{n,m}=\{T^{n,m}_{j}\}$ , where $T^{n,m}_{j}$ is a triangular cell with the barycenter $(x^{n,m}_{j},y^{n,m}_{j})$ , and index $m=0,1,2...$ is the level of refinement ( $\mathcal{T}^{n,0}\equiv\mathcal{T}^{0,0}$ for all $n$ and $\mathcal{T}^{0,0}$ represents the initial mesh with no refinement). The triangular cells in the mesh $\mathcal{T}^{n,m}$ are flagged for refinement/de-refinement based on the weak local residual (WLR) error estimator, see Section 4.4. On grid $\mathcal{T}^{n,m}$ , we apply “regular refinement” described in [13] on the triangles flagged for refinement to obtain a new mesh $\mathcal{T}^{n,m+1}$ with the refinement level $m+1$ . Namely, each flagged triangle (“parent” triangle) is split into four smaller triangles (“children” triangles) by inserting a new node at the mid-point of each edge of the “parent” triangle. Fig. 4.1 (a) is an illustration of the “regular refinement”, where a flagged cell $T^{n,m}_{j}$ is refined to obtain the “children” cells $T^{n,m+1}_{j_{s}},s=1,2,3,4$ by using the mid-points of the sides. In addition, due to the insertion of new nodes on the edges of the non-flagged triangles adjacent to refined triangles, we must also refine these neighboring cells by creating a new edge between the hanging node and the opposite corner as illustrated in Fig. 4.1 (b).

In practice, there may be some cells having very large WLR error (4.9), and we need to reach a higher level of refinement for those cells to improve the accuracy. This can be done by repeating the refinement for the flagged triangles in the refined mesh $\mathcal{T}^{n,m+1}$ to get the mesh with higher level $\mathcal{T}^{n,m+2},m=0,1,2,....$ . Fig. 4.2 is the illustration of the “regular refinement” procedure with two levels of refinement.

Note that in the numerical simulations of the wave phenomena, the regions of the domain flagged for refinement changes over time evolution and the refinement in some cells may become no longer needed. Hence, in [13], the de-refinement/coarsening procedure is performed to deactivate unnecessarily fine cells in the grid. Namely, at time $t^{n}$ , we deactivate “children” cells in the mesh $\mathcal{T}^{n,m+1},m=0,1,...,M_{n}-1$ based on the WLR and activate the corresponding “parent” cell from the mesh $\mathcal{T}^{n,m}$ back. More details of refinement/de-refinement process can be seen in [13].

Finally, at time $t^{n}$ , a hierarchical system of grids $\mathcal{S}^{n}=\{\mathcal{T}^{n,0},\mathcal{T}^{n,1},\mathcal{T}^{n,2},...,\mathcal{T}^{n,\mathcal{M}_{n}}\}$ is obtained, where $\mathcal{T}^{n,m}$ , $m=1,2,...,\mathcal{M}_{n}$ is the grid with the level of refinement $m$ reconstructed from the grid $\mathcal{T}^{n,m-1}$ . We will use the final mesh $\mathcal{T}^{n,\mathcal{M}_{n}}\in\mathcal{S}^{n}$ for adaptive central-upwind scheme (4.5a)-(4.5b) to evolve the numerical solution from time $t^{n}$ to time $t^{n+1}$ . Next, at $t^{n+1}$ , we generate a new adaptive grid $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}\in\mathcal{S}^{n+1}$ from the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ . After the mesh reconstruction at time $t^{n+1}$ , the obtained cell averages $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}^{n+1}$ on the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ need to be accurately projected on the new mesh $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ , using the ideas as summarized briefly below.

Case 1. At $t^{n+1}$ , a triangle $T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is the same cell as in the grid $\mathcal{T}^{n,\mathcal{M}_{n}}$ , we will maintain the cell averages for that triangle at $t^{n+1}$ .

Case 2. A cell $T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is obtained by de-refining some finer cells $T^{n,\mathcal{M}_{n}}_{j_{s}}\in\mathcal{T}^{n,\mathcal{M}_{n}},s=1,2,..,S$ . The cell average of solution, $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}$ in the cell $T^{n+1,\mathcal{M}_{n+1}}_{j}$ , is computed as

\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}

\displaystyle=\dfrac{1}{|T^{n+1,\mathcal{M}_{n+1}}_{j}|}\sum\limits_{s=1}^{S}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n}_{j_{s}}|T^{n,\mathcal{M}_{n}}_{j_{s}}|,

where $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n}_{j_{s}}$ is the solution in $T^{n,\mathcal{M}_{n}}_{j_{s}}$ .

Case 3. A triangle $T^{n+1,\mathcal{M}_{n+1}}_{j}\in\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ is obtained from the refinement of the cell $T^{n,\mathcal{M}_{n}}_{i}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ . If $T^{n,\mathcal{M}_{n}}_{i}$ is a single-fluid cell, the cell averages at $t^{n+1}$ in $\mathcal{T}^{n+1,\mathcal{M}_{n+1}}$ are approximated by using the the piecewise linear reconstruction (2.7) of the solution at $t^{n+1}$ in the triangle $T^{n,\mathcal{M}_{n}}_{i}$ . If $T^{n,\mathcal{M}_{n}}_{i}$ contains the density discontinuity, we will compute the cells averages in the “son” cell $T^{n+1,\mathcal{M}_{n+1}}_{j}$ based on the interface tracking, see Section 3.2, and the information of the nearby single-fluid cells. For example, suppose that from the reconstructed interface in “dad” cell $T^{n,\mathcal{M}_{n}}_{i}$ , the “child” cell $T^{n+1,\mathcal{M}_{n+1}}_{j}$ completely lies in one fluid, called fluid 1. Hence, triangle $T^{n+1,\mathcal{M}_{n+1}}_{j}$ is a single-fluid cell and the cell averages $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$U$\kern-0.50003pt}}}^{n+1}_{j}$ are set to be equal to the cell averages of another “reliable” cell which is in fluid 1 and is the closest cell to the “dad” cell $T^{n,\mathcal{M}_{n}}_{i}$ .

Remark 4.1

The update of cell averages in case 3 does not ensure the mass conservation which means the total mass of “son” cells in the adaptive grid is not equal to the mass of their “dad” cell. However, the volume of fluid and the exact location of the interface are conserved. In addition, the reconstructing method also maintains the steady state solution for “lake at rest” situations, see example 2 Section 5, since the values of the cell averages are obtained by using the information from nearby “reliable” cells.

4.3 Second-order Adaptive Time Evolution

Note that using a global time step in the adaptive algorithm may lead to a very small time step due to the presence of much finer cells in the mesh. To reduce the computational cost, we consider the approach based on the adaptive time step from [13, 11, 12, 30]. The main idea of this approach is to group cells into different levels based on the cell sizes and applying local time step on each level to evolve from $t^{n}$ to $t^{n+1}$ . Recently, in our work on the shallow water model [13], we have derived a simple and efficient adaptive time evolution algorithm based on the second-order strong stability preserving Runge-Kutta methods (SSPRK2) in [16, 11, 12, 30]. This method is also capable to perform on the multi-fluid system (1.1a)-(1.1b). The algorithm can be briefly described by an example as illustrated in Fig. 4.3.

First, we group all cells in the grid $\mathcal{T}^{n,\mathcal{M}_{n}}$ at time $t^{n}$ in cell levels $l=0,1,..,L$ based on their sizes. Namely, a cell $T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ belongs to the level $l$ , if $l$ is the smallest positive integer satisfying,

2^{l}\geq\dfrac{\displaystyle\max_{j}\left(\displaystyle\min_{k}\left(r_{jk}\right)\right)}{\displaystyle\min_{k}(r_{j,k})},

(4.1)

where $r_{jk},k=1,2,3$ are three altitudes of triangle $T^{n,\mathcal{M}_{n}}_{j}$ . Next, at $t^{n}$ , we define the reference time step $\Delta t$ as the local time step on the coarsest level $l=0$ of cells in the mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ by considering the CFL-type condition (2.20) locally on level $l=0$ .

\displaystyle\Delta t\equiv\Delta t^{n,0}_{0}=\frac{0.9\displaystyle\max_{j}\left(\displaystyle\min_{k}\left(r_{jk}\right)\right)}{18a_{max}},

(4.2)

where

a_{max}:=\displaystyle\max_{j,k}(a^{\rm in}_{jk},a^{\rm out}_{jk}),

(4.3)

and $(a^{\rm in}_{jk},a^{\rm out}_{jk})$ are the local one-sided speeds of propagation (2.8) at $t^{n}$ for sides $k=1,2,3$ in the triangle $T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ . We set, $t^{n+1}=t^{n}+\Delta t$ .

Next, assume that $\mathcal{P}_{l}$ is the number of steps needed for the higher levels $l=1,..L$ to evolve from $t^{n}$ to $t^{n+1},$ namely $[t^{n},t^{n+1}]=\cup[t_{l}^{n,p},t_{l}^{n,p+1}],p=0,...,\mathcal{P}_{l}-1$ with $t^{n,0}_{l}\equiv t^{n}$ , $t_{l}^{n,\mathcal{P}_{l}}\equiv t^{n+1}\quad\forall l$ . At $t^{n,p}_{l}$ , the local time step for cells on these levels $l=1,...,L$ is calculated by

\Delta t^{n,p}_{l}=\dfrac{2^{-l}\Delta t}{\max(\mu^{n,p}_{l},1)},

(4.4)

where parameter $\mu^{n,p}_{l}$ takes into account change in the local one-sided speeds of the propagation,

\mu^{n,p}_{l}=\dfrac{\displaystyle\max_{j,k}(a^{\rm in}_{jk},a^{\rm out}_{jk})^{n,p}_{l}}{a_{max}},

where $(a^{\rm in}_{jk},a^{\rm out}_{jk})^{n,p}_{l}$ are the local one-sided speeds of propagation of the cell $T^{n,\mathcal{M}_{n}}_{j}$ in the level $l$ at $t^{n,p}_{l}$ . Therefore, on each cell $T^{n,\mathcal{M}_{n}}_{j}$ of level $l$ , for each substep $[t_{l}^{n,p},t_{l}^{n,p+1}]\equiv[t_{l}^{n,p},t_{l}^{n,p}+\Delta t_{l}^{n,p}],p=0,1,2,3,...,\mathcal{P}_{l}-1$ of the evolution from $t^{n}$ to $t^{n+1}$ , we apply the following two adaptive steps of the SSPRK2 method, see [16, 11, 12, 30],


$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{(1)}$	$\displaystyle=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p}-\Delta t^{n,p}_{l}\left(\frac{1}{\|T^{n,\mathcal{M}_{n}}_{j}\|}\sum_{k=1}^{3}{\bm{H}}^{n,p}_{jk}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{S}$\kern-0.50003pt}}}^{n,p}_{j}\right):={\bm{R}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p},\Delta t_{l}^{n,p}),$	(4.5a)
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p+1}$	$\displaystyle=\dfrac{1}{2}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,p}+\dfrac{1}{2}{\bm{R}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{(1)},\Delta t_{l}^{n,p}).$	(4.5b)

The flux term ${\bm{H}}^{n,p}_{jk}$ in (4.5a) -(4.5b) is the flux (2.6) computed at $t=t^{n,p}_{l}$ . The source term $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt${\bm{S}}$\kern-0.50003pt}}}^{n,p}_{j}$ in (4.5a) -(4.5b) is the source (2.25) - (2.26) computed at $t=t^{n,p}_{l}$ with the time step $\Delta t_{l}^{n,p}$ . Note that, $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,0}\equiv\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n,\mathcal{P}_{l}}\equiv\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern-0.50003pt}}}_{j}^{n+1}$ .

Remark 4.2

If cells from different cell levels are neighbors, we use linear interpolation in time to match the time levels of such cells, see also Fig. 3.6, for the illustration of the interpolation.

4.4 A Posteriori Error Estimator

To create a robust indicator for the adaptive mesh refinement, in our prior work [13], we have derived local error estimator from the idea of Weak Local Residual (WLR) presented in [40, 18]. This error indicator has shown its advantages in accurately capturing the regions with large error in numerical simulation for Saint-Venant system of shallow water model. Hence, for the adaptive central-upwind scheme for SWEDs, we will extend the error estimator by applying the computation performed in [13] for the last equations in the system (1.1d).

Let us recall that from the weak form of the mass conservation equation (1.1a), in [13], the WLR errors $E^{w,n+\frac{1}{2}}_{i}$ at each node $N_{i}$ on mesh $\mathcal{T}^{n,\mathcal{M}_{n}}$ is given by the formula,

$\displaystyle E^{w,n+\frac{1}{2}}_{i}$	$\displaystyle=\frac{1}{\Delta}(\mathcal{U}^{w,n+\frac{1}{2}}_{i}+\mathcal{F}^{w,n+\frac{1}{2}}_{i}+\mathcal{G}^{w,n+\frac{1}{2}}_{i}),$	(4.6)
$\displaystyle\mathcal{\bm{U}}^{w,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}\frac{1}{3}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}^{n}_{j_{c}}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$w$\kern-0.50003pt}}}^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{F}^{w,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}a^{(i)}_{c}\frac{\Delta t}{2}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|((\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hu$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hu$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{G}^{w,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}b^{(i)}_{c}\frac{\Delta t}{2}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|((\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hv$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hv$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$

where $C_{i}$ is the number of triangles $T^{n,\mathcal{M}_{n}}_{j_{c}}$ having common vertex at node $N_{i}$ . Here, the quantity $(a^{(i)}_{c},b^{(i)}_{c})$ is the gradient of the linear piece restricted to $T^{n,\mathcal{M}_{n}}_{j_{c}}$ ,

	$\displaystyle a^{(i)}_{c}$	$\displaystyle=\dfrac{\widetilde{y}_{2}-\widetilde{y}_{3}}{(\widetilde{y}_{3}-\widetilde{y}_{i})(\widetilde{x}_{2}-\widetilde{x}_{i})-(\widetilde{y}_{2}-\widetilde{y}_{i})(\widetilde{x}_{3}-\widetilde{x}_{i})},$		(4.7)
	$\displaystyle b^{(i)}_{c}$	$\displaystyle=\dfrac{\widetilde{x}_{3}-\widetilde{x}_{2}}{(\widetilde{y}_{3}-\widetilde{y}_{i})(\widetilde{x}_{2}-\widetilde{x}_{i})-(\widetilde{y}_{2}-\widetilde{y}_{i})(\widetilde{x}_{3}-\widetilde{x}_{i})},$		(4.7)

where $N_{i}=(\widetilde{x}_{i},\widetilde{y}_{i}),(\widetilde{x}_{2},\widetilde{y}_{2}),$ and $(\widetilde{x}_{3},\widetilde{y}_{3})$ are the three vertices of triangle $T^{n,\mathcal{M}_{n}}_{j_{c}}$ .

Now, by applying the same calculation in [13] on the weak form of the last equations in the system (1.1d), we then define the WLR error of variable $h\rho$ at node $N_{i}$ in the grid as,

$\displaystyle E^{h\rho,n+\frac{1}{2}}_{i}$	$\displaystyle=\frac{1}{\Delta}(\mathcal{U}^{h\rho,n+\frac{1}{2}}_{i}+\mathcal{F}^{h\rho,n+\frac{1}{2}}_{i}+\mathcal{G}^{h\rho,n+\frac{1}{2}}_{i}),$	(4.8)
$\displaystyle\mathcal{\bm{U}}^{h\rho,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}\frac{1}{3}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{n}_{j_{c}}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}}^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{F}^{h\rho,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}a^{(i)}_{c}\frac{\Delta t}{2}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|((\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hu\rho$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hu\rho$\kern-0.50003pt}}})^{n+1}_{j_{c}}),$
$\displaystyle\mathcal{G}^{h\rho,n+\frac{1}{2}}_{i}$	$\displaystyle=\sum\limits_{c=1}^{C_{i}}b^{(i)}_{c}\frac{\Delta t}{2}\|T^{n,\mathcal{M}_{n}}_{j_{c}}\|((\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hv\rho$\kern-0.50003pt}}})^{n}_{j_{c}}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hv\rho$\kern-0.50003pt}}})^{n+1}_{j_{c}}).$

Hence, the error in a cell $T^{n,\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n.\mathcal{M}_{n}}$ takes into account both WLR errors of water surface $w$ and of variable $h\rho$ as,

e_{j}=\max_{k}\left(\left|E^{w,n+\frac{1}{2}}_{jk}\right|,\left|E^{h\rho,n+\frac{1}{2}}_{jk}\right|\right),\quad k=1,2,3,

(4.9)

where $E^{w,n+\frac{1}{2}}_{jk}$ and $E^{h\rho,n+\frac{1}{2}}_{jk}$ are the WLR errors computed in (4.6) and (4.8) at a node $k$ of triangle $T_{j}$ .

The error $e_{j}$ in each cell $T^{n.\mathcal{M}_{n}}_{j}\in\mathcal{T}^{n,\mathcal{M}_{n}}$ is compared to the error tolerance computed by

\omega=\sigma\max_{j}(e_{j}),

(4.10)

where $\sigma<1$ is a given problem-dependent constant. Based on the error comparison, the cell is then either “flagged” for refinement/de-refinement or “no-change”.

5 Numerical Experiments

In this section, we will verify the computational efficiency of the designed adaptive central-upwind scheme. We compare the results of the adaptive method with the results of the central-upwind scheme without the adaptivity, see Section 4, on uniform triangular meshes (example of such uniform triangular mesh is outlined in Fig. 5.1). To this end, in all examples, we calculate the $L^{1}$ -errors and a ratio, $\mathcal{R}_{CPU}=\frac{CPU_{uniform}}{CPU_{adaptive}}$ , which is the ratio of the CPU times of the central-upwind algorithm without the mesh reconstruction to the CPU time of the adaptive central-upwind algorithm. To compare $L^{1}$ -errors, as well as to compare the CPU times and to compute $\mathcal{R}_{CPU}$ , we consider uniform mesh and the adaptive mesh with the same size of the smallest cells. Namely, in Tables 5.1-5.8, $L^{1}$ -errors and $\mathcal{R}_{CPU}$ are computed by using the uniform meshes $2\times N\times N,N=100,200,400$ and the corresponding adaptive meshes which are reconstructed from the coarser uniform mesh $2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}$ ( $\mathcal{M}=1,2$ is the highest level of refinement in the adaptive mesh). To compute the $L^{1}$ -errors, the reference solution is obtained by applying the central-upwind method without implementing adaptivity techniques on the uniform mesh with $2\times 800\times 800$ triangles. In all experiments, we consider a zero-order extrapolation at all boundaries. In addition, we use the gravitational acceleration, $g=1.0$ and the reference density, $\rho_{0}=997$ [15] for all examples. The desingularization parameters $\tau$ and $\varepsilon$ for calculations of the velocity components $u$ and $v$ are set $\tau=\max_{j}\{|T^{n,\mathcal{M}_{n}}_{j}|^{2}\}$ and $\varepsilon=10^{-4}$ (see Section 2.1 formula (2.6) in [29]).

5.1 Example 1:

In the first example taken from [9], we will compare the performance of the adaptive central-upwind scheme and the central-upwind scheme without adaptivity on uniform triangular meshes. We also verify experimentally the advantages of the interface reconstruction in compressing the diffusion of variable density at the interface of fluids and improving the accuracy of computed solutions.

We consider the bottom topography $B(x,y,0):=0$ and the following initial condition,

(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(2,0,0,1.5\rho_{0}),\quad\mbox{if}\quad&x^{2}+y^{2}<0.5,\\ (1,0,0,\rho_{0}),\quad&\mbox{otherwise}.\end{cases}

(5.1)

The data is simulated in the domain $[-1,1]\times[-1,1]$ . The error tolerance (4.10) for the mesh refinement in this example is set to $\omega=0.01\max_{j}(e_{j})$ .

In Fig. 5.2 and Fig. 5.3, we show the numerical solution of the water surface (first column) and the density (second column) at time $t=0.15$ . The solution is calculated by using the central-upwind scheme on uniform meshes on Fig. 5.2 (a, b) and Fig. 5.3 (a,b) and using the adaptive central-upwind scheme on Fig. 5.2 (c, d) and Fig. 5.3 (c,d). The adaptive meshes in Fig. 5.2 (third column) are obtained from the uniform mesh $2\times 100\times 100$ , Fig. 5.2 (a). The adaptive mesh with one level of refinement $\mathcal{M}=1$ (as the highest level of refinement) is on Fig. 5.2 (c), and with two levels of refinement $\mathcal{M}=2$ (as the highest level of refinement) is on Fig. 5.2 (d). As can be observed in Fig. 5.2, both $w$ and $\rho$ are much sharper resolved by using the adaptive central-upwind scheme. We note also, that by increasing the level of refinement from $\mathcal{M}=1$ to $\mathcal{M}=2$ , the number of cells in the mesh increases from $40,292$ cells with $\mathcal{M}=1$ to $60,326$ cells with $\mathcal{M}=2$ , but the accuracy is clearly improved with higher resolution as seen in Fig. 5.2 (c) and Fig. 5.2 (d). Also from the adaptive meshes in Fig. 5.2 (third column), one can easily see that only cells in the region having steep gradients are refined. This means that the WLR error estimator accurately detects regions in the domain for adaptive refinement/coarsening.

Next, in Table 5.1 we calculate the $L^{1}$ -errors obtained on the adaptive grids and on the fixed uniform grids. The errors obtained in the uniform meshes are approximate to the errors calculated in the corresponding adaptive meshes (the adaptive meshes have the same size of the smallest cells with the uniform meshes). However, the adaptive scheme uses fewer cells than the central-upwind scheme which does not consider the mesh reconstruction. In addition, in Table 5.2, we also compute the $\mathcal{R}_{CPU}$ ratio to compare the computational cost of the two methods. The results in Table 5.1 and Table 5.2 show that the adaptive scheme produces similar accuracy as the scheme designed on fixed uniform triangular meshes, but at a less computational cost.

algorithm without adaptivity			adaptive algorithm
			one level $\mathcal{M}=1$			two levels $\mathcal{M}=2$
cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate
$2\times 100\times 100$	0.0256		13,516	0.0257		12,800	0.0265
$2\times 200\times 200$	0.0127	1.01	40,292	0.0128	1.00	38,284	0.0133	0.99
$2\times 400\times 400$	0.0047	1.43	153,152	0.0049	1.39	60,326	0.0050	1.41

Table 5.1: Example 1:

L^{1}

-errors of the water surface

w

of the IVP (5.1) at

t=0.15

, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh

2\times N\times N,N=100,200,400

) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh

2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}

uniform mesh	adaptive mesh $\mathcal{M}=1$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=1$		adaptive mesh $\mathcal{M}=2$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
(cells)	(cells)	total	without grid generation	(cells)	total	without grid generation
$2\times 100\times 100$	13,516	1.81	2.01	12,800	1.83	1.94
$2\times 200\times 200$	40,292	1.82	2.14	38,284	2.14	2.30
$2\times 400\times 400$	153,152	1.98	2.29	60,326	2.33	2.50
$\mathcal{R}_{CPU}$ average:		1.87	2.15		2.10	2.27

Table 5.2: Example 1:

\mathcal{R}_{CPU}

ratio for the IVP (5.1) at

t=0.15

, where for adaptive central-upwind scheme, we consider the total CPU times and CPU times without the grid generation.

In addition, we will use this example to show that interface reconstruction presented in Section 3.2 plays an important part in preserving the sharpness of the solution as well as improving the accuracy of the adaptive central-upwind scheme. On Fig. 5.4 (a), we plot the water surface $w$ (first) and density $\rho$ (second column) at $t=0.15$ by using the adaptive central-upwind method, but implemented without the interface tracking technique. We then compare the results shown in Fig. 5.4 (a) to the results calculated by using the adaptive scheme with the interface tracking, Fig. 5.4 (b). The adaptive meshes on Fig. 5.4 (third column) are reconstructed from the uniform mesh $2\times 100\times 100$ with one level of refinement. As can be seen in Fig. 5.4 (a), both $w$ and $\rho$ are very scattered around the contact wave when we do not track and reconstruct the interface. Meanwhile, in Fig. 5.4 (b), the proposed adaptive scheme, though using an adaptive mesh with fewer cells (10828 cells fewer), provides more accurate results.

In the next numerical test, we replace the flat bottom with the bottom topography that consists of two Gaussian shaped humps as

B(x,y,t)=\begin{cases}0.5e^{-100(x+0.5)^{2}+(y+0.5)^{2}},\quad&\mbox{if}\quad x<0,\\ 0.6e^{-100(x-0.5)^{2}+(y-0.5)^{2}},\quad&\mbox{if}\quad x\geq 0.\end{cases}

(5.2)

The purpose of this test is to illustrate the performance of the adaptive algorithm in situations having irregular bottom topography. In Fig. 5.5, we show the contour plots of the water surface $w$ (first column) and density $\rho$ (second column) obtained at $t=0.2$ by using the central-upwind scheme with and without adaptivity. The computed solutions of the water surface exhibit reflecting waves where the flow meets the submerged humps. Clearly, from the plots of the adaptive meshes in Fig. 5.5 (third column), the meshes are adapted to the behavior of the flow. Hence, the WLR error estimator is capable to exactly detect the location of the steep local gradients in the solution.

We then recompute the accuracy of the solution for this example 5.1-5.2, see Table 5.3, and the CPU time ratio, see Table 5.4. The results show that the adaptive scheme uses fewer cells and takes a smaller CPU time to achieve the approximately small $L^{1}$ -errors as computed by the scheme without adaptivity. Therefore, the advantages of the adaptive scheme is maintained in examples with irregular bottom level.

algorithm without adaptivity			adaptive algorithm
			one level $\mathcal{M}=1$			two levels $\mathcal{M}=2$
cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate
$2\times 100\times 100$	0.0271		16,252	0.0273		16,006	0.0282
$2\times 200\times 200$	0.0134	1.01	57,300	0.0136	1.00	42,602	0.0140	0.99
$2\times 400\times 400$	0.0050	1.43	213,268	0.0054	1.40	120,654	0.0055	1.41

Table 5.3: Example 1:

L^{1}

-errors of the water surface

w

of the IVP (5.1)-(5.2) at

t=0.2

, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh

2\times N\times N,N=100,200,400

) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh

2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}

uniform mesh	adaptive mesh $\mathcal{M}=1$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=1$		adaptive mesh $\mathcal{M}=2$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
(cells)	(cells)	total	without grid generation	(cells)	total	without grid generation
$2\times 100\times 100$	16,252	1.69	1.89	16,006	1.84	1.94
$2\times 200\times 200$	57,300	2.57	2.95	42,602	3.41	3.68
$2\times 400\times 400$	213,268	2.24	2.60	120,654	2.63	2.84
$\mathcal{R}_{CPU}$ average:		2.17	2.48		2.63	2.82

Table 5.4: Example 1:

\mathcal{R}_{CPU}

ratio for the IVP (5.1)-(5.2) at

t=0.2

, where for adaptive central-upwind scheme, we consider the total CPU times and CPU times without the grid generation.

5.2 Example 2:

The second numerical example here was proposed in [9] to verify the capability of the adaptive scheme in preserving the steady state solution in “lake at rest’ problems, (1.2) and (1.3). In particular, the initial data consists of two “lake at rest” states of type (1.2) connected through the density jump corresponding to the “lake at rest” state of type (1.3) as

(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(3,0,0,\dfrac{4}{3}\rho_{0}),\quad\mbox{if}\quad&x^{2}+y^{2}<0.25,\\ (2,0,0,3\rho_{0}),\quad&\mbox{otherwise}.\end{cases}

(5.3)

In this example, we consider the bottom topography (5.2) on a computational domain $[-1,1]\times[-1,1]$ . To reconstruct the adaptive meshes, the threshold is set, $\omega=0.1\max_{j}(e_{j})$ . In Fig. 5.6, we present the plots of the computed water surface (first column) and density (second column) at $t=0.15$ obtained by using the central-upwind scheme, but without adaptivity, Fig. 5.6 (a, b) and by using the adaptive algorithm Fig. 5.6 (c, d). The adaptive grids plotted on Fig. 5.6 (third column) are generated from the uniform mesh $2\times 100\times 100$ with one level of refinement $\mathcal{M}=1$ , Fig. 5.6 (c), and two levels of refinement $\mathcal{M}=2$ , Fig. 5.6 (d). Fig. 5.7 shows the 3D plots of the numerical solution computed by the two methods. As expected, in Fig. 5.6 and Fig. 5.7, the adaptive scheme with interface tracking exactly preserves the steady state. Hence, in Fig. 5.6 (third column), the WLR error only marks cells surrounding the circle of density jump for refinement. In addition, no pressure oscillations are observed at the interface.

Next, we will illustrate the advantages of the adaptive central-upwind scheme. We compute the $L^{1}$ -error, Table 5.5, and the CPU ratio, Table 5.6, by using the central-upwind method without adaptivity and using the adaptive algorithm presented in our work. In order to compare the computational costs and calculate $\mathcal{R}_{CPU}$ , we consider uniform and adaptive meshes with the same size of the smallest cells. Table 5.5 shows that in the adaptive meshes, we achieve $L^{1}$ -errors as small as the errors obtained in the corresponding uniform meshes. However, the adaptive algorithm uses fewer cells and reduces the CPU times up to eight times. The computational cost is remarkably cut down since as illustrated in Fig. 5.6 and Fig. 5.7, only a few cells in the neighborhood of the density discontinuity have large WLR errors and are therefore marked for refinement. This example has clearly show the efficiency of the proposed scheme for numerically solving the system of multi-fluid flow.

algorithm without adaptivity			adaptive algorithm
			one level $\mathcal{M}=1$			two levels $\mathcal{M}=2$
cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate
$2\times 100\times 100$	0.0045		6,284	0.0045		5,928	0.0045
$2\times 200\times 200$	0.0021	1.10	22,546	0.0021	1.10	11,260	0.0021	1.10
$2\times 400\times 400$	7.9252e-04	1.41	85,132	8.5079e-04	1.30	33,904	8.9818e-04	1.23

Table 5.5: Example 2:

L^{1}

-errors of the water surface

w

t=0.15

, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh

2\times N\times N,N=100,200,400

) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh

2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}

uniform mesh	adaptive mesh $\mathcal{M}=1$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=1$		adaptive mesh $\mathcal{M}=2$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
(cells)	(cells)	total	without grid generation	(cells)	total	without grid generation
$2\times 100\times 100$	6,284	3.01	3.54	5,928	3.48	3.75
$2\times 200\times 200$	22,546	3.78	4.76	11,260	7.76	8.66
$2\times 400\times 400$	85,132	4.46	5.57	33,904	11.13	12.86
$\mathcal{R}_{CPU}$ average:		3.75	4.62		7.46	8.42

Table 5.6: Example 2:

\mathcal{R}_{CPU}

ratio at

t=0.15

, where for adaptive central-upwind scheme, we consider the total CPU times and CPU times without the grid generation.

5.3 Example 3:

The last example is designed to illustrate the capability of the proposed adaptive algorithm to handle irregular density interfaces. Hence, in a domain $[-1,1]\times[-1,1]$ , the density jump at $t=0$ is given by a curve which consists of a horizontal segment, a vertical segment, and a quarter of a circle connected at their endpoints as illustrated in Fig. 5.8.

The initial condition is described by

(w,u,v,\rho)^{T}(x,y,0)=\begin{cases}(2,0,0,\rho_{0}),\quad\mbox{if}\quad&(x,y)\in\Omega,\\ (1,0,0,1.5\rho_{0}),\quad&\mbox{otherwise}.\end{cases}

(5.4)

where

\Omega:=\{x<-0.5,y<0\}\cup\{(x+0.5)^{2}+(y+0.5)^{2}<0.25\}\cup\{x<0,y<-0.5\}.

. We consider a bottom topography with a surbmerged hump as

B(x,y,t)=0.5e^{-100(x^{2}+y^{2})}.

In this example, we will also perform the same numerical tests which are done in previous examples. Namely, we first calculate the water surface and density at $t=0.15$ using central-upwind scheme without adaptivity and present the results in Fig. 5.9 (a, b) and Fig. 5.10 (a, b). In 5.9 (c, d) and 5.10 (c, d), we plot the results for $w$ (first column) and $\rho$ (second column) obtained by the adaptive scheme. The adaptive grids in 5.9 (third column) are generated from the uniform grid $2\times 100\times 100$ for one level of refinement $\mathcal{M}=1$ and $\mathcal{M}=2$ using the threshold $\omega=0.01\max_{j}(e_{j})$ . As expected, the solutions obtained in the adaptive meshes with high levels of refinement are much sharper than the solutions computed in fixed uniform meshes. The density jump moves Northeast and does not diffuse. There is no non-physical spurious waves generated at the interface. Also, as can be seen in the adaptive meshes, the WLR error indicator captures subtle features of the solution.

Finally, we compute the $L^{1}$ -errors in Table 5.7 and the CPU ratio in Table 5.8 by using the adaptive scheme and using the central-upwind method without the mesh generation. By comparing the results presented in Tables 5.7 and 5.8, one can easily see that at a reduced computational cost, the proposed adaptive central-upwind method is still able to obtain the accurate solutions for this example. In our experiments, we considered only $\mathcal{M}=1$ and $\mathcal{M}=2$ , but to further enhance the accuracy of the numerical solution at the lower computational cost, one can consider higher levels of refinement.

algorithm without adaptivity			adaptive algorithm
			one level $\mathcal{M}=1$			two levels $\mathcal{M}=2$
cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate	cells	$L^{1}$ -error	rate
$2\times 100\times 100$	0.0155		12,164	0.0145		8,994	0.0166
$2\times 200\times 200$	0.0074	1.07	40,814	0.0074	0.97	30,263	0.0071	1.23
$2\times 400\times 400$	0.0027	1.45	148,473	0.0028	1.40	87,214	0.0028	1.34

Table 5.7: Example 3:

L^{1}

-errors of the water surface

w

t=0.15

, and the convergence rates of the central-upwind scheme without adaptivity (uniform mesh

2\times N\times N,N=100,200,400

) and the adaptive scheme (the corresponding adaptive mesh is reconstructed from the uniform mesh

2\times N/2^{\mathcal{M}}\times N/2^{\mathcal{M}}

uniform mesh	adaptive mesh $\mathcal{M}=1$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=1$		adaptive mesh $\mathcal{M}=2$	$\mathcal{R}_{CPU}$ with $\mathcal{M}=2$
(cells)	(cells)	total	without grid generation	(cells)	total	without grid generation
$2\times 100\times 100$	12,164	1.75	1.98	8,994	2.69	2.87
$2\times 200\times 200$	40,814	2.71	3.14	30,263	3.57	3.86
$2\times 400\times 400$	148,473	2.86	3.38	87,214	4.60	5.03
$\mathcal{R}_{CPU}$ average:		2.44	2.83		3.62	3.92

Table 5.8: Example 3:

\mathcal{R}_{CPU}

ratio at

t=0.15

, where for adaptive central-upwind scheme, we consider the total CPU times and CPU times without the grid generation.

6 Conclusion

We have developed a new adaptive well-balanced and positivity preserving central-upwind scheme on unstructured traingular meshes for shallow water equations with variable density. The scheme is designed as an extension the scheme in [13] by utilizing the interface tracking method in [9] and the interface reconstruction in [15]. The proposed scheme is capable to preserve the steady state solutions (1.2) and (1.3) and prevent the oscillation at the density jumps. In addition, to achieve an efficient strategy for the adaptive mesh reconstruction, we also obtain a robust local error indicator. We performed several challenging numerical tests for multi-fluid models and we demonstrated that the new adaptive central-upwind scheme maintains well-balanced and positivity-preserving properties and obtains high-accuracy at a reduced computational cost.

Acknowledgements

I would like to sincerely thank my advisor, Prof. Yekaterina Epshteyn, for her suggestion of the problem to investigate, her encouragement and help. The proposed approach in this paper is an extension of the adaptive method developed in [13] which is my joint work with her. Without her support, this work would not have started, progressed, or ended.

References

[1] D. S. Balsara, M. Dumbser, and R. Abgrall, Multidimensional hllc riemann solver for unstructured meshes–with application to euler and mhd flows, Journal of Computational Physics, 261 (2014), pp. 172–208.
[2] A. Bollermann, G. Chen, A. Kurganov, and S. Noelle, A well-balanced reconstruction of wet/dry fronts for the shallow water equations, J. Sci. Comput., 56 (2013), pp. 267–290.
[3] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova, Central Upwind Scheme on Triangular Grids for the Saint Venant System of Shallow Water Equations, AIP Conference Proceedings, 1389 (2011), pp. 686–689.
[4] S. Bryson, Y. Epshteyn, A. Kurganov, and G. Petrova, Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system, ESAIM Math. Model. Numer. Anal., 45 (2011), pp. 423–446.
[5] S. Bryson and D. Levy, Balanced central schemes for the shallow water equations on unstructured grids, SIAM J. Sci. Comput., 27 (2005), pp. 532–552 (electronic).
[6] A. Chertock, Y. Epshteyn, H. Hu, and A. Kurganov, High-order positivity-preserving hybrid finite-volume-finite-difference methods for chemotaxis systems, Adv. Comput. Math., 44 (2018), pp. 327–350.
[7] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, and P. A. Markowich, Sinking, merging and stationary plumes in a coupled chemotaxis-fluid model: a high-resolution numerical approach, Journal of Fluid Mechanics, 694 (2012), pp. 155–190.
[8] A. Chertock, A. Kurganov, and Y. Liu, Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients, Numer. Math., 127 (2014), pp. 595–639.
[9] , Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients, Numerische Mathematik, 127 (2014), pp. 595–639.
[10] P. J. Dellar, Common hamiltonian structure of the shallow water equations with horizontal temperature gradients and magnetic fields, Physics of Fluids, 15 (2003), pp. 292–297.
[11] M. O. Domingues, S. M. Gomes, O. Roussel, and K. Schneider, An adaptive multiresolution scheme with local time stepping for evolutionary pdes, Journal of Computational Physics, 227 (2008), pp. 3758–3780.
[12] R. Donat, M. Martí, A. Martínez-Gavara, and P. Mulet, Well-balanced adaptive mesh refinement for shallow water flows, Journal of Computational Physics, 257 (2014), pp. 937–953.
[13] Y. Epshteyn and T. Nguyen, Adaptive central-upwind scheme on triangular grids for the saint-venant system, arXiv preprint arXiv:2011.06143, (2020).
[14] M. A. Ghazizadeh and A. Mohammadian, An adaptive central-upwind scheme on quadtree grids for variable density shallow water equations, arXiv preprint arXiv:2008.02111, (2020).
[15] M. A. Ghazizadeh, A. Mohammadian, and A. Kurganov, An adaptive well-balanced positivity preserving central-upwind scheme on quadtree grids for shallow water equations, Computers $\&$ Fluids, 208 (2020), p. 104633.
[16] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
[17] C. W. Hirt and B. D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, Journal of computational physics, 39 (1981), pp. 201–225.
[18] S. Karni and A. Kurganov, Local error analysis for approximate solutions of hyperbolic conservation laws, Adv. Comput. Math., 22 (2005), pp. 79–99.
[19] A. Kurganov, Finite-volume schemes for shallow-water equations, Acta Numerica, 27 (2018), pp. 289–351.
[20] A. Kurganov and D. Levy, Central-upwind schemes for the Saint-Venant system, M2AN Math. Model. Numer. Anal., 36 (2002), pp. 397–425.
[21] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys, 2 (2007), pp. 141–163.
[22] A. Kurganov and Y. Liu, New adaptive artificial viscosity method for hyperbolic systems of conservation laws, J. Comput. Phys., 231 (2012), pp. 8114–8132.
[23] A. Kurganov and J. Miller, Central-upwind scheme for Savage-Hutter type model of submarine landslides and generated tsunami waves, Comput. Methods Appl. Math., 14 (2014), pp. 177–201.
[24] A. Kurganov, S. Noelle, and G. Petrova, Semi-discrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740.
[25] A. Kurganov, S. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and hamilton–jacobi equations, SIAM Journal on Scientific Computing, 23 (2001), pp. 707–740.
[26] A. Kurganov and G. Petrova, Central-upwind schemes on triangular grids for hyperbolic systems of conservation laws, Numer. Methods Partial Differential Equations, 21 (2005), pp. 536–552.
[27] A. Kurganov and G. Petrova, A second-order well-balanced positivity preserving scheme for the Saint-Venant system, Commun. Math. Sci., 5 (2007), pp. 133–160.
[28] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations, J. Comput. Phys., 160 (2000), pp. 720–742.
[29] X. Liu, J. Albright, Y. Epshteyn, and A. Kurganov, Well-balanced positivity preserving central-upwind scheme with a novel wet/dry reconstruction on triangular grids for the Saint-Venant system, Journal of Computational Physics, 374 (2018), pp. 213 – 236.
[30] M. Moreira Lopes, M. O. Domingues, K. Schneider, and O. Mendes, Local time-stepping for adaptive multiresolution using natural extension of runge–kutta methods, Journal of Computational Physics, 382 (2019), pp. 291 – 318.
[31] S. Osher and R. P. Fedkiw, Level set methods: an overview and some recent results, Journal of Computational physics, 169 (2001), pp. 463–502.
[32] T. Qin, C.-W. Shu, and Y. Yang, Bound-preserving discontinuous galerkin methods for relativistic hydrodynamics, Journal of Computational Physics, 315 (2016), pp. 323–347.
[33] W. Rider and D. Kothe, Stretching and tearing interface tracking methods, in 12th computational fluid dynamics conference, 1995, p. 1717.
[34] W. J. Rider and D. B. Kothe, Reconstructing volume tracking, Journal of computational physics, 141 (1998), pp. 112–152.
[35] P. Ripa, Conservation laws for primitive equations models with inhomogeneous layers, Geophysical & Astrophysical Fluid Dynamics, 70 (1993), pp. 85–111.
[36] , On improving a one-layer ocean model with thermodynamics, Journal of Fluid Mechanics, 303 (1995), pp. 169–201.
[37] M. L. Sæ tra, A. R. Brodtkorb, and K.-A. Lie, Efficient GPU-implementation of adaptive mesh refinement for the shallow-water equations, J. Sci. Comput., 63 (2015), pp. 23–48.
[38] J. A. Sethian and P. Smereka, Level set methods for fluid interfaces, Annual review of fluid mechanics, 35 (2003), pp. 341–372.
[39] H. Shirkhani, A. Mohammadian, O. Seidou, and A. Kurganov, A well-balanced positivity-preserving central-upwind scheme for shallow water equations on unstructured quadrilateral grids, Comput. & Fluids, 126 (2016), pp. 25–40.
[40] E. Tadmor, Local error estimates for discontinuous solutions of nonlinear hyperbolic equations, SIAM J. Numer. Anal., 28 (1991), pp. 891–906.
[41] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y.-J. Jan, A front-tracking method for the computations of multiphase flow, Journal of computational physics, 169 (2001), pp. 708–759.
[42] S. O. Unverdi and G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of computational physics, 100 (1992), pp. 25–37.
[43] Y. Xing and X. Zhang, Positivity-preserving well-balanced discontinuous galerkin methods for the shallow water equations on unstructured triangular meshes, Journal of Scientific Computing, 57 (2013), pp. 19–41.
[44] X. Yang and A. J. James, Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids, Journal of computational physics, 214 (2006), pp. 41–54.
[45] X. Yang, A. J. James, J. Lowengrub, X. Zheng, and V. Cristini, An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids, Journal of Computational Physics, 217 (2006), pp. 364–394.

	$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h\rho$\kern-0.50003pt}}})^{n+1}_{j}=$	$\displaystyle\dfrac{1}{9}\sum_{m,s}h_{j}(M_{jm})\rho_{j}(M_{js})$
		$\displaystyle-\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\cos(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})u_{jk}(M_{jk})+a_{jk}^{\rm out}h_{j}(M_{jk})\rho_{j}(M_{jk})u_{j}(M_{jk})\Big{]}$
		$\displaystyle-\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\frac{\ell_{jk}\sin(\rho_{jk})}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\Big{[}a_{jk}^{\rm in}h_{jk}(M_{jk})\rho_{jk}(M_{jk})v_{jk}(M_{jk})+a_{jk}^{\rm out}h_{j}(M_{jk})\rho_{j}(M_{jk})v_{j}(M_{jk})\Big{]}$
		$\displaystyle+\dfrac{\Delta t}{\|T_{j}\|}\sum_{k=1}^{3}\ell_{jk}\frac{a_{jk}^{\rm in}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\big{[}h_{jk}(M_{jk})\rho_{jk}(M_{jk})-h_{j}(M_{jk})\rho_{j}(M_{jk})\big{]}$
	$\displaystyle=$	$\displaystyle\dfrac{1}{9}\sum_{m\neq s}h_{j}(M_{jm})\rho_{j}(M_{js})$
		$\displaystyle+\sum_{k=1}^{3}h_{j}(M_{jk})\rho_{j}(M_{jk})\Big{[}\dfrac{1}{9}-\dfrac{\Delta t}{\|T_{j}\|}\frac{\ell_{jk}a_{jk}^{\rm out}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm in}+u^{\perp}_{j}(M_{jk})\right)\Big{]}$
		$\displaystyle+\sum_{k=1}^{3}h_{jk}(M_{jk})\rho_{jk}(M_{jk})\dfrac{\Delta t}{\|T_{j}\|}\frac{\ell_{jk}a_{jk}^{\rm in}}{a_{jk}^{\rm in}+a_{jk}^{\rm out}}\left(a_{jk}^{\rm out}-u^{\perp}_{jk}(M_{jk})\right).$