Adaptive discontinuous Galerkin methods in multiwavelets bases

Applied Numerical Mathematics 61 (2011) 879–890

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Adaptive discontinuous Galerkin methods in multiwavelets bases Rick Archibald ∗ , George Fann, William Shelton Computer Science and Mathematics, Oak Ridge National Laboratory, Oak Ridge, TN 37831, United States

a r t i c l e

i n f o

Article history: Received 29 December 2009 Received in revised form 31 January 2011 Accepted 1 February 2011 Available online 9 March 2011 Keywords: Multiwavelets Discontinuous Galerkin

a b s t r a c t We use a multiwavelet basis with the Discontinuous Galerkin (DG) method to produce a multi-scale DG method. We apply this Multiwavelet DG method to convection and convection–diffusion problems in multiple dimensions. Merging the DG method with multiwavelets allows the adaptivity in the DG method to be resolved through manipulation of multiwavelet coefficients rather than grid manipulation. Additionally, the Multiwavelet DG method is tested on non-linear equations in one dimension and on the cubed sphere. Published by Elsevier B.V. on behalf of IMACS.

1. Introduction The discontinuous Galerkin (DG) method is a finite element method that is locally conservative and stable with highorder accuracy. The DG formulation utilizes an elementwise discontinuous approximation, where numerical information only passes locally through numerical fluxes, to represent the dynamics and structure of highly complex solutions. Additionally, the formulation of DG allows this method to handle complex geometries, irregular meshes with hanging nodes, hp-adaptation, scalability, and solution approximations that may have different basis functions of different order in different elements. DG has been applied to various engineering and scientific problems [19]. Wavelets in computational PDEs began by demonstrating wavelet numerical solution of PDEs with periodic boundary conditions [10,12,14,15]. Wavelets have proved to be an efficient tool in developing adaptive numerical methods which control the global approximation error [6,8,11]. This research focuses on multiwavelets, that significantly differ from wavelet in that the basis functions are discontinuous, organized at different scales in small groups of several functions, and generated by multiple scaling functions. The development of multiwavelets has been along a similar time-line as DG. Multiwavelets were first established as a multiresolution basis to generate effective sparse integro-differential operator representation [1], and then later were demonstrated to be well suited for high-order adaptive solvers of partial differential equations [2]. As a bases, multiwavelets are a discontinuous, orthogonal, compactly supported, multiscale set of functions with vanishing moments that yield highorder adaptive approximations of L 2 functions. With respect to integro-differential operators, the properties of multiwavelets has made them a bases that can provide an effective sparse representation of a wide class of integro-differential operators. This research combines multiwavelets with the DG method, and demonstrates how this union results in a multi-scale adaptive DG method. We refer to this new approach as multiwavelet DG. Merging the DG method with multiwavelets allows the adaptivity in the DG method to be resolved through manipulation of multiwavelet coefficients rather than grid manipulation. Wavelets have been used with the DG method as a post-processing step on cell averages for mesh adaptivity [9]. The work presented in this paper is unique in the fact that multiwavelets are fully merged into the DG method, and adaptivity

*

Corresponding author. E-mail addresses: [email protected] (R. Archibald), [email protected] (G. Fann), [email protected] (W. Shelton).

0168-9274/$30.00 Published by Elsevier B.V. on behalf of IMACS. doi:10.1016/j.apnum.2011.02.005

880

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

can be achieved directly without a post-processing step using all the information contained in each element as opposed to cell averages. This paper is organized as follows. In Section 2 we introduce the multiwavelet bases and its key features. Sections 3 and 4 describe the DG method for convection and convection–diffusion problems, respectively, and further demonstrate how multiwavelet bases can be incorporated into these methods. Section 5 describes the method of time stepping used in this research, which has been demonstrated to be particularly effective and efficient for multiwavelet based schemes [2,3]. Examples in multiple dimensions and different geometries are given in Section 6 for both convection and convection– diffusion problems. Section 7 finishes this paper with conclusions and a discussion of future directions. 2. Multiwavelet bases In this section we summarize some properties of the multiwavelet bases derived and developed in [1] and [2]. For k = 1, 2, . . . , and n = 0, 1, 2, . . . , we define Vnk as a space of piecewise polynomial functions,

 





Vkn = f : f ∈ Πk I n , for  = 0, . . . , 2n − 1, and supp( f ) = I n ,

(2.1)

where Πk ( I  ) is the space of all polynomials of degree less than k on the interval I  = n

n

(2−n , 2−n (

+ 1)]. In Sections 3

and 4 we strictly use the I n notation to describe all intervals. The space Vnk has dimension 2n k and has the following nested property,

Vk0 ⊂ Vk1 ⊂ · · · ⊂ Vkn ⊂ · · · . The multiwavelet subspace

Vkn ⊕ Wkn = Vkn+1 ,

Wnk ,

(2.2)

n = 0, 1, 2, . . . , is defined as the orthogonal complement of

Vnk

in

Vnk +1

or

Wkn ⊥Vkn ,

(2.3) 1

where the norm is defined as  f  =  f , f 2 , with inner product

1  f , g =

f (x) g (x) dx.

(2.4)

0

From the previous definition, we can split Vnk into n + 1 orthogonal subspaces as,

Vkn = Vk0 ⊕ Wk0 ⊕ Wk1 ⊕ · · · ⊕ Wkn−1 . Given an orthonormal basis φ0 , . . . , φk−1 of φ0 , . . . , φk−1 by dilation and translation,

  φ nj (x) = 2n/2 φ j 2n x −  ,

(2.5) Vk0 ,

the space

Vnk

n

is spanned by 2 k functions which are obtained from

j = 0, . . . , k − 1,  = 0, . . . , 2n − 1.

(2.6)

By construction similar properties hold for multiwavelets. If the piecewise polynomial functions, ψ0 , . . . , ψk−1 form an orthonormal basis for Wk0 , then by dilation and translation the space Wnk is spanned by 2n k functions

  ψ nj = 2n/2 ψ j 2n x −  ,

j = 0, . . . , k − 1,  = 0, . . . , 2n − 1.

(2.7)

The relations (2.2) and (2.3) between the subspaces may be expressed by the two-scale difference equations,

φi (x) =

k −1 √  

( 0)

(1 )



h i j φ j (2x) + h i j φ j (2x − 1) ,

2

i = 0, . . . , k − 1 ,

j =0

ψi (x) =

k −1 √  

2

( 0)

(1 )



g i j φ j (2x) + g i j φ j (2x − 1) ,

i = 0, . . . , k − 1 ,

(2.8)

j =0

where the coefficients can be determined through Gauss–Legendre quadrature, by multiplying and integrating Eqs. (2.8) by specific orthogonal basis functions to obtain, k −1



k −1





1  xm ( 0) φ j (xm ), hi j = √ w m φi 2 2 m =0



1  xm + 1 (1 ) φ j (xm ), hi j = √ w m φi 2 2 m =0

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

k −1 1 

( 0)

gi j = √

2 m =0

 w m ψi

k −1

xm

881



2

φ j (xm ),





1  xm + 1 (1 ) gi j = √ w m ψi φ j (xm ), 2 2 m =0

(2.9)

for the Gauss–Legendre nodes, x0 , . . . , xk−1 , and weights, w 0 , . . . , w k−1 . Using the two-scale difference equations and coefficients of Eqs. (2.8), the multiwavelet basis function can be generated. The Legendre polynomials P 0 , . . . , P k−1 are used on [−1, 1], to construct an orthonormal basis for Vk0 . Specifically, for j = 0, . . . , k − 1, we define the scaling functions as

√

2 j + 1P j (2x − 1),

φ j (x) =

if x ∈ [0, 1],

0,

(2.10)

otherwise.

We use the algorithm given in [1] to construct the wavelets basis with the choice of vanishing moments. Given a function, f ∈ Vnk , it can be represented by an expansion of scaling functions

f h (x) =

n k −1 2 −1 

snj φ nj (x),

(2.11)

=0 j =0

where the coefficients snj are n 2− (+1)

snj =

f (x)φ nj (x) dx.

(2.12)

2−n l

The function f (x) has a multiwavelet expansion, derived by multiple applications of the two scale relationship, given by

f h (x) =

k −1 



n−1 2 −1 



m

s0j0 φ j (x) +

m dm j  ψ j  (x)

(2.13)

,

m=0 =0

j =0

with the coefficients n 2− (+1)

dm j

f (x)ψ nj (x) dx.

=

(2.14)

2−n l

Using the two-scale difference Eqs. (2.8), the following relations between the coefficients on two consecutive levels m and m + 1 can be generated,

sm i =

k −1  



+1 1 , h i j sm + h i j smj ,+ j ,2 2+1 ( 0)

(1 )

j =0

dm i =

k −1  



+1 1 , g i j sm + g i j smj ,+ j ,2 2+1 ( 0)

(1 )

(2.15)

j =0

and +1 sm = i ,2

k −1  

( 0)

( 0)



m h ji sm j , + g ji d j , ,

j =0

+1 dm = i ,2+1

k −1  

(1 )

(1 )



m h ji sm j , + g ji d j , .

(2.16)

j =0

These equations are used to create fast transforms between Eqs. (2.11) and (2.13), and can be used recursively to refine the function approximation to arbitrary but finite resolution. 3. Multiwavelet discontinuous Galerkin for convection We review the DG method [5] and show how the multiwavelet scaling function can be merged into the DG method. This section will provide a brief introduction to the DG method for convection dominated problems and present the notation necessary to explain a multiwavelet DG method. Consider the one-dimensional scalar non-linear convection equation,

u t + f ( u ) x = 0.

(3.1)

When solutions of Eq. (3.1) become complicated or the time of evolution becomes large, it is necessary to utilize high-order numerical methods to obtain acceptable accuracy. Advances in finite difference, finite element, and finite volume methods

882

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

have produced various high-order schemes for each of these numerical techniques. There are pros and cons associated with these high-order methods that are well explained in [18]. This study focuses on the DG method, which is a type of finite element method that retains local conservation by incorporating the concepts of numerical fluxes and limiters developed for high-order finite difference and finite volume schemes. The DG method is flexible in the choice of mesh and basis functions used to approximate the solution of Eq. (3.1). Specifically setting the DG solution approximation space of Eq. (3.1) to the space of piecewise polynomials as defined in Eq. (2.1), it is possible to seamlessly incorporate Multiwavelets into the DG method. Given a fixed order k  0 and resolution n  0 and assuming a solution of the form u h (x, t ) ∈ Vnk , the variational formulation of the DG method is derived by multiplying Eq. (3.1) by the test functions φ j  ∈ Vnk and integrating by parts to obtain











u h (x, t ) t φ nj (x) dx =

I n







f u h (x, t ) φ nj (x)

x









dx − ˆf u h 2n ( + 1), t φ nj 2n ( + 1)

I n

     + ˆf uh 2n , t φ nj 2n  ,

(3.2)

for j = 0, . . . , k − 1 and  = 0, 1, . . . , 2n − 1, and ˆf (u h ) is a monotone numerical flux. Here we adopted notation standard to multiwavelets to express the numerical DG scheme. It should be noted that the numerical scheme as stated in (3.2) has both space and time dimension, but also added a scale dimension without increasing the degrees of freedom in the DG scheme. Using the two scale relationship Eq. (2.15) on the test functions in Eq. (3.2) we can re-write the numerical DG scheme as

1



1



u h (x, t ) t φ j (x) dx =

0

and











f u h (x, t ) φ j (x)

x



x=1 dx − ˆf u h (x, t ) φ nj (x) x=0 ,

(3.3)

0



 u h (x, t ) t ψ m j  (x) dx =





  f u h (x, t ) ψ m j  (x) x dx +

+1 Im 2

Im 









f u h (x, t ) ψ m j  (x)

x

dx

+1 Im 2+1

x=2m−1 (+1)

x=2m (+1)    

− ˆf uh (x, t ) ψ m − ˆf uh (x, t ) ψ m j  (x) x=2m  j  (x) x=2m−1 (+1) ,

(3.4)

for j = 0, . . . , k − 1 and  = 0, 1, . . . , 2 − 1 and m = 0, . . . , n − 1. Eq. (3.4) is broken up into the appropriate intervals where m the wavelet ψ m j  is smooth. Since multiwavelet functions are discontinuous, evaluation of the wavelet ψ j  for the integration by parts terms are to be the limit of the wavelet as taken inside each corresponding open interval. Transformation of the numerical DG scheme into an equivalent formulation using multiwavelets as the test functions has two immediate advantages. First, this transformation isolates a set of equations in Eq. (3.3), where the numerical flux is directly known from the multiwavelet coefficients. Second, this transformation displays, instead of the standard notion of numerical fluxes only interacting with its nearest neighbors, how the numerical flux can have a multiscale structure that has specific interactions at various scales. Finally, we remark that the multiwavelet functions are piecewise smooth functions and the integrals are evaluated on each piecewise smooth section. Next we assume that the approximate solution in Eqs. (3.3) and (3.4) at (x, t ) have the multiwavelet expansion m

u h (x, t ) =

 k −1  j =0

u h (x, t + t ) =

n−1 2 −1 



m

s0j0 φ j (x) +

 k −1 

m dm j  ψ j  (x)

m=0 =0 n−1 2 −1  m

s˜ 0j0 φ j (x) +

,

˜dm ψ m (x) j

(3.5)

j

m=0 =0

j =0

with,

u h (x, t )t = lim

u h (x, t + t ) − u h (x, t )

t →0

(3.6)

t

and as a direct result of multiwavelets orthogonality, upon substitution we have the following system of equations,

lim

t →0

s˜ 0j0 − s0j0

t

1 0

and







f u h (x, t ) φ nj (x)

=



x



x=1 dx − ˆf u h (x, t ) φ nj (x) x=0 ,

(3.7)

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

lim



d˜ m − dmj j

t →0



  f u h (x, t ) ψ m j  (x) x dx +

=

t

+1 Im 2









f u h (x, t ) ψ m j  (x)

x

883

dx

+1 Im 2+1

x=2m−1 (+1)

x=2m (+1)    

− ˆf uh (x, t ) ψ m − ˆf uh (x, t ) ψ m j  (x) x=2m  j  (x) x=2m−1 (+1) ,

(3.8)

for j = 0, . . . , k − 1 and  = 0, 1, . . . , 2n − 1 and m = 0, 1, . . . , n − 1. 4. Multiwavelet local discontinuous Galerkin for convection–diffusion problems We demonstrate in this section how, similar to the previous section on convection, incorporating multiwavelets into the local discontinuous Galerkin (LDG) method [4] for convection–diffusion problems yields an adaptive multi-scale structured formulation. This section reviews the LDG method for convection–diffusion problems and present the notation necessary to explain a multiwavelet LDG method. Consider the one-dimensional convection–diffusion equation,



u t + f (u ) x = a (u )u x



in [0, 1] × [0, T ],

x

(4.1)

and a(u )  0. The idea of LDG is to formulate Eq. (4.1) as a first order system, then use appropriate numerical fluxes that guarantee stability and local solvability of the solution approximation to the system. We begin by using the LDG approach to re-write the convection–diffusion equation as the first order system,



ut + f (u )x = b(u )q



and q − B (u )x = 0,

x

(4.2)

where

b (u ) =





a (u )

u B (u ) =

and

b(s) ds.

(4.3)

Using the same techniques as Section 3, for a fixed order k  0 and resolution n  0 the numerical LDG scheme is given by,









I n

and





u h (x, t ) t φ nj (x) dx =



g u h (x, t ), qh (x, t ) φ nj (x)



x













dx − gˆ u h 2n ( + 1), t , qh 2n ( + 1), t φ nj 2n ( + 1)

I n

       + gˆ uh 2n , t , qh 2n , t φ nj 2n  ,



 qh (x, t )φ nj (x) dx = I n







B qh (x, t ) φ nj (x)

 

x

(4.4)







dx − Bˆ qh 2n ( + 1), t φ nj 2n ( + 1)

I n

     + Bˆ qh 2n , t φ nj 2n  ,

(4.5)

for













g u h (x, t ), qh (x, t ) = f u h (x, t ) − b u h (x, t ) qh (x, t ),

(4.6)

with j = 0, . . . , k − 1 and  = 0, 1, . . . , 2 − 1. Here the approximate LDG solution is u h , qh ∈ with the test functions given in Eq. (2.6). The numerical flux design follows the sufficient conditions described in [4] in order to guarantee the stability of the scheme. Finally, using the two scale relationship of Eq. (2.15), multiwavelets orthogonality, and assuming u h follows the multiwavelet expansion given in Eq. (3.5), then qh can be written as the multiwavelet expansion, n

qh (x, t ) =

 k −1 

n−1 2 −1 

Vnk



m

s0j0 φ j (x) +

qh (x, t + t ) =

(4.7)

,

m=0 =0

j =0

and

m dm j  ψ j  (x)

 k −1  j =0

s0j0 φ j (x) +

m n −1 2 −1 

m dm j  ψ j  (x)

,

m=0 =0

the LDG solution to Eq. (4.2) is given by the following system of differential equations,

(4.8)

884

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

lim

s˜ 0j0 − s0j0

t →0

t

1









g u h (x, t ), qh (x, t ) φ nj (x)

=

x



dx − gˆ u h (1, t ), qh (1, t ) φ nj (1)

0

  + gˆ uh (0, t ), qh (0, t ) φ nj (0),

lim

d˜ m − dmj j

t →0

t



(4.9)







g u h (x, t ), qh (x, t ) ψ m j  (x)

=



+1 Im 2







g u h (x, t ), qh (x, t ) ψ m j  (x)

dx + x

x

dx

+1 Im 2+1

x=2m−1 (+1)

x=2m (+1)    

− gˆ uh (x, t ), qh (x, t ) ψ m − gˆ uh (x, t ), qh (x, t ) ψ m j  (x) x=2m  j  (x) x=2m−1 (+1) ,

lim

s˜ 0j0 − s0j0

t →0

t

1









B qh (x, t ) φ nj (x)

=

x







dx − Bˆ qh (1, t ) φ nj (1) + Bˆ qh (0, t ) φ nj (0),

(4.10)

(4.11)

0

and

lim

d˜ m − dmj j

t →0

t











B qh (x, t ) ψ m j  (x)

=

x



x=2m−1 (+1)

dx − Bˆ qh (x, t ) ψ m j  (x) x=2m 

I ml

x=2m (+1)  

+ Bˆ qh (x, t ) ψ m j  (x) x=2m−1 (+1) ,

(4.12)

for j = 0, . . . , k − 1 and  = 0, 1, . . . , 2n − 1 and m = 0, 1, . . . , n − 1. Just as in the previous section, evaluation of the wavelets for the integration by parts terms are to be the limit of the wavelets as taken from the inside of each corresponding open interval. 5. Time discretization We use a method of time stepping that has been used in multiwavelet based numerical methods [2,3]. Given a differential equations of the form,

ut = Lu + N (u ),

(5.1)

where the system is split into a linear operator L and a nonlinear operator N , the equivalent integral equation is given as,

u (t ) = e

tL

t u0 +

e (t −τ )L N (u ) dτ .

(5.2)

0

The discretization of Eq. (5.2) is given in [2,3] and uses scaling and squaring methods to produce approximations to the exponential linear operators in Eq. (5.2). This time-stepping schemes is called an exact linear part (ELP) schemes. 6. Examples To demonstrate the properties of multiwavelet DG and LDG we will begin by focusing on Examples 6.1 and 6.2, which deal with linear convection and convection–diffusion problems, respectively. Next we consider a non-linear problem in Example 6.3. We end off this section with Examples 6.4 and 6.5, which is standard climate models test. For all examples throughout this paper we use the well-known simple Lax–Friedrichs flux [13]. Additionally, we emphasize that when we mention DG we refer to the solution of Eq. (3.2), where as multiwavelet DG is the solution of Eqs. (3.7) and (3.8). Similarly, when we mention LDG we refer to the solution of Eqs. (4.4) and (4.5), where as multiwavelet LDG is the solution of Eqs. (4.9) to (4.12). Example 6.1. For arbitrary dimension, d, consider the linear convection problem

u t + ∇ · u = 0,

in (0, 1)d × (0, T ),

for periodic boundary conditions and initial conditions u 0 (x) = u (0, t ).

(6.1)

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

885

Table 1 Convergence rates for DG and multiwavelet DG for Example 6.1. Here the initial conditions are u 0 (x) = sin(2π x), with k = 2 and k = 3 and error tolerance = 10e−8 . Resolution (n)

DG

Multiwavelet DG

L 2 error

Order

L 2 error

Order

k=2 4 5 6 7

1.1954e–02 1.1186e–03 1.0268e–04 1.0098e–05

– 3.42 3.45 3.35

1.1954e–02 1.1186e–03 1.0268e–04 1.0111e–05

– 3.42 3.45 3.35

k=3 4 5 6 7

4.7328e–04 2.0856e–05 9.3186e–07 5.8985e–08

– 4.50 4.48 3.98

4.7328e–04 2.0835e–05 9.4271e–07 5.0292e–08

– 4.51 4.47 4.22

Table 2 Convergence rates for LDG and multiwavelet LDG for Example 6.2. Here the initial conditions are u 0 (x) = sin(2π x), with k = 2 and k = 3 and error tolerance = 10e−8 . Resolution (n)

LDG

Multiwavelet LDG

L 2 error

Order

L 2 error

Order

k=2 4 5 6 7

1.1529e–03 1.4610e–04 2.2534e–05 3.8269e–06

– 2.98 2.70 2.56

1.1529e–03 1.4610e–04 2.2534e–05 3.8241e–06

– 2.98 2.70 2.56

k=3 4 5 6 7

4.1736e–05 2.5938e–06 2.0724e–07 2.0666e–08

– 4.01 3.65 3.33

4.1739e–05 2.5920e–06 2.0067e–07 1.7749e–08

– 4.01 3.69 3.50

Example 6.2. For arbitrary dimension, d, consider the linear convection–diffusion problem

u t + ∇ · u − ∇ · ∇ u = 0,

in (0, 1)d × (0, T ),

(6.2)

for periodic boundary conditions and initial conditions u 0 (x) = u (0, t ). We begin by analyzing the convergence properties of the proposed numerical methods for Examples 6.1 and 6.2. It can be seen in Tables 1 and 2 that the expected convergence rate of k + 1 is demonstrated for Example 6.1 and to a lesser degree for Example 6.2. We note that in Table 2 the convergence rate for the convection–diffusion problem falls short of the rate of k + 1 when the error tolerance is set at = 10e−8 . However, this convergence rate is restored when the error tolerance is lowered. The table reports the number of levels, n, used in the multiwavelet expansion to express the level of resolution, it should be noted that the number of elements for each simulation is 2dn . Finally, we note that the multiwavelet LDG is more robust with respect to error tolerance in maintaining its convergence rate for higher resolution and polynomial basis order. Fig. 1 displays the property of similar sparsity pattern for the one-dimensional forward and inverse DG linear operator of Example 6.1 and LDG operator of Example 6.2. The forward operator represents an explicit time step and the inverse operator represents an implicit time step using the ELP scheme discussed in Section 5. The inverse operators of Fig. 1 display a generous time step, of the order of the CFL condition. It can be seen that sparsity pattern of multiwavelet DG/LDG inverse operators are self-similar to their forward operators. Predictive sparse structure is a highly desired property in high-performance computing, allowing for the potential development of effective scalable computational algorithms. Table 3 shows the ratio of CPU time in the one-dimensional case for DG as compared to multiwavelet DG for Example 6.1 and LDG as compared to multiwavelet DG for Example 6.2. It can be seen that the benefits of using multiwavelets in DG increases strongly with the order and resolution of the method. The underlying phenomenon behind the speed-up of the multiwavelet DG method is the property that given a specified error tolerance, the number of multiwavelet coefficients needed to represent a function, as expressed by Eq. (2.13), is less than is needed for the local DG basis, as expressed by Eq. (2.11). The improvement is slightly higher for the convection–diffusion since the inverse operator is denser for Example 6.2 as compared to Example 6.1, and hence the compactness of the multiwavelet representation is amplified. Table 4 extends the analysis of the ratio of CPU times to consider what occurs in higher dimensions. It is demonstrated that for both Examples 6.1 and 6.2 the benefits of using multiwavelets in DG and LDG improve with the dimension of the problem and this improvement is amplified for higher orders. Again, the underlying phenomena behind the speed-up of

886

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

Fig. 1. (a) Forward and (b) inverse multiwavelet DG operator for Example 6.1 with n = 8, k = 3, and d = 1 with 41 179 and 30 726 non-zero elements respectively. (c) Forward and (d) inverse LDG operator for Example 6.2 with n = 8, k = 3, and d = 1 with 47 772 and 37 920 non-zero elements respectively. The time step is t = 21n k . Table 3 Ratio of CPU time in the one-dimensional case for DG as compared to multiwavelet DG for Example 6.1 and LDG as compared to multiwavelet LDG for Example 6.2. Here the initial conditions are u 0 (x) = sin(2π x), with error tolerance = 10e−8 . Order (k)

Resolution (n) Example 6.1

1 2 3 4 5

Example 6.2

7

8

9

7

8

9

0.86 0.98 1.19 2.03 2.95

0.98 1.04 1.36 3.75 5.61

1.07 1.39 2.31 4.58 7.36

0.86 0.99 1.33 2.91 4.28

0.98 1.04 1.52 4.00 6.07

0.96 1.25 2.36 5.68 9.14

Table 4 Ratio of CPU time with increasing dimensions for DG as compared to multiwavelet DG for Example 6.1 and LDG as compared to multiwavelet LDG for d Example 6.2. Here the initial conditions are u 0 (x) = i =1 sin(2π xi ), with error tolerance = 10e−8 . Resolution (n)

Order (k)

Dimension (d) Example 6.1

5 3 2

1 2 3

Example 6.2

1

2

3

1

2

0.33 0.65 1.04

0.32 0.82 3.96

0.36 1.12 9.22

0.37 0.67 1.08

0.37 0.88 4.81

3 0.38 1.17 11.28

the multiwavelet DG and LDG method is because the compactness of the multiwavelet solution representation relative to the local DG basis. As seen by Tables 3 and 4, the benefit of using multiwavelets consistently occurs for k > 2 for all tested values of resolution. Next we consider the viscous Burgers’ equation. Example 6.3. The one-dimensional Burger’s equation given by,

ut − ν u xx + uu x = 0,

in (0, 1) × (0, T ),

(6.3)

for periodic boundary conditions has the exact solution given as,

u = −2ν for

φx (x − ct , t + τ ) φ(x − ct , t + τ )

(6.4)

τ > 0, with φ(x, t ) =

n =∞

2 e −(x−n) /4ν t .

(6.5)

n=−∞

This problem displays the natural adaptivity of multiwavelet DG. Fig. 2(b) displays the multiwavelet DG solution of Ex1 . For this solution, a third order explicit ELP scheme is used with c = 4, ν = 101π , τ = 21π , k = 3, and ample 6.3 for T = 16

error tolerance = 10e−8 . Fig. 2(c) displays the adaptivity pattern of the multiwavelet coefficients at T = 0, where the blue dot represents a group of k coefficients and the red line shows the support of the associated basis functions. For instance, at level zero there is one blue dot that represents the multiwavelet coefficients {d000 , d010 , d020 } that are used in the solution

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

887

1 Fig. 2. Multiwavelet DG solution of Example 6.3 for (a) T = 0 and (b) T = 16 . Adaptivity pattern of the multiwavelet coefficients at (c) T = 0 and (d) T = (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Table 5 Convergence rates for DG and multiwavelet DG for Example 6.3, along with ratio of CPU time. A third order explicit ELP scheme is used with c = 4,

τ=

1 2π

, final time t =

1 , 16

k = 3, and error tolerance

Resolution (n)

5 6 7

= 10e−8 .

LDG

Multiwavelet LDG

ν=

1 . 16

1 10π

,

Ratio

L 2 error

Order

L 2 error

Order

5.5290e–04 4.9922e–05 3.2172e–06

– 3.47 3.96

5.5290e–04 4.9922e–05 3.2168e–06

– 3.47 3.96

0.91 1.07 1.22

0 0 0 approximation, where the one red line represents the support of the associated basis functions, {ψ00 , ψ10 , ψ20 }. Skipping to level two, we can see that out of the four possible groups that can be used in reconstruction only the coefficients {d202 , d212 , d022 } and {d203 , d213 , d023 } are used. The two lines on level two, where red circles mark the starting and finishing positions, depict the domain of the basis functions associated with the coefficients on this level. The number of possible groups of coefficients doubles for each increase in level. However, it can be seen in Fig. 2(c) that more multiwavelet coefficients are focused near the discontinuity of the solution. Fig. 2(d) shows that the adaptivity changes with the movement of the solution discontinuity. Finally, Table 5 shows that favorably properties of multiwavelet DG demonstrated for the previous linear problems are maintained for this non-linear example. We end this section with a couple of problems that has specific importance to the development of climate models, which are convection and the shallow water equations on the sphere. The purpose of showing these examples is two fold. First, the multiwavelet DG method is capable of challenging non-linear problems. Second, the major constraint in the multiwavelet DG, which is the restriction to quadrilateral elements, can be overcome and the method expanded to different geometries by developing appropriate transforms.

Example 6.4. Given the advecting field, h, the equation for convection on a sphere in flux form is,

∂h + ∇ · (hv) = 0. ∂t

(6.6)

888

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

Fig. 3. (a) Red points on the sphere, transfered to the blue points on the cube forms the basis in developing the so-called ‘cubed-sphere’ model for problems such as Example 6.4. (b) The advecting wind on the sphere of Eq. (6.8) for α = π4 , and (c) the advecting wind on the sphere transforms to the cube surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. (a) Initial conditions (6.7) on the cube surface and (b) the relative error after one complete revolution, with k = 3 and 364 elements.

The first test in the standardized suit developed by the climate modeling community [20] is to solve Eq. (6.6) with initial conditions



h0 (λ, θ) =

1000 cos(3π r ),

if r < 13 ,

0,

else,

(6.7)

for r (λ, θ) = arccos(cos(θ) cos(λ − 32π )) and advecting wind



v=

cos(θ) cos(α ) + sin(θ) cos(λ) sin(α ) − sin(λ) sin(α )



.

(6.8)

The design of local methods to solve problems such as Example 6.4 have received much attention since they have the potential of exploiting high-performance computing with much more efficiency than current climate models. In this paper we use the DG method developed by [16] to solve Example 6.4. We describe the major features of the DG model and refer the reader to [16] for the extensive details. The goal of this example is to demonstrate how the incorporation of multiwavelet DG into this problem has a favorable impact. Fig. 3(a) depicts the geometry used to solve Example 6.4, an inscribed cube is used to develop cube-to-sphere and sphere-to-cube transformations. The convection problem is transformed onto the cube and [16] shows that the DG method used on each face of the cube with appropriate coupling are a very effective local method to solve Example 6.4. Fig. 3(b) depicts the most computationally challenging direction for the advecting wind on the sphere (see Eq. (6.8)), which is when α = π4 . Fig. 3(c) depicts how the advecting wind on the sphere transforms to the cube surface. Fig. 4(a) depicts the initial conditions (6.7) on the cube surface. The direction of the advection wind takes the solution h through four corners and along two cube edges for every rotation around the cube. Fig. 4(b) depicts the final multiwavelet DG solution after one rotation. In this simulation the order was k = 3, the number of elements was 364, and the error tolerance was 10e−8 . In comparison to the DG solution, the multiwavelet DG ran two times faster. The major contributor to

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

889

Fig. 5. Shallow water equations on the sphere at day 15, with conditions given in Example 6.5 for k = 3 and 364 elements. Table 6 Ratio of CPU time as compared to multiwavelet DG for Examples 6.4 and 6.5, with error tolerance

= 10e−8 .

Resolution (n)

Order (k)

Example 6.4

2 3

Example 6.5

2

3

4

2

3

4

0.43 0.45

0.42 0.52

0.46 0.69

0.97 1.07

1.03 1.28

1.38 1.75

the speed-up of multiwavelet DG stems from the fact that for the specified error tolerance, the initial conditions (6.7) can be expressed with three times fewer coefficients using multiwavelets over the local DG basis. Example 6.5. The shallow water equations on the sphere is given by,

∂ h∗ v + ∇ · (vhv) + f kˆ × h∗ v + gh∗ ∇ h = 0, ∂t   ∂ h∗ + ∇ · h∗ v = 0 ∂t

(6.9)

where f is the Coriolis parameter, h∗ is the depth of the fluid, kˆ is the unit vector in the z-direction, and h is the height of the free surface, or h = h∗ + h s with h s the height of the earth surface. The vector velocity is denoted as v with components u in the longitudinal (λ) direction and v in the latitudinal (θ) directions. Operators are given as,

∇( ) ≡ and

∇ ·v≡

i

∂

r cos θ ∂λ

+

j ∂

(6.10)

r ∂θ

 ∂ u ∂ v cos θ + , r cos θ ∂λ ∂θ 1



(6.11)

with r being the radius of the earth. The fifth test in the standardized suit developed by the climate modeling community [20] is to solve Eq. (6.9) with initial conditions

h = 5400 +

2r Ω + 1 2g

(− cos λ cos θ sin α + sin θ cos α )2 ,

(6.12)

and

 h s = 2000 1 −

2 9 min[ π89 , (λ − 32π )2 + (θ − π6 )2 ]

π

(6.13)

,

where Ω is the rotational rate of the earth and the vector velocity is given in Eq. (6.8) for

α = 0.

Fig. 5 depict the shallow water equations on the sphere at day 15, with conditions given in Example 6.5 for k = 3 and 364 elements. We used a third order explicit ELP scheme using a time step of 22.5 mins. The results compare favorable to other published results [17]. Table 6 gives the ratio of CPU time as compared to multiwavelet DG for Examples 6.4 and 6.5, with error tolerance = 10e−8 . Example 6.4 is representative of a special case where it is not advantageous to use multiwavelets. The solution of Example 6.4 is almost zero everywhere, and for this case the number of significant coefficients needed to represent this solution is almost the same for both DG and multiwavelet DG methods. The added overhead multi-scale structure of multiwavelets cannot be offset in this case. However, Example 6.5 is a more complex solution and it can be seen in Table 6 that multiwavelets provide a distinct advantage in computation.

890

R. Archibald et al. / Applied Numerical Mathematics 61 (2011) 879–890

7. Conclusions This research demonstrates that merging the DG method with the advances in multiwavelets provides an adaptive multiscale DG method. We use convection and convection–diffusion problems in multiple dimensions to highlight the partial decoupling across different length scales that occurs in multiwavelet DG. One current restriction of using multiwavelet DG is the method is only developed for the quad-mesh, extension to triangle-mesh should be possible using work done on the multiwavelet basis by [7]. Finally, all of the favorable properties of multiwavelet DG hold, and even are enhanced, for higher-dimensional problems with large degrees of freedom. Acknowledgements The submitted manuscript has been authored by a contractor [UT-Battelle, manager of Oak Ridge National Laboratory (ORNL)] of the U.S. Government under Contract No. DE-AC05-00OR22725. Accordingly, the U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes. References [1] B. Alpert, A class of bases in L 2 for the sparse representation of integral operators, SIAM J. Math. Anal. 24 (1) (1993) 246. [2] B. Alpert, G. Beylkin, D. Gines, L. Vozovoi, Adaptive solution of partial differential equations in multiwavelet bases, Journal of Computational Physics 182 (1) (2002) 149. [3] G. Beylkin, J.M. Keiser, L. Vozovoi, A new class of stable time discretization schemes for the solution of nonlinear PDEs, Journal of Computational Physics 147 (1998) 362. [4] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection diffusion systems, SIAM Journal on Numerical Analysis 35 (1998) 2440. [5] B. Cockburn, C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, Journal of Scientific Computing 16 (3) (2001) 173. [6] A. Cohen, R. Masson, Wavelet methods for second order elliptic problems preconditioning and adaptivity, preprint R 97036, Univ. P. et M. Curie, Paris, 1997. [7] N. Coult, Introduction to Discontinuous Wavelets, in: B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory Computation and Applications, Springer-Verlag, February 2000. [8] S. Dahlke, W. Dahmen, R. Hochmuth, R. Schneider, Stable multiscale bases and local error estimation for elliptic problems, Applied Numerical Mathematics: Transactions of IMACS 23 (1) (1997) 21–47. [9] J. Daz Calle, P. Devloo, M. Gomes, Wavelets and adaptive grids for the discontinuous Galerkin method, Numerical Algorithms 39 (1) (2005) 143. [10] R. Glowinski, W. Lawton, M. Ravachol, E. Tenenbaum, Wavelet solutions of linear and nonlinear elliptic, parabolic and hyperbolic problems in one space dimension, in: Computing Methods in Applied Sciences and Engineering, vol. 55, SIAM, Philadelphia, 1990. [11] A. Harten, Multiresolution algorithms for the numerical solution of hyperbolic conservation laws, Comm. Pure Appl. Math. 48 (1995) 1305. [12] W. Lawton, W. Morrell, E. Tenenbaum, J. Weiss, The wavelet-Galerkin method for partial differential equations, Technical Report, Aware, Inc., 1990. [13] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhauser Verlag, Basel, 1990. [14] S. Malla, W.L. Hwang, Singularity detection and processing with wavelets, Technical Report 549, Courant Institute of Mathematical Sciences, New York University, March 1991. [15] C. Meneveau, Analysis of turbulence in the orthonormal wavelet representation, J. Fluid Mech. 323 (1991) 469. [16] R. Nair, S. Thomas, R. Loft, A discontinuous Galerkin transport scheme on the cubed sphere, Monthly Weather Review 133 (4) (2005) 814. [17] R. Nair, S. Thomas, R. Loft, A discontinuous Galerkin global shallow water model, Monthly Weather Review 133 (4) (2005) 876. [18] 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 (2) (2003) 107–118. [19] C.-W. Shu, Discontinuous Galerkin methods: general approach and stability, in: Numerical Solutions of Partial Differential Equations, in: Adv. Courses Math. CRM Barcelona, Birkhäuser, Basel, 2009, pp. 149–201. [20] D. Williamson, J. Hack, R. Jakob, P. Swarztrauber, J. Drake, A standard test set for numerical approximations to the shallow water equations in spherical geometry, Journal of Computational Physics 102 (1992) 211.