A multiscale method applied to shallow water flow

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) A multiscale method applied to shallow water flow Anna Mart´ ınez Gavara1, Guillaume Chiavassa2, Rosa Donat1 1Dpt. de Matem`atica Aplicada, Univ. de Val`encia, Aptdo. 46100, Dr. Moliner,50 Burjassot (Val`encia). E-mails: [email protected], [email protected]. 2Ecole G´en´eraliste D’Ing´enieurs de Marseille, Technopˆole de Chˆateau-Gombert, 13383 Marseille Cedex, 13. E-mail: [email protected]. Palabras clave: 65M06, 65D05, 35L65 Resumen A flux-limited second order scheme with the C-property is used to solve the one dimensional or two dimensional Saint-Venant system for shallow water flows with non-flat bottom and friction terms, as is introduced in [7]. High resolution at low cost can be obtained by applying a point-value multiresolution transform [2, 3, 9] in order to detect regions with singularities. The above method is applied in these regions, while a cheap polynomial interpolation is used in the smooth zones, thus lowering the computational cost. 1. Introduction The numerical simulation of physical problems modeled by systems of conservation laws is difficult due to the presence of discontinuities in the solution. High-resolution shock capturing (HRSC) schemes succeed in computing highly accurate numerical solutions, typically second-or thirdorder in smooth regions, while maintaining sharp, oscillation-free numerical profiles at discontinuities. The power of a HRSC scheme lies usually in a computation of the numerical flux function of the scheme often expensive, which is the main drawback of these schemes, specially in multi-dimensional computations. However, the costly numerical flux function of a HRSC scheme is only necessary in a neighborhood of singularities, so, Harten in [9] proposed a scheme based on reducing the computational cost using the smoothness information of the data obtained from a multiresolution transform. The goal is to save time in the evaluation of numerical flux functions, replacing the expensive numerical flux evaluation with a cheap polynomial interpolation in the smooth regions. 1 A. Mart´ınez-Gavara, R. Donat, G. Chiavassa On the other hand, the shallow water system is used to model real-life applications in which the flow regime is steady or quasi-steady, and much effort has been devoted to design numerical techniques that are capable to preserve steady states at the discrete level as well as to accurately compute the evolution of small dynamical perturbations of these. The inclusion of the source term in a direct discretization of the system becomes a nontrivial issue, because many schemes do not respect stationary solutions. In [7], the authors seek to obtain an extension of the numerical scheme developed by Donat and Marquina in [4], that avoids the use of averaged quantities in computing the numerical flux function at cell interfaces, for non-homogeneous conservation laws by incorporating the idea of flux gradient and source term balancing in [6]. However, the extension based on the use of two spectral decompositions at each computational interface does not satisfy the exact C-property of [1], and a combined 1-Jacobian/2-Jacobian scheme is proposed. In this work, we apply the multilevel technique developed in [3] to the shallow water system. The basic underlying shock capturing scheme is the 1J−2Jscheme of [7]. The inclusion of the source term is, then, done directly through the numerical divergence operator, and the multilevel technique can be applied in a straightforward manner to the shallow water system. The paper is organized as follows: In section 2, we recall the basic ingredients of the multilevel algorithm. In section 3, we present the main steps of the method used to solve the 2D shallow water equations, as well as some remarks on the C-property for the scheme. Finally, we present several numerical experiments that validate the technique in section 4. 2. The multilevel algorithm The multilevel strategy has been described an analyzed in [3], here we only present the main steps. Let us consider a 2D system of hyperbolic conservation laws: Ut+F(U)x+G(U)y= 0 (1) where U is the vector of conserved quantities. We consider discretizations of this system on a Cartesian grid G0={(xi=i4x, yj=j4y), i = 0,· · · , Nx j = 0,· · · , Ny}using the semi-discrete formulation: dUij dt +Div(U)ij = 0,(2) with the numerical divergence computed as: Div(U)ij = Fi+1 2,j −Fi−1 2,j 4x+ Gi,j+1 2−Gi,j−1 2 4y,(3) where F(u1,· · · , uk+m) and G(u1,· · · , uk+m) are consistent numerical flux functions, which are the trademark of the scheme. The goal of the multilevel method is to reduce the cpu time associated to the underlying scheme by reducing the number of expensive flux evaluations. To undestand the basic mechanism, let us consider Euler’s method applied to (2), Un+1 ij =Un ij − 4tDiv(U)n ij (4) 2 A multiscale method applied to shallow water flow If both Unand Un+1 are smooth around (xi, yj) at time tn, then (4) implies that the numerical divergence is also smooth and we can avoid using the numerical flux functions of the scheme in its computation. On the other hand, if a discontinuity appears during the time evolution, the Riemann solver of the scheme has to be called to compute the numerical divergence. So, the smoothness analysis of Unto Un+1 and the computation of Div(U) using this information, are the most important part of the algorithm. As in [2], the computation of the numerical divergence Div(Un) on the finest grid is carried out in a sequence of steps. The numerical divergence is evaluated in all the points on the coarsest grid GLusing the numerical flux function of the scheme, and for the finer grids, the divergence is evaluated recursively, either by the same procedure or with a cheap interpolation using the values obtained on the coarser grids. The different resolution levels are specified by a set of nested grids {Gl, l = 1,· · · , L}given as follows, (xi, yj)∈ Gl⇐⇒ (x2li, y2lj)∈ G0(5) The information about the regularity of the data contained in the multiresolution transform of the numerical solution is used to determine a flag vector ((bl ij)l,ij), whose value (0 or 1) will determine the choice of the procedure to evaluate the divergence. Given a tolerance parameter ε, the value is obtained by applying two tests to the detail coefficients (or wavelet coefficients) (dl ij)l,i,j: if |dl ij| ≥ ε⇒bl i−k,j−m= 1 k, m =−2,· · · ,2 if |dl ij| ≥ 2rεand l > 1⇒bl−1 2i−k,2j−m= 1 k, m =−1,0,1 (6) where large values of the detail coefficients correspond to non-smooth zones of the solution like shocks or contact discontinuities, existing or in formation. The multilevel evaluation of the numerical divergence is made as follows: The divergence Div(U) is computed at the points of the coarsest grid using the scheme. Then for the finer grids, for example if the divergence is known on Gl, the values of Divl−1(U) on Gl−1are computed using the boolean flag: if bl ij = 1,compute Divl−1(U)ijdirectly with the scheme if bl ij = 0,compute Divl−1(U)ij =I[(xi, yj); Divl(U)],(7) where I[(xi, yj); Divl(U)] is a 2D polynomial interpolation of Div(U) around the point (xi, yj) using the values computed on Gl. The process is repeated from l=L, · · · ,1 and, once it is completed, we obtain the values of Div(U) on the finest grid G0, which are needed by the ODE solver. 3. The shallow water equations The shallow water equations is a hyperbolic system of conservation laws that approximately describes various geophysical flows. The source terms are due to topography, we do not consider wind effects and Coriolis force. The resulting system of equations becomes: Ut+F(U)x+G(U)y=S(8) 3 A. Mart´ınez-Gavara, R. Donat, G. Chiavassa   h q1 q2  t +   q1 q2 1 h+1 2gh2 q1q2 h   x +   q2 q1q2 h q2 2 h+1 2gh2   y =  0 −ghzx −ghzy   where his the water depth, q1and q2are the two components of the discharge, and zis the bottom topography. The numerical technique we use, follows [7]. We consider a method of lines approach in which the time integration is performed via a TVD-Runge-Kutta method (see [8]). The spatial terms are treated in a dimension by dimension fashion, thus, it is sufficient to give a description of the technique for the one-dimensional system. Let us consider the 1D system: Ut+F(U)x=S(x, U) (9) and rewrite in the form Ut+G(U)x= 0, with G(U, x) = F(U) + B(U, x) and B(U, x) = (0,Rx ¯xghzxds)T. As in [7], we use a semi-discrete formulation of the type Ut+ G+ i+1 2 −G− i−1 2 4x= 0 (10) where the computation of G± i+1 2 only involves integral terms over consecutive cell centers and follows the basic design strategy in Marquina’s flux formula: two states are computed at each side of a cell-interface, ULand UR, and the numerical flux functions are obtained by applying the scalar algorithm to “sided” local characteristic fluxes. The states ULand URat each side of a given interface are obtained by ENO interpolation of the physical variables as specified in [5]. Unless specifically stated, the order of the interpolation used to compute these states is the same as the order of the scheme. Given UL=UL i+1 2 and UR=UR i+1 2 , the left and right states at the i+1 2cell-interface, the flux functions G± i+1 2 shall be defined as G± i+1 2 = 2 X p=1 (ˆ G± i+1 2 )p,LRp(UL)+(ˆ G± i+1 2 )p,RRp(UR) (11) where Lp(UL), Rp(UL)(Lp(UR), Rp(UR)), p = 1,2, are the left and right eigenvectors of the Jacobian matrix J(U) = F0(U), associated to the eigenvalues λp(UL)(λp(UR)). ( ˆ G±,p i+1 2 )L,R are the local modified characteristic fluxes, whose high order terms involve only quantities of the form Bi,i+1 =Bi+1 −Bi= (0,Zxi+1 xi ghzxds)T(12) and the contribution of the source terms at first order depends only of the wind coming from the right (+) or left (−) at the interface (more details in [7]). 3.1. C-property Berm´udez and V´azquez [1] and V´azquez-Cend´on [10] discussed an approach for approximating source terms which is designed for quasi-steady and steady flow. Consider the shallow water equation for the quiescent flow case, u(x, t) = 0 and h(x, t) = D−z(x, t)∀(x, t) 4 A multiscale method applied to shallow water flow Figura 1: Exact water surface and topography at steady state. For this stationary case Ut= 0, if the source term approximation balances with the numerical fluxes, then the numerical scheme satisfies: the approximate C-property, if the numerical scheme is accurate to the order O(4x2) when applied to the quiescent flow case; the exact C-property, if the numerical scheme is exact when applied to the quiescent flow case. In this way, it is proven in [7] that if UL i+1 2 =UR i+1 2 (e.g. = Ui+Ui+1 2) (1J), the scheme verifies the exact C-property, and if UL i+1 2 6=UR i+1 2 (2J), the scheme verifies the approximate C- property, provided the order of accuracy is at least 2. Hence, the preferred option is to combine both, the 1J-2J scheme, to get the benefits of both alternatives. This is the scheme of our choice in the next section. 4. Numerical experiments This section is devoted to the presentation of the results obtained with the multilevel algorithm. In order to evaluate the quality and efficiency of the algorithm, there are some parameters to be tested. The quality is analyzed by measuring the difference between the multilevel solution Unand the reference one, Un ref , which corresponds to the scheme without multiresolution, in some appropriate norm (we choose the discrete l1-norm). The efficiency of the multilevel algorithm with respect to the reference simulation is controlled by two parameters, from one side, the percentage of numerical divergences computed directly per time step, %f, is an important quantity, but a more concret measure is given by θiter, the cpu gain for a given iteration, and θ, the gain for the global simulation. 5 A. Mart´ınez-Gavara, R. Donat, G. Chiavassa First of all, we are going to validate numerically the C-property in the multilevel algorithm, we follow [11] and consider a smooth topography given by z(x, y) = 0,8e−50((x−0,5)2+(y−0,5)2)(13) with (x, y)∈[0,1] ×[0,1], we see that in a quiescent state (q1=q2= 0), we maintain the steady flow (see figure 1). In table 1, we show the measurements of the L1-error and the percentage of numerical divergences computed directly per time step with multiresolution. We can see that the l1-error is of the order of roundoff error. The C-property is, thus, preserved. Grid size G0%fmin - %fmax l1-error 257×257 6,5784 - 6,5784 5,4674 ·10−15 513×513 1,6510 - 1,6510 1,1376 ·10−14 Tabla 1: Steady state with smooth topography. L1-error and percentage of divergence computed Next, we consider a test containing shocks and rarefaction waves (Leveque 2D test [11]). The bottom topography is given as z(x, y) = 0,5e−50((x−0,5)2+(y−0,5)2)(14) on [0,1] ×[0,1] with g= 1. The initial conditions are q1=q2= 0 and h(x, y) = ½1,01 −z(x, y),0,1< x < 0,2; 1−z(x, y), otherwise. (15) In figure 2, we display the level curves of the numerical solution obtained with and without the multilevel algorithm, and we can observe that the numerical simulation is of the same “quality” as the reference simulation. On Figure 3 we show the l1-error measured for variable h(x, y) when applying the multilevel algorithm with Nx =Ny = 64 and 2 levels of refinement, for different values of the tolerance ε. We can observe that, as in [3] the closeness to the reference simulation, can be controlled by adjusting the tolerance suitably. Grid size G0%fmin - %fmax cpu gain θ 64 ×64 62.77 - 81.17 1.2838 128 ×128 33.20 - 56.99 2.0174 256 ×256 15.49 - 30.76 3.5156 Tabla 2: Leveque 2D test at time t= 0,7. Percentage of resolved flux and cpu gain Finally, we show in Table 2 the global gain for each simulation and the maximum and the minimum values for %fin the simulation. We can observe that the finer the grid, the smaller the percentage of direct flux evaluations. In Figure 4 we represent θ(t) and %f(t). We can observe that, there are few non-smooth structures in the flow, the gain is quite large for fine grids. The behavior of θ(t) is roughly inversely proportional to that of %f(t). 6 A multiscale method applied to shallow water flow x y reference 129x129 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.999 1 1.001 1.002 1.003 1.004 1.005 1.006 x y multiresolution 129x129 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.999 1 1.001 1.002 1.003 1.004 1.005 1.006 Figura 2: Leveque 2D test at time t= 0,7.Left: reference simulation. Right: multilevel simulation with ε= 10−4 10−4 10−3 10−2 10−6 10−5 10−4 10−3 10−2 tolerance ε error e1 64 x 64 Figura 3: Leveque 2D test at time t= 0,7. Error between the multilevel algorithm and the reference one for different values of ε.Nx =Ny = 64 and L= 2 Referencias [1] A. Berm´udez, M.E. V´azquez, Upwind Methods for Hyperbolic Conservation Laws with Source Terms, Computers and Fluids, Vol. 23 No. 8, 1049-1071 (1994) [2] B. Bihari, A. Harten. Multiresolution shemes for the numerical solutions of 2D conservation laws. SIAM J. Sci. Comput., 18 (1997), 315-354. [3] G. Chiavassa, R. Donat Point value multiscale algorithms for 2D compressible flows, SIAM J. Sci. Compt., Vol. 23, No. 3, pp. 805-823. [4] R. Donat, A. Marquina. Capturing Shock Reflections: An improved Flux Formula, Journal of Computational Physics, 125:42-58, 1996. [5] R.P. Fedkiw, B. Merriman, R. Donat, S. Osher. The Penultimate Scheme for Systems of Conservation Laws: Finite Difference ENO with Marquina’s Flux Splitting. Progress in Numerical Solutions of Partial Differential Equations, Arcachon, France, edited by M.Hafez, July 1998. [6] Ll. Gasc´on, J.M. Corber´an, “ Construction of Second-Order TVD Schemes for Nonhomogeneous Hyperbolic Conservation Laws”, Journal of Computational Physics,172, 261–297 (2001). [7] G. Haro, Numerical simulation of shallow water equations amd some physical models in image processing. Ph.D.Thesis, Departament of Technologies, Universitat Pompeu Fabra, Barcelona, 2005. 7 A. Mart´ınez-Gavara, R. Donat, G. Chiavassa 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 time t cpu gain ε=1e−3 64 x 64 128 x 128 256 x 256 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 10 20 30 40 50 60 70 80 90 time t %f ε=1e−3 64 x 64 256 x 256 Figura 4: Leveque 2D test at time t= 0,7. Time evolution of θand %ffor different initial grid.Left: cpu gain. Right: percentage flux [8] C. W. Shu, S. Osher. Efficient Implementation of Essentially Non-Oscillatory Shock Capturing Schemes II. Journal of Computational Physics, 83:32-78,1989. [9] A. Harten. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Comm. Pure Appl. Math., 48 (1995), 1305-1342. [10] M.E. V´azquez-Cend´on, Improved Treatment of Source Terms in Upwind Schemes for the Shallow Water Equations in Channels with Irregular Geometry, J.Comput. Phys. 148:497-526 (1999). [11] Y. Xing, C.W. Shu, High order finite difference WENO schemes with the exact conservation property for the shallow water equations. Journal of Computational Physics, Vol. 208, No. 1, 206-227 (2005). 8