An adaptive mesh method for 1D hyperbolic conservation laws

Accepted Manuscript An adaptive mesh method for 1D hyperbolic conservation laws Fuxing Hu, Rong Wang, Xueyong Chen, Hui Feng

PII: DOI: Reference:

S0168-9274(14)00209-8 10.1016/j.apnum.2014.10.008 APNUM 2898

To appear in:

Applied Numerical Mathematics

Received date: 19 November 2013 Revised date: 1 September 2014 Accepted date: 3 October 2014

Please cite this article in press as: F. Hu et al., An adaptive mesh method for 1D hyperbolic conservation laws, Applied Numerical Mathematics (2015), http://dx.doi.org/10.1016/j.apnum.2014.10.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An adaptive mesh method for 1D hyperbolic conservation laws Fuxing Hua,∗, Rong Wangb , Xueyong chenc , Hui Fengd a Department

of Mathematics, Huizhou University, Huizhou, Guangdong, 516007, P.R.China b Department of General Education, South University of Science and Technology of China, Shenzhen, Guangdong, 518055, P.R. China c Faculty of Mathematics & Statistics, Xuchang University, Xuchang, Henan, 461000, P.R.China d Faculty of Mathematics & Statistics, Wuhan University, Wuhan, Hubei, 430072, P.R.China

Abstract An adaptive method is developed for solving one-dimensional systems of hyperbolic conservation laws, which combines the rezoning approach with the finite volume weighted essentially non-oscillatory (WENO) scheme. An a posteriori error estimate, used to equidistribute the mesh, is obtained from the differences between respective numerical solutions of 5th-order WENO (WENO5) and 3rdorder ENO (ENO3) schemes. The number of grids can be adaptively readjusted based on the solution structure. For higher efficiency, mesh readjustment is performed every few time steps rather than every time step. In addition, a high order conservative interpolation is used to compute the physical solutions on the new mesh from old mesh based on the finite volume ENO reconstruction. Extensive examples suggest that this adaptive method exhibits more accurate resolution of discontinuities for a similar level of computational time comparing with that on a uniform mesh. Keywords: adaptive method, finite volume WENO methods, hyperbolic conservation laws 2010 MSC: 65M50, 76M12

∗ Corresponding

author Email addresses: [email protected] (Fuxing Hu), [email protected] (Rong Wang), [email protected] (Xueyong chen), [email protected] (Hui Feng)

Preprint submitted to Elsevier

December 18, 2014

1. Introduction Solutions of partial differential equations arising in science and engineering frequently have large variations occurring over small portion of the physical domain, and a major challenge when solving such problems is to resolve these portions with sufficient accuracy while also keeping the computational work within acceptable limits. The vast majority of numerical methods for solving hyperbolic problems have been developed for fixed, uniform mesh. Using a uniform mesh throughout the physical domain can become a substantial and computational expense, especially in multidimensions where the number of mesh points required can be prohibitively large. In principle therefore, it is desirable to employ a non-uniform mesh that is sparse in regions where the solution is smooth and more concentrated near large-variation portions [1]. Since the solutions are not stationary, it is attractive to allow the grids to move in time and fine grid resolution can be maintained near the regions with large solution variation, thereby attaining a balance between accuracy and efficiency. The significance and advantages of using high order schemes in numerical simulation have been demonstrated by many applications [2]. WENO schemes, one of most popular high order methods, are designed based on the successful ENO schemes [3, 4]. WENO schemes divided into two categories, finite difference and finite volume ones, have been constructed in the past decade [5, 6, 7]. It is worth mentioning that the finite volume WENO schemes can be applied to arbitrary nonuniform grids, while high order finite difference WENO schemes can keep only conservative property on uniform or smoothly grids. During the past few decades, adaptive mesh methods have been become another useful tool used to solve partial differential equations with large solution variations, such as shock waves, boundary layers and so on [8]. Roughly speaking, there are three types of mesh adaptive methods. The first one is hrefinement methods (Adaptive Mesh Refinement [9, 10, 11, 12]), in which mesh are adjusted at fixed time levels, where points can also be added to or removed from the mesh in order to obtain an approximate solution with the desired level of accuracy. specifically, mesh points are added in the regions where the solution variation or error is large, and mesh points are removed in the regions where the solution is smooth. The second type is p-refinement methods [13, 14, 15, 16]

2

with which polynomial order adaptively varies from place to place according to a certain error estimate or indicator. The third type is r-refinement methods [17, 18, 1, 19, 20], also called moving mesh methods, which relocate grid points in a mesh having a fixed number of nodes in such a way that the nodes remain concentrated in regions of rapid variation of the solution. In this paper, we study rezoning approach [1] for the numerical solution of conservation laws together with finite volume WENO5 scheme (the quasiLagrange approach combine with finite difference WENO schemes has been studied in [21]). With the rezoning approach, the mesh is updated in an intermittent way using certain error equidistribution principle. An a posteriori error estimate is computed by the differences between the numerical solutions of WENO5 and ENO3. This error estimate possess two purposes. the first purpose is the application for equidistributing mesh, and the second one is to predict the number of cells under current solution structure. A conservative interpolation is served to interpolate physical solutions on the new mesh from the old mesh. Since the procedure of computing the coefficients in the WENO scheme and interpolation is very time-consuming, the restriction on the frequency of mesh adjustment is devised to improve the efficiency. Rather than adjusting the mesh every time step, we do the mesh adjustment while the estimated error exceed a prescribed threshold. However, we find that this intermittent mesh adjustment is not suitable for the nonlinear systems of conservation laws. As we know, the post-shock oscillations will occur when the shock pass through the interface in the discontinuous grids for the nonlinear systems of conservation laws. Since the mesh need to be readjusted every time step, the solution of auxiliary scheme is replaced by the solution of WENO5 scheme at the beginning of every time step. Hence the life cycle of the auxiliary solution is only one time step. In order to decrease the computational cost of the adaptive method, the auxiliary scheme is chosen to be the modified version of ENO3 scheme and coupled with the first order Euler forward time discretization. Specifically, the smoothest stencil in the modified ENO3 scheme is measured by the weights or smoothness indicators which have been computed in the procedure of implementing the WENO5 scheme. Another important component of adaptive mesh method is how to construct a high order conservative interpolation. The high order conservative interpola3

tion based on cell averages is relatively simpler than that on point values. The finite volume ENO reconstruction can be used to construct high order polynomial which is both globally conservative and essentially non-oscillatory. To keep the adaptive method maintain 5th-order accuracy, a 5th-order ENO conservative interpolation is used here to preserves the solution quantities. The existing adaptive mesh methods are almost integrated with fixed number of mesh points. Since the solution structure of hyperbolic conservation laws can largely change with the time, the number of mesh points also should be adaptively adjusted depend on the metamorphic solutions. Hence, a rough forecast algorithm in grid number is presented here, which is firstly devised by Wang et al. [22] for 1-D parabolic PDEs. In addition, it has advantage that no sensitive parameter need to be set in this adaptive method. We start in Sec. 2 with a description of the WENO methods on non-uniform mesh. Sec. 3 and 4 present the detailed adaptive algorithm to resolve the discontinuous solutions and gives a high order conservative interpolation and forecast system. In Sec. 5, numerical tests are performed on several problems encompassing both scalar conservation laws and systems. The accuracy, efficiency and robustness of the method are demonstrated here, and no numerical instability is found under abrupt mesh movement.

2. WENO formulations on non-uniform mesh To illustrate, we consider the 1-D, scalar hyperbolic conservation law ut + fx (u) = 0,

a ≤ x ≤ b,

t ≥ 0,

(2.1)

on a non-uniform mesh :

a = x1/2 < x3/2 < · · · < xN −1/2 < xN +1/2 = b.

(2.2)

The cells Ii , cell centers xi and length of cells Δxi are respectively defined by Ii = [xi−1/2 , xi+1/2 ],

xi = (xi−1/2 + xi+1/2 )/2,

4

Δxi = xi+1/2 − xi−1/2 ,

with i = 1, · · · , N . Integrating (2.1) over Ii we obtain  1  du(xi , t) =− f (u(xi+1/2 , t)) − f (u(xi−1/2 , t)) , dt Δxi

(2.3)

where ui (t) denotes the cell average of u(x, t) on Ii , i.e., 1 u(xi , t) = Δxi



xi+1/2

u(ξ, t)dξ. xi−1/2

We can approximate the conservative scheme (2.3) by 1 ˆ dui (t) =− (fi+1/2 − fˆi−1/2 ), dt Δxi

(2.4)

where ui (t) is the numerical approximation to u(xi , t), and the numerical flux + fˆi+1/2 = G(u− i+1/2 , ui+1/2 )

is the approximation to f (u(xi+1/2 , t)). In addition, as shown in Fig. 1, u± i+1/2 are obtained from the WENO5 reconstruction procedure,

(i)

xi+1/2 -

(i)

- 

(i)

- 

(i) P2 (xi+ 12 ), ω ˆ2

u− i+ 1 2

(i) P1 (xi+ 12 ), ω ˆ1 (i) P0 (xi+ 12 ), ω ˆ0



Fig. 1.

(i+1)

(xi+ 12 )

(i+1)

(xi+ 12 )

(i+1)

, P2

(i+1)

, P1

ω2 ω1

(i+1)

ω0

(i+1)

, P0

(xi+ 12 )

The procedure of WENO5 reconstruction

u+ i+1/2 =

2 

ωr(i+1) Pr(i+1) (xi+1/2 )

r=0

5

(2.5)

u+ i+ 1 2

and u− i+1/2 =

2 

ω ˆ r(i) Pr(i) (xi+1/2 ),

(2.6)

r=0

(i)

where Pr (x) is quadratic polynomial, defined in each stencil, which is third(i)

order accurate approximation to the exact solution u(x, t) on Ii . ω ˆr weight of polynomial (i) P2 (xi+1/2 ),

(i) Pr (x)

(i) P1 (xi+1/2 )

is the

at x = xi+1/2 . Through convex combination of

and P0 (xi+1/2 ), u− i+1/2 can achieve the 5th-order (i)

accuracy to u(xi+1/2 , t) in smooth regions while keep ENO property [23, 4] around large-variation regions. Similarly, u+ i+1/2 can be constructed by the (i+1)

similar way. We omit the detailed formulae of ωr (i) Pr (xi+1/2 )

(i)

(i+1)

,ω ˆ r , Pr

(xi+1/2 ) and

on non-uniform mesh for saving space, and the interested reader

refer to [24]. In the end, we emphasize that using WENO schemes on nonuniform mesh is computational cost, but they can give high-resolution solutions.

3. Moving mesh methods The moving mesh methods can be divided into five components: the error estimations, the strategy employed to move the mesh, the method used to estimate the number of new cells, the interpolation used to obtain the physical solutions on the new mesh, and the approach served to integrate the physical equations (which is given in Sec. 2). 3.1. Error estimate Denote the cell average of WENO5 and ENO3 scheme on cell Ii by Ui (t) = (u1i (t), · · ·

P , uN (t))T and Vi (t) = (vi1 (t), · · · , viN P (t))T , respectively, where N P i

is the number of PDEs. Two sets of error estimation are considered; the first set is for the error on each cell,

Ei = Δxi

NP  j=1

|uji (t) − vij (t)|

AT OL + RT OL|uji (t)|

,

i = 1, · · · , N,

(3.1)

and the second set is for the error of each equation on the physical domain, Ej =

N  i=1

Δxi

|uji (t) − vij (t)|

AT OL + RT OL|uji (t)| 6

,

j = 1, · · · , N P,

(3.2)

where AT OL and RT OL are the absolute and relative tolerances respectively. The absolute tolerance AT OL, a threshold error value, represents the acceptable error as the value of the measured state approaches zero. The relative tolerance RT OL measures the error relative to the size of each state and represents a percentage of the state’s value. These two error estimation will be used to distribute the grid points and forecast the number of cells. 3.2. Redistribution of mesh Here, the method used to relocate the mesh is similar to de Boor’s algorithm [25], which is a simple yet useful approximation method for finding an equidistributing mesh. Assume that the errors (3.1)(3.2) and the numerical solutions is available on an arbitrary background mesh (2.2). The idea of de Boor’s algorithm is to find a new mesh ˆ: 

ˆ3/2 < · · · < x ˆN −1/2 < x ˆN +1/2 = b, a=x ˆ1/2 < x

(3.3)

which satisfies the following conditions 

N P x ˆi+1/2 x ˆi−1/2

1 p

(Ei ) dx =

1

(E j ) p

j=1

= C,

N

i = 1, . . . , N

(3.4)

where p is chosen to be the order of ENO3 scheme. The formula (3.4) is often referred to equidistribution principle (EP). Since uji (t) is the cell average of uj (x, t) on Ii , we set uj (x, t) = uji (t),

f or x ∈ [xi−1/2 , xi+1/2 ],

E(x, t) = Ei (t),

f or x ∈ [xi−1/2 , xi+1/2 ],

and similarly,

i.e., both uj (x, t) and E(x, t) are piecewise constant functions for variable x. In order to find the new grid points x ˆk+1/2 , one first determines on which

7

cell x ˆk+1/2 should locate in the old mesh. If l−1  i=0

1

Eip ≤ k · C ≤

l  i=0

1

Eip ,

k = 1, · · · N − 1

then x ˆk+1/2 lies on the cell Il = [xl−1/2 , xl+1/2 ]. Since E(x, t) is piecewise constant, x ˆk+1/2 can be directly calculated by k·C − x ˆk+1/2 = xl−1/2 +

l−1  i=0

C

1

Eip .

ˆN +1/2 = b are fixed all the time. Note that the boundary points x ˆ1/2 = a and x 3.3. Estimation of cell number For boundary value ODE problems, one has the following result [26]. When the number of cells N → ∞, using a quasiuniform mesh, a numerical method of order p gives a O(N −p ) error, |U (x) − u(x)| = O(N −p ), where U (x) is the numerical solution of the pth-order methods and u(x) is the exact solution of ODE. Here, we assume that a similar property holds for the piecewise error E j in (3.2) of hyperbolic conservation laws (2.1), E≡

max E j ∝ N −p .

1≤j≤N P

(3.5)

The equation (3.5) will be directly used to calculate the number of cells, and numerical experiments suggest that the cross-regional application seems reasonable. 3.4. Interpolation procedure Conservative interpolation from old mesh to new mesh plays an important role in rezoning approach. We use Gaussian quadrature formulas to compute the cell averages on the new meshes. The fifth-order finite volume ENO (ENO5) reconstruction [3, 4] is used to calculate the values at the Gaussian points. Two reasons lead us to the ENO5 reconstruction instead of WENO5. Firstly, for

8

different Gaussian points, the formulas of the weights needed in the WENO5 reconstruction scheme are cumbersome and distinct from one another because it is a nonuniform mesh. On the other hand, the computation cost of the ENO5 reconstruction is much lower than the one of WENO5. Secondly, negative weights may appear in the interior Gaussian integral points [27, 28]. Thus spurious oscillations may appear near the shock if WENO5 reconstruction is used. Similar observations can be found in [27, 28]. We now give the ENO5 interpolation procedure for 1D system. ˆj on some new interval Iˆj = [ˆ xj−1/2 , x ˆj+1/2 ] Suppose that the cell average U is wanted, and Ui on the old mesh Ii = [xi−1/2 , xi+1/2 ],

1≤i≤N

is known. Fig. 2 gives the local mesh distribution.

xi−1/2

xi+1/2

x ˆj−1/2

Fig. 2.

xi+3/2

x ˆj+1/2

Local mesh distribution

Denote the fourth-order piecewise polynomial on cell Ii by P (i) (x) which is obtained by ENO reconstruction. With the idea of ENO reconstruction, there are five stencils can be chosen, i.e., Sir = {Ii+r−4 , Ii+r−3 , Ii+r−2 , Ii+r−1 , Ii+r },

0≤r≤4

and r means the number of cell at the right of Ii . Through divided differences in these five different stencils, we can determined the most smoothest stencil which is used to reconstruct the interpolation polynomial on the cell Ii . As shown in Fig. 2, suppose the piecewise polynomial P (i) (x) obtained by ENO ˆj on new cell reconstruction is available on each cell Ii , the new cell average U

9

Iˆj is the integral of the form ⎛ ⎞  x 1  xˆ 1 i+ j+ 1 2 2 (i) (i+1) ˆj = ⎝ P (x)dx + P (x)dx⎠ , U Δˆ xj x ˆ 1 x 1 j−

i+

2

(3.6)

2

where Δˆ xj = [ˆ xj−1/2 , x ˆj+1/2 ]. Since the integrants in (3.6) is polynomial, these two integrals can be calculated exactly by Gaussian integral formula. Note that ˆj− 12 ] has been computed if the integral of function P (i) (x) on interval [xi− 12 , x xj− 12 , xi+ 12 ] can be in advance, then the integral of function P (i) (x) on interval [ˆ achieved directly by 

xi+ 1 2

x ˆj− 1

 P (i) (x)dx = Δxi Ui −

2

x ˆj− 1

2

P (i) (x)dx.

xi− 1 2

Obviously, the interpolation procedure used here is globally conservative. 3.5. Summary In summary, our adaptive-mesh method has several major contributions as follows. • We use a high-order (actually 5th-order WENO) scheme for the spacial discretization coupled with adaptive-mesh approach. An a posteriori error estimate is obtained from the differences between respective numerical solutions of WENO5 and ENO3 schemes. • The number of cells are evaluated by the error estimate. Then the equidistribution principle is applied to obtain a new mesh. • A high order conservative interpolation (based on ENO5 reconstruction) is used to compute the cell averages on the new meshes from old meshes.

4. Algorithm procedure Assume at time tn the mesh n = {xni } and the physical solutions Uin = U (xi , tn ) on each cell Ii of n both are available. Then we need to redistribute the mesh with the errors Ei and E j . The algorithm procedure is presented detailed as follows.

10

step 1. Firstly, let auxiliary solutions (numerical solutions of ENO3) Vin := Uin for i = 1, · · · , N , then evolve the solutions Uin and Vin one time step on fixed mesh n with WENO5 and ENO3 methods, respectively. Denote these evolved solutions at tn+1 by Uin+1 and Vin+1 . Using formulae (3.1) and (3.2), the errors Ei for each cell and E j for each equation can be computed easily. step 2. If E ≡

max E j < 1, then tn := tn+1 , Uin := Uin+1 and repeat step

1≤j≤N P

1. Else go to step 3. step 3. Using the error Ei and E j obtained in step 1 to redistribute the mesh ˆ n has been gained by the redisn at time tn . Suppose the new mesh  ˆ (ˆ ˆn = U x i , tn ) tribution procedure in Sec. 3.2, then the new cell average U i can be get with the conservative interpolation in Sec. 3.4. Until now, the new mesh and cell averages have been acquired at time tn . ˆ n and repeat the step 1, we get the error Ei and E j . If the step 4. Let Uin = U i error E = 1. If E =

max E j < 1, then tn := tn+1 , Uin := Uin+1 and repeat step

1≤j≤N P

max E j ≥ 1 as before, then we say that the number of cells

1≤j≤N P

N is not enough to resolve the physical solution. A system of evaluation ˆ , which is given in Sec. 3.3. is used to predict the new number of cells N ˆ should be corresponding to step 5. Since E ∝ N −p , the new cell number N a relative small error. Generally, we require this small error computed ˆ cells lies in a suitable interval [e1 , e2 ]. Here, we choose e1 = 0.1 with N and e2 = 0.3, which are the empirical values according to the numerical ˆ satisfy experiments. So we expect the new errors E ˆ −p ∝ E ˆ = 0.2. N

(4.1)

Comparing the equation (3.5) and (4.1), we have  ˆ =N N

E 0.2

p1 .

(4.2)

ˆ obtained in step 5, we can get the cell step 6. Using the new cell number N ˆ n for i = 1, · · · , N ˆ with conservative interpolation as described average U i

11

ˆ , U n := U ˆ n , and repeat step 1. If E = in Sec. 3.4. Let N := N

max E j

1≤j≤N P

lies in the interval [e1 , e2 ], then we say that the mesh is equidistributed ˆ n+1 , and N := N ˆ and repeat step reasonably. And let tn := tn+1 , U n := U 1. Else if E =

max E j lies out of the interval [e1 , e2 ], we go to step 5

1≤j≤N P

unless E is inside this interval. The following issues should be noted: The errors Ei and E j present in step 1 is computed with the cell averages ˆ n+1 at time tn+1 , and so no time-lag is introduced in the mesh Uin+1 and U i movement. The empirical values e1 = 0.1 and e2 = 0.3 are specified for the following reason. If E < 0.1, we say that the reason for E achieving so small is that the number of cells is too many. In contrast, If E > 0.3, we say that the reason for E achieving so big is that the number of cells is too little. In addition, two restriction are placed on the mesh generation procedure. One of the restriction is that the ratio R of the length of the adjacent cells satisfies the condition 0.5 ≤ R ≤ 2. The object of this restriction is to smooth the mesh and keep the stability of the move mesh methods. Another restriction is to make the generated cell not too small. Small cell will lead to small time-step length. Generally, we don’t want the length of the smallest cell less than the absolute tolerance AT OL. A practical handling of this difficulty is specified as follows. If the length of the cell Ii is smaller than AT OL, this cell will be merged with one of its adjoined cells that is less than the other; i.e., Ii is merged into Ii−1 if Δxi−1 ≤ Δxi+1 , and vice versa.

5. Numerical results In this section some classical 1D test examples are used to examine the proposed moving mesh method. In the following experiments, the ODEs resulting from the semi-discretized PDEs are evolved in time by the third-order total

12

variation diminishing Runge-Kutta scheme [29]:

U1

=

Un + ΔtL(Un ),

U2

=

3 4 Un

+ 14 U1 + 14 ΔtL(U1 ),

Un+1

=

1 3 Un

+ 23 U2 + 23 ΔtL(U2 ),

(5.1)

where Δt = CF L · min{Δxi , i = 1, . . . , N }, and CF L number is set to be 0.1. The error shown in the following numerical results is the l1 -error,

l1 =

N 

|Ui − U (xi )|Δx.

i=1

For all the examples, we let the relative error tolerance RT OL = AT OL. Example 1. (Linear advection equation) In this subsection, we present numerical experiments of the advection equation, ut + ux = 0

(5.2)

with initial conditions ⎧ ⎪ ⎪ ⎪ 1, ⎪ ⎪ ⎨ u(x, 0) = 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1,

if 0 ≤ x < 0.5, if 0.5 ≤ x < 1.5, if 1.5 ≤ x < 2,

⎧ ⎪ ⎪ ⎪ 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ 1,

and

u(x, 0) =

(5.3)

if 0 ≤ x < 0.25, if 0.25 ≤ x < 0.75,

⎪ ⎪ ⎪ 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ sin(2πx),

(5.4)

if 0.75 ≤ x < 1, if 1 ≤ x < 2.

Here the user-specified number of initial cells is chose to be N = 100. For the initial-boundary-value problem (IBV) (5.2)(5.3), the initial condition contains two discontinuities and they propagate to the right with speed s = 1. Fig. 3 shows the numerical solution with AT OL = 10−3 at time t = 0.1 and we can find that the grid points concentrate around the discontinuities with high 13

1.2

1

0.8

0.6

0.4

0.2

0

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

The numerical solution of problem (5.2)(5.3) with AT OL = 10−3 at time t = 0.1. The result number of cells required for AT OL = 10−3 is about N = 110.

Fig. 3.

resolution. Table 1 gives the l1 -error and number of cells for different choice of AT OL. Note that the number of cells given in the table is just a approximate value. It means that the number of cells used at each time step is around this value. In this test, the main reason of the number of cells lingers around some constant is that the structure of solution is not change essentially. In the latter examples we will demonstrate the solution evolved dynamically and resolved with proper grid number adaptively. Table 1.

The l1 -error and number of cells required is shown with different AT OLs.

AT OL

l1 -error

number of cells

1.0e-03 1.0e-04 1.0e-05

2.92e-02 7.11e-03

138 151

3.59e-03 8.91e-04

156 205

1.96e-04

242

1.0e-06 1.0e-07

Fig. 4 gives the comparison of efficiency between adaptive method and WENO5 method on uniform mesh. Although these two methods have almost identical efficiency if relative large error is required, the adaptive method has 14

−1

10

data1 data2 data3

−2

error

10

−3

10

−4

10

1

2

10

3

10

4

10

5

10

6

10

10

CPU timing

Fig. 4. The l1 -error in the numerical solutions of (5.2)(5.3) at time t = 0.1 is plotted as a function of CPU timing (in millisecond).

higher efficiency if high-resolution solution is required. 1.5 1.3 1

1.2 1.1

0.5 1 0

0.9 0.8

−0.5 0.7

−1

0.6 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.1

1.2

1.3

1.4

1.5

1.6

Two numerical solutions are shown at t = 0.1. The solution denoted by ’∗’ is obtained with the restriction M LC = 0.05 and another solution denoted by ’+’ is obtained without the restriction.

Fig. 5.

The initial condition (5.4) has two discontinuities and a very smooth sine function sin(2πx). The IBV problem (5.2)(5.4) is used to introduce another restriction applied on mesh generation. From Fig. 5 we can find that the main difference between these two numerical solutions locates on the sine function. As shown in the right of the Fig. 5, the solution marked by ’+’ is not well fitting the sine curve. The presentation of the phenomenon is due to that too 15

little cells are distributed on this area. The solution marked by ’∗’ is computed with the restriction: the maximum length of the cell should not too large. Denote by M LC the maximum tolerant length of the cell. Fig. 5 shows the numerical solution with M LC = 0.05 can allocate adequate cells on the sine curve with higher resolution. The numerical solution of (5.2)(5.3) also applied this restriction. In the following tests M LC = 0.05 will be chosen to be a default value except specified. Example 2. (Burgers’ equation) The inviscid Burgers’ equation is an example of scalar nonlinear hyperbolic PDE, given by ut + (

u2 )x = 0 2

(5.5)

with initial conditions u(x, 0) = 0.5 + sin(πx) and

⎧ ⎪ ⎪ ⎪ 1, if 0 ≤ x < 0.5, ⎪ ⎪ ⎨ u(x, 0) = 0, if 0.5 ≤ x < 1, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ −4, if 1 ≤ x < 2.

(5.6)

(5.7)

The problem (5.5)(5.6) often serves as a standard test for adaptive moving methods. For comparison, we plot (the solid line in the Fig. 6) a reference solution using 1000 uniform cells. The initial condition (5.7) is devised to show that the number of cells can change adaptively according to the structure of solution. Initially, there are two shocks with speed s1 = 0.5 and s2 = −2 respectively, and they collide together at time t = 0.2. Fig. 7 gives the numerical solutions at several different times with AT OL = 10−3 . Fig. 8 shows the mesh distribution and the number of cells at some prescribed times. The amount of cells required for AT OL = 10−3 before two shocks colliding together is around N = 105. However, when they merge into one shock, the amount of cells decreases to around N = 75. The phenomenon occurred here is justified since the structure of solution become much simpler when the time t ≥ 0.2. This adjustment in the number of cells is owe to the prediction procedure in Sec. 3.3.

16

2

1.5

1

0.5

0

−0.5

−1

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 6. The numerical solution of problem (5.5)(5.6) with AT OL = 10−3 at times t = 0, 0.2, 0.4, 0.6, 0.8, 1.0. The result number of cells required for AT OL = 10−3 is about N = 85.

Example 3. (Euler equations) In this subsection we consider Euler equations ⎛

⎞

⎞

⎛

ρu ⎜ ρ ⎟ ⎜ ⎟ ⎜ ⎜ ⎟ ⎜ ⎜ ⎜ ρu ⎟ + ⎜ ρu2 + p ⎟ ⎜ ⎜ ⎠ ⎝ ⎝ E u(E + p) t

⎟ ⎟ ⎟ ⎟ = 0, ⎟ ⎠

(5.8)

x

with the following four initial conditions,

(ρ, u, p) =

⎧ ⎪ ⎨ (0.445, 0.698, 3.528),

if 0 ≤ x ≤ 0.5,

⎪ ⎩ (0.5, 0, 0.571),

if 0.5 ≤ x ≤ 1,

(ρ, v, p) =

⎧ ⎪ ⎨ (1, 0, 1),

if 0 ≤ x ≤ 0.5,

⎪ ⎩ (0.125, 0, 0.1),

if 0.5 ≤ x ≤ 1,

⎧ ⎪ ⎪ ⎪ ⎪ (1, 0, 1000), if 0 ≤ x < 0.1, ⎪ ⎨ (ρ, u, p) = (1, 0, 0.01), if 0.1 ≤ x < 0.9, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1, 0, 100), if 0.9 ≤ x ≤ 1,

17

(5.9)

(5.10)

(5.11)

2

2

1

1

0

0

−1

−1

−2

−2

−3

−3

−4

−4

−5

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

−5

2

2

2

1

1

0

0

−1

−1

−2

−2

−3

−3

−4

−4

−5

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

−5

2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 7. The numerical solutions of the initial-value problem of (5.5)(5.7) is shown at times t = 0, 0.1, 0.2, 0.3. Two shocks devised in the initialization collide together when t = 0.2. and (ρ, v, p) =

⎧   √ ⎪ ⎨ 27 , 4 35 , 31 ,

if 0 ≤ x < 0.1,

⎪ ⎩ (1 + 0.2 sin 50x, 0, 1),

if 0.1 ≤ x ≤ 1.

7

9

3

(5.12)

Firstly, we consider the Lax’s shock tube problem (5.8)(5.9) which results in a fast moving shock wave to the right ahead of a slower moving contact discontinuity and a rarefaction wave to the left of the initial discontinuity [30]. Fig. 9 gives the numerical solution of Lax problem computed with adaptive method in an intermittent way. Some serious postshock oscillations appear behind the shock and these oscillations cause that a great deal of grids are distributed here since oscillations necessarily introduce large numerical errors. The reason for these oscillations we think is the shock crossing the grid interfaces. For illustrative purpose, Fig. 10 shows how the postshock oscillations can form when the shock passes through the grid interface. To overcome this defect, we change the algorithm which readjusts the mesh only when the total error E > 1 into the

18

0.35 0.3

N=63

0.25

N=74

0.2

N=80

0.15

N=100

0.1

N=109

0.05

N=105 N=98

0 −0.05

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Fig. 8. The number of cells required for AT OL = 10−3 is shown at times t = 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3.

algorithm which readjusts the mesh every time step. 1.6 1.36 1.4 1.34 1.2 1.32 1 1.3

0.8

0.6

1.28

0.4

1.26

0.2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.68

0.7

0.72

0.74

0.76

0.78

0.8

0.82

0.84

Fig. 9.

The left is the numerical solution of Lax’s shock tube problem and the right is the zoomed version between the shock wave and contact discontinuity.

With the modified algorithm, we need redistribute the mesh every time step, which including distribution of grid points, interpolation of physical solutions on the new mesh and computation of the coefficients of WENO scheme. So the efficiency of the modified algorithm should have lower efficiency than original one. For comparison of efficiency for the original and modified algorithm, we still use the problem (5.2)(5.3) for testing. As shown in Fig. 11, it is surprised

19

1.6 1.33 1.4 1.32 1.2

1.31

1

1.3 1.29

0.8

1.28 0.6 1.27 0.4

0.2

1.26

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.72

0.74

0.76

0.78

0.8

0.82

0.84

0.86

0.88

0.9

Fig. 10. The physical domain [0, 1] is divided into two subdomain [0, 0.8] and [0.8, 1]. Then these two subdomain are both subdivided into uniform meshes but with different cell length. Two conspicuous oscillations can be observed behind the shock. that the modified algorithm has slightly higher efficiency than original algorithm. The reasons for this we think may be owe to that the modified algorithm can repress the error at each time step, yet the original algorithm introduces some errors when the mesh keep fixed. In addition, all the previous testes are computed with modified algorithm but the numerical results in Table 1 and Fig. 4. In the end, it is noted that two algorithm have similar results for scalar conservation laws and linear systems of conservation laws, but they have different performance for nonlinear systems of conservation laws. The second problem (5.8)(5.10) is Sod’s shock tube problem [31] which has similar solution with Sod’s shock tube problem. Fig. 12 shows these two Riemann problems with different absolute tolerance AT OL. The numerical solutions on the left are the Sod problem and the ones on the right are the Lax problem. The top two figures are computed with AT OL = 10−2 and the numerical solutions can not well enough to simulate the exact solutions. The medial two computed with AT OL = 10−3 arrive a significant improvement compared with the top two numerical solutions. And the bottom two solutions computed with AT OL = 10−4 are comparable with the medial two but with more resolution around the contact discontinuities. Generally, we except that the numerical solution will be more accurate when less AT OL is used. However, in some special cases, the excessively small AT OL will lead to bad solutions (refer to the next test). The problem (5.8)(5.11) is first used as a test problem by Woodward and

20

−1

10

original algorithm modified algorithm uniform mesh

−2

error

10

−3

10

−4

10

1

10

2

10

3

4

10

10

5

10

6

10

CPU timing

Fig. 11. The comparison of efficiency between the original and modified algorithm with the problem (5.2)(5.3).

Colella [32] and is often referred to as the Woodward-Colella blast-wave problem. The two discontinuities in the initial condition each have the form of a shocktube problem and yield strong shock wave and contact discontinuities going inwards and rarefaction waves going outwards. The boundaries at x = 0 and x = 1 are both reflecting walls and the reflected rarefaction waves interact with the other waves. The numerical results obtained with AT OL = 10−3 and AT OL = 10−4 , respectively, are presented in Fig. 13. The numerical solution computed with AT OL = 10−4 has some trouble at the tip around x = 0.78. It is unclear for the reason of this phenomenon and quite a number of interpolations from old mesh to new mesh may be the reason leading to the trouble. The last example (5.8)(5.12) is the Shu-Osher problem [33], which contains both shocks and complicated smooth regions structures. Because of the presentation of many complex smooth regions, two choice of M LC are shown in Fig. 14. It can be found that the numerical solution get some improvement through choosing less M LC, especially around the region [0.2, 0.3].

21

1

1.6

0.9

1.4

0.8 1.2 0.7 0.6

1

0.5

0.8

0.4 0.6 0.3 0.4

0.2 0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.2

1

1

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

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.6

0.9

1.4

0.8 1.2 0.7 0.6

1

0.5

0.8

0.4 0.6 0.3 0.4

0.2 0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.2

1

1

1.6

0.9

1.4

0.8 1.2 0.7 0.6

1

0.5

0.8

0.4 0.6 0.3 0.4

0.2 0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.2

1

Fig. 12. The top two numerical solutions are computed with AT OL = 10−2 . The medial two numerical solutions are computed with AT OL = 10−3 and the bottom two numerical solutions are computed with AT OL = 10−4 . 6. Conclusions An adaptive mesh method combined with the fifth-order WENO scheme is presented in this paper. Here, we directly applied the 5th-order WENO scheme on the non-uniform mesh. Since the solution structure of hyperbolic equations may evolve very complicated, the number of cells should be measured naturally by the solution structure. In order to change the number of cells adaptively 22

7

7

6

6

5

5

4

4

3

3

2

2

1

1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

the left numerical solution is computed with AT OL = 10−3 and the right is computed with AT OL = 10−4

Fig. 13.

5

5

4.5

4.5

4

4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

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

the left numerical solution is computed with M LC = 0.05 and the right is computed with M LC = 0.005

Fig. 14.

according to the structure of solutions, a prediction system is proposed based on an a posteriori error estimate. In addition, a high order conservative interpolation is used to compute the cell averages on the new meshes from old meshes. Numerical tests demonstrate that, when using the 5th-order WENO method, the adaptive-mesh method is much more efficient than the uniformmesh one. Finally, this method can be extended directly to two dimensions if a tensor-product mesh is used. Acknowledgment This paper is supported by National Natural Science Foundation of China(No. 11272277, No. 11301455, No. 11326096), Foundation for Distinguished Young Talents in Higher Education of Guangdong, China (No. 2013LYM 0089), and Doctoral Scientific Research Project of Huizhou University(C5120201).

23

References [1] W. Huang, R. D. Russell, Adaptive Moving Mesh Methods, SpringerVerlag, New York, 2011. [2] C. W. Shu, High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD, International Journal of Computational Fluid Dynamics 17 (2003) 107–118. [3] A. Harten, S. Osher, Uniformly high order Essentially Non-Oscillatory schemes I, SIAM J. Numer. Anal. 24 (1987) 279–309. [4] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order Essentially Non-Oscillatory schemes III, J. Comput. Phys. 71 (1987) 231– 303. [5] X.-D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [6] G.-S. Jiang, C.-W. Shu, Efficient implementation of Weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [7] C.-W. Shu, Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws, in: Advanced numerical approximation of nonlinear hyperbolic equations, Lecture Notes in Mathematics, Vol. 1697, Springer-Verlag, Berlin, 1998, pp. 325–432. [8] D. F. Hawken, J. J. Gottlieb, J. S. Hansen, Review of some adaptive nodemovement techniques in finite element and finite difference solutions of PDEs, J. Comput. Phys. 95 (1991) 254–302. [9] M. Berger, Adaptive mesh refinement for hyperbolic partial differential equations, Ph.D. thesis, Stanford university, stanford, CA, 1982. [10] M. J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (1984) 484–512. [11] M. J. Berger, P. Colella, Local adaptive mesh refinement for shock hydrodynamics, J. Comput. Phys. 82 (1989) 64–84.

24

[12] S. McCormick, Multilevel Adaptive Methods for Partial Differential Equations, Frontiers in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, 1989. [13] W. Gui, I. Babuˇ ska, The h, p and h-p versions of the finite element method in 1 dimension. Part I. the error analysis of the p-version, Numer. Math. 49 (1986) 577–612. [14] W. Gui, I. Babuˇ ska, The h, p and h-p versions of the finite element method in 1 dimension. Part II. the error analysis of the h-and h-p version, Numer. Math. 49 (1986) 613–657. [15] W. Gui, I. Babuˇ ska, The h, p and h-p versions of the finite element method in 1 dimension. Part III. the the adaptive h-p version, Numer. Math. 49 (1986) 577–612. [16] S. Adjerid, J. E. Flaherty, P. K. Moore, Y. Wang, High-order adaptive methods for parabolic systems, Physica D 60 (1992) 94–111. [17] K. Miller, R. N. Miller, Moving finite elements. I, SIAM J. Numer. Anal. 18 (1981) 1019–1032. [18] K. Miller, Moving finite elements. II, SIAM J. Numer. Anal. 18 (1981) 1033–1057. [19] W. Huang, Y. Ren, R. D. Russell, Moving mesh partial differential equations (MMPDES) based on the equidistribution principle, SIAM J. Numer. Anal. 31 (1994) 709–730. [20] H. Z. Tang, T. Tang, Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws, SIAM J. Numer. Anal 41 (2003) 487–515. [21] X. Yang, W. Huang, J. Qiu, A moving mesh WENO method for onedimensional conservation laws, SIAM J. Sci. Comput. 34 (2012) 2317–2343. [22] R. Wang, P. Keast, P. Muir, A high-order global spatially adaptive collocation method for 1-D parabolic PDEs, Appl. Numer. Math. 50 (2004) 239–260.

25

[23] A. Harten, S. Osher, B. Engquist, S. Chakravarthy, Some results on uniformly high order accurate essentially Non-Oscillatory schemes, Appl. Numer. Math. 2 (1986) 347–377. [24] R. Wang, H. Feng, R. J. Spiteri, Observations on the fifth-order WENO method with non-uniform meshes, Applied Math. and Comp. 196 (2008) 433–447. [25] C. de Boor, Good approximation by splines with variable knots, in: A. Meir and A. Sharma, editors, Spline Functions and Approximation Thoery, Birkh¨ auser verlag, Basel und Stuttgart, 1973, pp. 57–73. [26] U. Ascher, R. M. M. Mattheij, R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Diferential Equations, SIAM, Philadelphia, PA, 1995. [27] K. Sebastian, C.-W. Shu, Multidomain WENO finite difference method with interpolation at subdomain interfaces, J. Sci. Comput. 19 (2003) 405– 438. [28] J. shi, C. Hu, C.-W. Shu, A technique of treating negative weights in WENO schemes, J. Comput. Phys. 175 (2002) 108–127. [29] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comp. 67 (1998) 73–85. [30] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Commun. Pure. Appl. Math. 7 (1954) 159–193. [31] G. A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. comput. phys. 107 (1978) 1–31. [32] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. comput. phys. 54 (1984) 115–173. [33] C.-W. Shu, S. Osher, Efficient implementation of Essentially NonOscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78.

26