A Hierarchical Multiresolution Adaptive Mesh Refinement for ... - DCSL

c 2008 Society for Industrial and Applied Mathematics

SIAM J. SCI. COMPUT. Vol. 31, No. 2, pp. 1221–1248

A HIERARCHICAL MULTIRESOLUTION ADAPTIVE MESH REFINEMENT FOR THE SOLUTION OF EVOLUTION PDEs∗ S. JAIN† , P. TSIOTRAS† , AND H.-M. ZHOU‡ Abstract. In this paper, we propose a novel multiresolution adaptive mesh refinement algorithm for solving initial-boundary value problems (IBVP) for evolution PDEs. The proposed algorithm dynamically adapts the grid to any existing or emerging irregularities in the solution, by refining the grid only at those places where the solution exhibits sharp features. The main advantage of the proposed grid adaptation method is that it results in a grid with a fewer number of nodes when compared to adaptive grids generated by existing multiresolution-based mesh refinement techniques. Several examples show the robustness and stability of the proposed algorithm. Key words. mesh refinement, multiresolution, adaptive PDE methods, ENO, weighted ENO AMS subject classifications. 65M50, 65N50, 65M06, 65D05 DOI. 10.1137/070708329

1. Introduction. It is well known that the solution of evolution PDEs is often not smooth even if the initial data are smooth. For instance, shocks may develop in hyperbolic conservation laws. To capture discontinuities in the solution with high accuracy, one needs to use a fine resolution grid. The use of a uniformly fine grid requires a large amount of computational resources in terms of both CPU time and memory. Hence, in order to solve evolution equations in a computationally efficient manner, the grid should adapt dynamically to reflect local changes in the solution. Several adaptive gridding techniques for solving PDEs have been proposed in the literature. A nice survey of the early works on the subject can be found in [5, 47]. Currently, popular adaptive methods for solving PDEs are (i) moving mesh methods [2, 1, 4, 7, 6, 13, 14, 17, 29, 32, 35, 36, 45], in which an equation is derived that moves a grid of a fixed number of finite difference cells or finite elements so as to follow and resolve any local irregularities in the solution; (ii) the so-called adaptive mesh refinement method [8, 9, 10, 11], in which the mesh is refined locally based on the difference between the solutions computed on the coarser and the finer grids, and (iii) wavelet-based or multiresolution-based methods [3, 12, 15, 18, 19, 20, 23, 24, 26, 37, 38, 48, 49, 50], which take advantage of the fact that functions with localized regions of sharp transition can be compressed efficiently. Our proposed method falls under this latter category. Mallat [34] formulated the basic idea of multiresolution analysis for orthonormal wavelets in L2 (R). Harten [19, 20, 21] later proposed a general framework for multiresolution representation of data by integrating ideas from three different fields, namely, theory of wavelets, numerical solution of PDEs, and subdivision schemes. Recently, Alves et al. [3] proposed an adaptive multiresolution scheme, similar to the multiresolution approach proposed by Harten [19, 20] and Holmstrom [23] for solving hyperbolic ∗ Received by the editors November 15, 2007; accepted for publication (in revised form) August 18, 2008; published electronically December 31, 2008. This work was partially supported by NSF through awards CMS-0510259 and DMS-0410062. http://www.siam.org/journals/sisc/31-2/70832.html † School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0150 ([email protected], [email protected]). ‡ School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332-0160 (hmzhou@ math.gatech.edu).

1221

1222

S. JAIN, P. TSIOTRAS, AND H.-M. ZHOU

PDEs. These approaches share similar underlying ideas. Namely, the first step is to interpolate the function values at the points belonging to a particular resolution level from the corresponding points at the coarser level and find the interpolative error at the points of that particular resolution level. Once this step has been performed for all resolution levels, all the points that have an interpolative error greater than a prescribed threshold are added to the grid, along with their neighboring points at the same level and the neighboring points at the next finer level. The main difference between these approaches is that in Harten’s approach [19, 20], the solution for each time step is represented on the finest grid, and one calculates the interpolative errors at all the points of the finest grid at each mesh refinement step. On the other hand, Holmstrom [23] and Alves et al. [3], compute the interpolative error only at the points that are in the adaptive grid. If a value that does not exist is needed, Holmstrom interpolates the required function value recursively from coarser scales. Alternatively, Alves et al. [3] add to the grid the points that were used to predict the function values at all previously added points, in order to compute the interpolative error during the next mesh adaptation. Parallel to Harten’s original idea, a modified approach has also been developed by Cohen, M¨ uller, and coauthors for the solution of evolution PDEs [15, 18, 37, 38]. In this paper, we first propose a novel multiresolution scheme for data compression, which results in a higher compression rate compared to the multiresolution approach by Harten [19, 20, 21] for the same desired accuracy. Subsequently, we apply the proposed encoding scheme to solve initial-boundary value problems (IBVP) encountered in evolution PDEs. The proposed multiresolution scheme for data compression works with any of the interpolation techniques mentioned in [21]. One of the key features of our algorithm is that it is a “top-down” (from coarse to fine scale) approach, and we use the most recently updated information to make predictions. Moreover, our interpolations are not restricted to the use of only the retained points at the coarser level, but also use the retained points at the same level (and even the next finer level in the case of solving PDEs). This allows for a more accurate interpolation, which, in turn, leads to fewer points in the final grid. In the proposed algorithm, we continuously keep on updating the grid as we go from the coarsest level to finer levels. If the interpolative error at a point that belongs to a particular level is greater than the prescribed threshold, we add that point to the grid. In the case of solving PDEs, we also add the neighboring points at the same level and the neighboring points at the next level to the grid. We make use of the fact that the point at which the interpolative error is greater than a prescribed threshold, this point is added to the grid, and, in addition, it can be used to predict the remaining points at the same level and the levels below it. Moreover, for refining the mesh for solving evolution PDEs, we predict the function value at a particular point only from the points that are already present in the grid; hence we avoid recursive interpolations from the coarser scales as is done by Holmstrom [23]. At the same time, we do not need to add any extra points to the grid that are required just for computing the interpolative errors at the next mesh refinement step, as is done by Alves et al. [3]. The paper is organized as follows. In the first part of the paper, we start by formulating the problem, and we present the new multiresolution scheme for data compression, followed by the mesh refinement algorithm for the solution of evolution PDEs. Next, we compare the proposed data compression scheme with Harten’s data compression scheme (and the mesh refinement approach with the one of Alves et al. for the case of solution of hyperbolic PDEs). We show that the proposed algorithm

MESH REFINEMENT FOR EVOLUTION PDEs

1223

results, in general, in a fewer number of grid points compared to Harten’s approach (for the case of data compression) and the approach of Alves et al. [3] (for the case of the mesh refinement for the solution of PDEs). In the second part of the paper, we present an algorithm for solving the IBVP for evolution equations on an adaptive nonuniform grid, generated using the proposed mesh refinement technique. This analysis is followed by several numerical examples that show the robustness of the proposed approach and the advantages in terms of computational time compared to the uniform mesh case. 2. Problem statement. Many problems in engineering and physics can be written in the form of an IBVP for an evolution equation: ut + f (uxx , ux , u, x) = 0 in D × (0, ∞), (2.1) (IBVP) : u=g on D × {t = 0}, where D = D ∪ ∂D, with D ⊂ R bounded. The function f : Rm × Rm × Rm × D → R and the initial function g : D → Rm are given. The unknown is the function u : D × [0, ∞) → Rm . The algorithm proposed in this paper works for other boundary conditions as well, but, for simplicity in the analysis below, we only use periodic, Dirichlet, and Neumann boundary conditions. Without loss of generality, we will further assume that D = (0, 1). In (IBVP), the initial function g can be irregular. Even if g is smooth, discontinuities such as shocks (in hyperbolic conservation laws) and kinks (in Hamilton–Jacobi (HJ) equations) can develop in the solution u at some later time. Therefore, we would like to adapt the grid dynamically to any existing or emerging irregularities in the solution, instead of using a fine mesh over the whole spatial and temporal domain. In the next section, we propose a novel grid refinement technique for solving (IBVP) in a computationally efficient manner. 3. Adaptive gridding. First, we give a brief overview on dyadic grids, which are used in the proposed multiresolution mesh refinement scheme for solving (IBVP). 3.1. Dyadic grids. Since D = [0, 1], we consider dyadic grids of the form (3.1)

Vj = {xj,k ∈ [0, 1] : xj,k = k/2j , 0 ≤ k ≤ 2j },

Jmin ≤ j ≤ Jmax ,

where j denotes the resolution level, k the spatial location, and Jmin , Jmax ∈ Z+ 0 . We denote by Wj the set of grid points belonging to Vj+1 \ Vj . Therefore, (3.2) Wj = {yj,k ∈ [0, 1] : yj,k = (2k + 1)/2j+1 , 0 ≤ k ≤ 2j − 1},

Jmin ≤ j ≤ Jmax − 1.

Throughout the rest of this paper, we will use the notations xj,k and yj,k to represent elements of Vj and Wj , respectively. Hence, xj+1,k ∈ Vj+1 is given by xj,k/2 , if k is even, (3.3) xj+1,k = yj,(k−1)/2 , otherwise. An example of a dyadic grid with Jmin = 0 and Jmax = 5 is shown in Figure 3.1. With a slight abuse of notation, we write Vj+1 = Vj ⊕ Wj , although Wj is not an orthogonal complement of Vj in Vj+1 . The subspaces Vj are nested, VJmin ⊂ VJmin +1 · · · ⊂ VJmax , with limJmax →∞ VJmax = D. The sequence of subspaces Wj satisfy Wj ∩ W = ∅ for j = .

1224


V0

k=0

k=1

W0

k=0

W1

k=0

k=1

W2 W3 W4

k=0

k=3 k=7

k=0

k=15

k=0

Fig. 3.1. Example of a dyadic grid.

Before presenting the multiresolution grid adaptation scheme for solving IBVP for evolution PDEs, we propose a novel multiresolution scheme for data compression that will form the basis of our grid adaptation scheme for solving (IBVP). 3.2. Data compression. Suppose g : D → R is specified on a grid VJmax , UJmax = {gj,k : xj,k ∈ VJmax },

(3.4)

where gj,k = g(xj,k ). Let I(x; XGrid ) denote any pth order interpolation of U = {gj,k : xj,k ∈ XGrid }, where XGrid = {xj ,k }i+p =i ⊂ Grid, where Grid ={xji ,ki : xji ,ki ∈ [0, 1], 0 ≤ ki ≤ 2ji , Jmin ≤ ji ≤ Jmax , for i = 0 . . . N, (3.5)

and xji ,ki < xji+1 ,ki+1 , for i = 0 . . . N − 1} ⊂ VJmax

and x ∈ [xji ,ki , xji+p ,ki+p ]. In (3.5), Grid can be uniform or nonuniform. Then the encoding algorithm for compressing the signal g is as follows. Step 1. Initialize an intermediate grid Gridint = VJmin , with function values Uint = Umin , where Umin = {gJmin,k : 0 ≤ k ≤ 2Jmin }. Set j = Jmin . Step 2. DO for k = 0, . . . , 2j . (a) Compute the interpolated function value gˆ(yj,k ) = I(yj,k , XGridint ). (b) If the interpolative error coefficient at the point yj,k , (3.6)

dj,k = |g(yj,k ) − gˆ(yj,k )| > ,

where is the prescribed threshold, then add yj,k to the intermediate grid Gridint and the corresponding function value g(yj,k ) to Uint . Step 3. Increment j by 1. If j < Jmax , go to Step 2, otherwise move on to the next step. Step 4. Terminate the algorithm. The final nonuniform grid representing the compressed information is GridM = Gridint , and the corresponding function values are the set UM = Uint . If we represent the above nonlinear encoding procedure by an operator M , then we can write (3.7)

UM = M UJmax .

One should note that in Harten’s approach [19, 20, 21], the points of a particular resolution level Wj are interpolated only from the corresponding points belonging to Vj . In our approach, instead, we continuously keep on updating the grid, and the 2j −1 points {g(yj,k )}k=0 of level Wj are interpolated from the function values at the points

1225


in Vj ⊕ Wj . Hence, by making use of the extra information from levels Wj —which, in any case, will be added to the adaptive grid—we are able to reduce the number of grid points in the final grid. This process results in higher compression factors, as will be shown in section 3.4 via several examples. The interpolation operator I can be easily constructed using piecewise-polynomial interpolation, essentially nonoscillatory (ENO) interpolation, cubic-spline interpolation, trigonometric interpolation, etc., [21]. For piecewise-polynomial interpolation, the stencil XGridint consists of the p + 1 nearest points to x in Gridint . By p + 1 nearest points, here we mean one neighboring point on the left of x, one neighboring point on the right of x, and the remaining p − 1 points are the points nearest to x in the set Gridint . In case two points are at the same distance from x, that is, if a point on the left and a point on the right are equidistant to x, then we choose a point so as to equalize the number of points on both sides. One may use Neville’s algorithm to construct the respective interpolating polynomials on the fly. For ENO interpolation, the stencil XGridint consists of one neighboring point on the left and one neighboring point on the right of x in the set Gridint , and the remaining p − 1 points are selected from the set Gridint that give the least oscillatory polynomial. For more details on ENO interpolation, the reader is referred to [21, 39]. Next, we give the decoding algorithm, that is, the algorithm for computing (3.8)

UˆJmax = M −1 UM .

One way of decoding the information back from the compressed signal in nonlinear schemes is to keep track of the stencils that were used for interpolating the function values at a particular point while encoding the information and use the same stencils to decode the information from the compressed signal. An alternative way (described below) is to follow the same approach as in the encoding algorithm. Step 1. Initialize Gridint = VJmin , with function values UˆJmax = Uint = Umin , where Umin = {gJmin,k : 0 ≤ k ≤ 2Jmin }. Set j = Jmin . Step 2. DO for k = 0, . . . , 2j . If g(yj,k ) ∈ UM , then add g(yj,k ) to UˆJmax , Uint and yj,k to Gridint , otherwise add gˆ(yj,k ) = I(yj,k , XGridint ) to UˆJmax . Step 3. Increment j by 1. If j < Jmax , go to Step 2, otherwise move on to the next step. Step 4. Terminate the algorithm. Next, we derive an error estimate between the original signal UJmax and the decoded signal UˆJmax . Proposition 1. Let UJmax be defined as in (3.4), and let UˆJmax = M −1 UM , where −1 denotes the decoding algorithm described above. Then for 1 ≤ m < ∞, M ⎛ (3.9) Em (g) = UJmax − UˆJmax m = ⎝ ≤

Jmax 2

1 2Jmax + 1

⎞ m1 |gJmax ,k − gˆJmax ,k |m ⎠

k=0

2Jmax − 2Jmin 2Jmax + 1

m1

and (3.10)

E∞ (g) = UJmax − UˆJmax ∞ =

max

0≤k≤2Jmax

|gJmax ,k − gˆJmax ,k | ≤ .

1226


Proof. First, we note that k = 0, . . . , 2Jmin ,

|gJmin ,k − gˆJmin ,k | = 0,

(3.11)

since gJmin ,k ∈ UM for all k = 0, . . . , 2Jmin . Next, since the function values in the set UˆJmax are interpolated directly only from the function values in UM , we have direct control over the error. Therefore, k = 0, . . . , 2j − 1

|g(yj,k ) − gˆ(yj,k )| ≤ ,

(3.12)

for j = Jmin , . . . , Jmax − 1. Hence, we have UJmax − UˆJmax ∞ =

(3.13)

max

0≤k≤2Jmax

|gj,k − gˆj,k | ≤ .

For 1 ≤ m < ∞, we have (3.14) (3.15)

UJmax − UˆJmax m ⎞ ⎛J m j min Jmax 2 −1 −1 2 1 ⎝ = Jmax |gJmin ,k − gˆJmin ,k |m + |g(yj,k ) − gˆ(yj,k )|m ⎠ 2 +1 k=0

(3.16)

≤

m +1

2Jmax

j=Jmin k=0

j Jmax −1 −1 2

1 = m

j=Jmin k=0

2Jmax − 2Jmin . 2Jmax + 1

Consequently, for 1 ≤ m < ∞, (3.17)

UJmax

− UˆJmax m ≤

which completes the proof. Example 1. Consider g1 : D → R, (3.18) and g2 : D → R,

(3.19)

g1 (x) =

2Jmax − 2Jmin 2Jmax + 1

m1 ,

1, 13 ≤ x ≤ 23 , 0, otherwise,

⎧ ⎪ ⎪0, ⎪ ⎪ ⎪ ⎪1, ⎪ ⎪ ⎨0, g2 (x) = ⎪ sin(πx), ⎪ ⎪ ⎪ ⎪ ⎪ 0, ⎪ ⎪ ⎩ x,

0 ≤ x < 16 , 1 1 6 ≤ x < 3, 1 1 3 ≤ x < 2, 1 2 2 ≤ x < 3, 2 5 3 ≤ x < 6, 5 6 ≤ x ≤ 1.

For this example, we consider a grid with Jmin = 2 and Jmax = 10 and use ENO interpolation for the encoding and the decoding algorithms described above. For both g1 and g2 , the data compression factor (3.20)

C=

2Jmax + 1 − Ngp × 100%, 2Jmax + 1

where Ngp denotes the number of grid points, along with the decoding errors Em , m = 1, 2, ∞, with an interpolating polynomial of degree p = 3 and different thresholds

1227

MESH REFINEMENT FOR EVOLUTION PDEs Table 3.1 Example 1. Data compression along with the decoding errors for the proposed approach.

g1 g2 g2

10−3 10−3 10−7

C 97.56 94.93 93.85

E∞ 0 2.2741 × 10−5 9.0426 × 10−8

E1 0 8.2103 × 10−7 3.7983 × 10−9

10

9

9

8

8

7

7

j

j

10

E2 0 2.8111 × 10−6 1.2662 × 10−8

6

6

5

5

4

4

3

0

0.2

0.4

0.6

0.8

1

3

0

0.2

x

(a) g1 (x).

0.4

0.6

0.8

1

x

(b) g2 (x).

Fig. 3.2. Example 1. Grid point distribution for = 1 × 10−3 .

, are summarized in Table 3.1. First, we consider = 10−3 for both the functions g1 and g2 . The proposed algorithm compressed the given signals g1 and g2 using only 25 and 52 points, respectively. The decoding errors Em (m = 1, 2, ∞) are well below the threshold for both these functions. The grid point distributions for both g1 and g2 are shown in Figure 3.2. Next, for g2 , we decrease the threshold to 10−7 . It is observed that the proposed encoding algorithm increased the number of points used for compressing the signal; once again, the decoding errors are below the prescribed threshold. Next, we present the grid adaptation algorithm for solving (IBVP) using the decoding algorithm mentioned above. 3.3. Grid adaptation for the solution of (IBVP). Assume we are given a nonuniform grid of the form (3.5). For simplicity, we denote by unj,k the value of u(x, t) evaluated at x = xj,k and t = tn , where 0 ≤ k ≤ 2j , Jmin ≤ j ≤ Jmax , n ∈ Z+ 0, t0 = 0, tn = tn−1 + Δtn for n > 0, and Δtn is the time step based on the Courant– Friedrichs–Levy (CFL) condition [46] for hyperbolic equations and the von Neumann condition [46] for all other evolution equations. The “top-down” approach of our algorithm allows one to add and remove points using the most recently updated information. To this end, suppose u(x, tn ) is specified on the grid Gridold , with corresponding solution values Uold = {unj,k : xj,k ∈ Gridold }, where Gridold can be either regular or irregular.1 We assume Gridold ⊇ VJmin . Our aim is to find a new grid Gridnew , by adding or removing points from Gridold , reflecting local changes in the solution. To this end, we initialize an intermediate grid Gridint = VJmin , with the function values Uint = {unJmin,k : unJmin ,k ∈ Uold , 0 ≤ k ≤ 2Jmin }, and we set j = Jmin . The mesh refinement algorithm proceeds as follows. 1 Typically,

Gridold at time t = 0 is regular with Gridold = VJmax .

1228


Step 1. Find the points belonging to the intersection of Wj and Gridold , that is, (3.21)

Y = {yj,ki : yj,ki ∈ Wj ∩ Gridold , for i = 1, . . . , M, 1 ≤ M ≤ 2j }.

Step 2. Set i = 1. (a) Compute the interpolated function values at point yj,ki ∈ Y , uˆ(yj,ki ), that is, u ˆ (yj,ki ) = I(yj,ki , XGridint ), where u ˆ is the th element of u ˆ for = 1, . . . , m.2 (b) If at the point yj,ki ,3 (3.22)

dj,ki (un ) = max |u (yj,ki , tn ) − uˆ (yj,ki )| < , =1,...,m

go to Step 2(f), otherwise add yj,ki to the intermediate grid Gridint and move on to the next step. (c) Add to Gridint N1 points on the left and N1 points on the right neighboring to the point yj,ki in Wj . This step accounts for the possible displacement of any sharp features of the solution during the next time integration step. The value of N1 dictates the frequency of mesh adaptation and is provided by the user. The larger the N1 , the smaller the frequency of mesh adaptation will be, at the expense of a larger number of grid points in the adaptive grid. Hence, there is a trade-off between the frequency of mesh adaptation and the number of grid points. (d) Add to Gridint 2N2 neighboring points at the next finer level 2 {yj+1,2ki + }N =−N2 +1 , where 1 ≤ N2 ≤ 2N1 . This step accounts for the possibility that the solution becomes steeper in this region. Our experience has shown that N2 = N1 is a good choice. (e) Add the function values at all of the newly added points to Uint . If the function value at any of the newly added points is not known, we interpolate the function value at that point from the points in Gridold and their function values in Uold using I(·, XGridold ). (f) Increment i by 1. If i ≤ M , go to Step 2(a), otherwise move on to the next step. Step 3. Increment j by 1. If j < Jmax , go to Step 1, otherwise move on to the next step. Step 4. Terminate the algorithm. The final nonuniform grid is Gridnew = Gridint , and their corresponding function values are the set Unew = Uint . Remark 1. Although the proposed grid adaptation algorithm will work for any interpolation technique, in this paper, we use ENO interpolation to avoid any unphysical interpolation of the data. Remark 2. For sake of brevity, in this paper, we work only with the point-value discretization of data, but the proposed encoding and decoding algorithms (or the grid adaptation algorithm for solving PDEs) will also work for discretizations based on cell-averages. Next, we explain the proposed grid adaptation algorithm with the help of a simple example. Example 2. Consider a dyadic grid V4 and the function 1, x = x4,k , (3.23) g(x) = 0, otherwise, 2 We

would like to remind the reader that u ∈ Rm . that u(yj,k , tn ) ∈ Uold for all yj,k ∈ Y .

3 Note


1229

with k = 6, so that g denotes an impulse located at x = x4,6 = 0.375. Let Jmin = 0, Jmax = 4, p = 1, = 0.1, N1 = N2 = 1, and consider Gridold = VJmax . For this example, the proposed grid adaptation algorithm is illustrated in Figure 3.3. In Figure 3.3, the solid circles show the points belonging to the intermediate grid Gridint and those belonging to Gridnew . The empty squares show the points that are being tested or have been tested for inclusion in Gridint . If the interpolative error coefficient at a point is greater than the prescribed threshold, then we show that point by a solid square. The left and the right neighbors are shown by left and right triangles, respectively. For reference, all points at that particular level are shown by empty circles. Next, we give a point-by-point illustration of the proposed algorithm for Example 2. V0 . Initialize Gridint with VJmin = V0 . W0 . Check if the function value at the point y0,0 ∈ W0 can be interpolated from the nearest p + 1 = 2 points in Gridint , which, in this case, are the points x0,0 and x0,1 . Since, for this example, g(y0,0 ) can be interpolated from the points in Gridint , we do not include y0,0 in Gridint . Next, we move on to level W1 . W1 . y1,0 . Since g(y1,0 ) can again be interpolated from the function values at points x0,0 , x0,1 ∈ Gridint , we do not include y1,0 in Gridint and move on to the next point y1,1 . y1,1 . For the same reason as before, we do not include this point and move on to the next level. W2 . y2,0 . For the same reason as before, we do not include this point. y2,1 . Moving further to y2,1 , we find that g(y2,1 ) cannot be interpolated from the neighboring two points x0,0 , x0,1 ∈ Gridint . Hence, we include y2,1 in Gridint along with points y2,0 , y2,2 ∈ W2 and y3,2 , y3,3 ∈ W3 . y2,2 . Next, we check point y2,2 . The nearest points to y2,2 in Gridint are y3,3 and x0,1 . Since g(y2,2 ) can be interpolated from y3,3 and x0,1 , we move to the next point y2,3 belonging to the level W2 . y2,3 . For the same reason as for point y2,2 , we do not include y2,3 . W3 . y3,0 . Since y3,0 can be interpolated from the existing points in Gridint , we do not include y3,0 in the grid. y3,1 . For the same reason as for y3,0 , we do not include y3,1 in the grid. y3,2 . Subsequently, we check y3,2 . Since g(y3,2 ) cannot be interpolated from the nearest two points y2,0 , y2,1 ∈ Gridint , we include y3,2 (which, in any case, is already present in Gridint ) along with points y3,1 and y3,3 (already present in Gridint ) in Gridint . y3,3 . Moving on to the next point y3,3 , we see again that g(y3,3 ) cannot be interpolated from the nearest two points y2,1 , y2,2 ∈ Gridint . Hence, we include y3,3 (already present in Gridint ) along with y3,2 (already present in Gridint ) and y3,4 in Gridint . y3,4 . Since g(y3,4 ) can be interpolated from the two nearest points y3,3 , y2,2 ∈ Gridint , we move on to the next point y3,5 . y3,5 The nearest two points to y3,5 in Gridint are y2,2 , x0,1 , and since g(y3,5 ) can be interpolated from these two points, we do not include y3,5 in Gridint . y3,6 . For the same reason as for y3,5 we do not add point y3,6 to the grid. y3,7 . For the same reason as for y3,5 we do not add point y3,7 to the grid. Gridnew . The final adaptive grid Gridnew is shown by the solid circles in Figure 3.3. The adaptive grid generated using the previous algorithm depends on how we select points along the grid, that is, whether we move from left to right or from right

1230


g k=0

k=6

k=16

Grid int = VJmin = V0 W0 Grid int W1 Grid int W2 W2 W3 Grid int W2 Grid int W3 W3 Grid int W3 W3 Grid int W3 Grid new Fig. 3.3. Demonstration of the grid adaptation algorithm using Example 2.


1231

to left across each level. It also depends on the location of the singularity. If the singularity is located in the middle, then it does not matter whether we move from left to right or from right to left. The result will be the same nonuniform grid. If, on the other hand, the singularity is not in the middle, then the final grid depends on the way in which we traverse across each level. To illustrate this fact, we again consider Example 2, but this time with k = 1. Hence, the impulse is located at x = x4,1 . If we go from left to right, then the adaptive grid consists of the points x4,0 , x4,1 , x4,3 , x4,5 , x4,16 . If we go from right to left, then the grid consists of the points x4,0 , x4,1 , x4,3 , x4,16 . Let k = 15, which implies that the impulse is located at x = x4,15 . If we go from left to right, then the grid consists of the points x4,0 , x4,13 , x4,15 , x4,16 , and if we go from right to left, then the grid consists of the points x4,0 , x4,11 , x4,13 , x4,15 , x4,16 . Note that, in the proposed algorithm, it is not mandatory to traverse across a level only from the leftmost point or from the rightmost point. We can instead start from any point at that level, each time resulting in a different grid. This suggests that by using a suitable probability distribution function to choose the order in which the points at each particular level are selected, one may be able to further optimize the grid. We will not elaborate on this observation in this paper. 3.4. Comparison with existing multiresolution-based approaches. The proposed encoding and grid adaptation algorithms result, in general, in a fewer number of grid points when compared to Harten’s multiresolution scheme [19, 20, 21], the grid adaptation algorithms of Harten [19, 20], Holmstrom [23], and Alves et al. [3]. First, we explain why this is so and then we give several examples to demonstrate this fact. In the encoding algorithms (grid adaptation algorithms for the solution of IBVP) 2j −1 only from the of existing approaches [19, 20, 23, 3], one interpolates {g(yj,k )}k=0 function values at the points belonging to Vj for j = Jmin . . . Jmax − 1, and only then, 1 one adds to the adaptive grid, the points yj,k (along with the points {yj,k+ }N =−N1 and 2 {yj+1,2k+ }N =−N2 +1 ) for all the pairs (j, k) such that dj,k > . In the proposed method, we continuously keep on updating the adaptive grid instead. If the interpolative error coefficient at yj,k , where 0 ≤ k ≤ 2j − 1 and Jmin ≤ j ≤ Jmax − 1, is greater than the prescribed threshold, we add yj,k to the adaptive grid, while at the same time we 1 also add to the adaptive grid the neighboring points at the same level {yj,k+ }N =−N1 , 2 as well as the neighboring points at the next level {yj+1,2k+ }N =−N2 +1 . We use the newly added point(s) also for interpolating the remaining points at level Wj and the 2j −1 levels below it. In other words, in the proposed approach, the values {g(yj,k )}k=0 are interpolated from the function values at the points in Vj ⊕Wj for j = Jmin , . . . , Jmax −1 (in Vj ⊕ Wj ⊕ Wj+1 for Jmin ≤ j ≤ Jmax − 2 and in Vj ⊕ Wj for j = Jmax − 1). Hence, by making use of the extra information from levels Wj (and Wj+1 ), which, in any case, will be added to the adaptive grid, we are able to reduce the number of grid points in the final grid. We illustrate this fact with the help of several examples. Example 3. We again consider the functions g1 and g2 given by (3.18) and (3.19), respectively, and a grid with Jmin = 2 and Jmax = 10. In order to compare the proposed encoding algorithm (using ENO interpolation) with Harten’s encoding algorithm (using ENO interpolation) and Harten’s encoding algorithm (using central interpolation), we set N1 = N2 = 0 in the grid adaptation algorithm of section 3.3. One should note that, for N1 = N2 = 0, the grid adaptation algorithm given in section 3.3 is simply the decoding algorithm described in section 3.2. The number

1232

S. JAIN, P. TSIOTRAS, AND H.-M. ZHOU Table 3.2 Example 3. Comparison of the proposed decoding approach with Harten’s approach.

g1 g2 g2

10−3 10−3 10−7

Ngp 25 52 63

NgHeno 29 58 73

NgHc 53 108 119

Ngp /NgHeno 0.86 0.90 0.86

Ngp /NgHc 0.47 0.48 0.53

of grid points used by the proposed algorithm Ngp , Harten’s algorithm (using ENO interpolation) NgHeno , and Harten’s algorithm (using central interpolation) NgHc for different thresholds using interpolating polynomial of degree p = 3 are summarized in Table 3.2. For both functions g1 and g2 , we found that the proposed algorithm results in up to 14% fewer number of points in the compressed data compared to Harten’s approach (using ENO interpolation) and up to 53% fewer points compared to Harten’s approach (using central interpolation). It is reminded that Harten’s encoding algorithm (using central interpolation) forms the basis of the multiresolution scheme of Holmstrom [23] and Alves et al. [3] for solving PDEs. In Harten’s approach, the solution at each time step is represented on the finest grid, and one encodes and decodes the solution at each time step in order to calculate the interpolative errors. In other words, the interpolative errors are computed at all points of the fine grid at each mesh refinement step. Holmstrom [23], on the other hand, calculates the interpolative error coefficients only at the points that are in the adaptive grid; if a function value is needed that does not exist in the present grid, the function value is interpolated recursively from coarser scales. In the algorithm of Alves et al. [3], one also adds to the grid the points that were used to predict the function values at all the previously added points in order to compute the interpolative error during the next mesh adaptation. Therefore, in the approach of Alves et al. when a point yj,k (0 ≤ k ≤ 2j − 1 and Jmin ≤ j ≤ Jmax − 1) is added to the grid, one also include its parents, which were used to predict the function value at that point. The parents are not needed for approximating the given function to the prescribed accuracy but are included just for calculating the interpolative error coefficient at the point yj,k during the next mesh adaptation. In the proposed algorithm, on the other hand, whenever a point is being checked for inclusion in the adaptive grid, we predict the function value at that point only from the points already existing in the adaptive grid. Hence, if that point is inserted in the grid, we do not need to add any extra points (i.e., its parents). This also alleviates the task of keeping track of the parents from the rest of the points, as in the approach of Alves et al., and the task of recursively calculating the function values from the coarser resolution levels as is done in the approach of Holmstrom. Next, we give several examples to compare the proposed grid adaptation approach with the algorithm of Alves et al. [3] for solving evolution PDEs. Example 4. First we consider a very simple example. For this example, we consider a dyadic grid V4 and the function 1, x = x4,k , (3.24) g(x) = 0, otherwise, with an impulse located at x = k/24 , where 0 ≤ k ≤ 16. Let Jmin = 0, Jmax = 4, p = 1, = 0.1, and N1 = N2 = 1. Table 3.3 shows the number of grid points used by the proposed grid adaptation algorithm Ngp and the number of points NgA used

1233

MESH REFINEMENT FOR EVOLUTION PDEs Table 3.3 Example 4. Comparison of the proposed algorithm with the algorithm of Alves et al. k 0 1 2 3 4 5 6 7 8

Ngp 9 5 7 6 11 6 9 6 13

NgA 9 6 9 8 12 9 12 9 13

Ngp /NgA 1 0.83 0.78 0.75 0.92 0.67 0.75 0.67 1

k 9 10 11 12 13 14 15 16

Ngp 6 9 6 11 5 7 4 9

NgA 9 12 9 12 7 9 6 9

Ngp /NgA 0.67 0.75 0.67 0.92 0.71 0.78 0.67 1

Table 3.4 Example 5. Comparison of the proposed algorithm with the algorithm of Alves et al.

g1 g2

10−3 10−3

Ngp 53 108

NgA 93 185

Ngp /NgA 0.57 0.58

by the grid adaptation scheme of Alves et al. [3] for k = 0, . . . , 16. We found that when the impulse is located at either the left boundary (k = 0) or the right boundary (k = 16) or in the middle of the domain (k = 8), both the traditional approach and the proposed approach result in the same grid. For all other cases, the grids generated are different. Moreover, we see that the proposed algorithm results in a fewer number of grid points. For this example, the proposed algorithm outperforms the algorithm of Alves et al. [3] by up to 33%. Example 5. We again consider the functions g1 and g2 given by (3.18) and (3.19), respectively, and a grid with Jmin = 2 and Jmax = 10. This time we set N1 = N2 = 1 in the proposed grid adaptation algorithm. Table 3.4 gives the number of points used by the proposed grid adaptation algorithm Ngp and the number of points NgA used by the grid adaptation scheme of Alves et al. [3]. For this example, we observe that the proposed grid adaptation algorithm outperforms the algorithm of Alves et al. [3] by up to 43%. We are now ready to present the algorithm for solving the (IBVP) on an adaptive, nonuniform grid. 4. Numerical solution of the IBVP for evolution equations. The numerical scheme for discretizing (IBVP) depends on f (uxx , ux , u, x). The proposed grid adaptation algorithm will work for many numerically stable discretization schemes for (IBVP). We use different schemes for the numerical examples discussed in this paper, depending on the problem. Hence, in the next section, we describe only the techniques we use for calculating the spatial derivatives ux and uxx on the nonuniform grid, and we state the numerical schemes in the examples themselves. 4.1. Calculation of spatial derivatives. To calculate the derivative ux on the adaptive nonuniform grid Gridnew , we use the weighted ENO (WENO) scheme [27, 28, 33, 39] on nonuniform grids. To this end, let the nonuniform grid be given as in (3.5). Now define (4.1)

D+ unji ,ki =

unji+1 ,ki+1 − unji ,ki xji+1 ,ki+1 − xji ,ki

,

D− unji ,ki =

unji ,ki − unji−1 ,ki−1 xji ,ki − xji−1 ,ki−1

.

1234


n ± A third-order ENO approximation [22, 42, 43] to (u± x )ji ,ki = ux (xji ,ki , tn ) is given by one of the following expressions: n ((u± x )ji ,ki )1 =

(4.2a)

7v2 11v3 v1 − + , 3 6 6

or n ((u± x )ji ,ki )2 = −

(4.2b)

5v3 v4 v2 + + , 6 6 3

or n ((u± x )ji ,ki )3 =

(4.2c)

5v4 v5 v3 + − , 3 6 6

n − n − n where, for calculating (u− x )ji ,ki , we use v1 = D uji−2 ,ki−2 , v2 = D uji−1 ,ki−1 , v3 = − n − n − n n D uji ,ki , v4 = D uji+1 ,ki+1 , v5 = D uji+2 ,ki+2 , and, for calculating (u+ x )ji ,ki , we use v1 = D+ unji+2 ,ki+2 , v2 = D+ unji+1 ,ki+1 , v3 = D+ unji ,ki , v4 = D+ unji−1 ,ki−1 , v5 = D+ unji−2 ,ki−2 . The basic idea behind a third-order ENO scheme is to choose either n ± n ± n ± n ((u± x )ji ,ki )1 or ((ux )ji ,ki )2 or ((ux )ji ,ki )3 for approximating (ux )ji ,ki by choosing the smoothest possible polynomial interpolation of u. n It is reminded that a WENO approximation of (u± x )ji ,ki is a convex combination of the approximations in (4.2a), (4.2b), and (4.2c), that is,

n (u± x )ji ,ki =

(4.3)

3 =1

n ω ((u± x )ji ,ki ) ,

where 0 ≤ ω ≤ 1 for = 1, 2, 3 and ω1 + ω2 + ω3 = 1. The weights for fifth-order accuracy are given by [27, 39] (4.4)

ω =

α , α1 + α2 + α3

= 1, 2, 3,

where (4.5) (4.6) (4.7) (4.8)

α =

α ¯ , (S + δ)2

= 1, 2, 3,

13 1 (v1 − 2v2 + v3 )2 + (v1 − 4v2 + 3v3 )2 , 12 4 13 1 (v2 − 2v3 + v4 )2 + (v2 − v4 )2 , S2 = 12 4 13 1 S3 = (v3 − 2v4 + v5 )2 + (3v3 − 4v4 + v5 )2 , 12 4 S1 =

and (4.9)

α ¯1 = 0.1,

α ¯2 = 0.6,

α ¯3 = 0.3.

In (4.5), δ is used to prevent the denominator from becoming zero. In our computations, we have used δ = 10−6 . For the sake of brevity, we denote the cell walls by (4.10)

xji−1/2 ,ki−1/2 =

xji−1 ,ki−1 + xji ,ki , 2

xji+1/2 ,ki+1/2 =

xji ,ki + xji+1 ,ki+1 . 2


1235

In order to calculate (uxx )nji ,ki = uxx (xji ,ki , tn ) on a nonuniform grid (3.5), we use the centered second difference scheme [46] n n uj −un un ji ,ki ji ,ki −uji−1 ,ki−1 i+1 ,ki+1 − xji+1 ,ki+1 −xji ,ki xji ,ki −xji−1 ,ki−1 (4.11) (uxx )nji ,ki = . xji+1/2 ,ki+1/2 − xji−1/2 ,ki−1/2 We are now ready to give the algorithm for solving the (IBVP) for evolution equation (2.1). 4.2. Solution of the (IBVP) for evolution PDEs. Based on the problem, the desired accuracy, and the computational hardware, we choose the minimum resolution level Jmin , the maximum resolution level Jmax , the threshold , the order of the interpolating polynomial p, and the parameters N1 , N2 required for the grid adaptation algorithm described in section 3.3. The final time tf is assumed to be given. To solve (IBVP) on an adaptive grid, we first initialize Gridold = VJmax , Uold = Jmax {g(xJmax ,k )}2k=0 and set t = 0, n = 0.4 Then the algorithm proceeds as follows. Step 1. Given Gridold , Uold , find the new grid Gridnew and the function values at all the points in Gridnew , Unew = {unj,k : xj,k ∈ Gridnew } using the grid adaptation algorithm given in section 3.3. The new grid Gridnew is the grid on which we will propagate the solution from time t to time t + Δtadapt , where Δtadapt = kΔtn (k ∈ N) is the time after which the grid should be adapted again. For hyperbolic equations, N1 · Δxmin (4.12) k= , Δtn · wave speed where Δxmin = minGridnew (xji+1 ,ki+1 −xji ,ki ). The reader is referred to [46, 39] for details on computing the wave speed. For all other evolution equations, use k = 1. Step 2. Compute the solution at time t = tn+1 at all the points belonging to Gridnew , Unew = {un+1 j,k : xj,k ∈ Gridnew } using any numerically stable scheme and increment n by 1. Keep on repeating this step while t < Δtadapt . If t ≥ tf , terminate the algorithm. Step 3. Reassign the sets: Gridold ← Gridnew , Uold ← Unew . It should be noted that we do not interpolate the function values at the finest level during the mesh refinement process. In the proposed mesh refinement algorithm, we check only the retained points in Gridold to further add and remove points in the grid. The interpolative error coefficients are computed only at the points yj,k ∈ Gridold , and the solution Uold for all yj,k ∈ Gridold is known from the previous step. Step 4. Calculate the new time at which the next mesh adaptation should take place tadapt = t + Δtadapt . Go to Step 1. Remark 3. As pointed out earlier, Δtn is computed based on the CFL condition [46] for hyperbolic equations and the von Neumann condition [46] for all other evolution equations. For both the CFL condition and the von Neumann condition, Δtn depends on Δxmin . Hence, in the proposed algorithm, Δtn changes adaptively depending on Δxmin , which also changes adaptively. J

int case of hardware limitations, we suggest using Gridold = VJint , Uold = {g(xJint ,k )}2k=0 , where Jmin < Jint < Jmax is chosen based on the specific computational hardware used.

4 In

1236


5. Numerical examples. In this section, we present several examples to demonstrate the stability and robustness of our algorithm. These examples also illustrate the ability of the algorithm to automatically capture and follow any existing or selfsharpening features of the solution that develop in time. Example 6. Consider a nonlinear conservation law (5.1)

ut + (F (u))x = 0.

To ensure that shocks and other steep gradients move at the right speed, (5.1) should be written in a discrete conservation form [31, 39]: (5.2)

n un+1 ji ,ki = uji ,ki − Δtn+1

Fjni+1/2 ,ki+1/2 − Fjni−1/2 ,ki−1/2 xji+1/2 ,ki+1/2 − xji−1/2 ,ki−1/2

,

where Fjni±1/2 ,ki±1/2 = F (xji±1/2 ,ki±1/2 , tn ). For a specific example, we consider the inviscid Burgers’ equation 1 2 u (5.3) ut + = 0. 2 x We use the same smooth initial condition and the Dirichlet boundary condition as in [3], that is, (5.4)

g(x) = sin(2πx) +

1 sin(πx), 2

u(0, t) = u(1, t) = 0,

to check the ability of the proposed algorithm to capture the shock. The solution is a wave that develops a very steep gradient and subsequently moves towards x = 1. Because of the zero boundary values, the wave amplitude diminishes with increasing time. For solving (5.3)–(5.4), we use (5.2) along with the ENO–Roe scheme proposed by Shu and Osher [43] on a nonuniform grid for calculating the numerical flux functions Fjni±1/2 ,ki±1/2 . For temporal integration, we use a third-order total variation diminishing (TVD) Runge–Kutta (RK) scheme [42]. The numerical solution at times t = 0 seconds (s), t = 0.158 s, t = 0.5 s, and t = 1 s using a grid with Jmin = 4 and Jmax = 12 are shown in Figure 5.1a. The other parameters used in the grid adaptation procedure are p = 3, = 0.01, and N1 = N2 = 1. Figure 5.1 also shows the grid point distribution in the adaptive mesh at times t = 0 s, t = 0.1 s, t = 0.158 s, and t = 1 s. We see that, as the shock continues to develop, the algorithm adds points at the finer levels of resolution in the region where the shock is developing and removes points from the regions where the solution is getting smoother. Similar conclusions can be drawn by looking at the time evolution of the number of grid points (Figure 5.1f). We observe that the number of grid points increases as the shock continues to develop, and once the solution is smooth everywhere except for the region of the shock, the number of grid points remains pretty much steady, with the number of grid points oscillating about a mean value of 43 points. This shows that the proposed strategy uses only the grid points that are actually necessary to attain a given precision, and the algorithm is able to add and remove points when and where is needed. A comparison of CPU times for the uniform and adaptive grids, along with the L1 error (E1 (u)) between the solution of the proposed multiresolution algorithm and the fine grid solution evaluated at grid VJmax , and the number of grid points used by

1237

MESH REFINEMENT FOR EVOLUTION PDEs 12

1.6 t=0.158

11

t=0 1.2 t=0.5

10

0.8

9

u(x)

t=1

j

0.4

8 7

0 6 −0.4

V12

5

Adaptive −0.8

0

0.2

0.4

0.6

0.8

4

1

0

0.2

0.4

12

12

11

11

10

10

9

9

8

8

7

7

6

6

5

5

0

0.2

0.4

0.8

1

(b) Grid point distribution at t = 0 s.

j

j

(a) Solution u(x, t).

4

0.6

x

x

0.6

0.8

4

1

0

0.2

0.4

x

0.6

0.8

1

x

(c) Grid point distribution at t = 0.1 s.

(d) Grid point distribution at t = 0.158 s.

12

60

11 10 45

Ng

j

9 8 7 30 6 5 4

0

0.2

0.4

0.6

0.8

x

(e) Grid point distribution at t = 1 s.

1

15

0

0.2

0.4

0.6

0.8

1

t

(f) Time evolution of the number of grid points.

Fig. 5.1. Example 6. Parameters used in the simulation are Jmin = 4, Jmax = 12, p = 1, = 0.01, and N1 = N2 = 1.

the proposed algorithm at the final time step for different Jmax = 8, 9, 10, 11, 12 are summarized in Table 5.1. We observe a major speedup in the computational time compared to the uniform mesh, and the speedup factors increase at an approximate rate of two. The proposed approach results in speedup factors that are higher than those reported in [3]. For scale Jmax = 12, the speedup factor using the proposed approach is 63.7, which is about 27% higher than the one reported in [3]. It is reminded that we chose N1 = 1 and Alves et al. [3] chose N1 = 2, which implies that,

1238

S. JAIN, P. TSIOTRAS, AND H.-M. ZHOU Table 5.1 Example 6. L1 error and computational times for uniform versus adaptive mesh.

Jmax 8 9 10 11 12

Uniform mesh Ng in VJmax tcpu (s) 28

+ 1 = 257 29 + 1 = 513 210 + 1 = 1025 211 + 1 = 2049 212 + 1 = 4097

Ng at tf in Gridnew

2.7106 9.3851 36.6631 223.8606 804.9415

Adaptive mesh E1 (u) 10−3

7.1991 × 7.1717 × 10−3 7.7397 × 10−3 7.7220 × 10−3 7.8012 × 10−3

31 34 37 40 43

tcpu (s)

Speedup

0.5835 1.2737 2.6622 6.0399 12.6301

4.6454 7.3684 13.7717 37.0636 63.7320

Table 5.2 Example 6. L1 errors at different times for Jmax = 12. t 0.158 0.5 1

Ng in Gridnew 49 42 43

C (%) 98.80 98.97 98.95

E1 (u) 8.2093 × 10−3 9.3552 × 10−3 7.8012 × 10−3

in our case, mesh refinement was performed twice as many times as was performed in [3] for the same problem, and, even then, the speedup factor is 27% higher than the one reported in [3]. The L1 errors (E1 (u)) along with the number of grid points used by the proposed algorithm at times t = 0.158 s, t = 0.5 s, and t = 1 s for Jmax = 12 have been summarized in Table 5.2. In the next two examples, we consider HJ equations, that is, evolution equations as in (2.1), where (5.5)

f (uxx , ux , u, x) = f (ux ).

For discretizing f (ux ), we use the Lax–Friedrich’s (LF) scheme [16, 39, 40, 43] + ux + u− 1 x + − , u ) = f (5.6) f (ux ) = fˆLF (u− − αx (u+ x x x − ux ), 2 2 where, αx = maxux ∈I x |f1 (ux )|, f1 is the partial derivative of f with respect to ux , max ], and the minimum and the maximum values of ux are identified by I x = [umin x , ux + considering all the values of u− x and ux on the nonuniform grid. Example 7. First, we consider the HJ equation with convex f (ux ), taken from [40]: (5.7)

ut +

(ux + 1)2 = 0, 2

with the initial condition and the periodic boundary condition as in [40], that is, (5.8)

g(x) = − cos πx,

u(−1, t) = u(1, t),

− 1 ≤ x < 1.

For solving the problem using the proposed algorithm, we first convert the above mentioned problem from x ∈ [−1, 1] to x ˆ ∈ [0, 1] by using a simple change of variables x = 2ˆ x − 1. With a slight abuse of notation, we denote x ˆ by x, and hence the problem(5.7)–(5.8) transforms to (5.9)

ut +

( 12 ux + 1)2 = 0, 2


1239

with the initial condition and the periodic boundary condition (5.10)

g(x) = − cos π(2x − 1),

u(0, t) = u(1, t).

− The derivatives u+ x , ux in the LF discretization are approximated using a WENO scheme, and the temporal integration is performed using a third-order TVD RK scheme. The numerical solution at times t = 0 s, t = 1.5/π 2 s, t = 3.5/π 2 s, t = 7/π 2 s, t = 10/π 2 s, and t = 14/π 2 s using a grid with Jmin = 4 and Jmax = 12 are shown in Figure 5.2a. The other parameters used in the grid adaptation algorithm are p = 3, = 0.001, and N1 = N2 = 2. Figure 5.2 also shows the grid point distribution in the adaptive mesh at times t = 0 s, t = 3.5/π 2 s, t = 7/π 2 s, and t = 14/π 2 s. We see that as the kink continues to develop, the algorithm adds points at the finer levels of resolution in the region where the kink is developing and removes points from the regions where the solution is getting smoother and smoother. As HJ equation (5.9) continues to evolve further in time, the discontinuity in the first derivative of the solution is smoothing out, and, as a result, the algorithm starts removing points from the finer levels of resolution. This, again, demonstrates that the proposed strategy uses only the grid points that are actually necessary to attain a given precision, and the algorithm is able to add and remove points when and where needed. The L1 errors (E1 (u)) along with the number of grid points used by the proposed algorithm at times t = 1.5/π 2 s, t = 3.5/π 2 s, t = 7/π 2 s, t = 10/π 2 s, and t = 14/π 2 s for Jmax = 12 have been summarized in Table 5.3. Example 8. Next, we consider the HJ equation with nonconvex f (ux ),

(5.11)

ut − cos(αux + 1) = 0,

with (5.12)

g(x) = − cos π(2x − 1),

u(0, t) = u(1, t),

where α is a constant. We again use an LF scheme (5.6) for solving the IBVP (5.11)– − (5.12). The derivatives u+ x , ux in the LF discretization are approximated using a WENO scheme, and the temporal integration is performed using a third-order TVD RK scheme. The choice of α = 0.5 results in the commonly used test problem for onedimensional (1-D) HJ equations given in [40]. In order to make the problem more interesting and challenging, in this work, we consider two more choices for α, namely, α = 1 and α = 1.5. The choices α = 1 and α = 1.5 result in more kinks in the solution at time t = 1.5/π 2 s. The numerical solutions for all the cases at time t = 1.5/π 2 s using a grid with Jmin = 4 and Jmax = 12 along with the corresponding grid point distributions are shown in Figure 5.3. The other parameters used in the grid adaptation algorithm are p = 3, = 0.001, and N1 = N2 = 2. The solutions at t = 1.5/π 2 s for α = 0.5, 0.1, and 1.5 have two, four, and six kinks, respectively. We once again observe that the proposed algorithm is able to capture all the kinks in the solutions accurately and efficiently by adding points at the finer resolution levels in the region of kinks, while resolving the smoother regions using only the points at the coarse resolution levels. The L1 errors (E1 (u)) along with the number of grid points used by the proposed algorithm at time t = 1.5/π 2 s for α = 0.5, 1, 1.5, and Jmax = 12 have been summarized in Table 5.4.

1240

S. JAIN, P. TSIOTRAS, AND H.-M. ZHOU 1

12

V12

t=0

Adaptive

11

0.5 2

t=1.5/π

10

0

2

t=3.5/π

−1

j

u(x)

9

−0.5

8 7

2

t=7/π

6

2

t=10/π −1.5

2

5

t=14/π −2 −1

−0.5

0

0.5

4 −1

1

−0.5

0

x

12

12

11

11

10

10

9

9

8

8

7

7

6

6

5

5

−0.5

0

1


j

j


4 −1

0.5

x

0.5

4 −1

1

−0.5

0

x

0.5

1

x

(c) Grid point distribution at t = 3.5/π 2 s.

(d) Grid point distribution at t = 7/π 2 s. 80

12 11 10

60

Ng

j

9 8 7

40

6 5 4 −1

−0.5

0

0.5

1

20

0

0.5

1

1.5

t

x

(e) Grid point distribution at t = 14/π 2 s.

(f) Time evolution of the grid points.

Fig. 5.2. Example 7. Parameters used in the simulation are Jmin = 4, Jmax = 12, p = 3, = 0.001, and N1 = N2 = 2. Table 5.3 Example 7. L1 errors at different times for Jmax = 12. t 1.5/π 2 3.5/π 2 7/π 2 10/π 2 14/π 2

Ng in Gridnew 52 50 39 44 31

C (%) 98.73 98.78 99.05 98.93 99.24

E1 (u) 1.7603 × 10−3 4.2828 × 10−3 5.4194 × 10−3 6.0113 × 10−3 6.5118 × 10−3

1241


1.2

12 11

0.8 10 0.4

0

j

u(x)

9 8 7 −0.4 2

t=1.5/π

6

−0.8 5

V12

t=0

Adaptive −1.2

0

0.2

0.4

0.6

0.8

4

1

0

0.2

0.4

x

0.6

0.8

1

x

(a) Solution u(x, 1.5/π 2 ) for α = 0.5.

(b) Grid point distribution at t = 1.5/π 2 s for α = 0.5.

1.2

12 11

0.8 10 0.4

0

j

u(x)

9 8 7 −0.4 t=1.5/π2

6

−0.8 5

V12

t=0

Adaptive −1.2

0

0.2

0.4

0.6

0.8

4

1

0

0.2

0.4

x

0.6

0.8

1

x

(c) Solution u(x, 1.5/π 2 ) for α = 1.

(d) Grid point distribution at t = 1.5/π 2 s for α = 1.

1.2

12 11

0.8 10 0.4

j

u(x)

9 0

8 7

−0.4 t=1.5/π2

6

−0.8 5

V12

t=0

Adaptive −1.2

0

0.2

0.4

0.6

0.8

1

4

0

0.2

0.4

x

0.6

0.8

1

x

(e) Solution u(x, 1.5/π 2 ) for α = 1.5.

(f) Grid point distribution at t = 1.5/π 2 s for α = 1.5.

Fig. 5.3. Example 8. Parameters used in the simulation are Jmin = 4, Jmax = 12, p = 3, = 0.001, and N1 = N2 = 2. Table 5.4 Example 8. L1 errors at different times for Jmax = 12. α 0.5 1 1.5


C (%) 98.93 97.07 96.95

E1 (u) 1.1259 × 10−3 8.8733 × 10−4 5.6106 × 10−4

1242


Example 9. Consider the scalar reaction-diffusion problem that appears in combustion problems [25, 41]: (5.13)

ut − uxx −

(5.14)

ux (0, t) = 0,

Reδ (1 + a − u)e−δ/u = 0, aδ u(1, t) = 1,

u(x, 0) = 1.

The solution u represents the temperature of a reactant in a chemical system, a is the heat release, δ is the activation energy, and R is the reaction rate. For small times, the temperature gradually increases from unity with a “hot spot” forming at x = 0. After some finite time, ignition occurs, and the temperature at x = 0 jumps rapidly from near unity to near 1 + a. A flame front then forms and propagates towards x = 1 with a speed proportional to eaδ /2(1 + a). In real problems, a is close to unity and δ is large, thus the flame front moves exponentially fast after the ignition. We use the same problem parameters as in [25, 41], namely, a = 1, R = 5, and δ = 30. This is the same problem as the one given in [2], except for the value of the parameter δ, which in [2] is taken to be 20. Instead, we consider δ = 30 as in [25, 41], since the flame layer in this case is much thinner, and higher mesh adaptation is required. We use (4.11) to discretize uxx and use a third-order TVD RK scheme for temporal integration. To illustrate how we apply the Neumann boundary condition on a nonuniform grid, we again consider a grid of the form (3.5). To apply the Neumann boundary condition ux (0, t) = 0, we introduce a fictitious node xj−1 ,k−1 = −xj1 ,k1 , which lies outside the physical domain,5 and approximate the boundary condition by (ux )nj0 ,k0 =

(5.15)

unj1 ,k1 − unj−1 ,k−1 xj1 ,k1 − xj−1 ,k−1

= 0,

which implies unj−1 ,k−1 = unj1 ,k1 . Hence, at the boundary x = 0, (4.11) reduces to n n un u −un j ,k0 −uj−1 ,k−1 j0 ,k0 − 00−x 2 j1x,kj 1,k −0 j−1 ,k−1 2(unj1 ,k1 − unj0 ,k0 ) 1 1 (5.16) (uxx )nj0 ,k0 = = . xj1 ,k1 − xj−1 ,k−1 (xj1 ,k1 )2 The numerical solutions at times t = 0 s, t = 0.24 s, t = 0.241 s, and t = 0.244 s using a grid with Jmin = 4 and Jmax = 12 are shown in Figure 5.4a. The other parameters used in the grid adaptation algorithm are p = 3, = 10−5 /2Jmax−j , and N1 = N2 = 1. One of the main challenges of this problem is the fact that one needs to use a very small time step to capture the transition layer during the time of ignition. This is achieved automatically by the proposed algorithm, since the algorithm is adaptive both in time and space. As the mesh gets refined, the Δtn in the proposed algorithm for the solution of evolution PDEs given in section 4.2 also decreases. We see from Figure 5.4f that for time t < 0.195 s, the proposed algorithm found the solution using only about 50 to 65 points. Starting from t = 0.195 s to t = 0.2385 s, the algorithm slowly increased the number of points to around 95 points and thereafter added points at finer levels starting at t = 0.2385 s. As points from finer grid levels are being added, the algorithm automatically decreases the time step and is able to capture the transition layer during the time of ignition. The L1 errors (E1 (u)) along with the number of grid points used by the proposed algorithm at times t = 0.24 s, t = 0.241 s, and t = 0.244 s for Jmax = 12 have been summarized in Table 5.5. 5 Note

that xj0 ,k0 = 0.

1243

MESH REFINEMENT FOR EVOLUTION PDEs 2

12

1.9

V12

11

Adaptive

1.8

10 1.7 9 t=0.244

t=0.241

1.5

j

u(x)

1.6

1.4

8 7

1.3 6 t=0.24

1.2

5

1.1 1

0

0.2

0.4

0.6

0.8

4

1

0

0.2

0.4

12

12

11

11

10

10

9

9

8

8

7

7

6

6

5

5

0

0.2

0.4

0.8

1


j

j


4

0.6

x

x

0.6

0.8

4

1

0

0.2

0.4

0.6

0.8

1

x

x

(c) Grid point distribution at t = 0.24 s.

(d) Grid point distribution at t = 0.241 s.

12 300

11 250

10 200

8

Ng

j

9

150

7 100

6 50

5 4

0

0.2

0.4

0.6

0.8

1

0

0

x

0.05

0.1

0.15

0.2

0.25

t

(e) Grid point distribution at t = 0.244 s.

(f) Time evolution of the number of grid points.

Fig. 5.4. Example 9. Parameters used in the simulation are Jmin = 4, Jmax = 12, p = 3, = 10−5 /2Jmax −j , and N1 = N2 = 1. Table 5.5 Example 9. L1 errors at different times for Jmax = 12. t 0.24 0.241 0.244


C (%) 96.66 94.46 95.61

E1 (u) 4.6714 × 10−4 3.4834 × 10−3 3.2098 × 10−3

1244


Example 10. Finally, we consider a Riemann initial value problem (shock tube) for the Euler equations of gas dynamics as follows: (5.17)

ut + f (u)x = 0,

(5.18)

u(x, 0) =

uL , x < 0.5, uR , x > 0.5,

where (5.19)

u = [ρ m E]T ,

(5.20)

f (u) = νu + [0 p pν]T ,

ρ, m, E are the gas density, momentum, total energy per unit volume, respectively, ν = m/ρ is the velocity, and (5.21)

ρν 2 p = (γ − 1) E − 2

is the pressure. In (5.21), γ is the ratio of specific heat, which takes the usual value of 1.4 (for air). We consider the two well-known problems, namely, Sod’s problem [44], the initial data for which is given by (5.22)

uL = [1 0 2.5]T ,

uR = [0.125 0 0.25]T ,

and Lax’s problem [30], the initial data for which is given by (5.23)

uL = [0.445 0.698 8.82]T ,

uR = [0.5 0 1.4275]T.

We use the characteristic numerical scheme given in [43, 39] for solving this problem. The basic idea behind the characteristic scheme is to transform the nonlinear system (5.17) to a system of (nearly) independent scalar conservation laws and discretize each scalar conservation law independently in an upwind biased fashion. Then we transform the discretized system back to the original variables. We use the ENO– Roe fix scheme [43] on a nonuniform grid for obtaining the numerical flux function Fjni±1/2 ,ki±1/2 in the scalar field, and we use a third-order TVD RK scheme for temporal integration. The numerical solution of the density ρ(x, t) along with the grid point distributions in the adaptive mesh for Sod’s problem and Lax’s problem at times t = 0.2 s, t = 0.13 s, respectively, using a grid with Jmin = 4 and Jmax = 12 are shown in Figure 5.5. The other parameters used in the grid adaptation procedure are p = 3, = 0.001, and N1 = N2 = 2. Figure 5.5 also shows the time evolution of the number of grid points for both Sod’s and Lax’s problems. The L1 errors (E1 (ρ), E1 (m), E1 (E)) along with the number of grid points used by the proposed algorithm for solving Sod’s problem at times t = 0.05 s, t = 0.1 s, t = 0.15 s, and t = 0.2 s and Lax’s problem at times t = 0.05 s, t = 0.1 s, and t = 0.13 s for Jmax = 12 are summarized in Table 5.6 and Table 5.7, respectively.

1245


1

V12

1.3

Adaptive

0.7 ρ

1.1

ρ

0.9

0.4 0.7

0.5

V12

Adaptive 0.1

0

0.2

0.4

0.6

0.8

0.3

1

x

11

11

10

10

9

9

8

8

j

12

7

7

6

6

5

5

0

0.2

0.4

0.4

0.6

0.8

1

(b) Lax’s problem. Solution ρ(x, 0.13).

12

4

0.2

x

(a) Sod’s problem. Solution ρ(x, 0.2).

j

0

0.6

0.8

4

1

0

0.2

0.4

x

0.6

0.8

1

x

(c) Sod’s problem. Grid point distribution at (d) Lax’s problem. Grid point distribution at t = 0.2,s. t = 0.13,s. 240

350

300

200

250

g

N

N

g

160 200

120 150 80

40

100

0

0.05

(e) Sod’s problem. points.

0.1 t

0.15

0.2

50

0

0.05

Time evolution of grid (f) Lax’s problem. points.

0.1

0.14

t

Time evolution of grid

Fig. 5.5. Example 10. Parameters used in the simulation are Jmin = 4, Jmax = 12, p = 3, = 0.001, and N1 = N2 = 2.

1246

S. JAIN, P. TSIOTRAS, AND H.-M. ZHOU Table 5.6 Example 10. Sod’s problem. L1 errors at different times for Jmax = 12. t 0.05 0.1 0.15 0.2

Ng in Gridnew 212 189 173 195

C (%) 94.83 95.39 95.78 95.24

E1 (ρ) 1.0300 × 10−4 2.8712 × 10−4 4.9362 × 10−4 7.8443 × 10−4

E1 (m) 1.1859 × 10−4 3.2164 × 10−4 5.4437 × 10−4 8.1571 × 10−4

E1 (E) 2.9885 × 10−4 8.3684 × 10−4 1.4215 × 10−3 2.1954 × 10−3

Table 5.7 Example 10. Lax’s problem. L1 errors at different times for Jmax = 12. t 0.05 0.1 0.13


C (%) 93.36 93.41 93.48

E1 (ρ) 5.6092 × 10−5 1.8641 × 10−4 2.7005 × 10−4

E1 (m) 1.2380 × 10−4 4.1361 × 10−4 5.9612 × 10−4

E1 (E) 7.7312 × 10−4 3.3487 × 10−3 4.9735 × 10−3

6. Conclusions. In this paper, we have proposed a novel multiresolution scheme for data compression along with a new grid adaptation algorithm for solving evolution equations. Both are shown to outperform similar data compression and grid adaptation schemes in the literature. Specifically, we have shown that the proposed approach results in fewer grid points when compared to a common adaptive grid approach. For the examples considered here, we have observed savings of up to 43% in terms of the number of grid points and a gain of about 27% in terms of speedup factors. Several examples have demonstrated the stability and robustness of the proposed algorithm. In all examples considered, the algorithm adapted dynamically to any existing or emerging irregularities in the solution by automatically allocating more grid points to the region where the solution exhibited sharp features and fewer points to the region where the solution was smooth. As a result, the computational time and memory usage can be reduced, while maintaining an accuracy equivalent to the one obtained using a fine uniform mesh. For the solution of evolution equations in which the mesh refinement is performed at each time step, that is, Δtadapt = Δt, choosing a value of N1 > 1 will only result in more points in the grid; it will not decrease the frequency of mesh adaptation. For all other cases, the larger the N1 , the smaller the frequency of mesh adaptation will be, at the expense of a larger number of grid points in the adaptive grid. Hence, there is a trade-off between the frequency of mesh adaptation and the number of grid points. Our practical experience has shown that N1 = N2 = 2 is a good choice for solving hyperbolic equations. Future work should investigate the issue of how to choose the parameters N1 , N2 in an optimal way, including the adaptivity of these parameters across the spatial domain. Followup work should also concentrate on extending the proposed approach to multiple dimensions. It is expected that the savings in terms of CPU time and the number of grid points observed for the single spatial dimension case will be greater in multiple dimensions. One approach in this direction is to work directly with interpolating functions in higher dimensions and then follow the same approach as for the 1-D case. That is, use the error between the actual and interpolated values from neighboring points to determine which points to retain in the grid and which to remove. The challenge is to find a consistent way of selecting neighboring points. Another idea is to proceed in a dimension-by-dimension fashion. That is, to compute the interpola-


1247

tive error coefficients at a particular point using an interpolation operator based on function values in the intermediate grid along one direction, while keeping the other coordinates fixed and repeating the same for all directions. These issues are left for future investigation. Acknowledgments. We would like to thank the anonymous reviewers for their insightful suggestions. REFERENCES [1] S. Adjerid and J. Flaherty, A moving mesh finite element method with local refinement for parabolic partial differential equations, Comput. Methods Appl. Mech. Engrg., 56 (1986), pp. 3–26. [2] S. Adjerid and J. E. Flaherty, A moving finite element method with error estimation and refinement for one-dimensional time dependent partial differential equations, SIAM J. Numer. Anal., 23 (1986), pp. 778–796. ˜ es, F. T. Pinho, and P. J. Oliveira, [3] M. A. Alves, P. Cruz, A. Mendes, F. D. Magalha Adaptive multiresolution approach for solution of hyperbolic PDEs, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 3909–3928. [4] D. C. Arney and J. E. Flaherty, A two-dimensional mesh moving technique for time dependent partial differential equations, J. Comput. Phys., 67 (1986), pp. 124–144. [5] I. Babuˇska, J. Chandra, and J. E. Flaherty, eds., Adaptive Computational Methods for Partial Differential Equations, SIAM, Philadelphia, 1983. [6] M. Baines and A. Wathen, Moving finite element methods for evolutionary problems. I. Theory, J. Comput. Phys., 79 (1988), pp. 245–269. [7] M. J. Baines, Moving Finite Elements, Oxford University Press, New York, 1994. [8] J. Bell, M. Berger, J. Saltzman, and M. Welcome, Three-dimensional adaptive mesh refinement for hyperbolic conservation laws, SIAM J. Sci. Comput., 15 (1994), pp. 127– 138. [9] M. J. Berger and P. Colella, Local adaptive mesh refinement for shock hydrodynamics, J. Comput. Phys., 82 (1989), pp. 64–84. [10] M. J. Berger and R. J. Leveque, Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems, SIAM J. Numer. Anal., 35 (1998), pp. 2298–2316. [11] M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys., 53 (1984), pp. 484–512. [12] S. Bertoluzza, Adaptive wavelet collocation method for the solution of Burger’s equation, Transport Theory Statist. Phys., 25 (1996), pp. 339–352. [13] N. N. Carlson and K. Miller, Design and application of a gradient-weighted moving finite element code I: In one dimension, SIAM J. Sci. Comput., 19 (1998), pp. 728–765. [14] H. D. Ceniceros and T. Y. Hou, An efficient dynamically adaptive mesh for potentially singular solutions, J. Comput. Phys., 172 (2001), pp. 609–639. ¨ ller, and M. Postel, Fully adaptive multiresolution finite [15] A. Cohen, S. M. Kaber, S. Mu volume schemes for conservation laws, Math. Comp., 72 (2003), pp. 183–225. [16] M. Crandall and P. Lions, Two approximations of solutions of Hamilton-Jacobi equations, Math. Comp., 43 (1984), pp. 1–19. [17] E. A. Dorfi and L. O. Drury, Simple adaptive grids for 1-d initial value problems, J. Comput. Phys., 69 (1987), pp. 175–195. ¨ ller and S. Mu ¨ ller, Adaptive finite volume schemes for conservation [18] B. Gottschlich-Mu laws based on local multiresolution techniques, in Hyperbolic Problems: Theory, Numerics, Applications, M. Fey and R. Jeltsch, eds., Birkh¨ auser, Basel, 1999, pp. 385–394. [19] A. Harten, Adaptive multiresolution schemes for shock computations, J. Comput. Phys., 115 (1994), pp. 319–338. [20] A. Harten, Multiresolution algorithms for the numerical solution of hyperbolic conservation laws, Comm. Pure Appl. Math., 48 (1995), pp. 1305–1342. [21] A. Harten, Multiresolution representation of data: A general framework, SIAM J. Numer. Anal., 33 (1996), pp. 1205–1256. [22] A. Harten, B. Enquist, S. Osher, and S. Chakravarthy, Uniformly high-order accurate essentially non-oscillatory schemes III, J. Comput. Phys., 71 (1987), pp. 231–303. ¨ m, Solving hyperbolic PDEs using interpolating wavelets, SIAM J. Sci. Comput., [23] M. Holmstro 21 (1999), pp. 405–420.

1248


¨ m and J. Wald´ [24] M. Holmstro en, Adaptive wavelet methods for hyperbolic PDEs, J. Sci. Comput., 13 (1998), pp. 19–49. [25] W. Huang, Y. Ren, and R. D. Russell, Moving mesh methods based on moving mesh partial differential equations, J. Comput. Phys., 113 (1994), pp. 279–290. [26] L. Jameson, A wavelet-optimized, very high order adaptive grid and order numerical method, SIAM J. Sci. Comput., 19 (1998), pp. 1980–2013. [27] G. S. Jiang and D. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput., 21 (2000), pp. 2126–2143. [28] G. S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys., 126 (1996), pp. 202–228. [29] I. W. Johnson, A. J. Wathen, and M. J. Wathen, Moving finite element methods for evolutionary problems. II. Applications, J. Comput. Phys., 79 (1988), pp. 270–297. [30] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Commun. Pure Appl. Math., 7 (1954), pp. 195–206. [31] R. LeVeque, Numerical Methods for Conservation Laws, Birkh¨ auser, Boston, 1992. [32] R. Li and T. Tang, Moving mesh discontinuous Galerkin method for hyperbolic conservation laws, J. Sci. Comput., 27 (2006), pp. 347–363. [33] X. D. Liu, S. Osher, and T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys., 115 (1994), pp. 200–212. [34] S. G. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2 (R), Trans. Amer. Math. Soc., 315 (1989), pp. 69–87. [35] K. Miller, Moving finite elements II, SIAM J. Numer. Anal., 18 (1981), pp. 1033–1057. [36] K. Miller and R. N. Miller, Moving finite elements I, SIAM J. Numer. Anal., 18 (1981), pp. 1019–1032. ¨ ller, Adaptive Multiscale Schemes for Conservation Laws, Lect. Notes Comput. Sci. [37] S. Mu Eng. 27, Springer, New York, 2003. ¨ ller and Y. Stiriba, Fully adaptive multiscale schemes for conservation laws employing [38] S. Mu locally varying time stepping, J. Sci. Comput., 30 (2007), pp. 493–531. [39] S. Osher and R. P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, New York, 2003. [40] S. Osher and C.-W. Shu, High-order essentially nonoscillatory schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal., 28 (1991), pp. 902–922. [41] L. R. Petzold, Observations on an adaptive moving grid method for one-dimensional systems for partial differential equations, Appl. Numer. Math., 3 (1987), pp. 347–360. [42] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys., 77 (1988), pp. 439–471. [43] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes II, J. Comput. Phys., 83 (1989), pp. 32–78. [44] G. A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys., 27 (1978), pp. 1–31. [45] H. Tang and T. Tang, Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws, SIAM J. Numer. Anal., 41 (2003), pp. 487–515. [46] J. W. Thomas, Numerical Partial Differential Equations, Springer, New York, 1995. [47] J. F. Thompson, A survey of dynamically-adaptive grids in the numerical solution of partial differential equations, Appl. Numer. Math., 1 (1985), pp. 3–27. [48] O. V. Vasilyev, Solving multi-dimensional evolution problems with localized structures using second generation wavelets, Int. J. Comput. Fluid Dyn., 17 (2003), pp. 151–168. [49] O. V. Vasilyev and C. Bowman, Second-generation wavelet collocation method for the solution of partial differential equations, J. Comput. Phys., 165 (2000), pp. 660–693. [50] J. Wald´ en, Filter bank methods for hyperbolic PDEs, J. Numer. Anal., 36 (1999), pp. 1183– 1233.