Dynamical multilevel schemes for the solution of evolution equations by hierarchical finite element discretization

APPLIED ELSEVIER

NUMERICAL MATHEMATICS Applied Numerical Mathematics 23 (1997) 403--442

Dynamical multilevel schemes for the solution of evolution equations by hierarchical finite element discretization C. C a l g a r o ,

J. L a m i n i e ,

R. T e m a m * , l

Lab©rat©ire d'Analyse Numdrique et Equations aux Ddriv~es Partielles, Universit# Paris-Sud et CNRS URA 760, Bdtiment 425, 91405 Orsay Cedex, France

Abstract The full numerical simulation of turbulent flows in the context of industrial applications remains a very challenging problem. One of the difficulties is the presence of a large number of interacting scales of different orders of magnitude ranging from the macroscopic scale to the Kolmogorov dissipation scale. In order to better understand and simulate these interactions, new algorithms of incremental type have been recently introduced, stemming from the theory of infinite dimensional dynamical systems, see, e.g., the algorithms of Foias, Jolly et al. (1988), Foias, Manley and Temam (1988), Laminie et al. (1993, 1994), Marion and Temam (1989, 1990), Temam (1990). These algorithms are based on decompositions of the unknown functions into a large and a small scale component, one of the underlying ideas being to approximate the corresponding attractor. In the context of finite elements methods, the decomposition of solution into small and large scale components appears when we consider hierarchical bases (Yserentant, 1986; Zienkiewicz et al., 1982). In the present article we derive several estimates of the size of the structures for the linear and nonlinear terms which correspond to interactions of different hierarchical components of the velocity field, and also their time variation. The one-step time variation of these terms can be smaller than the expected accuracy of the computation. Using this remark, we implement an adaptive spatial and temporal multilevel algorithm which treats differently the small and large scale components of the flow. We derive several a priori estimates in order to study the perturbation introduced into the approximated equations. All the interaction terms between small and large structures are frozen during several time steps. Finally we implement the multilevel method in order to simulate a bidimensional flow described by the Burgers' equations. We perform a parametric study of the procedure and its efficiency. The gain on CPU time is significant and the trajectories computed by our multi-scale method remain close to the trajectories obtained with the classical Galerkin method. @ 1997 Elsevier Science B.V. Keywords: Hierarchical finite elements; 2D Burgers' equations; Auto-adaptive and multi-scale solvers; Long

time integration

* Corresponding author. i The numerical simulations of the work was supported by the Institut du Drveloppement et des Ressources en Informatique Scientifique (IDRIS), Bat. 506, B.P, 167, 91403 Orsay Cedex, France. 0168-9274/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII SO 168-9274(96)00074-8

404

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403~142

1. Introduction

The integration of evolution equations on large intervals of time remains a difficult problem, despite the significant increase in the computing power. It is well known that the large time behaviour of solutions of dissipative evolution equations depends on certain non-dimensional parameters such as the Reynolds or Grashof numbers or other similar numbers. When the values of those physical parameters are small the system evolves, as t --+ c~, towards a stationary solution. On the contrary, for large values of such parameters the system exhibits nontrivial dynamics related to turbulence. For dissipative evolution equations, such as the Navier-Stokes equations, it was shown that the longtime behaviour is encompassed by a compact attractor and that such an attractor has a finite fractal dimension (see [8]). Hence, although the attractor can be complex, or even fractal, the long-time dynamics of the system are finite dimensional. In practice, Constantin et al. in [2] estimate the number of degrees of freedom, needed to insure that a numerical solution of the Navier-Stokes equations is at least qualitatively correct, to be of order of the fractal dimension of the global attractor. To solve dissipative evolution equations in dynamically nontrivial situations, Foias, Sell and Temam have introduced in [9] the concept of inertial manifolds, which are finite dimensional manifolds containing the global attractor. The approximate inertial manifolds, introduced by Foias, Manley and Temam in [11], are more flexible and a less restrictive idea. They are explicit smooth manifolds approximating the attractor up to a certain level of accuracy and they allow to build numerical schemes providing orbits lying on the approximate set. From this Marion and Temam have subsequently developed new algorithms of incremental type, called nonlinear Galerkin methods (NLG methods), in the context of spectral and finite element discretizations (cf. [ 16,17]) and incremental unknowns methods in the context of finite differences [22]. These innovative methods are based on a decomposition of the unknowns, such as the velocity field, into a small scale component and a large scale one. Essential in the NLG methods is a differentiated treatment of small and large scales and in that respect the nonlinear Galerkin methods and their variations appear as dynamical multi-resolutions methods. A multilevel scheme, which develops a self-adaptive procedure, is implemented for spectral discretizations by Dubois, Jauberteau and Temam (cf. [5] and the reference therein; see also [6]); moreover, Debussche, Dubois and Temam develop further the study of the numerical implementation of the NLG method for turbulent flows (cf. [3]). In the context of the finite elements discretizations, implementations of the NLG method have been done by Laminie, Pascal and Temam [14,15]. The classical Galerkin approximation of Burgers' and Navier-Stokes problems have been studied using the hierarchical basis decomposition (with a coarse initial triangulation and a finer one obtained by subdividing any triangle into four similar subtriangles). The approximations made for the NLG method are justified when the mesh size h converges to zero. Although the quantities which represent the small scale components of the flow are small, numerical experiments done by Pascal [18] show that they cannot be neglected, especially in the case of higher Reynolds numbers. So we introduce a new idea which is to consider the time variations of terms depending upon the small scale components of the flow, with the intention of freezing them when their time variation is small. In order to split the solution into its small and large scale components, we use the hierarchical finite element bases, defined by considering several levels of nested triangulations. The study is performed for the two-dimensional Burgers' problem, considering the evolution of all linear and

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

405

nonlinear terms. The first results show that the time variation, over one time step, of the coupled terms can be small. Therefore these terms can be frozen during a short period of time representing a certain number of time steps. Based on this results, we propose a spatial and temporal multilevel self-adaptive method. The level of refinement, which defines the separation of the solution into its small and large scales, evolves in time as a multilevel V-cycle. Moreover the small components of the solution are frozen and, besides, we freeze the suitable terms in the equation during a suitable short time period. The structure of a series of V-cycles is dynamically determined using criteria which guarantee that all frozen quantities are smaller than a given parameter, representing the expected accuracy of the computation. Significant savings of computing time are obtained as compared to the classical finite element Galerkin implementation. The promises of the multi-scale algorithm are important: most often only one coarse grid nonlinear equation is solved instead of two nonlinear equations as in the previous implementations [14,15]. The article is organized as follows. In Section 2 we introduce the nonlinear Burgers' problem and its weak formulation. In Section 3 we recall briefly the hierarchical finite element basis and we split the solution into its components of different scales. The coarsest mesh gives the large structures of the flow, whereas on the finer grids one derives the corrections of the solution, which are smaller and smaller according to the level of refinement. In Section 4 we recall the classical Galerkin (FEM) implementation: the approximate solution computed by the standard Galerkin method will be compared with that provided by our multilevel algorithm. Moreover this solution will be the object, in Section 5, of an a posteriori estimate performed on the different terms of the equations. These quantities and the time variation of the coupled terms are analyzed in order to emphasize the numerical justifications which are applied in the development of the algorithm proposed in Section 6. The level of refinement evolves in time undergoing successive multilevel V-cycles. Moreover, the corrections to the approximate solution are frozen (indeed only the nonlinear equation corresponding to the large structures is solved) and also we freeze the coupled interaction terms during a suitable short period of time. After a complete description of this multilevel procedure, we derive important estimates of the characteristic numbers which determine the structure of the V-cycle series. The last section is devoted to the description and the analysis of numerical results performed with our multilevel method. This simulation is compared with the approximate solution determined by the classical Galerkin method (FEM).

2. Test problem and weak formulation Let f2 be a bounded domain of ~2 with boundary 0f2. We consider, as a model of nonlinear evolution equations, the two-dimensional Burgers' problem. The Burgers' equations are a simplified formulation of the incompressible Navier-Stokes equations. They constitute a model of viscous flows which incorporates the diffusive viscous terms and the nonlinear convection processes but the pressure gradient terms and the incompressibility condition are removed. These equations encompass some qualitative aspects of the Navier-Stokes equations and they produce a suitable model for the development of algorithms for Navier-Stokes equations themselves. In the two-dimensional Burgers' problem the unknown function u maps $2 x [0, T] into R 2 and it satisfies the following equations:

C. Calgaro et al. /Applied Numerical Mathematics 23 (1997) 403-442

406

0u Ot

A u + ( u . V ) u = y in ~ x [0, T],

Rl-e

u=g

(1)

on 0~2, V x • ~.

~(x, o) = uo(x)

We denote by (.,.) the scalar product in L2(I2) and by ((.,.)) that of H01(~) defined by ((u,v)) = f o V u . Vv. We use the following function spaces: 7-/= L2(f2) × L2(~2),

V = Hi(/2) x H I ( o ) ,

Vo = Hd(Y2) x H I ( O ) ,

Vg = {u • V; u = g on BY2}.

Let B(., .) be a bilinear operator from V x V into V ~. We denote by b the trilinear continuous form on V given by b(~,v,w) = (B(~,~),w>v,v '

V u , ~ , w • V.

For the Burgers' equations 2

f

b(u, v, w) = ~ ] i,j=l j

0vj

(2)

ui-~---wj dx. oxi

For the finite element treatment, we write (1) in weak form. Let f • 7-/and u0 • Vg; for all t > 0, the unknown function u = u(t) • L2(O, T; Vg) N C(0, T; 7-/) is solution of

-~,u

+

1

((u, ~)) + b(u, u, "~) = (f, "5), V~ • Vo,

(~(o), ~) = (uo, ~),

(3)

v~ • Vo.

3. The hierarchical finite element basis In order to split the solution into its small and large scale components, we use the hierarchical finite element basis. This concept has been introduced by Zienkiewicz et at. in [24] and has been developed by Yserentant [23] in the context of elliptic linear problems. From an initial finite element mesh ~ of the computational domain Y2, we construct a family of nested meshes

%cTc.--c% by subdividing each triangle into four similar subtriangles. Here, 7~ represents the finest triangulation, called fine mesh or nodal mesh. In Fig. 1 we represent three nested meshes (d = 2); for each level of triangulation, we have considered a uniform grid. Let r k be the set of vertices of triangulation 7~ and let Vk be the finite element space which approximates H 1(12):

vk = {vk • c ° ( o ) I w • 7~, vkI~ • Pl(V)},

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

407

level 0 __

level 1

.......... level 2

Fig. 1. Hierarchical triangulation of domain/2:3"0 C 7~ C "-/~.

where P1 (T) is the set of polynomials of degree less than or equal to one on each triangle of 7~. Obviously we have ~ k C ~kq_l

and

Vk C Vk+l

(0<~k~
The finite element space Vd associated with the fine grid 7~ is split into

Vd= Vd- + Wd, and recursively Vd ~-- Vo-I- W1 q-...-~- W d.

The usual nodal basis of Vk is denoted by ¢k,i. It is the canonical basis consisting of functions which take the value 1 at the node i of Tk and 0 at all other nodal points. The functions in the hierarchical basis of Wk (k = 1,... ,d) are denoted by ~bk,j, where j is a nodal point of "Fk \ Tk-1; Ck,j is a function of the nodal basis of Vk which vanishes on the nodes of triangulation 7k-1 (i.e., ~bk,j = ¢k,i at the node j of 7~ \ "Fk-X). Then we denote Vk -----Span{¢k,1,..., Ck,nk},

nk = card(Zk) -----dim(Vk),

Wk = Span{¢k,nk_l+l,... ,¢k,nk},

nk -- nk-1 = card(Z?k \ Sk-1) = dim(Wk).

If Ud is an element of Vd, then the corresponding decomposition reads Ud = Yd-1 + Zd,

where Yd-1 E Vd-1 and Zd E Wd.

Recursively we find ud=YO+ Zl+'"+

Zd,

wherey0cV0andzkEWk

For each level k, 1 ~< k ~< d, we write nk--1

nk

Yk ~-Yk-1 q-Zk : E yk-l'iCk-l'iq'Z Zk'i~Pk'i" i=l i=nk_lq-1

(k= l,...,d).

408

C. Calgaro et al. / A p p l i e d Numerical Mathematics 23 (1997) 403-442

The values of Yk-1,i and Zk, i are uniquely determined by Vi E [1,nk-1],

Yk-l,i = yk(Ai),

Vi E [nk-i -% 1, nk],

Zk,i = y k ( B i ) -- l ( y k - i ( A i l )

(4) q- Y k - l ( A i 2 ) ) ,

where Ai, A~l, Ai2 E "l-k-1 and Bi E 7~ \ 7k-1, Bi being the midpoint of the edge [Ail, Ai2] (with generalized notation, we set Yd = Ud). The large structures defined by Yk-1 are of the same order as Ud. From Taylor's formula, we easily infer the order of the corrections zk over the grid 7k, for each k ~< d. Let Ail,Ai2 E Tk-1 and Bi E Tk \ 7k-l, Bi the midpoint of the edge [Ail, Ai2]. We obtain It

, ___, 2 ud(Ai2) -~ ud(Bi)q- BiAi2 Ud(Bi)+ BiAi2

2

---+

3

--"+

3

+ O(IIBiA,211 ),

tl

ud(Ail) = Ud(Bi)-- A i--' l B i Ud(Bi)+ A i--' lBi

,.,d(B,) + O(IIA,,B,I I ). 2

From (4) we find It

zk(Bi)=ud(Bi)--

---+ 2 Ud(Bi) ~ 3 Ud(Ail) + ua(Ai2) =-ailBi 2 +O(llZilBill )" 2

Hence the incremental components zk can be of order h 2, where hk is the mesh width of the triangulation 7~ (cf. [22]); but in fact, if the solution has strong gradients, the size of the corrections zk can grow considerably. We will keeps this difficulty in mind when we build our algorithm. Our approach is based on differentiated treatments of the small scale and the large scale components. This approach utilizes the hierarchization of finite element spaces introduced above. We will consider the size of different structures and especially the time variation over one time step of the small scale components as well as the time variation of the linear and nonlinear terms which depend on the small scales. When these variations are small, we will be able to freeze the corresponding terms during a short period of time representing a few time steps. Then we will freeze the small corrections Zk computed over one to several finer grids, and we will approximate the solution only over the coarser grids. The solution computed by our procedure will be compared with that given by the classical Galerkin method applied to the nodal basis which we now present in the following section.

4. Discretization on the nodal basis

Let Vd be a finite element space which approximates H 1 (/2). We set Vd = Vd x V~,

Vo, = (v,, n ];g,d = (Vd E )24, Vd = gd on i)/2},

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403~142

409

with gd an appropriated approximation of g, and we equip "1,)9,d with the usual nodal basis. The discrete form of (3) consists in finding ud(t) C Vg,d such that

laUd _ "~ 1 ~-~-,UdJ + ~ee((Ud,~d)) + b(Ud, Ud,'~d) = (f, Ud), V~d e VO,d,

(ud(0), ~) = (u0,~,~d),

(5)

v~d e v0,~.

The time derivative is approximated by a finite difference scheme. We use a second-order implicit Crank-Nicholson scheme for the linear terms and a first-order explicit Euler scheme for the nonlinear terms. Once we known the solution u~ at time nAt, then the unknown u~ +l E Vg,d at time (n + 1)At is solution of the following linear system: \(

)

))

u~+l-u'~ ~i7 ' Ud + Ree

,Ud

= ( f , ud) -

b(u~,u"d, Ud),

V~d E ];0,d.

(6)

The scheme given by (6) is found to be conditionally stable: it requires the ratio lUdlooAt/hd not to exceed a critical value (CFL condition). Moreover the linear system for u~ +l is well conditioned.

5. Analysis of the decomposition based on the hierarchical basis Using the hierarchical basis decomposition, we split the solution into its small scale and its large scale components. In this section we give some numerical (empirical) estimates of the terms which appear in the equations projected on the different meshes. More precisely we will analyze the norm L 2 of all linear and nonlinear terms and also the time variation of the terms which depend on the small scale z. We aim at estimating the size of these terms and especially at checking if their variation over one time step can be small. The hierarchical basis for the finite element spaces considered, gives us the decomposition

Vg,d = Vg,~_~ + Wg,~, and recursively

Vg,~_l = Vg,d-2 + Wg,d-1,

...,

Vg,l = Vg,o + Wg,1.

From the variational formulation (5), we infer a spatial discretization corresponding to the decomposition of Vg,d into ];a,d-1 + Wg,d. For all t > 0, we look for Yd-l(t) E ];a,a-I and Zd(t) E kYa,d solutions of the following system:

0

1 + '-~-,Yd-1

+

((Yd-1

+Zd, Yd-1))

+b(yd-l,Yd-l,'Yd-l) + bint (Yd-I,Zd,'Yd-1) = (f,Y'd-1), VYd-1 E "l)O,d-1, OZd l((y~_l + zd, ~d))

(7)

+ b(yd-1, Yd-1, "Zd) ~- bint (Yd-1, Zd, "Zd) : (f, Z'd),

(8)

(~t-1 + -~-,~) +

VZ"d E "l/~O,d .

C. Calgaroet al. / Applied NumericalMathematics23 (1997) 403-442

410

We recall that b stands for the trilinear form defined in (2). Here, we have split the nonlinear term b (u, u, w) into b(y, y, w) + bint(Y, z, w), where

hint(Y, z, w) : b(y, z, w) -~-b(z, y, w) -]- b(z, z, w) = b(u,

w ) - b(y, y,

This term represents the interaction between the small and large structures and the interaction of the

small structures among themselves. Due to the decomposition of 12g,d_l into ~g,d-2 q- ]'~g,d-1, we can replace (7) by another system equivalent to that consisting of Eqs. (7)-(8) where we replace d - 1 and d by d - 2 and d - 1. Recursively, the same decomposition is made for each level k replacing d by k into Eqs. (7)-(8). The process can be iterated till we reach the decomposition 12g,0 + 142g,1. For the sake of generality, we rewrite the system (7)-(8) for the decomposition Vg,k = Vg,k-1 +l/Yg,k, for each level k:

(ayk-1 ~,Yk-1

)

(aZk ) 1 1 + ~ Ot ,Yk-1 + ~ee((Yk-i,Yk-1)) + ~ee((Zk,Yk-1))

q- b(yk-l, Yk-1, Yk-1) + bint (Yk-1, Zk, Yk-1) ---=(f, Yk-1), ( Oyk-1)

---~,"Zk

[" ~Zk ~ "~

1

VYk-1 E ~O,k-1, (9)

1

q- ~--~,Zk; -]- ~ee ((Yk-l,'Zk)) Jr- ~((Zk,~k))

+ b(Yk-1, Yk-1, Zk) q- bint (Yk-1, Zk, "Zk) = (f, Zk),

V'Zk E ]/YO,k.

(10)

During a simulation of the Burgers' problem solved by the classical Galerkin method, we perform several tests in order to evaluate the terms which appear in Eqs. (9)-(10), for each level k, 1 ~< k ~< d. Moreover we estimate the time variation of the terms depending upon the corrections Zk in Eq. (9) relative to Yk-1- The time variation is given by the difference of values over two successive time steps. This means that, here, the hierarchical basis decomposition of the nodal solution is made a posteriori, just for evaluation purposes.

5.1. Small and large scale structures In our simulation we will consider d = 4 hierarchical levels. Fig. 2 shows the evolution of the 12norm a n d / a - n o r m of the ratio zk/Yk-1 (for k = 1 , . . . , d) and also the norm of the nodal solution Ud. The/2-norm of a vector vk = (Ve,l,..., Vk,Nk) C Vk is given by the following relation:

This norm is of the same order as the L2-norm of vk:

IVklL2 = ~/vT M(k)Vk, where M (k) is the mass matrix at level k:

M(,~ ) = (Ca,i, ek,i)"

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

Hier Sol • Z(n) / Y(n)

411

Hier Sol • Z(n) / Y(n)

101

10 ~

10 °

10 ° E o Z

E o Z 10-+ CXl

/

10-1

cO

10 - 2

10 - 2

10 -3

I

i

,

1

2

3

10-3 4

,

i

i

J

o

1

2

3

Time

Time

Fig. 2. 12- a n d / ~ - n o r m s

4

-Nodal U ....... level 1 : z l / y O - - level 2 : z 2 / y l level 3 : z3/y2 level 4 : z 4 / y 3

o f ua and z k / y k - i (k ~ d).

Similarly, the l~-norm of Vk is given by

Ivkl~= where vk,i

sup Ivk,il,

l~i~Nk

1 Vk,i) 2 = (Vk,i, • R 2, and Ivk,il 2 = ((V]~,i)2 + (Vk,2 i) 2 ).

We recall that h k is the width of triangulation Tk, and that h k = h o / 2 k (Fig. 2 corresponds to a coarsest mesh of size h0 = 1/16). We can see that the incremental values z k a r e not of the order of h 2. The estimation of the order of small structures was

Izkl~2

tt

~ h~

ludlt2,

(11)

and this relation shows that the increments zk are not small if the solution displays strong gradients. We know of course that these gradients grow up when the Reynolds number increases. In Fig. 3 we plotted the 12- and/°C-norms of the quantities n--1 =

z~ - z k

At Oz~ ~

= At

;~n

k,

as well as the norms of y~ - y ~ - i and u ~ - u dn-1 . These quantities represent the time variation over one time step of the hierarchical and nodal solutions. Let us now introduce e, the desired accuracy of the computation. We assume that e is a given parameter, which can be chosen bases on a number of practical considerations. When we freeze the small structures zk during a few time iterations, we introduce in the solution an error which depends on the time variation of the corrections Zk. This error can be evaluated to be of order

Ozk(t) 12 cat

Ot

~ e,

412

C. Calgaro et al. /Applied

1o-2

E

Hier Sol Diff : U-h(n) LJ

Numerical Mathematics

23 (1997) 403442

Hier Sol Diff : U-h(n)

- U_h(n-1)

- U_h(n-1)

10-l

1o-3

5 z cu

1o-4

1o-5 L 0

0

1

2

Time

Fig. 3. Z*-and Im-noms

of (us - uz-‘),

(?/on- ytl”‘) and (z; - z;-‘)

3

4

Nodal U __ level 0 : y0 --.-- level 1 : zl - -

level2:z2

_

level3:z3 level4:z4

(k < d).

where c is a positive constant. We will justify the above estimate in Section 6 devoted to the description of the space-time multilevel algorithm. 5.2. The time derivative terms In Fig. 4 we can see the time evolution of the y and z components of the variational term (aud, /at, &). For each time iteration n and for all hierarchical levels k, 1 < k < d, we plotted the Z2-norm of the following quantities:

(%~_l/%

+k-1,j)

(respectively

(%E_,/%

on the coarse grid k - 1 (respectively

re p resents the time derivative of Y&l projected $k,j)) on the fine grid k), whereas (a.~$/%, &_l,j) (respectively

413

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

12 norm of (dy/dt, y-)

12 norm of (dz/dt, y-)

10- 2

10-3

1 0 -3

lo'

10-4

x\/

10-5

lo" 10-7

i

0

i

1

L

2 Time

1 0 -7

3

4

i

i

k

1

2

3

Time

12 norm of ( dy/dt, z-)

12 norm of (dz/dt, z-)

- ...... - - -- --

level level level level

1 2 3 4

10-2 10 -3

10-~

10-4

~'~

~

"

--.._

......,.,,

lo-'

..........

10-5

lo-" 10-~ 0

i

i

i

1

2

3

10-7 4

0

Time

i

i

i

1

2

3

Time

Fig. 4./2-norm of (Oyk-i/Ot, Ck-,,j), (Ozk/Ot, Ck-,,j), (Oyk-,/Ot, Ck,j) and (Ozk/Ot, Ck,j).

(Oz,~St, Ck,j)) is the time derivative of Zk projected on the grid k - 1 (respectively grid k). Obviously, these quantities are vectors. Hereafter the subscripts j will be omitted. Considering Fig. 4, we observe that the quantities (12) verify the following inequalities: ( ~zk

\

Nevertheless, the difference of size between these terms we do not take into consideration the idea of neglecting others. Such approximations done by Marion and Temam mesh width converged to zero and thus could be assumed

is not significant and consequently, certain terms by comparison with in [17] were justified because the to be much smaller than in actual

C. Calgaroet al. / Applied NumericalMathematics23 (1997) 403~442

414

12 norm of (dz/dt, y-)

Variation of (dz/dt, y-)

1 0-2

10 -2

10-3

10 -a

lo-'

10-"

10 -5

10-s 1 0 -6

10-~

1 0 -7

o

i

i

i

1

2

3

I

4

0

1

Time

2 Time

I

3

4

- .... - -

level 1 level 2 level 3 level 4

Fig. 5. Time evolution of

(Ozk/Ot,qSk-l,j) and At(O2zk/Ot2, Ok-l,j).

computations with realistic mesh sizes. Similarly the first numerical implementations done by Laminie and Pascal (cf. [14,15,18]), developed using only two hierarchical levels, showed that the contributions of evolutive term cannot be neglected. Therefore the schemes introduced in [17] are not practical, especially when the velocity vector field presents steep gradients. To overcome this difficulty we aim to define a multilevel procedure which, during a suitable short time period, solves the Burgers' equations (9) on a grid coarser than the nodal grid 7~. Hence we analyze the time variation of (Ozk/Ot, Ck-l,j) which gives a linear interaction between the small scale and the large scale structures in the y-equation on the grid Tk-1. For each level k <~ d and for each time iteration n ~> 2, we evaluate

(.i~Z~ ~ i~t

0z~ -1

()t ' (9k-l'J

) 12 At ( -~

~

zk,i_2z~,Tl+n

At 2

zn-2 k,i ~)k,i, *k-l,j ) 12.

\ i=nk_l+l

In Fig. 5 we can observe that the time variation of (Ozk/i~t, Ck-l,j) over one time step is much smaller than the term itself. More precisely the variation is locally small even though the term itself is not so small. The above result means that the time variation of (i~zk/Ot, Y'k-1) is locally negligible, and this will allow us to freeze it during a few steps of the iterations.

5.3. The diffusive viscous terms As for the time derivative terms, we now analyze the evolution during actual simulations of the Burgers' problem of the different parts of (1/Re)((gd, ua)) corresponding to the splitting associated with the hierarchical basis. For each level 1 ~< k ~< d and for each time iteration n, we display the /2-norm of the following terms (see Fig. 6):

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

12 norm of l / R e * ( ( y , y - ) )

415

12 norm of l / R e * ( ( z , y - ) )

10-2

10 -2

10 -3

10-3

: !/ /

10 ~

lo-' /

10 -s

10 -s

," if " \

f ~ _

~ J

/ /

1 0 -6

0

i

i

i

1

2

3

10 -6 4

0

J

i

i

I

2 Time

3

Time

4 - -

12 norm of 1/Re * (( y , z - ) )

12 norm of 1/Re * ( ( z , z - ) ) - -

10-2

level level level level

1 2 3 4

10 -2

10-3 1

J

10 -3

.....-

lo-'

/

lo-'

i

/

,//

/

/

10-~

10 -s

10 -6

' 1

0

' 2

' 3

10 4 4

' 1

0

' 2

Time

' 3

4

Time

Fig. 6. /Z-norm of the viscous terms (1/re)((yk_l,¢k

l,j)), (l/re)((Zk,qSk-Uj)),

(l/re)((yk-,,¢k,j)),

and

(1/Re)((zk,¢k,j)).

1

1

R ~ ( ( Y ~ - I , ¢ k - I , j ) ) ~--- ~

R---e ((z~, Ck-l,j))

=

( ~ Y~-l,i" V C k - l , i , V C k - l , j \ i=1

R-e

k# k

1 n -R-e ((Yk-1,2/~k,j))

,

Re ((z~,¢k,j))

-

-

1

= -Re -

~

)

Ink-!

Re

-

-

l <~j<~ n k - i ,

,

l <~j<~nk-l,

1+1

(13)

n

i=,

,

Yk-l,i" V C k - l , i , VCk,j

Z z k,i n • VCk,i, VCk,j i~-nk_l+l

,

)

,

n k _ l + l <<.j <~ nk.

416

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

Variation of 1/Re * ( ( z , y~))

12 norm of 1/Re * ( ( z , y - ) )

10-2

10-6 10-6 10 -7

1 0 -3

lo-6

lo-'

10-9 10 -1°

I

10-11

1 0 -5

/

10-12

/ / (

10 -6

i

i

1

0

i

2 Time

,

i

10-13

,

3

4

0

1

2 Time

3 - -

level 1

....

level 2

-

level 3 level 4

Fig. 7. Time evolutionof (1/Re)((zk, ek-l,j)) and At((Ozk/~t,

~bk-l,j)).

The different components of the Laplacian term are of the same order; consequently we cannot neglect certain terms by comparison with other, a sharp difference with spectral methods where some terms automatically vanish. Nevertheless, considering the projection of Laplacian term on the coarse grid, we obtain:

I((zk,¢k-l,j))l,2

<

Analyzing the interactions between Yk-1 and Zk in the diffusive terms we see that their variation over one time step is very small. In Fig. 7 we plotted

\~ \Z ~

1

( ~ k - l , j ) ) l 2 = Ree

z~_l,

-~ - -

Re

n

_ Zk,i ) , V ~)k,i ' V •k_ l,j

\ i=nk_l-F1

The time variation of the vector (1/Re)((Zk,¢k-l,j)) is much smaller than that of vector (Ozk/Ot,¢k-],j). Consequently, as for the time derivative term, we will propose to freeze, during a few time steps, the term (1/Re)((Zk, Y'k-1)) in the y-equation (9).

5.4. The nonlinear terms Concluding the analysis of the hierarchical terms, we consider for each level 1 ~< k ~< d, the time evolution of the/2-norm of the following quantities:

(2:. Calgaro et al./ Applied Numerical Mathematics 23 (1997) 403-442 12 norm of b (u,u,y~)

417

12 norm of b (u,u,z~)

10 -2

10 -2

f ....-

I 0 -3

10 -3 f J

lO-4

10 ~ f J

10- 5

10- 5

, 1

10~

0

, 2

, 3

10~

4

1

Time 12 norm of b (y,y,y~)

4

3 - ....... --

level level level level

12 norm of b_int (y,z,y~)

10-2

1: y 0 _ z l 2: y l z 2 3:y2 z3 4:y3z4

10-2

J 10- 3

J

10 -3 f

lo-4

2

Time

J

J

/ / f J /

10- s

10- s

10 -s

i

i

i

1

2

3

10 -6 4

0

i

i

i

1

2

3

Time

4

Time

Fig. 8. /2-norm of b(Ud,Ud,¢k--hj), b(Ud, Ud,¢kd), b(yk--l,yk--l,¢k--l,j) and bi,t (yk--l,Zk,¢k--l,j).

b (u'~, u d, n ~)k-l,j) = ((ur~ • V)It~, ~k-l,j),

1 <. j <. nk-1,

b(u ,

nk-1 + 1 < i

Ck,j)

=

V)

(14)

nk.

We can establish from the numerical graphics on top of Fig. 8 the following relation:

Ib(u'~, u~, Ck_l,j)ll 2 ~ Ib(u'~, urn, Ck,j) lt2. Hence we direct our attention towards the splitting of the nonlinear term on the coarse grid k - 1. The splitting consists, for 1 <~ j <~ nk-b of b(y'~_l, yr,_l, Ck-l,j) and

bint(YLl, Zrk, ~)k-l,j)~-b(~rd, Urd, C k - l , j ) -

b ( y L l , y L 1 , ~gk-l,j)

=b(y~_ l, z'~, Ck-l,j) + b(z'~, Y~-l, Ck-l,j) + b(z'~, z'~, Ck-l,j).

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

418

12 norm of b_int (y,z,y-)

Variation of b_int (y,z,y-)

10-2

I 0 -2

10 -3

/

10-~

10-s

10-6 10-r 0

i

i

i

1

2

3

10 -7 4

1

0

Time

2 Time

3

4

- ........

Fig. 9. Time evolution of the nonlinear interaction At(abint/OQ(yk-1, zk, ek-l,j).

term

l e v e l 1: y O _ z l

level 2:yl_z2

- -

level 3:y2_z3

.........

level 4; y3_z4

bint(Yk--l,Zk,q~k--l,j) and its variation with time

We can observe that the/2-norm of the vectors b(ud, Ud, ~)k-l,j) and b(yk_l, Yk-1, e k - l , j ) are almost the same (see the left-hand side of Fig. 8). On the other hand, we establish numerically the following inequality: n n n Ib(Yk-,, Yk-1, ek-l,j)l,2 ~ [bint(Yk-1, znk, ~)k-l,j)ll2, which confirms that the contribution of the nonlinear interaction terms bint is smaller than the component of the nonlinear term involving only the large scale Yk-1 (see the bottom of Fig. 8). However, the nonlinear interaction bint (Yk-1, Zk, Yk-1) cannot be neglected in the y-equation (9): indeed their effects on long time computations can modify the behaviour of the large scale structures. 2 Their time variation, represented in Fig. 9, is the following: At ()bint ()t ( Y nk - l ' Z k n' Ok-l'J)

12

n n [ n--1 zn--1 -= Ibint(Yk-l'Zk' ¢k-l,J) -- bint~Yk-l' k ' Ok-l,J) "~ll2"

As for the time derivative and viscous terms, we can consider to freeze the bint-term during some parts of the time iterations when particular relations are verified (we will present several inequalities in Section 6.2, devoted to estimating some parameters of our space-time multilevel procedure). Since the calculation of nonlinear terms are very expensive, we intend to take their effects into account in a simplified manner. Hence, we can expect a significant gain on CPU time. 2This is similar to the famous butterfly effect, well known in meteorology.

419

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

6. The s p a c e - t i m e multilevel algorithm In the previous section we have numerically shown that the time variation of certain terms related to the small structures zk can be considered as small during appropriate time intervals. Let us recall that e is the desired overall accuracy of the computation. From the previous results, we observe that for some hierarchical levels, the variation of the small scale Zk over one time step can be smaller than the tolerance e; so that we can neglect the time variation of the incremental components zk by freezing them to their last value, in the terms corresponding to the interactions between Yk-1 and zk in the y-equation (9) (i.e., the time derivative t e r m (Ozk/Oy, Yk-1), the viscous term (1/Re)((zk, Yk-1)) and the nonlinear convection term bint (Yk-1, Zk, Yk-l)). More precisely as long as their time variation remains smaller than a characteristic quantity depending on e, these terms can be frozen for a while (may be several time steps), without introducing an error larger than the accuracy ~. In this case, these terms become source terms in the y-equation, and during that period of time, there is no need to solve the z-equation (10). In this way we obtain a multilevel algorithm which realizes a completely self-adaptive procedure for solving our problem on Vd. The procedure depends on only one parameter, namely c, and it produces dynamical space-time V-cycles. The structure of a series of V-cycles is presented in Fig. 10. Let us now assume that the approximation Ud(x,t) is known at the time to (to = ( n - 1)At in Fig. 10). In order to characterize a series of V-cycles, we introduce the following numbers: P v and n V. These quantities are determined according to the tolerance ~: (i) P v represents the depth of a V-cycle and it is determined by estimating the coarsest acceptable hierarchical level l; consequently P v = d - I. Till the coarsest level l, the time variation of the small scale structures zk (l < k <~ d) is smaller than the tolerance ~. (ii) n v represents the global number of V-cycles which can be made with the same frozen interaction terms. Let "r(t0) be the length of the time interval during which we retain the depth P v and we freeze the nonlinear terms hint ( y k - l ( t ) , z k ( t ) , ~ k - ~ , j ) at their last value, at time to. After test = ok

n-I

x= ?

n+p

n

d d-I

--n+2

"

d-3 dt dt

I

F

n+p+l [] I p_V 3

d-2

I'~ f t o

test = ?

I x = (n_V * p_V* 2 + 1) * dt Fig. 10. The V-cycles structure.

420

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403~142 determination of Pv we can deduce the time interval tv, during which we perform only one V-cycle

t v = 2AtPv and the global interval T(t0): "r(t0) = n v t v + At. The nonlinear t e r m s hint can be frozen during the interval "r(t0) without losing the order of approximation on the large scales. The estimates for the depth Pv and the characteristic time "r(t0), from which we deduce n v, will be derived in Section 6.2. We remark that a series of V-cycles starts at the nodal level d (see Fig. 10). Indeed, when we verify the estimate for the time variation of hint (Yk-l(tO), zk(tO), q~k-l,j), we must compute the value of b(Ud(tO),Ud(tO), Cd,j). After this estimation, it is not expensive to solve the equation on the nodal grid, at time to + At. Finally, at the end of a complete series of cycles, one to several time steps of integration are made on the nodal mesh 7~. This procedure allows to compute the characteristic numbers for the next series. Certainly, it also plays an important role in our approach, like the projection on a Approximate Inertial Manifold in the original version of the nonlinear Galerkin method, which can be viewed as a filtering procedure. 3 We note here that the linear terms (~)zk/at(t), q~k-l,j) and (1/Re)((zk(t), Ck-l,j)), unlike the nonlinear term bint, are frozen at their last value during an interval less than or equal to t v. Indeed, the small scales zk, as well as the two linear z-terms, are updated during the upward phase of each V-cycle. It must be underlined that the fundamental role of the V-cycles is to allow a constant update of the z contributions. A choice obviously inappropriate is to realize a procedure which, during several time steps of integration, solves the equation on a mesh coarser than the nodal mesh. In this case we could produce a large accumulation of errors in the small structures. Let us consider in more details the integration process for the downward and the upward phases of a single V-cycle. Let us assume that the approximation ua(z, t) is known at time tn, on the nodal mesh 7~. For the downward phase, when we work on the mesh of level k, l < k <~ d, we determine at time tn + (d - k)At the approximation

yn+d-kk ~ yn+d-kk_l _~_Zkn+d-k as the solution of the equation projected on the grid 7~ (see Fig. 11). The incremental components zk are frozen to their value computed at tn + (d - k)At during the process of the downward phase and for the upward phase till level k. Hence they are not updated on the time interval of length 2(k - l)At (if k = d, then 2(k - 1)At = tv). As for the small structures zk, the linear terms ((azk/Ot)(t), Ck-l,j) and (1/Re) ((zk(t), Ck-l,j)) are kept to the values assumed at time tn + (d - k)At and frozen during the time interval 2(k - 1)At. At the end of each V-cycle, we derive an a posteriori estimate of the quantities Pv and r(to). Only if these estimates remain smaller than the tolerance ¢, can we execute the next downward phase; again the value of bint (Yk-l(t), Zk(t), Ck-l,j) is replaced by that at time to. 3 This question will be discussed elsewhere.

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

421

tv j4 t I !

d

D, i i i

I

i

,

!

;

I

Pv = d - l

I I

t,~

tn +Pv A t

t.-~ ( d - k ) A t

tn +' pv + ( k - l ) ) A t i

i

i

i pi

2(k--l)At Fig. 11. Description of a single V-cycle: estimate of the time interval for which the small structures zk are frozen. 6.1. The multilevel procedure

We have determined Ud(X , t) at the time n A t as the solution of the discretized equation corresponding to the nodal level d (see (6)). Let us consider now the split system at the levels d - 1 and d. The second step of a V-cycle, starting at level d, consists in solving the y-equation on the mesh 7d-l, putting each term which depends on Zd on the right hand side. We recall that the z-equation is not solved. Let us enter in a more detailed description of the algorithm. The linear part of the y-equation is discretized with the Crank-Nicholson scheme and we use the explicit Euler method for the nonlinear part. - 1 -I- Zdn - 1 and u~ = Yd-1 n d-Z~d given, w e c o m p u t e at t = (n + 1)At, u~ +1 = For Und - 1 = y nd-1 y~_+~ + z~ +1 by the two following steps: • Step 1. z~ +1 is frozen: z~ +1 = z~.

(15)

n+l is solution, for all Y'd-] in V0,d-1, of • Step 2. Yd-1

( 1 ~-

,

(f, Y'd-1)( --

) + 1 ((Yd-1 Ree\\

z~-

z~-I

At

,Yd-1

)

2

'

l

Yd-1 n

-}- b ( Y~d - l , Y d~ - l , Y d - ' )

n - i 1 ~Zdn - , -- ~.._ee((Zd,~-d_l) ) _ bint (Yd-

Yd-,)-

(16)

422

C. Calgaro et al. /Applied

Numerical Mathematics

23 (1997) 403442

In Eq. (16) some approximations have been made for the z-terms. More precisely: l For the time derivative term: as the time variation of ((azd/iX)(t), &__I) is small, it is approximated by (17) l

l

For the viscous term: due to Step 1, the Crank-Nicholson scheme can be reduced to an explicit Euler scheme on the term (l/Re)((Zd(t),&-I)); For the nonlinear interaction: bi,t(yk_t (t), zk (t), yd-1) is computed for t = to = (n - l)At and stored before the beginning of the series of V-cycle for each mesh k, 1 < k < d.

Remark 1. Due to Step 1, if the evolutive z-term is evaluated at (n + l)At, one obtains

The numerical schemes proposed by Marion and Temam in [17] neglected this evolutive term, by analogy with the schemes performed on the pseudo-spectral case (where the corresponding yk-t and zk are orthogonal for the L2- and H’-norms). The implementation of the nonlinear Gale&in methods, developed by Laminie and Pascal in the context of the finite elements discretization (cf. [14,15]), consisted in finding an approximation uh of the solution, decomposed into ~2h + zh, where y2h and zh are solutions of a highly coupled system (the coarse grid has mesh of width 2h and h is the mesh size of the finer grid). The schemes performed in the finite element case neglected some terms of the coupled system in y2h and zh, especially the evolutive terms

(here, 5& and ?& are functions test defined on the corresponding hierarchical spaces). Moreover, further approximations were performed on the z-equation which became a linear problem in the zh unknown. The authors deduced that these algorithms are performing in the case of sufficiently small Reynolds numbers of else for very small mesh sizes. Otherwise, steep gradients involve some instability and the approximations made on the z-terms are not justified. In the implementation of our multilevel procedure we do not enforce the approximation (18) which would introduce further significant error. For this reason the evolutive term is approximated at its last value, namely we enforce (17). The process of the downward phase stops on the coarsest level of the V-cycle, namely I = d - p,. We determine y,“” at time t = (m + l)At, where m = n + p, more precisely, for all Fl E VO,J, m+l

91

- Yin At

yy+'

+

Yp

2

7iii

+ b(YF,Yln,E)

-

1, by an equation similar to (16);

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

423

Table 1 Description of the multilevel algorithm • Initialization: * construct the family of nested meshes, . assemble the mass and stiffness matrices and the source terms, . initialize the solution, * solve the equation on the nodal mesh 7d determining ua, * set v c y c l e = O. • Time iterations: if (vcycle = 0) then . solve the equation on the nodal mesh 7~ determining Ud, * compute the time variation of the zk-terms, * if ( v a r i a t i o n < e) then evaluate the depth of a V-cycle P v = d - l, evaluate the maximum number of V-cycles n v, set vcycle = l; else

* evaluate and store bi.t(yk-1, zk,yk 1) (for k = d till l + 1), * solve the equation on the nodal mesh Ta determining Ud, * V-cycle series: execute 2 n v P V time steps: for i v ---- 1 to n v o downward phase: for k = d down to 1 + 1 compute Yk-1 solution of (9) on the grid "]~_ t, o upward phase: for k----l + 2 to d - 1

compute Yk-i solution of (9) on the grid 7~-i, for k = d determine ua on the nodal mesh 7~,

o evaluate Az = 2At ~~zt+l 12 and Ab =

Obi.t(yt,zz+l, el,j) z: ~t

o if (az > e) or (hb > e . (~-(t0)) -2) then set vcycle = O, get out of the V-cycle series end if next i v end if.

w h e r e ( F , Yl) is t h e p r o j e c t i o n o n the g r i d 7~ o f all the z - t e r m s f r o z e n o n the p r e v i o u s l e v e l s k = l + 2 , . . . , d. C e r t a i n l y , this t e r m w a s a l s o in the y - e q u a t i o n s r e l a t e d to the a b o v e l e v e l s , till l e v e l d - 2 . I f w e f u l l y n e g l e c t this c o n t r i b u t i o n , the d i f f e r e n c e b e t w e e n o u r a p p r o x i m a t i o n a n d the c l a s s i c a l s o l u t i o n i n c r e a s e s ; m o r e o v e r this t e r m h a s b e e n c o m p u t e d o n the a b o v e l e v e l s a n d o n l y a p r o j e c t i o n o n the present level of solution must be made. F o r the u p w a r d p h a s e , an e q u a t i o n s i m i l a r to (16) is s o l v e d in o r d e r to o b t a i n "Yk-1 - m + l , f o r l + 1 < k ~< d a n d m = n + P v + k - l - 1. W e n o t e that d u r i n g the u p w a r d p r o c e s s the e v o l u t i v e z - t e r m as w e l l as the v i s c o u s z - t e r m w e r e f r o z e n s i n c e 2 ( k - l) t i m e steps. W e c o n c l u d e a V - c y c l e w i t h the i t e r a t i o n o n the n o d a l m e s h 7d, w h i c h d e t e r m i n e s U d ( X , t ) , s o l u t i o n o f (6), at t = n A t + t v . A t the e n d o f e a c h V - c y c l e s o m e i n e x p e n s i v e tests are m a d e in o r d e r to v e r i f y the v a l i d i t y o f the c u r r e n t V - c y c l e series. H e n c e , w e t a k e into a c c o u n t the e v o l u t i o n o f the s m a l l s c a l e s as w e l l as that o f the n o n l i n e a r

424

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

term (these tests verify the estimates of the time variation relatively to zk and bint). If all estimates are smaller than the tolerance e, we can proceed and perform another V-cycle. Otherwise the procedure is stopped and, after one to several time steps of integration on the nodal (fine) mesh, we compute the new characteristic numbers Pv and n v for the next series. The complete algorithm is described in Table 1.

6.2. Time scale estimate for the small structures and for the linear and nonlinear interaction terms After the description of the numerical scheme realizing the space-time V-cycles, we direct our attention towards the criteria determining the depth of a V-cycle as well as the a priori estimates establishing the global time length T(to) during which a series of V-cycles can be executed. Debussche, in a personal communication, provides us the fundamental idea for obtaining the criteria performing the characteristic quantities of a series of V-cycles. Their theoretical justification is lightly different from that presented in the pseudo-spectral case (see [3]).

6.2.1. Estimate of p V, the depth of a V-cycle In the case of a pseudo-spectral discretization, one chooses an integer N which represents the total number of modes retained, relatively to each space direction, in the truncation of the Fourier expansion of the solution. In this case one can define two levels of discretization Nil(t) and Ni2(t), changing during the time, such that Nil(t) < Niz(t) < N. Then, a V-cycle consists of 2(i2 - i l ) + 1 time iterations between the levels Nil and Ni2. These levels are chosen so that suitable estimates of the L 2 and H 1 norms of the ratio ZN~p/yN~p (P = 1 or 2) are verified. Mainly, one requires Izu~2 IL2 <, O, lYN~2 IL 2

where 0 is a constant of the order of e 2. In the finite elements case, [15,17] and our numerical experiments show that the incremental terms Zk can be much larger than h~, where hk is the width of the mesh Tk. This is one of the most important differences with the spectral case: the ratio between the small scale and the large scale components is not as small as the spectral ratio. During a V-cycle, we restrict the solution of the Burgers' equation to the yk_l-equation, where l < k ~< d (see steps (15) and (16)). The corrections Zk are frozen during the time interval 2(k - l)At, where I is the coarsest level of solution (as shown in Fig. 11). Here, we impose that their time variation be smaller than the parameter e. For the upward phase, in order to get Yk (t), t = tn + (Pv + k - 1 + 1)At, we solve the y-equation on the grid 7~, but we commit an error on the initial data. Actually we have replaced the initial data

yk(t) = y k - l ( t ) + zk(t),

where t = tn + (Pv + k - 1)At,

by

+

+ k - 1)At) +

+

k)At).

C. Calgaro et al. /Applied Numerical Mathematics 23 (1997) 403-442

425

Level control "z' *2*dt < e p s 1 0 -2

L

-

level = 1

-

10-3

~

"'~

1/

- -

~ v "-" ',:-',.,,",=,..

level = 2 level = 3 level = 4

j'

1!'-

1 0 -4

0

1

2

3

4

Time Fig. 12. Evolution o f the

12-norm o f 2At(Oz'~/Ot)', here e = 10 - 3 .

From this, we infer an estimate of the error at the level k: tn+(Pv+k--l)At g~

e(k) =

/

[zk(s) - zk(tn + ( d - k)At)]

ds.

tn+(d-k)At

Thus we obtain

e(k) < 2(k - 1)/Xt • suPt -~0Zk(t) t2' with t E [tn + (d - k)At, tn + (py + k -/)At]. Hence, we impose, for k = l + 1 , . . . , d,

e(k) <. 2 ( k - l)Atsup Oz-~k(t) t

<

C.

12

W e observe in Fig. 12 that, for each time iteration m, the quantity

Iz• - zT-11l= -- At ~tk (mAt)z2 increases with decreasing values of the level k. In this case we choose the level l as the coarsest level for which the error e(l + 1) introduced on the initial data is smaller than e. Then we obtain the following estimate: 2At Oz~+l at t2 = 2lz?+l - z?+?l ],2 < e. This test is not expensive and it inhibits the beginning o f a series o f V - c y c l e s w h e n (see Fig. 12).

(20) m _m-- 1 112/>~/2 z~-~d

426

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

In order to control the time variation of zk on the different levels, a sufficient condition is to estimate (20) on the lower level l + 1, at the end of each V-cycle. This condition can be written in the form

2At ~i~Zl+l / ~ o /` + qtv) t2 < e,

for q = 0 , . . . ,n v - 1.

(21)

If (21) is violated, then the V-cycle process is stopped.

6.2.2. Estimate oft(to), the time interval for the transfer terms Now we derive the estimate of the global time length 7-(to) of the whole series of V-cycles on It0, to + T(t0)] during which the coupled z-terms are frozen in the y-equation (9), namely

(ozi~t rtv e, Yk-1 ),

-~ee((Zk(t),Yk_l))

and

bint(Yk_l(t),zk(t),yk_l)

(22)

We recall that we do not freeze these terms in the same manner. Indeed, the transfer terms bint are frozen during several V-cycles and kept at their value at time to (independently of the level k); whereas the two linear terms are updated inside a V-cycle and they are frozen during the time interval 2(k - 1)At, less than or equal to t V. We start by considering the error introduced in a single V-cycle, at level k - 1 (1 < k ~< d), when the computation proceeds from to + mAt to to + (m + 1)At, where m = 0 , . . . , 2Pv - 1. We approximate the exact equation 1

1

(~--l(t),Yk--1) -I- (-~k(t),Yk--1) -~-~ee((Yk-l(t),Yk-1)) -[- ~ee((Zk(t),Y'k-1)) +b(yk-l(t),yk-l(t),~k-1) + bint(Yk-l(t),zk(t),~k-1) = (f, Yk-1) VYk-1 E "~)O,k-l, (23) by the same equation where the coupled terms (22) are replaced by those frozen at their last value, namely

(~zk (tA 0t , w , Y k - l )

R---~((zk(ti),Yk-1))

and

bint(Yk-l(to),Zk(to),Yk-1).

(24)

We note to the time corresponding to the beginning of the V-cycle series and to + mat

ti=

to+(m-2(k-l)-l)At

for the downward phase, for the upward phase.

Hence, the error made on the equation used to compute y k - l ( t ) is exactly given by the integral, over one time step, of the difference between the exact equation (23) and its approximation, namely, for all Yk-1 c V0,k-1,

to+(m+l)At e(m) =

f

to+mAt

( zk zk., , ot (8)-

ds

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

427

to+(m+l)At 1

f to+rnAt

+ R---~

((Zk(8)-- zk(ti),Yk--1))d8

to+(rn+l)At

+

f

[bint(Yk-l(S),Zk(S),~k-1) --bint(Yk-l(to),zk(to),~k-1)] ds.

(25)

to+mAt R e m a r k 2. On the right-hand side of the equation approximating (23), we have not considered the possible z-terms frozen on the previous levels and projected on the current level k - 1. This choice involves an incomplete estimation of the error. Nevertheless (25) and the successive estimates are sufficient to insure the accuracy of the algorithm. Let Y'k-1 : ¢/¢-1,j (1 ~< j ~< n k _ 1); by the mean value theorem

to+(m+l)At

e(rn) <.

sup

tE [to+mAt,to+(rn+1)At]

+

f

Ot2

(s - ti) ds

to+mAt

5zk

1

to+(rn+l )At

tE[to+rnAtS,tUoP(rn+l)At] Ree ( (-~-(t), ~k_l,J) ) l2

f

(s - ti) ds

to+mAt to+(rn+ 1)At

+

sup

0bint -

(Yk-1 (t), Zk (t), ~gk_m,j)

t C [to+mAt, t0+ (m+ 1)At]

f

(s - to) ds.

(26)

to+mAt

Obviously,

f (s- ti)ds <<.f (s - to)ds, and observing that

to+(rn+l)At (8

to)

to+mAt

ds

((rn + 1)At) 2 - (mAt) 2 2

we finally have e(m) ~ (M1 + M2 + M 3 ) ( ( m + 1)At) 2 - (mAt) 2 2 where

02Zk .,

MI=

sup

tC[to,to+tv]

(~-(~),~)k_l,J)

lz ,

(27)

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

428

(28)

tC[to, to+tv] M3=

~,-~'

"

sup -0bint - ~ (Yk-1 (t)," Zk(t),¢k-l,j) tE [to,to+tV]

12.

We find the sum of the inequality (27) from rn = 0 to 2Pv - 1. The sum of these errors, during a single V-cycle, is lower than (M1 + M2 + M3) t y 2"

(29)

In conclusion, if we prescribe (M1 + M2 + M3) < etv2,

(30)

we introduce in the Yk_l-equation, during a single V-cycle, an error lower than e / 2 (this is true for each hierarchical level between the coarsest level I and the nodal level d). Henceforth, we consider the global error introduced when executing the whole series of V-cycles. The error over one time steps is given by inequality (25), but hence-forward m changes from 0 to r(to)/At -- 1. We deduce the following inequality:

e(m) ~< ( M1 + M2 + M3) ((m + 1)At) 2 - (mAt) 2 2

'

(31)

where we have set

M1 =

sup ( ~2Zk (t'l ) 12 te[to+qtv,to+(q+l)tv] ~ ~t 2 t /,ff)k-l,j , O ~ q ~ n v - - l , l
M2 =

O<~q<~nV--1 , l
M3 =

,

sup -tE[to+qtV,tO+(q+l)tv] Re

I Obint

sup

tC[to,to+'r(to)]

--~-(yk-l(t),Zk(t),Ok-l,j)

,,)),2

(32)

12.

l
To find the sum of (31) for m varying from 0 to r(to)/At - 1, we split the right hand side of this estimate in two parts: the two linear contributions and the nonlinear one. For the two quantities deriving from linear terms, the error over the whole series is given by n V times the error committed during one cycle (see the estimate (29)), whereas for the quantity which derives from the nonlinear interaction terms we compute

~'(to)/At-1

~3

K-"

m=0

((m + 1)At) 2 - (mAt) 2 : ~ 3 (r(to))_______~2 2 2

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

Evolution of (T_V)_k (eps=l .e-3) 101

429

Evolution of Tau_k (eps=l .e-3) 101

10°

t

~

,~ , , ~

:

10 °

F

t~\

10-I ~ 0

~/

~,

" 1

:

.........:.:.,.,...---:',"

2

.:" ..

3

,

10 -I

4

0

~lJ

.

.

.

.

.

...--

i

i

,

1

2

3

Time

Time

.-.",..-"

....,

4

--

level 1

....... level 2 - - - level 3 level 4 Fig. 13. Estimate of (tv)k and Tk for k = 0 , . . . , d.

Then, to insure that the global error introduced on the y-equation be lower than e, we prescribe 2

tv

(T(t°))2 < z.

(33)

Consequently, w e will enforce the two following inequalities:

sup { 'c ['o+q'v,'o+(.+l)'vl

e

\~--

Ozk < -

(34)

O<~q<~nV--I , l
and I -i)bint -~-(Yk-l(t),zk(t),¢k-l,j)

sup

l 2 < C(T(to)) -2.

(35)

tC[to,to+~'(to)] l
For the sake of simplicity, we note with the dot the differentiation with respect to t. We estimate numerically the two following quantities:

rk =

Ibint(Yk-l(t),Zk(t),¢k-l,j)ll2

and (tv)k

=

1(Sk(t),ek-l,j)lz2

+

"

In Fig. 13 we can see that 7-k, as well as (tv)k, decreases with decreasing values of the level k (l < k ~< d). Moreover, we note that (tv)k is of the order of rk. From (35) we deduce, obviously,

~-(to) <, rk,

fork=l+l,...,d.

C. Calgaro et al. / Applied NumericalMathematics 23 (1997) 403-442

430

Variation "level 1

Variation "level 2

1 0 -1

10-I

10-2

10-2

1 0 -3

1 0 -3

o

1

2 Time

3

4

o

1

Variation "level 3

2 Time

3 --

4 (var. of b_int) / b int (Var. of z) / z

Variation :level 4

10-1

10-1

,02

10 -2

1 0 -3

1 0 -3

0

1

2

3

4

0

1

Time

2 Time

3

4

Fig. 14. Correlation between mtlbi.t(yz, zz+., Cz,j)[t2/Ibi,t(yg, zl+,, Cz,~)l~ and Atl~klz2/Izk I~=Instead, as the time T(tO) is adjusted to be a multiple of t V, the period of a single V-cycle, we deduce from (34)

-C(to)~_nvtv= v/-fi-f ( v ~

tv) < x/-fi~ (tv)k,

fork=l+l,...,d.

In practice it appears that estimate (35) is generally more restrictive than estimate (34). Henceforth, the global time length of the whole series of V-cycles is written

1( r(to) =

e

Obim/Ot)(yl(to), ZIH-1(tO), (al,j)ll2

)1/2

A sufficient condition is to verify the estimate (35) (but not necessarily (34)) at the end of each V-cycle, only on the coarsest level l. The estimate (35) can be written under the form

~Dbint - (Yl(to + qtv), Zt+I (to + qtv), (at,j) 12 < ¢7(t0) -2,

for q = 0 , . . . , n v -- 1.

Once more, the series of V-cycles is stopped when (36) is violated.

(36)

431

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403~142

Remark 3.

The time variation of the nonlinear interaction terms bint is never computed inside a series of V-cycles, because its estimate is very expensive. In our numerical implementations, inequality (36) is exactly verified for q = 0, i.e., only before the beginning of the V-cycle series. In order to verify the constraint on the nonlinear term I(i3bint/Ot)(yl, Zl+l, ¢l,j)it 2, we put some terms, easily calculable, in the place of the time variation of bint. Naturally these terms must have the same size. The numerical experiments show, as in the pseudo-spectral case, a certain correlation between

Ati(Obint/i~t)(yl,Zt+l,¢l,j)ll 2 Ibint(y~, zt+l, Ct,y)lz~

and

AtlOzk/Otil2

(37)

Izk Iz2

This kind of correlation is accurate on the finer hierarchical levels (see Fig. 14), but we lose the accuracy on the coarsest levels. In the near future our efforts will be concentrated on finding others correlations better than (37), especially in the case of the Navier-Stokes problem.

7. Numerical results In this section we report on the numerical results obtained by using the space-time multilevel algorithm presented in the previous sections, and we compare these results with those obtained by using the classical Galerkin method. Let Y2 be the square domain [0, 1] x [0, 1]. We consider a nodal mesh of size hd = 1/256 and we follow the evolution of the Burgers' flow from T = 0 to T = 4. The initial condition is given by

COS(TrXl) COS(Trx2)

\u2(xl,x2,ol)

(xl + x2)/4

)

The boundary conditions are the following nonhomogeneous Dirichlet conditions:

g2(Xl,X2)

t'lZ2(Xl,X2,0) li~Fl)

and we force the Burgers' flow by a time independent external source f term which acts on only some low frequency components of the velocity field:

(Si(xi,x2)/ = (sin(31rx,) sin(47rx:)+ 3 sin(16~-Xl)sin(87rx2) )

f(xl,xa) = t f2(x,,x2) ]

t

sin(47rxl) sin(27rx2) + 2 sin(lOTrxl) sin(167rx2)

"

The results are presented for the Reynolds number Re = 500, considering at most d = 4 hierarchical levels (h0 = 1/16 is the width of the coarsest mesh 7~). We choose a time step At equal to 10 -3, so that the numerical scheme is found to verify the CFL stability condition At l U l l e dd < C, where c is a constant of the order of the unity.

432

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

Fig. 15. Localization of the nodes /91, ...,/96. Time evolution of CPU time 3000 --

/ --

=

"

Classical method Multilevel method

2000

1000

0

2

3

4

Time

Fig. 16. Time evolution of CPU time. We recall that we actually use the Crank-Nicholson scheme for the linear terms and the Euler forward scheme for the nonlinear terms (this scheme is globally a first-order scheme). Finally we choose the parameter c of the same order as the time step discretization, namely c = 10 -3. R e m a r k 4, The nonlinear Galerkin method was introduced in the context of numerical simulations of turbulent flows on large intervals of time. Our spatial and temporal multilevel algorithm is implemented with success in the case of solutions whose long term behaviour is not merely the convergence to a steady state. But the solution of the Burgers' problem, when we take the right-hand side f equal to zero, converges very fast to a stationary solution depending on the initial condition u0(zl, z2) (see [1,12]). In order to obtain a Burgers' flow sufficiently perturbed, we have imposed a time independent source f . The external force f puts sufficient energy to keep an almost turbulent flow, during a time interval long enough. This forcing is a substitute for the energy transfer phenomena, displayed into the energy spectrum, from the low frequencies towards the high frequencies. These phenomena are an important characteristic of the Navier-Stokes flow, but they are probably absent in the Burgers' flow.

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442 B u r g e r s vector tJekl ." t - 1.0

B u r g e r s vector f i e l d : t = ~.0

_--IAAf,5-,.-e_

~~_./7_/_:_.:..,.../___..~_.... - " - " ~ "¢ : ". g- -. x"f-.-. .-5".-.". /-~' -~' ~ -~

433

. -

~ --.--.-.~-~ "---"

~

-"~\\\"~..~>::~, ~ v -,/

.....

~~ ~ ~ ~~_ . _~~: > .~- _ ~ - - - - - _~L Z -

- -

-'-~_"-~.~'<

~

~Z

- - -

. . ~

Fig. 17. Velocity fields at t = 1.0 a n d t = 2.0; Re = 500 and h a = 1 / 2 5 6 .

Burgers vector field,"

t =3.0

Burger~ vectorfield : t =4.0

LZ'---._--,,7-~'---"---~"-\\"-"-. . . . . .

~ - _ _ _ : ~ - _ ~ _

_ - 1_.

-z-__-.: - -

--.-

_

•

..,'"-"

.Q,.,.,~, " -

.....---~-.-"-~

- _..~-

-

~

. ; . _ Y ~ ~

. . . . . . . . . . . . . . .

....-..._

_

-.,.,,..

~

,

~

'

Z

<

~

,

<

-

,

-

~

~

Fig. 18. Velocity fields at t = 3.0 and t = 4.0; Re = 500 and ha = 1 / 2 5 6 .

Figs. 17 and 18 show the evolution of the velocity field u(xl, x2, t) for t = 1 . 0 , . . . , 4.0. We notice a deep gradient, relatively to the first component of the solution, into the lower part of the domain $2. At t = 4.0 we stop the simulation; the flow has not yet converged to the steady state, but it has almost totally lost its turbulent behaviour.

434

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

411

Evolution of the V-cycles

3 2

E ol

= 0 ~'4 (1) >3 O ----2 t~ -9 1 t" 90 •~ 2000

i

i

500

1000

1500

2000

2500

3000 Iteration number

3500

4000

¢..

,

,

,

i

I

,

,

,

,

I

Fig. 19. Time evolution of the levels in the V-cycles.

Characteristic time 5

~4 "~3

~

>,2 0

~ !

0

500

1000

1

1500

I

2000

N

"~ 3 >l 2

o _,1, 2000

,I,

,I , I, 2500

,1,

,I . . . .

3000 Iteration number

,

,

3500

4000

Fig. 20. Characteristictime r(to) of the V-cycles.

Figs. 19 and 20 represent the time evolution of the characteristic numbers determining the structure of the V-cycles realized by our space-time multilevel algorithm. Fig. 19 exhibits the time evolution of the levels in the V-cycles (we recall that the upward phase of a V-cycle go up until the nodal level that is d = 4 in our simulation). From Fig. 20 we can determine the characteristic time length r ( t o ) of a V-cycle series. More precisely, the time iterations are represented on the x-axis and the depth P v of a V-cycle on the y-axis. The quantity P v is computed before the beginning of the series and it is the same during the whole interval r(to), as far as the estimates (21) or (36) are verified. We can observe that, except during a transient time interval after the initial period, "r(to), and then n v , are rather large. Hence this graph confirms that several V-cycles are realized with no update of the

435

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

Table 2 Coordinates of nodes retained in the computational domain /)1

(0.320, 0.394)

/92

(0.605, 0.375)

t93

(0.929, 0.745)

t94

(0.488, 0.492)

P5

(0.746, 0.488)

P6

(0.683, 0.242)

Table 3 12- and 1~ - n o r m of the vector components of u cG

-

UNLG

lu cG -- uNLGll2

lU CG -- uNLGI/2

lUlCG -- uNLGIcx~

lU CG -- uNLGIcx~

1.0

0.00304

0.00187

0.07633

0.02263

2.0

0.00239

0.00159

0.05672

0.01489

3.0

0.00264

0.00158

0.05556

0.02255

4.0

0.00387

0.00299

0.09615

0.09304

Time

Evolution of U1 1.0

,,s"

......

:'

0.5

"-. ,,,.,

0.0 -0.5 " --,,

-1.0

........................................ J

500

0

1000

1500

2000

2500

3000

3500

4000 P1 P2 P3

Difference o f s o l u t i o n s ( e p s V-cycle = 1.e-3) ._1 Z

(.9 v !

1 0 -1

1 0-2

- ~=~"..;': ...~.

.~ . . . . .

..-..,

10 -3 10-4

!

0 10-s (.9 10-6

,!

. :t

v ,i.==

1 0-~

0

500

1000

1500 2000 2500 Iteration number

Fig. 21. Time evolution of the difference u cC

-

3000

3500

UNLG at nodes P1, P2, P3.

4000

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

436

Evolution of U1 1.0

0.5 0.0 -0.5 -1.0

0

500

1000

1500

2000

2500

3000

3500

4000

Difference of solutions (eps V-cycle = 1 . e - 3 )

P4 ...... P5 P6

~-, 10 -1

!

1 0 -2

10 -a 10 -4

~__ 10 -5 10 -6 10 -7 0

500

1000

1500 2000 2500 Iteration number

Fig. 22. Time evolution of the difference u c°

-

3000

3500

4000

UlNLG at nodes P4, /:'5, Pr.

nonlinear terms bint. Concluding, at the end of a whole series the process begins again on the nodal grid, with Pv = O. R e m a r k 5. We described briefly in Section 6.2.1 the structure of a V-cycle in the case of a pseudospectral discretization. A V-cycle consists of 2(/2 - i l ) + 1 time iterations between two levels of discretization Nil and Ni2, Nil < N e < N. In this case one can choose several levels between Nil and Ni2. This is an important difference with the finite elements case where, for a quasi-structured mesh, at each projection on a coarser mesh we neglect 3 / 4 of values of the solution. To perform the comparison between the multilevel algorithm and the classical Galerkin method, we have retained six different points of the computational domain. Their coordinates are written in Table 2 and their position is plotted in Fig. 15. They are vertices of the nodal mesh 7"d (i.e., Pi E Ea for i ---- 1 , . . . , 6) and some of them are close to the regions where the strong gradient of the solution appears. We have stored the value of the horizontal component of the velocity Ul(Xa, x2, t) at these points during the time evolution (see the top of Figs. 21 and 22). At the bottom of the same figures we have plotted the quantity uC1G(Pi) --ulNLr(Pi) which represents the difference between the first component of the velocity computed with the classical method (marked CG) and the component approximated with the multilevel method (marked NLG). The difference between the trajectories obtained with the two different algorithms exceeds sometimes the tolerance e. However, when we perform the series of V-cycles, the error accumulated along the

C. Calgaro et al. /Applied Numerical Mathematics 23 (1997) 403-442

437

Evolution of U1 1.0

:,. ........

,........

0.5 0.0 -0.5 .,

-1.0

0

.,

500

1000

1500

2000

2500

3000

3500

4000

Difference of solutions .... 1 1 0-2 5

i

m

,,. . . . . . '.

• •

,,' :i

!

ii ,,

.--~~

P1(256) P2(256) P3(256) P2(128)

. ~ .. :,.~. . . . . . . . . . . .

0--3

10_4

0_s

10_8 110-7

i I

I

,

0

500

Fig. 23. T i m e e v o l u t i o n

~ i

1000

' ,

I

L

i

i i

I

1500 2000 2500 Iteration number

of the difference between

uco(ha

= 1/256)

] ,

3000

and

uCC(ha --

I

3500

1/128)

II

~

4000

at n o d e s P I ,

P2, P3.

iterations is not so large as cause the blow up or the divergence of the approximate solution. Indeed the simulation of the flow obtained with the multi-scale process is satisfactory: the trajectories of the multilevel method remain close to those obtained with the classical Galerkin method. Table 3 gives, for each component of the velocity, the 12- and l~-norms of the difference between the solutions computed with the CG and NLG methods. We can see that the vector field u(z, t) is globally well approximated, indeed the/2-norm of the difference is only lightly greater than e. On the contrary t h e / ~ - n o r m of this difference is large. Considering this result, we contemplate to study in the future the local behaviour of the linear and nonlinear terms that we would like to freeze during the time iterations. Up to now we have kept into account only the/2-norm of these terms, but this norm is an average measure of these terms and it tends to mask the steep gradients localized upon a small size of the computational domain. We can take into account this local analysis by using some adaptive hierarchical meshes in order to sufficiently discretized the domain ~2 in regions of steep gradients. Finally, we show in Fig. 16 the time evolution of the CPU time required by CG and NLG methods. To achieve the comparison between the different algorithms, we want to note the significant fact that the multilevel method requires nearly 30% less CPU time than the classical Galerkin method. The saving in given by the ratio Tcc - TNLC

TCG

'

438

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403--442

Evolution of 01 1.0 H

0.5

0.0 -0.5 -1.0

J

I

0

500

,

I

1000

,

I

1500

,

J

,

I

2000

,

2500

I

3000

~

I

L

3500

4000 -

Difference of solutions .... 10 -1

~

D '

10-a

~

0-5

e~

0-6

.

.

.

.

.

. . . .

P5 P6

•",, ,-.,

10 -2 v .li.

P4

.

1 0 -4 II

~=

'

I

I

(

'

i\

it

\

II

I

:D 10-7 0

500

1000

1500 2000 2500 Iteration number

3000

3500

4000

Fig. 24. Time evolution of the difference between uC6(ha = 1/256) and uCG(ha = 1/128) at nodes P4, Ps, P6.

where TcG is the CPU time required by the classical Galerkin method and TNLG the CPU time required by the multilevel method. R e m a r k 6. In order to justify the mesh size chosen for the nodal grid, we compare two approximate solutions of the Burgers' problem determined by the classical Galerkin method; the first solution is obtained using a space step hd = 1/256 and the second solution is computed using hd = 1/128 (obviously the Reynolds number and the time step remain unchanged). The approximate solution with mesh size hd ---- 1/128 is not totally satisfactory because of the presence of some fluctuations near the zone where the solution exhibits steep gradients. These fluctuations appear when we consider a too large mesh size and they disappear when we consider hd = 1/256. Moreover we compute the difference, at time t -- 2.0 and t = 4.0, between Ul(Xl, x2, t) approximated with CG and a mesh size h d = 1/256 and U l (X l, X2: l~) approximated with the same method but with a mesh size twice greater: hd = 1/128. We show this difference on the top of Figs. 25 and 26. On the bottom of the same figures we represent the difference between ul computed with the CG method and Ul computed with the multilevel method, again at t = 2.0 and t = 4.0, using the same space step hd = 1/256. We want to note the significant fact that the difference between Ul computed with CG and NLG method is smaller than the difference between ul computed with the CG method and different mesh sizes. One must observe that these quantities are plotted with different extremes on the z-axis.

439

C. Calgaroetal./AppliedNumericalMathe~tics23(1997) 403~42 Di

Di

f

f

f

f

.

.

So

So

I

I

.

.

G C ( 1 2 9 )

G C ( 2 5 7 )

~-~-

--

--

~

G C ( 2 5 7 )

;

G N L ( 2 5 7 )

Y

;

=

T

2.

=

0

2 . 0

--4

Fig. 25. On top: difference between u~ (t = 2) computed with CG method and ha = 1/256 or ha = 1/128. At the bottom: difference between u~ (t = 2) computed with CG and NLG method; ha = 1/256 is the same.

Finally, Figs. 23 and 24 exhibit the evolution o f the difference b e t w e e n u cG (ha = 1 / 2 5 6 ) and u cG (ha = 1 / 1 2 8 ) , difference established at the points P 1 , . . . , / 9 6 (see Table 2). T h e s e figures can be c o m p a r e d with the Figs. 21 and 22, w h i c h give the difference UlcG (Pi) - ulNLC(P~) c o m p u t e d with the s a m e m e s h size ha = 1 / 2 5 6 .

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

440

Di

f

f

.

So

I

.

G C ( 1 2 9 )

--

G C ( 2 5 7 )

;

T

=

4 . 0

o

ib

o

Di

f

f

.

So

I

.

G C ( 2 5 7 )

--

G N L ( 2 5 7 )

;

T

=

4 . 0

Lq

0 0 O

~- ~

~, ~ ~ > ~ . . ~ o .

-~'~o-

~

Fig. 26. On top: difference between us (t = 4) computed with CG method and ha = 1/256 or ha = 1/128. At the bottom: difference between ul (t = 4) computed with CG and NLG method; ha = 1/256 is the same.

8. Concluding remarks In this article we have proposed a space-time multilevel method in the context of finite elements. This method treats differently the large and small components of the flow which is solution of the

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403~t42

441

two-dimensional Burgers' problem. This decomposition is induced by the utilization of the hierarchical bases. We have obtained with this multilevel algorithm the first results of an application of the nonlinear Galerkin method in the context of a finite element discretization using several hierarchical levels. In the previous numerical applications of the original schemes proposed by Marion and Temam, Laminie et al. [14,15] observe that the approximations made in the y-equations and in the z-equation are not always justified from a numerical point of view. Although the size of certain coupled terms can be small, they cannot be neglected, especially when the solution exhibits steep gradients in the computational domain. The new suggestion here is to consider instead the time variation of the linear and nonlinear terms. The subsequent multilevel algorithm, which freezes the coupled terms during suitable time intervals, converges with a good efficiency, producing a significant saving of computing time. Several numerical estimates of terms which appear in the equations have been performed before the conception and the realization of the multilevel algorithm. Moreover we have derived mathematical estimates for the parameters involved in this dynamical procedure. In the case of flows which present deep gradients, it will be possible to envisage a local analysis of suitably frozen terms; this analysis can be obtained using adaptive hierarchical meshes. In the future our efforts will be concentrated on the criteria for determining the structure of the V-cycles as well as finding better correlations concerning the nonlinear interaction terms, in order to decrease the perturbations introduced by the multi-scale procedure. Other efforts will be concentrated on the implementation of a second-order scheme for the time discretization of the linear and nonlinear terms, and also on the implementation of our multilevel algorithm in the case of the Navier-Stokes problem.

References [1] C. Calgaro, Mrthodes multi-rrsolution auto-adaptatives en 616ments finis: application aux 6quations de la m6canique des fluides, Th~se Univ. Paris-Orsay (1996). [2] P. Constantin, C. Foias, O. Manley and R. Temam, Determining modes and fractal dimension of turbulent flows, J. Fluid Mech. 150 (1985) 427~,40. [3] A. Debussche, T. Dubois and R. Temam, The nonlinear Galerkin method: a multi-scale method applied to the simulation of homogeneous turbulent flow, Theoret. Comput. Fluid Dynamics 7 (1995) 279-315. [4] T. Dubois, Simulation numrrique d'rcoulements homog~nes et non-homog~nes par des mEthodes multirrsolution, Thrse Univ. Paris-Orsay (1993). [5] T. Dubois, E Jauberteau and R. Temam, Solution of the incompressible Navier-Stokes equations by the nonlinear Galerkin method, J. Sci. Comput. 8 (1993) 167-194. [6] T. Dubois, E Jauberteau and R. Temam, Dynamical multilevel methods in turbulence simulation, Computational Fluid Dynamics Review, to appear. [7] C. Foias, M.S. Jolly, I.G. Kevrekidis, G.R. Sell and E.S. Titi, On the computation of inertial manifolds, Phys. Lett. A 131 (1988) 433-436. [8] C. Foias and R. Temam, Some analytic and geometric properties of the solutions of the evolution NavierStokes equations, J. Math. Pures Appl. 58 (1979) 339-368. [9] C. Foias, G.R. Sell and R. Temam, Inertial manifolds for nonlinear evolutionary equations, J. Differential Equations 73 (1988) 309-353. [10] C. Foias, O. Manley and R. Temam, Sur l'interaction des petits et grands tourbillons dans les 6coulements turbulents, C. R. Acad. Sci. Paris S~r. I Math. 305 (1987) 497-500. [11] C. Foias, O. Manley and R. Temam, On the interaction of small and large eddies in two dimensional turbulent flows, Math. Model. Numer. Anal. 22 (1988) 93-114.

442

C. Calgaro et al. / Applied Numerical Mathematics 23 (1997) 403-442

[12] EC. Jain and D.N. Holla, Numerical solutions of coupled Burgers' equations, Internat. J. Non-Linear Mech. 13 (1978) 213-222. [ 13] E Jauberteau, R6solution num6rique des 6quations de Navier-Stokes instationnaires par m6thodes spectrales. M6thode de Galerkin non lin6aire, Th~se Univ. Paris-Orsay (1990). [14] J. Laminie, E Pascal and R. Temam, Implementation of finite element nonlinear Galerkin methods using hierarchical bases, J. Comput. Mech. 11 (1993) 384-407. [15] J. Laminie, E Pascal and R. Temam, Implementation and numerical analysis of the nonlinear Galerkin methods with finite elements discretization, Appl. Numer. Math. 15 (1994) 219-246. [16] M. Marion and R. Temam, Nonlinear Galerkin methods, SIAM J. Math. Anal. 26 (1989) 1139-1157. [17] M. Marion and R. Temam, Nonlinear Galerkin methods: the finite elements case, Numer. Math. 57 (1990) 205-226. [18] E Pascal, M6thodes de Galerkin non lin6aires en discr6tisation par 616ments finis et pseudo-spectrale. Application h la m6canique des fluides, Th~se Univ. Paris-Orsay (1992). [19] R. Temam, Sur l'approximation des 6quations de Navier-Stokes, C. R. Acad. Sci. Paris Ser. A 262 (1966) 219-221. [20] R. Temam, Navier-Stokes Equations, (North-Holland, Amsterdam, 1984). [21] R. Temam, Infinite Dimensional Dynamical Systems in Mechanics and Physics, Applied Mathematical Sciences 68 (Springer, New York, 1988). [22] R. Temam, Inertial manifolds and multigrid methods, SIAM J. Math. Anal. 21 (1990) 154-178. [23] H. Yserentant, On the multi-level splitting of finite element spaces, Numer. Math. 49 (1986) 379-412. [24] O.C. Zienkiewicz, D.W. Kelly, J. Gago and I. Babugka, Hierarchical finite element approaches, error estimates and adaptive refinement, in: J.R. Whiteman, ed., The Mathematics of Finite Elements and Applications IV (Academic Press, London, 1982).