- Email: [email protected]
Two adaptive wavelet algorithms for non-linear parabolic partial differential equations
Computers & Fluids 31 (2002) 467–480 www.elsevier.com/locate/compfluid
Two adaptive wavelet algorithms for non-linear parabolic partial differential equations q G. Chiavassa a, M. Guichaoua b, J. Liandrat a,* a b
LATP and ESM2, Technop^ ole de Ch^ ateau-Gombert, 13451 Marseille Cedex 20, France IRPHE and ESM2, Technop^ ole de Ch^ ateau-Gombert, 13451 Marseille Cedex 20, France Received 1 August 2001
Abstract This paper is concerned with the construction of wavelet based adaptive algorithms for the numerical resolution of evolution equations. The adaptivity is applied into two complementary directions. The first direction shares the approaches involved in classical adaptive finite element methods and is related to a solution dependent definition of spaces of approximation. The second direction is related to the approximation of evolution operators that is made solution dependent following the philosophy of essentially nonoscillatory schemes. After the construction of the schemes, numerical tests are provided. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Wavelets; Adaptive methods; Non-linear PDEs
1. Introduction Adaptive methods and related algorithms have been a subject of research for many years in various domains of numerical analysis. In this paper, we are concerned by adaptivity considered in two complementary meanings applied to the numerical resolution of the following equation: 8 ot U þ GðU Þ ¼ mDU ; > > < ðx; tÞ 2 0; 1½ ½0; T ; ð1Þ boundary conditions at x ¼ 0 and x ¼ 1; > > : U ðx; 0Þ ¼ gðxÞ; q
This work has been supported in part by TMR Research Network Contract FMRX-CT98-0184. Corresponding author. E-mail addresses: [email protected] (G. Chiavassa), [email protected] (M. Guichaoua), liandrat@ esm2.imt-mrs.fr (J. Liandrat). *
0045-7930/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 6 1 - 5
468
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
~ U , g is a given function and m is a real constant. where G is a non-linear operator of type U 7! G As usual, t stands for the time variable and x for the spatial variable. The first meaning we give to adaptivity is the one usually given in finite elements methods where the space of approximation into which the solution is computed is solution dependent. The space is therefore said to be adapted. This leads to non-linear methods in contrast with linear methods in which the space of approximation is set in advance. Such non-linear methods allow to improve the approximation efficiency since a finer resolution can be involved only near singularities or strong gradients of the solution. These methods are efficient as soon as fast algorithms are available on the adapted spaces of approximation. The second meaning we give to adaptivity in this paper is the one used in the framework of finite volume methods for conservation law problems. Adaptivity there, leads to solution dependent schemes also called non-linear schemes. Such methods allow the approximation of singular solutions with a significant delay in the occurrence of Gibbs-like oscillations. We present in this paper, two algorithms based on wavelets. On one hand, wavelet spaces offer the possibility to compress very efficiently smooth functions with isolated singularities or strong gradients regions. If some specific property is imposed to the compressed space of approximation (see Definition 2), the fast algorithms related to the pyramidal structure connected to wavelets (see Section 2) remain available. This leads to the first adaptive algorithm presented in Section 3.1. This algorithm belongs to the family of wavelet algorithms developed in Marseille (France) since 1990 [5,14,15]. Full adaptivity and the capability of tracking singularities are emphasized in this paper. This method shares various properties with some other wavelet algorithms of the literature, in particular the ones developed by Fr€ ohlich and Schneider [9,10], by Charton and Perrier [4] or by Gomes [11]. A review on recently developed approaches using wavelets has been performed by Dahmen [6]. On the other hand, recent works by Arandiga et al. [1,2] have introduced the notion of nonlinear multi-resolution and corresponding algorithms. The spirit of these works is a solution dependent definition of details, according to a solution dependent interpolation procedure belonging to the family of essentially non-oscillatory (ENO) schemes [13]. Following the same kind of ideas [12], we present in Section 4 a second adaptive algorithm where the non-linear operators are approximated using a solution dependent approach. This makes the algorithm solution dependent in a new way in comparison to the first algorithm and to other existing adaptive wavelet algorithms. In Section 5 we present numerical tests for the univariate flame (Section 5.1) and Burgers (Section 5.2) equations. 2. General setting In this section, we introduce the general setting in which our algorithms are developed and we recall the main results of wavelet theory necessary for the understanding of the sequel of the paper. 2.1. Wavelet setting We focus on the theory on the line since, by periodization, it can be easily transposed to the periodic framework.
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
469
Following Ref. [7] or Ref. [17] we call a biorthogonal wavelet basis of L2 ðRÞ any set of functions fwj;k ¼ 2j=2 wð2j x kÞ; k 2 Z; j 2 Zg;
ð2Þ
forming a Riesz basis. The dual basis is given by fw~j;k ¼ 2j=2 w~ð2j x kÞ; k 2 Z; j 2 Zg;
ð3Þ
and satisfies hwj;k ; w~j0 k 0 i ¼ dðj;kÞ;ðj0 ;k 0 Þ where h; i stands for the L2 scalar product. When w ¼ w~, the basis is orthonormal. Classically, we note Wj (resp. W~j ) the detail spaces with Wj ¼ spanfwj;k ; k 2 Zg (resp. W~j ¼ spanfw~j;k ; k 2 Zg) and Vj (resp. V~j ) the approximation spaces of the multi-resolution analysis associated to w (resp. w~). We have Vj ¼ spanf/j;k ; k 2 Zg (resp. V~j ¼ spanf/~j;k ; k 2 Zg) where / (resp. /~) is the corresponding scaling function. The following properties should also be translated with : (1) We assume that the spaces Vj offer an order m approximation, i.e. kf PVj f k2 ¼ Oð2 jm Þ;
8 f 2 H m;
ð4Þ
where PVj denotes the orthogonal projection onto Vj and H m stands for the Sobolev space of order m 2 N. The order of approximation (m) is connected to the number of vanishing moments of the dual wavelet. More precisely we have: Z ð5Þ xa w~ðxÞ dx ¼ 0; a ¼ 1; . . . ; m 1: (2) The space L2 ðRÞ can be expanded as an approximation space plus a sum of detail spaces, i.e. L2 ðRÞ ¼ Vj0 þ
þ1 X
ð6Þ
Wj ;
j¼j0
where j0 is any fixed level. Therefore, the family f/j0 ;k ; wj;k ; j P j0 ; k 2 Zg is a Riesz basis of L2 ðRÞ. (3) The embedding relations Vj Vjþ1 and Wj Vjþ1 induce the following scaling and detail relations: X X hk /ð2x kÞ and wðxÞ ¼ gk /ð2x kÞ ð7Þ /ðxÞ ¼ k
k 2
with fgk gk2Z , fhk gk2Z 2 l ðRÞ. The sequences fhk gk2Z , fgk gk2Z are called the filters of the multiresolution. They are the key ingredients of the so-called wavelet decomposition, i.e. the mapping from the set of scaling coefficients fcj0 ;k ¼ hf ; /j0 ;k i; k 2 Zg to the set of wavelet coefficients fdj;k ¼ hf ; wj;k i; j 6 j0 1; k 2 Zg. (4) We also require that the scaling function / gets vanishing moments up to order n/ , more precisely: Z ð8Þ xa /ðxÞ dx ¼ 0; a ¼ 1; . . . ; n/ : (5) In all the paper, we will call dyadic point any real of the form k ¼ k2 j , k 2 Z, j 2 Z and we will identify k to the unique couple ðj; kÞ satisfying k2 j ¼ k with minimal value of j and with the convention 0 ¼ ð0; 0Þ.
470
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Remark 1. We always consider that the number of non-zero coefficients of the filters is finite and equal to lg. This may be true only numerically, i.e. up to a given precision. In the sequel of the paper, lg will be referred as the symbol length.
3. Adaptive space of approximation and corresponding method In this section we present an adaptive resolution method for Eq. (1) using the compression properties of wavelet bases. The resulting algorithm is fully non-linear in the sense that the corresponding spaces of approximation are completely solution dependent and fully adaptive in the sense that its complexity, in term of arithmetic operations and also of memory storage, is proportional to the number of significant wavelet coefficients of the approximation. These properties are not shared by other existing and implemented wavelet algorithms. 3.1. Definition of adapted spaces We begin by introducing the notion of cone-structured spaces of approximation. Definition 2. The space C is said to be cone-structured if 9c; j0 and I 2 N such that C ¼ Vj0 þ
I 1 X
Xj ;
j¼j0
where Xj Wj and satisfies the following property: 8 j; j0 þ 1 6 j 6 I 1; wj;k 2 Xj ) fwj 1;k 0 ; with j2k 0 kj 6 cg Xj 1 : The space C will be a simple cone-structured space if moreover 8 j; j0 þ 1 6 j 6 I 1; 9k1 ; k2 2 N with Xj ¼ fwj;k ; k1 6 k 6 k2 g; where k1 and k2 are functions of j. It will become clear (see Section 3.2) that such specific geometrical conditions imposed to the space of approximation lead to efficient numerical algorithms. From now we describe our method in the case of a simple cone-structured space. The final Remark 5 is devoted to the case of more general cone-structured spaces. We now define the notion of e-adapted space. Definition 3. Given a function U and e > 0, a cone-structured space C is e-adapted to U if: 1. 8 j P I kPWj ðUÞk2 6 ekU k2 , 2. 8 wl;k 2 Vj and wl;k 62 C : jhU ; wl;k ij 6 ekU k2 . Following Ref. [5], we use an algorithm that enforces the space of approximation to remain cone-structured and e-adapted to U assuming some regularity conditions on U. This algorithm
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
471
never uses the approximation of U on the regular space VI but constructs adaptively the conestructured space according to the previous definition. Our goal is now to define an efficient numerical approximation of Eq. (1) on cone-structured spaces. 3.2. Discretization of the equation The discretization sequence we use is a discretization in time followed by a space approximation of each quantities (see Remark 4). 3.2.1. Time discretization Involving the semi-group generated by the Laplacian D, the exact solution of Eq. (1) is Z t mtD emðt sÞD GðU; sÞ ds: U ðx; tÞ ¼ e gðxÞ
ð9Þ
0
Classically, we define an increasing sequence ðtn Þ0 6 n 6 N of ½0; T , with t0 ¼ 0, tN ¼ T and the related time step dtn defined by tnþ1 ¼ tn þ dtn . For notation convenience we will drop the index n for the time step, but, in the numerical experiments, the time step may be n dependent. Using Eq. (9) between tn and tnþ1 leads to: Z dt m dtD U ðx; tnþ1 Þ ¼ e U ðx; tn Þ emðdt sÞD GðU; tn þ sÞ ds: ð10Þ 0
The goal of the time discretization procedure is to give a scheme which allows to compute an approximation of U ðx; tnþ1 Þ, denoted by U nþ1 ðxÞ, knowing an approximation of Uðx; tj Þ, j 6 n. The homogeneous part of Eq. (10) is computed directly without approximation. In Section 3.3 we show how to apply the heat operator to a function expanded on a cone-structured space. For clarity reasons, we simply write SðdtÞ ¼ em dtD . To compute the integral part of Eq. (10), a classical quadrature formula of order g is used and is defined by the weight coefficients fqi ; 0 6 i 6 gg and the quadrature points ftni ; 0 6 i 6 gg. We then obtain, for tni ¼ tn i dt, the semi-discretized form of Eq. (1): U
nþ1
n
ðxÞ ¼ SðdtÞU ðxÞ dt
g X
qi Sðði þ 1Þ dtÞGðU n i ðxÞÞ:
ð11Þ
i¼0
3.2.2. Spatial approximation In this section, the different terms of Eq. (11) are approximated on a finite dimensional space. We introduce the time dependent cone-structured space V n , that is supposed to be e-adapted to U n ðxÞ, and we note Uen the representation of U n ðxÞ on this space. Defining Gnn i as an approximation of G from V n i to V n , the fully discretized form of Eq. (1) writes: Uenþ1 ¼ PV nþ1 SðdtÞUen dt
g X i¼0
n i qi PV nþ1 Sðði þ 1Þ dtÞGnþ1 n i ðUe Þ:
ð12Þ
472
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Remark 4. Thanks to the specific discretization of the equation, the time evolution and projection operators are decoupled. This allows to adapt the space of approximation at each time step after applying the evolution operators. This property is important with regards to the spatial adaptivity of our algorithm and is used in practical computations. 3.3. Application of the linear operator From relation (12) and definition of the cone-structured space V n , we have to compute Sði dtÞUen ¼
j0 1 2X
cj0 ;k Sði dtÞ/j0 ;k þ
I 1 X X
dj;k Sði dtÞwj;k
ð13Þ
j¼j0 wjk 2Xj
k¼0
with cj0 ;k ¼ hUen ; /j0 ;k i and dj;k ¼ hUen ; wj;k i. Denoting sij0 ;k ¼ Sði dtÞ/j0 ;k and hij;k ¼ Sði dtÞwj;k , we know from Ref. [17] that these functions (called vaguelettes) and their Fourier transforms are uniformly localized with regards to j. This implies (for instance for hij;k ) that, on one hand, 8 e; 9ðr; r0 Þ 2 N2
such that sij0 ;k 2e Vj0 þr and hij;k 2e Vjþr0 ;
ð14Þ
where f 2e F means that 9 g 2 F such that kf gkL2 6 e and that, on the other hand, 9d 2 N such that jhhij;0 ; wj0 ;k ij 6 e if jkj > d:
ð15Þ
The link between r, r0 and e can be evaluated numerically. For instance for spline wavelets of order m ¼ 6 and e ¼ 10 6 , we use in our applications r ¼ 2, r0 ¼ 3 and d 6 50. Similar relations are valid for the other scalar products. It follows that the estimate of PV nþ1 Sði dtÞ can be performed efficiently using fast numerical convolutions. Using the properties of wavelets and of the heat operator, an efficient algorithm for the computation of the scalar products satisfying relation (15) has been derived in Ref. [5]. This leads to the computation of the homogeneous part of Eq. (12) with a method of complexity proportional to #V nþ1 , the dimension of V nþ1 . 3.4. Evaluation of the non-linear term n i We recall in this section the method we use to evaluate the quantities Gnþ1 n i ðUe Þ involved in relation (12). It has been derived in Ref. [5] and can be considered as a variation around the algorithm constructed in Ref. [15]. To simplify the notations, we only detail the case i ¼ 0. We first associate to V n (supposed to be a simple cone-structured space), a non-regular grid of points G, defined by (see Section 2 for the notations):
G ¼ fk; such that j ¼ j0 or j0 < j < I 1 and 2k1 ðj 1Þ lg =2 6 k 6 2k2 ðj 1Þ þ lg =2g; where lg is the symbol length (see Remark 1). It must be emphazised that there exists an algorithm for the computation of the points values of Uen on G from its representation on V n (and reversely) which complexity is proportional to #V n . This algorithm, developed and detailed in Ref. [5], is
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
473
similar to the classical tree-algorithm of Mallat [16], but acts only on the space V n and never requires the regular space VI . The evaluation of the non-linear terms can be described as follows, assuming that GðUÞ ¼ ~ U: G fhU ; wj;k ig 7! fhU ; /k ik2G g; fhU ; /k ig 7! fU ðkÞg; ~ðU ðkÞÞg; fU ðkÞg 7! fG ~ðU ðkÞÞg 7! fhGðU Þ; /k ig; fG fhGðUÞ; /k ig 7! fhGðU Þ; wj;k ig; where 7! stands for a mapping or its numerical approximation. Remark 5. We have presented all the steps of our adaptive algorithm for a simple cone-structured space. For practical computations a global space of approximation is generally built as the union of different simple cone-structured spaces (each of them being, for instance, associated to a specific strong gradient region of the solution). The application of the linear operators (heat operator and adaptive tree algorithm) is then performed independently on each simple space. The evaluation of non-linear term is slightly modified computing the values of Gnþ1 on the union of all the grid points G. Nevertheless, the complexity of the global algorithm remains proportional to #V n .
4. Adaptive approximation of evolution operator and corresponding method This section is devoted to the description of the second algorithm where adaptivity is related to a solution dependent approximation of the operators. The approach developed in this section is then complementary to the one developed in Section 3 where adaption deals with the spaces of approximation. Since the motivation is essentially given by the approximation of the non-linear part of Eq. (1) we focus on an equation of type: 8 ot U þ GðU Þ ¼ 0; > > < ðx; tÞ 2 0; 1½ ½0; T ; ð16Þ boundary conditions at x ¼ 0 and x ¼ 1; > > : U ðx; 0Þ ¼ gðxÞ: 4.1. General discretization of the equation As in Section 3.2, and for the same reasons (see Remark 4), the discretization of the Eq. (16) is first performed on the time variable and secondly on the spatial variable. Since the goal, in this section is to derive a solution dependent approximation of the non-linear term, the time discretization we use is explicit. We first describe the general discretization and describe in Section 4.2 the solution dependent discretization.
474
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
4.1.1. Time discretization The discretization in time is performed using a finite difference approximation of order r of Taylor type that writes, once a segmentation (see Section 3.2) has been performed: dtn ½1 n dtns 1 ½s 1 n nþ1 n n ¼ U dtn GðU Þ þ ð17Þ G ðU Þ þ þ G ðU Þ ; U 2 r! where G½0 ¼ G and G½lþ1 ¼ Dx G½l G. Remark 6. Such Taylor schemes are not very popular since they involve the computation of extra terms (fG½l g0
2 r!
ð18Þ
The remaining difficulty is to compute the scaling coefficients of G½l ðUn Þ, 0 6 l 6 s 1 from those of U n . The general framework described in Section 3.4 is suitable but the occurrence of derivatives of high order of non-linear terms (see the right-hand side term of Eq. (18)) imposes to pay special attention to their approximation when the solution exhibits locally strong gradients. That explains why a solution dependent approximation has been developed. 4.2. Solution dependent approximation We come back to the evaluation of hG½l ðU n Þ; /j0 ;k i. If we assume that GðU Þ is a combination of ~ U and of derivatives, we have to evaluate the scaling coefficients of derivatives of terms of type G non-linear forms of U. This is known to be a critical issue for which various solutions have been proposed in particular in the framework of finite volume methods for conservation laws [13]. We ~ U by introducing data propose to mimic the so-called ENO scheme for the evaluation of G dependent scaling functions generating Vj0 . Depending on U, we define specific functions /\j0 ;k , introduced in the algorithm (18) as follows: fhU ; /j0 ;k ig 7! fhU; /\j0 ;k ig; fhU ; /\j0 ;k ig 7! fU ðk2 j0 Þg; ~ðUðk2 j0 ÞÞg; fU ðk2 j0 Þg 7! fG ~ðUðk2 j0 ÞÞ; /\ ig; ~ðU ðk2 j0 ÞÞg 7! fhG fG j0 ;k ~ðU ðk2 j0 ÞÞ; /j ;k ig: ~ðU ðk2 j0 ÞÞ; /\ ig 7! fhG fhG j0 ;k 0
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
475
The construction of /\j0 ;k is made as follows: Assuming that U is smooth but in the vicinity of a point a ¼ ka 2 j0 , we write f/\j0 ;k g ¼ Ma f/j0 ;k g;
ð19Þ
where, each line of the infinite and invertible matrix Ma is an interpolation stencil of finite length. The matrix Ma (see Eq. (20) for an example) is constructed such that a minimal number of stencil intersect the point a in order to have /\j0 ;k ðaÞ ¼ 0 for a maximum of values of k. The construction of Ma is not unique and the analysis of an optimal construction is under study. Fig. 1(c) and (d) exhibits two functions /\0;k , corresponding to the matrix Ma presented in Eq. (20), for different values of k. These functions, in comparison to / are less oscillating and not symmetric. ka 0
1 0 B0 1 B B0 0 B B0 : B B0 1 10 B B0 0 B B0 : B B 0 : Ma ¼ B B0 : B B0 : B B0 : B B0 : B B0 : B B0 : B B0 : B @0 :
#
: 0 1 0 1 2
: 0 : : : : : : : : :
: : 0 1
1 :
101 0 : : : : : : : :
: : : 0 2 : 1 2
: : : :
12 :
1
101
1 2
: 0 : : : : : :
: : : : : : : :
: : : :
: : : : 1 0 10 : : 2 12
1 2 1
12 10 1 0 10 : 0 : : : : : : : : : :
: : : : : 0 1 10
12
2
12 : 0 : : : :
: : : : : : 0 1 10
1 2 : 1 10
: : : :
: : : : : : : 0 1 2
1 :
12 : : : :
: : : : : : : :
101 1 2
: 2 0 : : :
: : : : : : : : :
101 :
1 1 0 : :
: : : : : : : : : 0 :
: : : : : : : : : : 0 1 1
2 10 0 : 1 0 0 1 : 0
1 0 0C C 0C C 0C C 0C C 0C C 0C C 0C C 0C C: 0C C 0C C 0C C 0C C 0C C 0C C 1A
ð20Þ To end the estimate of hG ðU Þ; /j0 ;k i, the computation of the derivatives is performed using a method developed in Ref. [19] based on biorthogonal multi-resolution analyses associated to ow=ox. One of these scaling functions is plotted on Fig. 1(b). ½l
n
5. Numerical tests 5.1. The 1D thermo-diffusive flame equation In this section we apply the method described in Section 3 to the following system of equations: 8 < ot h þ v ox h X ¼ Dh; o C þ v ox C þ X ¼ Le1 DC; ð21Þ : t ðx; tÞ 2 0; L½ 0; T ½
476
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Fig. 1. (a) Initial spline scaling function / of order m ¼ 6, (b) scaling function associated to the operator o=ox, (c, d) scaling functions modified by the matrix Ma : /\0;k and /\0;k 0 for k 6¼ k 0 .
with X ¼ ðb2 =2LeÞC expðbðh 1Þ=ð1 þ cðh 1ÞÞÞ. This system stands for the propagation of a planar flame in the case of two reactive chemical components (hydrogen–air). The temperature hðx; tÞ, the hydrogen concentration Cðx; tÞ and the flame velocity v are the unknowns. The parameters Le, b, c are respectively the Lewis and Zeldovich numbers and the temperature ratio. They are set to the following classical values [8]: Le ¼ 1, b ¼ 10, c ¼ 0:8. When Le ¼ 1, the relation between h and C is linear and the system (21) reduces only to the first equation. The boundary conditions are hðx ¼ 0; tÞ ¼ 0 and hðx ¼ L; tÞ ¼ 1 and L ¼ 60 in our simulations. Following Ref. [8], it is possible to recover periodic boundary conditions, introducing the function sðxÞ ¼ 12ð1 þ tanhð10 tanðpðx 12ÞÞÞÞ satisfying sð0Þ ¼ 0 and sðþ1Þ ¼ 1 and solving the equation for the unknown h~ ¼ h s (we suppose that L is large enough so that sðLÞ ’ 1 at the precision of the computer). Writing this equation for x 2 0; 1½ gives: v 1 ot h~ þ ox h~ X ¼ 2 Dh~ þ F ðs; v; LÞ; L L
ð22Þ
where b2 ð1 h~ sÞ exp X¼ 2Le
! bðh~ þ s 1Þ ; 1 þ cðh~ þ s 1Þ
F ðs; v; LÞ ¼
1 v Ds ox s: 2 L L
The algorithm of Section 3 is applied using an orthonormal basis of spline wavelets of order 6, the thresholding parameter e set to 5 10 7 and the lower level is j0 ¼ 7. The velocity v is evaluated every five iterations and the simulation ends when the value of v remains constant. In
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
477
Fig. 2. (a) Production rate X for the stationary state (t ¼ 180). (b) Representation of the non-regular grid of points G associated to the solution at this time. If xjk ¼ k=2j 2 G a line of length j is plotted at this point.
that case, we obtain a stationary state in which the flame does not move due to the presence of the convective term ðv=LÞ ox h~. On Fig. 2, we have plotted the computed reaction rate X and the non-regular grid of point G corresponding to the cone-structured space used to compute the solution at time t ¼ 180. The maximum number of degrees of freedom used during the computation is 232 and has to be compared to 211 , the dimension of the regular space (V11 ) we would have used without adaption (the compression factor is then 211 =232 ’ 8:8). This number of degrees of freedom is comparable to the ones obtained by Fr€ ohlich and Schneider [9], 180, using a different wavelet algorithm, and by Pignol [18], 256, using an adaptive multi-grid method. The precision of the algorithm could be evaluated comparing the asymptotical numerical value of v to its theoretical value, v1 ¼ 0:908. From Fig. 3, reporting the time evolution of v, we get the value 0.9083, while the values obtained in Refs. [9,18] were respectively 0.9173 and 0.9087. This results prove that this algorithm is precise and allows us to reduce significantly the dimension of numerical problems using the compression properties of wavelets even when strong gradient regions are moving in the computational domain. 5.2. The Burgers equation In this section, we apply the two adaptive algorithms described in Sections 3 and 4 to the system (1) with GðU Þ ¼ 12ox U 2 , gðxÞ ¼ sinð2pxÞ, v ¼ 10 4 and periodic boundary conditions. It is known (see e.g. Ref. [3] for a complete description of the solution) that the solution exhibits a sharp gradient localized at point x ¼ 0:5 that get its maximal value dmax ¼ jox U ð1=2; tmax Þj 1=ð2vÞ 4 at time tmax ¼ 0:25. For both algorithms, we used the m-spline multiresolution analysis periodized on [0,1] with m ¼ 6 and j0 ¼ 7. In the second algorithm we used a ¼ 0:5 and s ¼ 2.
478
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Fig. 3. Time evolution of the computed flame velocity v.
Fig. 4 shows the numerical results obtained using the different algorithms. Clearly, the nonadapted algorithm is not suitable since Fig. 4(a) exhibits unacceptable oscillations. With at most six levels of adaption (the maximum value of I is 13 from Fig. 4(d)), a very nice approximation of the solution is provided using the algorithm of Section 3 (Fig. 4(b)). It appears from Fig. 4(c) that even with a low fixed resolution (j0 ¼ 7), the algorithm of Section 4 leads to a numerically acceptable solution. Comparing Figs. 4(c) and 5 (the projection of the exact solution in V7 ), reveals
Fig. 4. Solution of the viscous Burgers equation at tmax ¼ 0:25 with v ¼ 10 4 obtained by (a) the non-adapted algorithm in V7 , (b) the first adaptive algorithm of Section 3, (c) the second adaptive algorithm of Section 4 applied in V7 , (d) space scale parameter I (see Definition 2) as a function of time when using the first adaptive algorithm of Section 3.
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
479
Fig. 5. Projection in the space V7 of the exact solution of the Burgers equation at time tmax ¼ 0:25.
that the oscillations of the numerical solution are then comparable to the ones of the orthogonal projection of the exact solution in V7 . That shows that the algorithm of Section 4 damps the Gibbs oscillations connected to the non-linear term without smoothing the solution.
6. Conclusion In this paper, we have presented two numerical schemes based on wavelets and devoted to the adaptive numerical resolution of evolution equations. The key ingredients for the construction of these schemes are the natural properties of wavelets and scaling functions combined, on one hand with a specific structure of the spaces of approximation that induce fast and stable algorithms and, on the other hand, with solution dependent approximations that improve the accuracy. The numerical tests show that the two adaptivity strategy and the related algorithms provide significant improvements. The simulation of the 1D thermo-diffusive flame equation using the first algorithm shows that it is able to track strong gradient regions thanks to simple adaptivity criteria, while retaining a high precision comparable to the one of spectral methods. Even when a large number of refinement levels is required (as in the case of the Burgers equation), the first algorithm provides excellent results. Simulations with even more refinement levels (see Ref. [5]) have shown that in highly adapted configurations (I ¼ 23), the algorithm remains stable, precise and fast. This property has not been proved for other existing adaptive methods. The second algorithm, implemented on a non-adaptive space of approximation, provides, thanks to an efficient estimate of non-linear terms, results close to the projection of the solution on the (low) resolution space. This algorithm opens an interesting direction of research at the intersection between wavelet methods and non-linear finite volume schemes for conservation laws. Even if the two algorithms does not follow the same track (the first provides a high precision estimate of the solution on a highly adapted space while the second provides a stable estimate on a low resolution space), a combination of the two algorithms can be of interest for instance to simulate chemistry equations where the non-linear term is really sharp.
480
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
References [1] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function I: linear reconstruction techniques. SIAM J Numer Anal 1998;36:160–203. [2] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function II: non-linear reconstruction techniques. SIAM J Sci Comput 1999;20:1053–93. [3] Basdevant C, Deville M, Haldenwang P, Lacroix JM, Ouazzani J, Peyret R. Spectral and finite difference solution of the Burgers equation. Comput Fluids 1986;14(1):23–41. [4] Charton P, Perrier V. A pseudo-wavelet scheme for the two-dimensional Navier–Stokes equation. Matem atica Aplicada e Computacional 1996;15. [5] Chiavassa G, Liandrat J. A fully adaptive wavelet algorithm for parabolic partial differential equations. Appl Numer Math 2001;36:333–58. [6] Dahmen W. Wavelet and multiscale methods for operator equations. Acta Numerica 1997;6:55–228. [7] Daubechies I. Ten lectures on wavelets. Philadelphia: SIAM; 1992. [8] Denet B, Haldenwang P. Numerical study of thermal-diffusive instability of premixed flames. Combust Sci Tech 1992;86:199–221. [9] Fr€ ohlich J, Schneider K. An adaptative wavelet Galerkin algorithm for one and two dimensional flame computation. Eur J Mech B––Fluids 1994;13(4):439–71. [10] Fr€ ohlich J, Schneider K. Numerical simulation of decaying turbulence in an adaptive wavelet basis. Appl Comput Harmon A 1996;3:393–7. [11] Gomes S. Convergence estimates for the wavelet-Galerkin method: superconvergence at the node points. Campinas, Brasil: Preprint IMEC-UNICAMP; 1995. quations aux [12] Guichaoua M. Analyses multiresolutions biorthogonales associees a la resolution numerique d’E Derivees Partielles. In: PhD thesis. IRPHE: Universite d’Axi-Marseille II; 1999. [13] Harten A. High resolution schemes for hyperbolic conservation laws. J Comput Phys 1983;49:357–93. [14] Liandrat J, Tchamitchian P. Resolution of the 1d regularized Burgers equation using a spatial wavelet approximation: algorithms and numerical results. ICASE report no. 90-83,1990. [15] Liandrat J, Tchamitchian Ph. On the fast approximation of some non linear operators in adapted wavelet spaces. Adv Comput Math 1998;8:179–92. [16] Mallat S. Multiresolution approximation and wavelets. Trans Am Math Soc 1989;315:69–88. [17] Meyer Y. Ondelettes et operateurs I: ondelettes. Paris: Hermann; 1990. [18] Pignol D. Etude de quelques methodes numeriques de raffinement local. Application a la dynamique des fronts de flamme. In: PhD thesis. Juillet: Universite de Provence; 1995. [19] Ponenti Pj, Liandrat J. Numerical algorithms based on biorthogonal wavelets. ICASE report no. 96-13, 1996.
Two adaptive wavelet algorithms for non-linear parabolic partial differential equations q G. Chiavassa a, M. Guichaoua b, J. Liandrat a,* a b
LATP and ESM2, Technop^ ole de Ch^ ateau-Gombert, 13451 Marseille Cedex 20, France IRPHE and ESM2, Technop^ ole de Ch^ ateau-Gombert, 13451 Marseille Cedex 20, France Received 1 August 2001
Abstract This paper is concerned with the construction of wavelet based adaptive algorithms for the numerical resolution of evolution equations. The adaptivity is applied into two complementary directions. The first direction shares the approaches involved in classical adaptive finite element methods and is related to a solution dependent definition of spaces of approximation. The second direction is related to the approximation of evolution operators that is made solution dependent following the philosophy of essentially nonoscillatory schemes. After the construction of the schemes, numerical tests are provided. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Wavelets; Adaptive methods; Non-linear PDEs
1. Introduction Adaptive methods and related algorithms have been a subject of research for many years in various domains of numerical analysis. In this paper, we are concerned by adaptivity considered in two complementary meanings applied to the numerical resolution of the following equation: 8 ot U þ GðU Þ ¼ mDU ; > > < ðx; tÞ 2 0; 1½ ½0; T ; ð1Þ boundary conditions at x ¼ 0 and x ¼ 1; > > : U ðx; 0Þ ¼ gðxÞ; q
This work has been supported in part by TMR Research Network Contract FMRX-CT98-0184. Corresponding author. E-mail addresses: [email protected] (G. Chiavassa), [email protected] (M. Guichaoua), liandrat@ esm2.imt-mrs.fr (J. Liandrat). *
0045-7930/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 6 1 - 5
468
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
~ U , g is a given function and m is a real constant. where G is a non-linear operator of type U 7! G As usual, t stands for the time variable and x for the spatial variable. The first meaning we give to adaptivity is the one usually given in finite elements methods where the space of approximation into which the solution is computed is solution dependent. The space is therefore said to be adapted. This leads to non-linear methods in contrast with linear methods in which the space of approximation is set in advance. Such non-linear methods allow to improve the approximation efficiency since a finer resolution can be involved only near singularities or strong gradients of the solution. These methods are efficient as soon as fast algorithms are available on the adapted spaces of approximation. The second meaning we give to adaptivity in this paper is the one used in the framework of finite volume methods for conservation law problems. Adaptivity there, leads to solution dependent schemes also called non-linear schemes. Such methods allow the approximation of singular solutions with a significant delay in the occurrence of Gibbs-like oscillations. We present in this paper, two algorithms based on wavelets. On one hand, wavelet spaces offer the possibility to compress very efficiently smooth functions with isolated singularities or strong gradients regions. If some specific property is imposed to the compressed space of approximation (see Definition 2), the fast algorithms related to the pyramidal structure connected to wavelets (see Section 2) remain available. This leads to the first adaptive algorithm presented in Section 3.1. This algorithm belongs to the family of wavelet algorithms developed in Marseille (France) since 1990 [5,14,15]. Full adaptivity and the capability of tracking singularities are emphasized in this paper. This method shares various properties with some other wavelet algorithms of the literature, in particular the ones developed by Fr€ ohlich and Schneider [9,10], by Charton and Perrier [4] or by Gomes [11]. A review on recently developed approaches using wavelets has been performed by Dahmen [6]. On the other hand, recent works by Arandiga et al. [1,2] have introduced the notion of nonlinear multi-resolution and corresponding algorithms. The spirit of these works is a solution dependent definition of details, according to a solution dependent interpolation procedure belonging to the family of essentially non-oscillatory (ENO) schemes [13]. Following the same kind of ideas [12], we present in Section 4 a second adaptive algorithm where the non-linear operators are approximated using a solution dependent approach. This makes the algorithm solution dependent in a new way in comparison to the first algorithm and to other existing adaptive wavelet algorithms. In Section 5 we present numerical tests for the univariate flame (Section 5.1) and Burgers (Section 5.2) equations. 2. General setting In this section, we introduce the general setting in which our algorithms are developed and we recall the main results of wavelet theory necessary for the understanding of the sequel of the paper. 2.1. Wavelet setting We focus on the theory on the line since, by periodization, it can be easily transposed to the periodic framework.
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
469
Following Ref. [7] or Ref. [17] we call a biorthogonal wavelet basis of L2 ðRÞ any set of functions fwj;k ¼ 2j=2 wð2j x kÞ; k 2 Z; j 2 Zg;
ð2Þ
forming a Riesz basis. The dual basis is given by fw~j;k ¼ 2j=2 w~ð2j x kÞ; k 2 Z; j 2 Zg;
ð3Þ
and satisfies hwj;k ; w~j0 k 0 i ¼ dðj;kÞ;ðj0 ;k 0 Þ where h; i stands for the L2 scalar product. When w ¼ w~, the basis is orthonormal. Classically, we note Wj (resp. W~j ) the detail spaces with Wj ¼ spanfwj;k ; k 2 Zg (resp. W~j ¼ spanfw~j;k ; k 2 Zg) and Vj (resp. V~j ) the approximation spaces of the multi-resolution analysis associated to w (resp. w~). We have Vj ¼ spanf/j;k ; k 2 Zg (resp. V~j ¼ spanf/~j;k ; k 2 Zg) where / (resp. /~) is the corresponding scaling function. The following properties should also be translated with : (1) We assume that the spaces Vj offer an order m approximation, i.e. kf PVj f k2 ¼ Oð2 jm Þ;
8 f 2 H m;
ð4Þ
where PVj denotes the orthogonal projection onto Vj and H m stands for the Sobolev space of order m 2 N. The order of approximation (m) is connected to the number of vanishing moments of the dual wavelet. More precisely we have: Z ð5Þ xa w~ðxÞ dx ¼ 0; a ¼ 1; . . . ; m 1: (2) The space L2 ðRÞ can be expanded as an approximation space plus a sum of detail spaces, i.e. L2 ðRÞ ¼ Vj0 þ
þ1 X
ð6Þ
Wj ;
j¼j0
where j0 is any fixed level. Therefore, the family f/j0 ;k ; wj;k ; j P j0 ; k 2 Zg is a Riesz basis of L2 ðRÞ. (3) The embedding relations Vj Vjþ1 and Wj Vjþ1 induce the following scaling and detail relations: X X hk /ð2x kÞ and wðxÞ ¼ gk /ð2x kÞ ð7Þ /ðxÞ ¼ k
k 2
with fgk gk2Z , fhk gk2Z 2 l ðRÞ. The sequences fhk gk2Z , fgk gk2Z are called the filters of the multiresolution. They are the key ingredients of the so-called wavelet decomposition, i.e. the mapping from the set of scaling coefficients fcj0 ;k ¼ hf ; /j0 ;k i; k 2 Zg to the set of wavelet coefficients fdj;k ¼ hf ; wj;k i; j 6 j0 1; k 2 Zg. (4) We also require that the scaling function / gets vanishing moments up to order n/ , more precisely: Z ð8Þ xa /ðxÞ dx ¼ 0; a ¼ 1; . . . ; n/ : (5) In all the paper, we will call dyadic point any real of the form k ¼ k2 j , k 2 Z, j 2 Z and we will identify k to the unique couple ðj; kÞ satisfying k2 j ¼ k with minimal value of j and with the convention 0 ¼ ð0; 0Þ.
470
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Remark 1. We always consider that the number of non-zero coefficients of the filters is finite and equal to lg. This may be true only numerically, i.e. up to a given precision. In the sequel of the paper, lg will be referred as the symbol length.
3. Adaptive space of approximation and corresponding method In this section we present an adaptive resolution method for Eq. (1) using the compression properties of wavelet bases. The resulting algorithm is fully non-linear in the sense that the corresponding spaces of approximation are completely solution dependent and fully adaptive in the sense that its complexity, in term of arithmetic operations and also of memory storage, is proportional to the number of significant wavelet coefficients of the approximation. These properties are not shared by other existing and implemented wavelet algorithms. 3.1. Definition of adapted spaces We begin by introducing the notion of cone-structured spaces of approximation. Definition 2. The space C is said to be cone-structured if 9c; j0 and I 2 N such that C ¼ Vj0 þ
I 1 X
Xj ;
j¼j0
where Xj Wj and satisfies the following property: 8 j; j0 þ 1 6 j 6 I 1; wj;k 2 Xj ) fwj 1;k 0 ; with j2k 0 kj 6 cg Xj 1 : The space C will be a simple cone-structured space if moreover 8 j; j0 þ 1 6 j 6 I 1; 9k1 ; k2 2 N with Xj ¼ fwj;k ; k1 6 k 6 k2 g; where k1 and k2 are functions of j. It will become clear (see Section 3.2) that such specific geometrical conditions imposed to the space of approximation lead to efficient numerical algorithms. From now we describe our method in the case of a simple cone-structured space. The final Remark 5 is devoted to the case of more general cone-structured spaces. We now define the notion of e-adapted space. Definition 3. Given a function U and e > 0, a cone-structured space C is e-adapted to U if: 1. 8 j P I kPWj ðUÞk2 6 ekU k2 , 2. 8 wl;k 2 Vj and wl;k 62 C : jhU ; wl;k ij 6 ekU k2 . Following Ref. [5], we use an algorithm that enforces the space of approximation to remain cone-structured and e-adapted to U assuming some regularity conditions on U. This algorithm
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
471
never uses the approximation of U on the regular space VI but constructs adaptively the conestructured space according to the previous definition. Our goal is now to define an efficient numerical approximation of Eq. (1) on cone-structured spaces. 3.2. Discretization of the equation The discretization sequence we use is a discretization in time followed by a space approximation of each quantities (see Remark 4). 3.2.1. Time discretization Involving the semi-group generated by the Laplacian D, the exact solution of Eq. (1) is Z t mtD emðt sÞD GðU; sÞ ds: U ðx; tÞ ¼ e gðxÞ
ð9Þ
0
Classically, we define an increasing sequence ðtn Þ0 6 n 6 N of ½0; T , with t0 ¼ 0, tN ¼ T and the related time step dtn defined by tnþ1 ¼ tn þ dtn . For notation convenience we will drop the index n for the time step, but, in the numerical experiments, the time step may be n dependent. Using Eq. (9) between tn and tnþ1 leads to: Z dt m dtD U ðx; tnþ1 Þ ¼ e U ðx; tn Þ emðdt sÞD GðU; tn þ sÞ ds: ð10Þ 0
The goal of the time discretization procedure is to give a scheme which allows to compute an approximation of U ðx; tnþ1 Þ, denoted by U nþ1 ðxÞ, knowing an approximation of Uðx; tj Þ, j 6 n. The homogeneous part of Eq. (10) is computed directly without approximation. In Section 3.3 we show how to apply the heat operator to a function expanded on a cone-structured space. For clarity reasons, we simply write SðdtÞ ¼ em dtD . To compute the integral part of Eq. (10), a classical quadrature formula of order g is used and is defined by the weight coefficients fqi ; 0 6 i 6 gg and the quadrature points ftni ; 0 6 i 6 gg. We then obtain, for tni ¼ tn i dt, the semi-discretized form of Eq. (1): U
nþ1
n
ðxÞ ¼ SðdtÞU ðxÞ dt
g X
qi Sðði þ 1Þ dtÞGðU n i ðxÞÞ:
ð11Þ
i¼0
3.2.2. Spatial approximation In this section, the different terms of Eq. (11) are approximated on a finite dimensional space. We introduce the time dependent cone-structured space V n , that is supposed to be e-adapted to U n ðxÞ, and we note Uen the representation of U n ðxÞ on this space. Defining Gnn i as an approximation of G from V n i to V n , the fully discretized form of Eq. (1) writes: Uenþ1 ¼ PV nþ1 SðdtÞUen dt
g X i¼0
n i qi PV nþ1 Sðði þ 1Þ dtÞGnþ1 n i ðUe Þ:
ð12Þ
472
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Remark 4. Thanks to the specific discretization of the equation, the time evolution and projection operators are decoupled. This allows to adapt the space of approximation at each time step after applying the evolution operators. This property is important with regards to the spatial adaptivity of our algorithm and is used in practical computations. 3.3. Application of the linear operator From relation (12) and definition of the cone-structured space V n , we have to compute Sði dtÞUen ¼
j0 1 2X
cj0 ;k Sði dtÞ/j0 ;k þ
I 1 X X
dj;k Sði dtÞwj;k
ð13Þ
j¼j0 wjk 2Xj
k¼0
with cj0 ;k ¼ hUen ; /j0 ;k i and dj;k ¼ hUen ; wj;k i. Denoting sij0 ;k ¼ Sði dtÞ/j0 ;k and hij;k ¼ Sði dtÞwj;k , we know from Ref. [17] that these functions (called vaguelettes) and their Fourier transforms are uniformly localized with regards to j. This implies (for instance for hij;k ) that, on one hand, 8 e; 9ðr; r0 Þ 2 N2
such that sij0 ;k 2e Vj0 þr and hij;k 2e Vjþr0 ;
ð14Þ
where f 2e F means that 9 g 2 F such that kf gkL2 6 e and that, on the other hand, 9d 2 N such that jhhij;0 ; wj0 ;k ij 6 e if jkj > d:
ð15Þ
The link between r, r0 and e can be evaluated numerically. For instance for spline wavelets of order m ¼ 6 and e ¼ 10 6 , we use in our applications r ¼ 2, r0 ¼ 3 and d 6 50. Similar relations are valid for the other scalar products. It follows that the estimate of PV nþ1 Sði dtÞ can be performed efficiently using fast numerical convolutions. Using the properties of wavelets and of the heat operator, an efficient algorithm for the computation of the scalar products satisfying relation (15) has been derived in Ref. [5]. This leads to the computation of the homogeneous part of Eq. (12) with a method of complexity proportional to #V nþ1 , the dimension of V nþ1 . 3.4. Evaluation of the non-linear term n i We recall in this section the method we use to evaluate the quantities Gnþ1 n i ðUe Þ involved in relation (12). It has been derived in Ref. [5] and can be considered as a variation around the algorithm constructed in Ref. [15]. To simplify the notations, we only detail the case i ¼ 0. We first associate to V n (supposed to be a simple cone-structured space), a non-regular grid of points G, defined by (see Section 2 for the notations):
G ¼ fk; such that j ¼ j0 or j0 < j < I 1 and 2k1 ðj 1Þ lg =2 6 k 6 2k2 ðj 1Þ þ lg =2g; where lg is the symbol length (see Remark 1). It must be emphazised that there exists an algorithm for the computation of the points values of Uen on G from its representation on V n (and reversely) which complexity is proportional to #V n . This algorithm, developed and detailed in Ref. [5], is
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
473
similar to the classical tree-algorithm of Mallat [16], but acts only on the space V n and never requires the regular space VI . The evaluation of the non-linear terms can be described as follows, assuming that GðUÞ ¼ ~ U: G fhU ; wj;k ig 7! fhU ; /k ik2G g; fhU ; /k ig 7! fU ðkÞg; ~ðU ðkÞÞg; fU ðkÞg 7! fG ~ðU ðkÞÞg 7! fhGðU Þ; /k ig; fG fhGðUÞ; /k ig 7! fhGðU Þ; wj;k ig; where 7! stands for a mapping or its numerical approximation. Remark 5. We have presented all the steps of our adaptive algorithm for a simple cone-structured space. For practical computations a global space of approximation is generally built as the union of different simple cone-structured spaces (each of them being, for instance, associated to a specific strong gradient region of the solution). The application of the linear operators (heat operator and adaptive tree algorithm) is then performed independently on each simple space. The evaluation of non-linear term is slightly modified computing the values of Gnþ1 on the union of all the grid points G. Nevertheless, the complexity of the global algorithm remains proportional to #V n .
4. Adaptive approximation of evolution operator and corresponding method This section is devoted to the description of the second algorithm where adaptivity is related to a solution dependent approximation of the operators. The approach developed in this section is then complementary to the one developed in Section 3 where adaption deals with the spaces of approximation. Since the motivation is essentially given by the approximation of the non-linear part of Eq. (1) we focus on an equation of type: 8 ot U þ GðU Þ ¼ 0; > > < ðx; tÞ 2 0; 1½ ½0; T ; ð16Þ boundary conditions at x ¼ 0 and x ¼ 1; > > : U ðx; 0Þ ¼ gðxÞ: 4.1. General discretization of the equation As in Section 3.2, and for the same reasons (see Remark 4), the discretization of the Eq. (16) is first performed on the time variable and secondly on the spatial variable. Since the goal, in this section is to derive a solution dependent approximation of the non-linear term, the time discretization we use is explicit. We first describe the general discretization and describe in Section 4.2 the solution dependent discretization.
474
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
4.1.1. Time discretization The discretization in time is performed using a finite difference approximation of order r of Taylor type that writes, once a segmentation (see Section 3.2) has been performed: dtn ½1 n dtns 1 ½s 1 n nþ1 n n ¼ U dtn GðU Þ þ ð17Þ G ðU Þ þ þ G ðU Þ ; U 2 r! where G½0 ¼ G and G½lþ1 ¼ Dx G½l G. Remark 6. Such Taylor schemes are not very popular since they involve the computation of extra terms (fG½l g0
2 r!
ð18Þ
The remaining difficulty is to compute the scaling coefficients of G½l ðUn Þ, 0 6 l 6 s 1 from those of U n . The general framework described in Section 3.4 is suitable but the occurrence of derivatives of high order of non-linear terms (see the right-hand side term of Eq. (18)) imposes to pay special attention to their approximation when the solution exhibits locally strong gradients. That explains why a solution dependent approximation has been developed. 4.2. Solution dependent approximation We come back to the evaluation of hG½l ðU n Þ; /j0 ;k i. If we assume that GðU Þ is a combination of ~ U and of derivatives, we have to evaluate the scaling coefficients of derivatives of terms of type G non-linear forms of U. This is known to be a critical issue for which various solutions have been proposed in particular in the framework of finite volume methods for conservation laws [13]. We ~ U by introducing data propose to mimic the so-called ENO scheme for the evaluation of G dependent scaling functions generating Vj0 . Depending on U, we define specific functions /\j0 ;k , introduced in the algorithm (18) as follows: fhU ; /j0 ;k ig 7! fhU; /\j0 ;k ig; fhU ; /\j0 ;k ig 7! fU ðk2 j0 Þg; ~ðUðk2 j0 ÞÞg; fU ðk2 j0 Þg 7! fG ~ðUðk2 j0 ÞÞ; /\ ig; ~ðU ðk2 j0 ÞÞg 7! fhG fG j0 ;k ~ðU ðk2 j0 ÞÞ; /j ;k ig: ~ðU ðk2 j0 ÞÞ; /\ ig 7! fhG fhG j0 ;k 0
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
475
The construction of /\j0 ;k is made as follows: Assuming that U is smooth but in the vicinity of a point a ¼ ka 2 j0 , we write f/\j0 ;k g ¼ Ma f/j0 ;k g;
ð19Þ
where, each line of the infinite and invertible matrix Ma is an interpolation stencil of finite length. The matrix Ma (see Eq. (20) for an example) is constructed such that a minimal number of stencil intersect the point a in order to have /\j0 ;k ðaÞ ¼ 0 for a maximum of values of k. The construction of Ma is not unique and the analysis of an optimal construction is under study. Fig. 1(c) and (d) exhibits two functions /\0;k , corresponding to the matrix Ma presented in Eq. (20), for different values of k. These functions, in comparison to / are less oscillating and not symmetric. ka 0
1 0 B0 1 B B0 0 B B0 : B B0 1 10 B B0 0 B B0 : B B 0 : Ma ¼ B B0 : B B0 : B B0 : B B0 : B B0 : B B0 : B B0 : B @0 :
#
: 0 1 0 1 2
: 0 : : : : : : : : :
: : 0 1
1 :
101 0 : : : : : : : :
: : : 0 2 : 1 2
: : : :
12 :
1
101
1 2
: 0 : : : : : :
: : : : : : : :
: : : :
: : : : 1 0 10 : : 2 12
1 2 1
12 10 1 0 10 : 0 : : : : : : : : : :
: : : : : 0 1 10
12
2
12 : 0 : : : :
: : : : : : 0 1 10
1 2 : 1 10
: : : :
: : : : : : : 0 1 2
1 :
12 : : : :
: : : : : : : :
101 1 2
: 2 0 : : :
: : : : : : : : :
101 :
1 1 0 : :
: : : : : : : : : 0 :
: : : : : : : : : : 0 1 1
2 10 0 : 1 0 0 1 : 0
1 0 0C C 0C C 0C C 0C C 0C C 0C C 0C C 0C C: 0C C 0C C 0C C 0C C 0C C 0C C 1A
ð20Þ To end the estimate of hG ðU Þ; /j0 ;k i, the computation of the derivatives is performed using a method developed in Ref. [19] based on biorthogonal multi-resolution analyses associated to ow=ox. One of these scaling functions is plotted on Fig. 1(b). ½l
n
5. Numerical tests 5.1. The 1D thermo-diffusive flame equation In this section we apply the method described in Section 3 to the following system of equations: 8 < ot h þ v ox h X ¼ Dh; o C þ v ox C þ X ¼ Le1 DC; ð21Þ : t ðx; tÞ 2 0; L½ 0; T ½
476
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Fig. 1. (a) Initial spline scaling function / of order m ¼ 6, (b) scaling function associated to the operator o=ox, (c, d) scaling functions modified by the matrix Ma : /\0;k and /\0;k 0 for k 6¼ k 0 .
with X ¼ ðb2 =2LeÞC expðbðh 1Þ=ð1 þ cðh 1ÞÞÞ. This system stands for the propagation of a planar flame in the case of two reactive chemical components (hydrogen–air). The temperature hðx; tÞ, the hydrogen concentration Cðx; tÞ and the flame velocity v are the unknowns. The parameters Le, b, c are respectively the Lewis and Zeldovich numbers and the temperature ratio. They are set to the following classical values [8]: Le ¼ 1, b ¼ 10, c ¼ 0:8. When Le ¼ 1, the relation between h and C is linear and the system (21) reduces only to the first equation. The boundary conditions are hðx ¼ 0; tÞ ¼ 0 and hðx ¼ L; tÞ ¼ 1 and L ¼ 60 in our simulations. Following Ref. [8], it is possible to recover periodic boundary conditions, introducing the function sðxÞ ¼ 12ð1 þ tanhð10 tanðpðx 12ÞÞÞÞ satisfying sð0Þ ¼ 0 and sðþ1Þ ¼ 1 and solving the equation for the unknown h~ ¼ h s (we suppose that L is large enough so that sðLÞ ’ 1 at the precision of the computer). Writing this equation for x 2 0; 1½ gives: v 1 ot h~ þ ox h~ X ¼ 2 Dh~ þ F ðs; v; LÞ; L L
ð22Þ
where b2 ð1 h~ sÞ exp X¼ 2Le
! bðh~ þ s 1Þ ; 1 þ cðh~ þ s 1Þ
F ðs; v; LÞ ¼
1 v Ds ox s: 2 L L
The algorithm of Section 3 is applied using an orthonormal basis of spline wavelets of order 6, the thresholding parameter e set to 5 10 7 and the lower level is j0 ¼ 7. The velocity v is evaluated every five iterations and the simulation ends when the value of v remains constant. In
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
477
Fig. 2. (a) Production rate X for the stationary state (t ¼ 180). (b) Representation of the non-regular grid of points G associated to the solution at this time. If xjk ¼ k=2j 2 G a line of length j is plotted at this point.
that case, we obtain a stationary state in which the flame does not move due to the presence of the convective term ðv=LÞ ox h~. On Fig. 2, we have plotted the computed reaction rate X and the non-regular grid of point G corresponding to the cone-structured space used to compute the solution at time t ¼ 180. The maximum number of degrees of freedom used during the computation is 232 and has to be compared to 211 , the dimension of the regular space (V11 ) we would have used without adaption (the compression factor is then 211 =232 ’ 8:8). This number of degrees of freedom is comparable to the ones obtained by Fr€ ohlich and Schneider [9], 180, using a different wavelet algorithm, and by Pignol [18], 256, using an adaptive multi-grid method. The precision of the algorithm could be evaluated comparing the asymptotical numerical value of v to its theoretical value, v1 ¼ 0:908. From Fig. 3, reporting the time evolution of v, we get the value 0.9083, while the values obtained in Refs. [9,18] were respectively 0.9173 and 0.9087. This results prove that this algorithm is precise and allows us to reduce significantly the dimension of numerical problems using the compression properties of wavelets even when strong gradient regions are moving in the computational domain. 5.2. The Burgers equation In this section, we apply the two adaptive algorithms described in Sections 3 and 4 to the system (1) with GðU Þ ¼ 12ox U 2 , gðxÞ ¼ sinð2pxÞ, v ¼ 10 4 and periodic boundary conditions. It is known (see e.g. Ref. [3] for a complete description of the solution) that the solution exhibits a sharp gradient localized at point x ¼ 0:5 that get its maximal value dmax ¼ jox U ð1=2; tmax Þj 1=ð2vÞ 4 at time tmax ¼ 0:25. For both algorithms, we used the m-spline multiresolution analysis periodized on [0,1] with m ¼ 6 and j0 ¼ 7. In the second algorithm we used a ¼ 0:5 and s ¼ 2.
478
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
Fig. 3. Time evolution of the computed flame velocity v.
Fig. 4 shows the numerical results obtained using the different algorithms. Clearly, the nonadapted algorithm is not suitable since Fig. 4(a) exhibits unacceptable oscillations. With at most six levels of adaption (the maximum value of I is 13 from Fig. 4(d)), a very nice approximation of the solution is provided using the algorithm of Section 3 (Fig. 4(b)). It appears from Fig. 4(c) that even with a low fixed resolution (j0 ¼ 7), the algorithm of Section 4 leads to a numerically acceptable solution. Comparing Figs. 4(c) and 5 (the projection of the exact solution in V7 ), reveals
Fig. 4. Solution of the viscous Burgers equation at tmax ¼ 0:25 with v ¼ 10 4 obtained by (a) the non-adapted algorithm in V7 , (b) the first adaptive algorithm of Section 3, (c) the second adaptive algorithm of Section 4 applied in V7 , (d) space scale parameter I (see Definition 2) as a function of time when using the first adaptive algorithm of Section 3.
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
479
Fig. 5. Projection in the space V7 of the exact solution of the Burgers equation at time tmax ¼ 0:25.
that the oscillations of the numerical solution are then comparable to the ones of the orthogonal projection of the exact solution in V7 . That shows that the algorithm of Section 4 damps the Gibbs oscillations connected to the non-linear term without smoothing the solution.
6. Conclusion In this paper, we have presented two numerical schemes based on wavelets and devoted to the adaptive numerical resolution of evolution equations. The key ingredients for the construction of these schemes are the natural properties of wavelets and scaling functions combined, on one hand with a specific structure of the spaces of approximation that induce fast and stable algorithms and, on the other hand, with solution dependent approximations that improve the accuracy. The numerical tests show that the two adaptivity strategy and the related algorithms provide significant improvements. The simulation of the 1D thermo-diffusive flame equation using the first algorithm shows that it is able to track strong gradient regions thanks to simple adaptivity criteria, while retaining a high precision comparable to the one of spectral methods. Even when a large number of refinement levels is required (as in the case of the Burgers equation), the first algorithm provides excellent results. Simulations with even more refinement levels (see Ref. [5]) have shown that in highly adapted configurations (I ¼ 23), the algorithm remains stable, precise and fast. This property has not been proved for other existing adaptive methods. The second algorithm, implemented on a non-adaptive space of approximation, provides, thanks to an efficient estimate of non-linear terms, results close to the projection of the solution on the (low) resolution space. This algorithm opens an interesting direction of research at the intersection between wavelet methods and non-linear finite volume schemes for conservation laws. Even if the two algorithms does not follow the same track (the first provides a high precision estimate of the solution on a highly adapted space while the second provides a stable estimate on a low resolution space), a combination of the two algorithms can be of interest for instance to simulate chemistry equations where the non-linear term is really sharp.
480
G. Chiavassa et al. / Computers & Fluids 31 (2002) 467–480
References [1] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function I: linear reconstruction techniques. SIAM J Numer Anal 1998;36:160–203. [2] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function II: non-linear reconstruction techniques. SIAM J Sci Comput 1999;20:1053–93. [3] Basdevant C, Deville M, Haldenwang P, Lacroix JM, Ouazzani J, Peyret R. Spectral and finite difference solution of the Burgers equation. Comput Fluids 1986;14(1):23–41. [4] Charton P, Perrier V. A pseudo-wavelet scheme for the two-dimensional Navier–Stokes equation. Matem atica Aplicada e Computacional 1996;15. [5] Chiavassa G, Liandrat J. A fully adaptive wavelet algorithm for parabolic partial differential equations. Appl Numer Math 2001;36:333–58. [6] Dahmen W. Wavelet and multiscale methods for operator equations. Acta Numerica 1997;6:55–228. [7] Daubechies I. Ten lectures on wavelets. Philadelphia: SIAM; 1992. [8] Denet B, Haldenwang P. Numerical study of thermal-diffusive instability of premixed flames. Combust Sci Tech 1992;86:199–221. [9] Fr€ ohlich J, Schneider K. An adaptative wavelet Galerkin algorithm for one and two dimensional flame computation. Eur J Mech B––Fluids 1994;13(4):439–71. [10] Fr€ ohlich J, Schneider K. Numerical simulation of decaying turbulence in an adaptive wavelet basis. Appl Comput Harmon A 1996;3:393–7. [11] Gomes S. Convergence estimates for the wavelet-Galerkin method: superconvergence at the node points. Campinas, Brasil: Preprint IMEC-UNICAMP; 1995. quations aux [12] Guichaoua M. Analyses multiresolutions biorthogonales associees a la resolution numerique d’E Derivees Partielles. In: PhD thesis. IRPHE: Universite d’Axi-Marseille II; 1999. [13] Harten A. High resolution schemes for hyperbolic conservation laws. J Comput Phys 1983;49:357–93. [14] Liandrat J, Tchamitchian P. Resolution of the 1d regularized Burgers equation using a spatial wavelet approximation: algorithms and numerical results. ICASE report no. 90-83,1990. [15] Liandrat J, Tchamitchian Ph. On the fast approximation of some non linear operators in adapted wavelet spaces. Adv Comput Math 1998;8:179–92. [16] Mallat S. Multiresolution approximation and wavelets. Trans Am Math Soc 1989;315:69–88. [17] Meyer Y. Ondelettes et operateurs I: ondelettes. Paris: Hermann; 1990. [18] Pignol D. Etude de quelques methodes numeriques de raffinement local. Application a la dynamique des fronts de flamme. In: PhD thesis. Juillet: Universite de Provence; 1995. [19] Ponenti Pj, Liandrat J. Numerical algorithms based on biorthogonal wavelets. ICASE report no. 96-13, 1996.