A second order pvm flux limiter method. Application to magnetohydrodynamics and shallow stratified flows

A second order PVM flux limiter method. Application to magnetohydrodynamics and shallow stratified flows M.J. Castro-D´ıaz∗, E. D. Fern´andez-Nieto†, G. Narbona-Reina† , M. de la Asunci´on‡ June 23, 2014 Abstract In this work we propose a second order flux limiter finite volume method, named PVM-2U-FL, that only uses information of the two external waves of the hyperbolic system. This method could be seen as a natural extension of the well known WAF method introduced by Prof. Toro in [21]. We prove that independently of the number of unknowns of the 1D system, it recovers the second order accuracy at regular zones, while in presence of discontinuities, the scheme degenerates to PVM-2U method, which can be seen as an improvement of the HLL method (see [4], [8]). Another interesting property of the method is that it does not need any spectral decomposition of the Jacobian or Roe matrix associated to the flux function. Therefore, it can be easily applied to systems with a large number of unknowns or in situations where no analytical expression of the eigenvalues or eigenvectors are known. In this work, we apply the proposed method to Magnetohydrodynamics and to stratified multilayer flows. Comparison with the twowaves WAF and HLL-MUSCL methods are also presented. The numerical results show that PVM-2U-FL is the most efficient and accurate among them. Key words: Finite Volumes, flux limiters, riemann solver, second order, magnetohydrodinamic, multilayer, stratified flows 1 Introduction The goal of this article is to design a robust, simple and fast second order flux limiter numerical scheme for solving one dimensional hyperbolic systems. An interesting technique to obtain second order accurate and robust schemes is to use a non-linear combination of first and second order methods in terms of flux limiters functions. An example of this type scheme can be defined ∗Departamento de An´alisis Matem´atico. Universidad de M´alaga.([email protected]) †Departamento de Matem´atica Aplicada I, E.T.S. Arquitectura. Universidad de Sevilla. Avda. Reina Mercedes 2. 41012 Sevilla, Spain ([email protected], [email protected]). ‡Dpto. Lenguajes y Sistemas Inform´aticos, Universidad de Granada ([email protected]) 1 by a suitable combination of Roe method (which is only first order) near discontinuities, and the Lax-Wendroff method (which is second order in space and time) in regular areas. Note that the previous scheme requires the explicit knowledge of the eigenstructure of the system, which is not straight forward for some hyperbolic systems, making this scheme computationally expensive in those cases. It is also well known that the use of incomplete Riemann solvers as Rusanov, Lax-Friedrichs, HLL, among others (see [11], [25], [5], [8], [29]) allows one to reduce the computing time required by a Roe solver (see, for instance, [9]). Although Roe scheme gives, in general, a better resolution of the discontinuities than incomplete Riemann solvers, when combined with high order methods may be indistinguishable. In [4] Castro and Fern´andez-Nieto introduce a family of incomplete simple Riemann solvers named as PVM (Polynomial Viscosity Matrix), for conservative and nonconservative hyperbolic systems, defined in terms of viscosity matrices computed by a suitable polynomial evaluation of a Roe linearization, that overcome the difficulty of the computation of the spectral decomposition of Roe matrices. PVM schemes can be seen as the natural extension of the one proposed in [8] for balance laws, and, more generally, for nonconservative systems. An interesting numerical scheme that uses flux limiters functions is the WAF (Weighted Average Flux) method, introduced by Toro in [21]. It is a one-step Godunov-type method to solve hyperbolic conservation laws that achieves second order accuracy by averaging the solution of the Riemann problem with piecewise constant initial data. As it is well known, due to Godunov’s theorem, linear schemes with high order accuracy generate spurious oscillations near discontinuities. To avoid this problem, WAF method uses flux limiter functions. The resulting scheme is a non-linear TVD (Total Variation Diminishing) scheme with second order accuracy. WAF scheme has been extensively used to approximate hyperbolic systems, see for example [22], [23], [10], [24], [28]. It has been also used as the base of higher order numerical solvers (see [27]). Nevertheless, WAF method needs the explicit knowledge of the structure of the approximated Riemann problem to achieve second order accuracy. For example, if we only consider the information of the two external waves and we use the HLL intermediate flux we obtain a WAF method – that we will name in what follows HLL-WAF method– that has second order accuracy for 1D 2x2 hyperbolic conservative systems. The main objective of this paper is to obtain a new flux limiter scheme that only uses the information of the two external waves, like the HLL-WAF scheme, and that achieves second order accuracy for 1D N×Nhyperbolic systems with N≥2. The resulting scheme can be seen as a natural extension of the original HLL-WAF scheme and it is defined in terms of a non-linear combination of a suitable PVM scheme, that is first order, with the second order Lax-Wendroff scheme. The paper is organized as follows: in Section 2, first we summarize how WAF and, in particular, HLL-WAF methods are derived. Next, HLL-WAF method is rewritten as a nonlinear combination of two PVM schemes. In section 3, the new flux limiter scheme is defined and finally, some numerical tests for the 1D ideal magnetohydrodynamics and the multilayer shallow-water systems are presented. Comparison with HLL-WAF and the second order HLL methods with MUSCL (see [12], [13], [32]) state reconstruction are also provided. 2 2 Preliminaries In this section we summarize the derivation of WAF method introduced by Prof. E.F. Toro in [21]. Let us consider the conservative hyperbolic system wt+F(w)x= 0, x ∈[0, L], t ∈[0, T],(2.1) where w(x, t) takes values on an open convex set O⊂RN,Fis a regular function from Oto RN. Let us consider a partition of the domain {xi}i={i∆x}iwhere, by simplicity, ∆xis supposed to be constant, and we denote tn=n∆twhere ∆tis the time step. A finite volume method in conservative form to approximate (2.1) can be written as wn+1 i=wn i−∆t ∆x(Fn i+1/2−Fn i−1/2),(2.2) where wn idenotes an approximation of the mean value of the solution on the control volume (xi−1/2, xi+1/2) at time t=tn: wn i≈1 ∆xZxi+1/2 xi−1/2 w(x, tn)dx, and Fn i+1/2=F(wn i, wn i+1) denotes the numerical flux function that characterizes each method. Let us consider a Riemann problem associated to (2.1) with initial data wn iand wn i+1:        wt+F(w)x= 0, w(x, 0) = wix < 0; wi+1 x > 0, (2.3) where we have removed superindex nfor sake of simplicity. In what follows, the dependency of the intercell i+ 1/2 will be dropped for clarity if there is no ambiguity. Le us denote Slfor l= 1,··· , N the approximation of the characteristic velocities and let us consider the computational cell V= [−∆x/2,∆x/2] ×[0,∆t]. Then, the WAF numerical flux is obtained by integrating the physical flux in Vusing the midpoint rule for the time integral: FWAF i+1/2=1 ∆xZ∆x/2 −∆x/2 F( ˜w(x, ∆t 2))dx, (2.4) where ˜wis an approximated solution of the Riemann problem (2.3), composed by N+1 constant states. If we define ωk,k= 0,··· , N + 1 (see Figure 1 for the case N= 2) as ωk=1 2(ck−ck−1); with c0=−1, cN+1 = 1,and cl=∆t ∆xSl,for 1 ≤l≤N, (2.5) 3 0 t x S !t !t / 2 S 1 ""3 1 2 "2 #! ! x /2 x /2 Figure 1: Computational grid and waves to compute the WAF method for problem (2.1) (N= 2) then, FWAF i+1/2can be rewritten as FWAF i+1/2= N+1 X k=1 ωkF(k) i+1/2,(2.6) where F(k) i+1/2is the value of the flux function in the interval k. Finally, taking into account the definition of ωk, we have that: FWAF i+1/2=1 2(Fi+Fi+1)−1 2 N X k=1 ck∆F(k) i+1/2,(2.7) where ∆F(k) i+1/2=F(k+1) i+1/2−F(k) i+1/2, and Fi=F(wi). The WAF scheme is second order accurate in time and space, therefore according to Godunov’s theorem, it produces spurious oscillations for non-smooth solutions. To overcome this fact, a TVD stabilization must be performed. If we denote by χ(v) a flux limiter function, then a limiter function can be defined by Ψ(v, c)=1−(1 −|c|)χ(v), and the TVD-WAF flux function becomes as follows: FWAF i+1/2=1 2(Fi+Fi+1)−1 2 N X k=1 sign(ck)Ψk∆F(k) i+1/2,(2.8) where Ψk= Ψ(v(k), ck) = 1 −(1 −|ck|)χ(v(k)).(2.9) Some suitable choices for χcan be found in [26]. In this work we consider the Beam-Warming limiter: χ(v(k)) = min(max(0, v(k)),1). 4 For v(k)=v(k)(Sk) at the interface xi+1/2we consider the following definition: v(k)(Sk) =                    ¯m((pi+1 −pi−1)/2, pi+1 −pi, pi−pi−1) pi+1 −pi ,if Sk>0, ¯m((pi+2 −pi)/2, pi+1 −pi, pi+2 −pi+1) pi+1 −pi ,if Sk<0, 1 if |pi+1 −pi| ≤ ε 1≤k≤N. (2.10) In this definition ¯mis the minmod limiter: ¯m(a, b, c) = sgn(a) + sgn(b) 2 sgn(b) + sgn(c) 2min(|a|,|b|,|c|); {pj}j=i+2 j=i−1is a set of scalar values that depend on the problem and εis a small parameter (in the numerical tests we will consider ε= ∆x3). Finally, using the former definition of Ψkand ck, TVD-WAF flux can be written as follows: FWAF i+1/2=1 2(Fi+Fi+1)−1 2 N X k=1 sgn(Sk)(1 −χk) ∆F(k) i+1/2−1 2 ∆t ∆x N X k=1 Skχk∆F(k) i+1/2,(2.11) where χk=χ(v(k)). 2.1 Two-waves WAF method In this section we consider the TVD-WAF method resulting when only the fastest (SR) and the slowest (SL) wave of the Riemann problem are used. Let us denote by χLand χRthe corresponding flux limiter evaluations. Now, as in the original paper of Prof. Toro (see [21]), using the HLL flux to evaluate the intermediate flux F(2) i+1/2, F(2) i+1/2=SRFi−SLFi+1 +SRSL(wi+1 −wi) SR−SL and taking into account that F(1) i+1/2=Fi=F(wi) and F(3) i+1/2=Fi+1 =F(wi+1), the two-waves WAF scheme (HLL-WAF in what follows) can be written: FHLL-WAF i+1/2=1 2(Fi+Fi+1)−1 2(ν1(χL, χR)(wi+1 −wi) + ν2(χL, χR)(Fi+1 −Fi)) −1 2 ∆t ∆x(µ1(χL, χR)(wi+1 −wi) + µ2(χL, χR)(Fi+1 −Fi)) ,(2.12) where 5 ν1(χL, χR) = SLSR((1 −χL)sgn(SL)−(1 −χR)sgn(SR)) SR−SL , ν2(χL, χR) = (1 −χR)|SR|−(1 −χL)|SL| SR−SL , µ1(χL, χR) = SLSR(SLχL−SRχR) SR−SL ,(2.13) µ2(χL, χR) = S2 RχR−S2 LχL SR−SL . In [4] a family of first order finite volume methods named PVM-lmethod is proposed. For the case of conservative systems in the form of (2.3), they can be defined in terms of a numerical flux function written as follows Fi+1/2=1 2(Fi+1 +Fi)−1 2Qi+1/2(wi+1 −wi).(2.14) The numerical viscosity matrix Qi+1/2, is defined in terms of a Roe Matrix Ai+1/2associated to F(w), that is, Ai+1/2verifies Fi+1 −Fi=Ai+1/2(wi+1 −wi).(2.15) In particular Qi+1/2is given by a polynomial evaluation of this Roe Matrix as Qi+1/2=Pi+1/2 l(Ai+1/2),(2.16) where Pi+1/2 l(x) is a polynomial of degree lverifying Pi+1/2 l(x) = l X j=0 αi+1/2 jxj,such that Pi+1/2 l(x)≥ |x| ∀x∈[SL, SR].(2.17) Taking into account the Roe property (2.15), we can write (2.12) under the structure of (2.14) by simply defining QHLL−W AF i+1/2(χL, χR) = QHLL−W AF o1,i+1/2(χL, χR) + ∆t ∆xQHLL−W AF o2,i+1/2(χL, χR),(2.18) with QHLL−W AF o1,i+1/2(χL, χR) = ν1(χL, χR)I+ν2(χL, χR)Ai+1/2 QHLL−W AF o2,i+1/2(χL, χR) = µ1(χL, χR)I+µ2(χL, χR)Ai+1/2.(2.19) That is, the usual two-waves HLL-WAF method (2.12) can be seen as a combination of two PVM schemes whose viscosity matrices are Qo1,i+1/2and Qo2,i+1/2, respectively, associated to the first degree polynomials: Po1 1(x) = ν1(χL, χR) + ν2(χL, χR)xand Po2 1(x) = µ1(χL, χR) + µ2(χL, χR)x. 6 The HLL scheme can also be interpreted as a PVM method (see [4] for details) for a viscosity matrix Qi+1/2=P1U(Ai+1/2) and the polynomial P1U(x) = SR|SL|−SL|SR| SR−SL +|SR|−|SL| SR−SL x. (2.20) Indeed, since HLL-WAF method is based on the HLL method, we can find a relation between their definitions as PVM schemes through the polynomials that define respectively their viscosity matrices. Thus, we can check that QHLL−W AF o1,i+1/2(χL, χR) = sgn(SL)(1 −χL) + sgn(SR)(1 −χR) 2Ai+1/2(2.21) +sgn(SR)(1 −χR) 2P1M(Ai+1/2)−sgn(SL)(1 −χL) 2P1M(Ai+1/2) QHLL−W AF o2,i+1/2(χL, χR) = SLχL+SRχR 2Ai+1/2+SRχR 2P1M(Ai+1/2)−SLχL 2P1M(Ai+1/2). where P1M(x) = −2SRSL SR−SL +SR+SL SR−SL x. (2.22) Note that for the case SL<0< SR, we have that P1M(x) = P1U(x), the polynomial associate to the HLL method. Moreover, it can be seen that P1U(x) = sgn(SL) + sgn(SR) 2x+sgn(SR)−sgn(SL) 2P1M(x).(2.23) Remark 2.1. Notice that: •If χL=χR= 0, then QHLL−WAF i+1/2(χL, χR) = QHLL−W AF o1,i+1/2(χL= 0, χR= 0). Then, FHLL-WAF i+1/2 reduces to the usual HLL flux. •If χL=χR= 1, then QHLL−WAF i+1/2(χL, χR) = ∆t ∆xQHLL−W AF o2,i+1/2(χL= 1, χR= 1). Then, for the case N= 2,FHLL-WAF i+1/2reduces to the usual Lax-Wendroff method. Therefore, HLLWAF method achieves second order accuracy in space and time for 1D, 2×2systems, but it is not true for N > 2. 3 A two-waves PVM flux limiter method (PVM-2U-FL) As we have seen in previous section HLL-WAF method can be interpreted as an improvement of the HLL method to achieve a second order accuracy scheme satisfying a TVD property through a flux limiter function for 2 ×2 1D systems. In this section we propose a new method based on the same idea, where the first order method is now replace by a suitable PVM scheme such 7 as the resulting method achieves second order accuracy at regular areas independently of the dimension of the system. In particular we present a new two-waves PVM flux limiter scheme with the following properties: •if χL=χR= 0 the method reduces to the first order PVM-2U method introduced in [8] and extended in [4]. •if χL=χR= 1 the method reduces to the usual Lax-Wendroff method for N≥2. Let us recall that the PVM-2U method is defined by the second degree polynomial, P2U(x) verifying that: P2U(SL) = |SL|, P2U(SR) = |SR|and P0 2U(SM) = sgn(SM),(3.1) where SM=SLif |SL|≥|SR|, SRif |SL|<|SR|. Moreover, the following relation can be derived: P2U(x) = sgn(SL) + sgn(SR) 2x+sgn(SR)−sgn(SL) 2P2,¯α(x).(3.2) where P2,¯α(x) = ¯αP2M(x) + (1 −¯α)P1M(x),(3.3) with ¯α=(SR−SL)sgn(SM)−(SR+SL) 4SM−2(SL+SR),(3.4) the polynomial P1M(x) is defined by (2.22) and P2M(x) = −SR+SL SR−SL x+2 SR−SL x2. Remark 3.1. The PVM-2U method can be seen as a generalization of HLL method, in the sense that the PVM-2U method can be obtained from the HLL method by replacing the polynomial P1M by the polynomial P2,¯α(see equations (2.23) and (3.2)). Moreover, if ¯α= 0 then P2,¯α=P1M. Remember that HLL-WAF method has been defined by (2.21) in terms of the polynomial P1M(x). We propose to define a new flux-limiter type scheme in terms of the polynomial P2,¯α(x) as follows: F2U-FL i+1/2=Fi+Fi+1 2−1 2Q2U−FL i+1/2(χL, χR)(wi+1 −wi),(3.5) where Q2U−FL i+1/2(χL, χR) is defined by: Q2U−FL i+1/2(χL, χR) = Q2U−FL o1,i+1/2(χL, χR) + ∆t ∆xQ2U−FL o2,i+1/2(χL, χR),(3.6) 8 with Q2U−FL o1,i+1/2(χL, χR) = sgn(SL)(1 −χL) + sgn(SR)(1 −χR) 2Ai+1/2 +sgn(SR)(1 −χR) 2P2,αR(Ai+1/2)−sgn(SL)(1 −χL) 2P2,αL(Ai+1/2), (3.7) and Q2U−FL o2,i+1/2(χL, χR) = SLχL+SRχR 2Ai+1/2+SRχR 2P2,αR(Ai+1/2)−SLχL 2P2,αL(Ai+1/2),(3.8) where αK= 1 −(1 −χK)(1 −¯α), K =L, R, with ¯αdefined by (3.4). Note that Q2U−F L o1,i+1/2(χL, χR) (respectively Q2U−F L o2,i+1/2(χL, χR)) can be obtained from QHLL−W AF o1,i+1/2(χL, χR) (respectively QHLL−W AF o2,i+1/2(χL, χR), by replacing the polynomial P1Mby the polynomials P2,αRor P2,αLif the right or left limiters are involved, respectively. Proposition 3.1. We have the following results: a)If χL=χR= 0 then, Q2U−FL i+1/2(χL= 0, χR= 0) = Q2U−FL o1,i+1/2(χL= 0, χR= 0) = P2U(Ai+1/2). Therefore, the two-waves PVM Flux limiter method reduces to the first order PVM-2U method. b)If χL=χR= 1 then, Q2U−F L i+1/2(χL= 1, χR= 1) = ∆t ∆xQ2U−FL o2,i+1/2(χL= 1, χR= 1) = ∆t ∆xA2 i+1/2. That is, the two-waves PVM Flux limiter method reduces to the Lax-Wendroff method. Proof: Let us suppose that χR=χL= 0, then αL=αR= ¯α. Therefore, using equations (3.3) and (3.4) we have that P2,αL(x) = P2,αR(x) = ¯αP2M(x) + (1 −¯α)P1M(x) = P2,¯α(x). Then, by (3.2) we obtain that Q2U−FL o1,i+1/2(χL= 0, χR= 0) = P2U(Ai+1/2). If χR=χL= 1, then αR=αL= 1, therefore P2,αL=P2,αR=P2M(x). And P2M(x) is the polynomial such that Q2U−FL o2,i+1/2(χL= 1, χR= 1) = A2 i+1/2.  Finally, Q2U−F L i+1/2can be written in a more compact form as follows: Q2U−FL i+1/2(χL, χR) = γ0,i+1/2Id +γ1,i+1/2Ai+1/2+γ2,i+1/2A2 i+1/2,(3.9) where γj,i+1/2=γj(χL, χR)∈R,j= 1,2,3, are given by: 9 HLL-MUSCL provide similar results. Observe that for a fixed mesh, PVM-2U-FL is also the most accurate, being HLL-MUSCL more accurate than HLL-WAF in this case. −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x vx Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x vx Ref. sol HLL−MUSCL PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 x vx Ref. sol HLL−WAF PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 x vx Ref. sol HLL−MUSCL PVM−2U−FL Figure 3: Velocity vxfor the Brio-Wu shock tube problem 4.1: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 4.2 High Mach shock tube problem This problem was presented also in [1] with the aim of testing the robustness of the numerical schemes for high Mach number flows. The initial conditions are (ρ, vx, vy, vz, Bx, By, Bz, P) = ((1,0,0,0,0,1,0,1000) for x≤0, (0.125,0,0,0,0,−1,0,0.1) for x > 0, and we take γ= 2. The Mach number of the right-moving shock is 15.5. The problem has been solved in [−1,1] using 200 grid points, CFL coefficient 0.8 and final time t= 0.012. For this test we also found that PVM-2U-FL provide the best results.They are plotted in Figures 7-10. A reference solution computed using HLL method with 25600 points has been also considered. 16 −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x By Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x By Ref. sol HLL−MUSCL PVM−2U−FL −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x By Ref. sol HLL−WAF PVM−2U−FL −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x By Ref. sol HLL−MUSCL PVM−2U−FL Figure 4: Magnetic field Byfor the Brio-Wu shock tube problem 4.1: HLL-WAF and PVM2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 4.3 Non-planar Riemann problem A non-planar Riemann problem with solution containing two strong rotational waves was proposed in [30]. The initial conditions are given by (ρ, vx, vy, vz, Bx, By, Bz, P) = ((1.7,0,0,0,1.1,1,0,1.7) for x≤0, (0.2,0,0,1.4968909,1.1,cos β, sin β, 0.2) for x > 0, where β= 2.3. Notice that although the problem has an unique solution, the initial conditions are close to initial conditions for which the problem admits non-unique solutions (see [30]). Figures 11-15 show the solution computed in the interval [−1,1.5] with 800 grid points, CFL number 0.8, γ= 5/3 and final time t= 0.4. Finally, an efficiency curve is shown in Figure 16, where CPU time vs error is shown in log scale for different mesh sizes from 100 up to 1600 grid points. Similar results to those in Test 1 are obtained: PVM-2U-FL is the most efficient among them and HLL-WAF and HLL-MUSCL provide similar results. Observe that PVM-2U-FL is also the most accurate among them for a fixed mesh, being HLL-MUSCL more accurate than HLL-WAF for a given mesh. 17 −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x P Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x P Ref. sol HLL−MUSCL PVM−2U−FL −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x P Ref. sol HLL−WAF PVM−2U−FL −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x P Ref. sol HLL−MUSCL PVM−2U−FL Figure 5: Hydrostatic pressure Pfor the Brio-Wu shock tube problem 4.1: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 10−210−1100101102 10−3 10−2 10−1 100 CPU time L1-Error HLL−MUSCL HLL−WAF PVM−2U−FL Figure 6: Brio-Wu shock tube problem 4.1: efficiency curve for HLL-MUSCL, HLL-WAF and PVM-2U-FL schemes. 18 −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x ρ Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x ρ Ref. sol HLL−MUSCL PVM−2U−FL 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.1 0.15 0.2 0.25 0.3 0.35 0.4 x ρ Ref. sol HLL−WAF PVM−2U−FL 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.1 0.15 0.2 0.25 0.3 0.35 0.4 x ρ Ref. sol HLL−MUSCL PVM−2U−FL Figure 7: Mass density ρfor the high Mach shock tube problem 4.2: HLL-WAF and PVM-2UFL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 x vx Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 x vx Ref. sol HLL−MUSCL PVM−2U−FL Figure 8: Velocity vxfor the high Mach shock tube problem 4.2: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). 19 −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 x By Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 x By Ref. sol HLL−MUSCL PVM−2U−FL 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 x By Ref. sol HLL−WAF PVM−2U−FL 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 x By Ref. sol HLL−MUSCL PVM−2U−FL Figure 9: Magnetic field Byfor the high Mach shock tube problem 4.2: HLL-WAF and PVM2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 2 2.5 3 x P Ref. sol HLL−WAF PVM−2U−FL −1−0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 2 2.5 3 x P Ref. sol HLL−MUSCL PVM−2U−FL Figure 10: Hydrostatic pressure Pfor the high Mach shock tube problem 4.2: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). 20 −1−0.5 0 0.5 1 1.5 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x ρ Ref. sol HLL−WAF PVM−2U−FL −1−0.5 0 0.5 1 1.5 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x ρ Ref. sol HLL−MUSCL PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.5 1 1.5 x ρ Ref. sol HLL−WAF PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.5 1 1.5 x ρ Ref. sol HLL−MUSCL PVM−2U−FL Figure 11: Mass density ρfor the non-planar Riemann problem 4.3: HLL-WAF and PVM-2UFL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 21 −1−0.5 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x vx Ref. sol HLL−WAF PVM−2U−FL −1−0.5 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x vx Ref. sol HLL−MUSCL PVM−2U−FL −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 x vx Ref. sol HLL−WAF PVM−2U−FL −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 x vx Ref. sol HLL−MUSCL PVM−2U−FL Figure 12: Velocity vxfor the non-planar Riemann problem 4.3: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 22 −1−0.5 0 0.5 1 1.5 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 x vy Ref. sol HLL−WAF PVM−2U−FL −1−0.5 0 0.5 1 1.5 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 x vy Ref. sol HLL−MUSCL PVM−2U−FL 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 x vy Ref. sol HLL−WAF PVM−2U−FL 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 x vy Ref. sol HLL−MUSCL PVM−2U−FL Figure 13: Velocity vyfor the non-planar Riemann problem 4.3: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 23 −1−0.5 0 0.5 1 1.5 0 0.5 1 1.5 x vy Ref. sol HLL−WAF PVM−2U−FL −1−0.5 0 0.5 1 1.5 0 0.5 1 1.5 x vy Ref. sol HLL−MUSCL PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 x vy Ref. sol HLL−WAF PVM−2U−FL −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 x vy Ref. sol HLL−MUSCL PVM−2U−FL Figure 14: Velocity vzfor the non-planar Riemann problem 4.3: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 24 −1−0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 x By Ref. sol HLL−WAF PVM−2U−FL −1−0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 x By Ref. sol HLL−MUSCL PVM−2U−FL −0.4 −0.3 −0.2 −0.1 0 0.1 0 0.2 0.4 0.6 0.8 1 x By Ref. sol HLL−WAF PVM−2U−FL −0.4 −0.3 −0.2 −0.1 0 0.1 0 0.2 0.4 0.6 0.8 1 x By Ref. sol HLL−MUSCL PVM−2U−FL Figure 15: Magnetic field Bzfor the non-planar Riemann problem 4.3: HLL-WAF and PVM2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). General view (top) and a zoom (down). 10−210−1100101102 10−3 10−2 10−1 100 CPU time L1-Error HLL−MUSCL HLL−WAF PVM−2U−FL Figure 16: Non-planar Riemann problem 4.3: efficiency curve for HLL-MUSCL, HLL-WAF and PVM-2U-FL schemes. 25 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x z Water surface Ref. inteface HLL−MUSCL PVM−2U−FL Bottom 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x z Water surface Ref. inteface HLL−WAF PVM−2U−FL Bottom Figure 21: Free surface and comparison of the interface for problem 5.3: HLL-WAF and PVM2U-FL schemes (top) and HLL-MUSCL and PVM-2U-FL schemes (down). 32 3 3.5 4 4.5 5 5.5 6 6.5 7 0 0.5 1 1.5 2 x z Water surface Ref. inteface HLL−MUSCL PVM−2U−FL Bottom 3 3.5 4 4.5 5 5.5 6 6.5 7 0 0.5 1 1.5 2 x z Water surface Ref. inteface HLL−WAF PVM−2U−FL Bottom Figure 22: Zoom of Figure 21: HLL-WAF and PVM-2U-FL schemes (left) and HLL-MUSCL and PVM-2U-FL schemes (right). 100101102103104105 10−2 10−1 100 101 CPU Error HLL−MUSCL HLL−WAF PVM−2U−FL Figure 23: Internal dam break problem 5.3: efficiency curve for HLL-MUSCL, HLL-WAF and PVM-2U-FL schemes. 33 6 Conclusions In this paper we present a computationally fast and efficient second order flux limiter finite volume method. It can be seen as a generalization of the HLL-WAF method, in the sense that it uses flux limiter functions to combine an incomplete Rieman solver with another one that allows to recover the second order accuracy in smooth regions. Moreover, both methods: HLL-WAF and PVM-2U-FL, only uses information of the two external waves. But: (i) the HLL-WAF method only recovers second order accuracy for 1D systems with two unknowns, while PVM-2U-FL recovers second order for arbitrary 1D systems; (ii) HLL-WAF degenerates to HLL method near discontinuities, while the PVM-2U-FL method degenerates to the PVM2U method. Let us remark that the PVM-2U can be seen as a generalization of the HLL method. Application to conservative and nonconservative systems are provided. In both cases, the numerical tests show that PVM-2U-FL is the most efficient among the compared methods and HLL-WAF and HLL-MUSCL provide similar results for MHD, but HLL-MUSCL provides better results than HLL-WAF for the multilayer shallow water system. Moreover, PVM-2U-FL is also the most accurate among them for a fixed mesh, being HLL-MUSCL more accurate than HLL-WAF for a given mesh. The extension to higher dimensions is not straightforward and it will be considered in future works. References [1] M. Brio, C. C. Wu. An upwind differencing scheme for the equations of ideal magnetohydrodynamics. J. Comput. Phys. 75: 400–422, (1988). [2] P. Cargo, G. Gallice. Roe matrices for ideal MHD and systematic construction of Roe matrices for systems of conservation laws. J. Comput. Phys. 136: 446–466, (1997). [3] M.J. Castro, P.G. LeFloch, M.L. Mu˜noz, and C. Par´es. Why many theories of shock waves are necessary: Convergence error in formally path-consistent schemes. Jour. Comp. Phys. 3227: 8107–8129, (2008). [4] M.J. Castro, E.D. Fern´andez-Nieto. A class of computationally fast first order finite volume solvers: PVM methods. SIAM J. Sci. Comput. 34(4): 2173–2196, (2012). [5] M.J. Castro, A. Pardo, C. Par´es, and E.F. Toro. On some fast well-balanced first order solvers for nonconservative systems. Math. Comp, 79 (271): 1427–1472, (2010). [6] S. F. Davis. Simplified Second-Order Godunov-Type Methods. SIAM J. Sci. Stat. Comput., 9: 445–473, (1988). [7] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74: 483–548, (1995). 34 [8] P. Degond, P.F. Peyrard, G. Russo, Ph. Villedieu. Polynomial upwind schemes for hyperbolic systems. C. R. Acad. Sci. Paris 328: 479–483, (1999). [9] B. Einfeldt. On Godunov-type methods for gas dynamics. SIAM J. Numer. Anal. 25: 294– 318, (1988). [10] E.D. Fern´andez-Nieto, G. Narbona-Reina. Extension of waf type methods to nonhomogeneous shallow water equations with pollutant, J. Sci. Comput. 36(2): 193–217, (2008). [11] A. Harten, P.D. Lax, B. van Leer B. On Upstream Differencing and Godunov Type Schemes for Hyperbolic Conservation Laws. SIAM Review, 25(1): 35–61, (1983). [12] N.E. Kolgan. Application of the minimum-derivative principle in the construction of nitedi?erence schemes for numerical analysis of discontinuous solutions in gas dynamics. Uchenye Zapiski TsaGI [Sci. Notes Central Inst. Aerodyn], 3(6): 68–77. 137, 138, (1972). [13] N.E. Kolgan. Finite-diference schemes for computation of three dimensional solutions of gas dynamics and calculation of a ow over a body under an angle of attack. Uchenye Zapiski TsaGI [Sci. Notes Central Inst. Aerodyn], 6(2):1-6. 137, 138, (1975). [14] C. Par´es, M.J. Castro. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to Shallow-Water Systems. M2AN, Vol. 38(5): 821–852, (2004). [15] C. Par´es. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J. Num. Anal. 44(1): 300–321, (2006). [16] C. Par´es and M.L. Mu˜noz Ru´ız. On some difficulties of the numerical approximation of nonconservative hyperbolic systems. Bolet´ın SEMA, 47:23–52, (2009). [17] Y. Loukili, A. Soulaymani, Numerical Tracking of Shallow Water Waves by the Unstructured Finite Volume WAF Approximation, Int. J. Comput. Meth. Eng. Sci. Mech. 8(2):75– 88, (2007). [18] P.L. Roe. Approximate Riemann solvers, parameter vectors and difference schemes. J. Comp. Phys., 43:357–371, (1981). [19] S. Serna. A characteristic-based nonconvex entropy-fix upwind scheme for the ideal magnetohydrodynamics equations. J. Comput. Phys. 228, 4232–4247, (2009). [20] J.B. Schijf and J.C. Schonfeld. Theoretical considerations on the motion of salt and fresh water. In Proc. of the Minn. Int. Hydraulics Conv., 321–333. Joint meeting IAHR and Hyd. Div. ASCE. (1953). [21] E.F. Toro. A Weighted Average Flux Method for Hyperbolic Conservation Laws. Proceedings of the Royal Society of London, Series A, Mathematical and Physical Sciences, 423:401–418, (1989). 35 [22] E.F. Toro, Riemann problems and the the WAF method for solving two-dimensional shallow water equations, Phil. Trans. Roy. Soc. London, A338, pp.43-68, (1992). [23] E.F. Toro, The weighted average flux method applied to the time dependent euler equations, Phil. Trans. Roy. Soc. London, A341: 499–530, (1992). [24] E. F. Toro, R. C. Millington, L.A.M. Nejad, Primitive upwind numerical methods for hyperbolic partial differential equations, in Bruneau, C. H. (ed.), Sixteenth International Conference on Numerical Methods for Fluids Dynamics. Lecture notes in physics, SpringerVerlag, Berlin, 421–426, (1998). [25] E.F. Toro, S.J. Billett. Centred TVD schemes for hyperbolic conservation laws. IMA Journal of Numerical Analysis (20): 47–79, (2000). [26] E.F. Toro: Shock-Capturing Methods for Free-Surface Shallow Flows, Wiley, England, (2001). [27] E.F. Toro, V.A. Titarev, TVD Fluxes for the High-Order ADER Schemes, J. Sci. Computing, (24).3: 285–309, (2005). [28] E.F. Toro, V.A. Titarev, ADER schemes for scalar hyperbolic conservation laws with source terms in three space dimensions, J. Com. Phys (202): 196–215, (2005). [29] E.F. Toro, V.A. Titarev. MUSTA fluxes for systems of conservation laws. J. Comput. Phys. 216(2): 403–429, (2006). [30] M. Torrilhon, D. S. Balsara. High order WENO schemes: investigations on non-uniform convergence for MHD Riemann problems. J. Comput. Phys. (201): 586–600, (2004). [31] I. Toumi. A weak formulation of Roe approximate Riemann solver. J. Comp. Phys. 102(2): 360–373, (1992). [32] Van Leer B. Toward the ultimate conservative difference scheme. III. Upstream-centered finite-difference schemes for ideal compressible flows. Journal of Computational Physics (7).23: 263–275, (1977). 36