- Email: [email protected]
Multilevel stabilization devices in CFD
Computers & Fluids 31 (2002) 421–436 www.elsevier.com/locate/compfluid
Multilevel stabilization devices in CFD q Claudio Canuto Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Turin, Italy Received 1 August 2001 Dedicated to Professor Roger Peyret
Abstract The numerical discretization of such fundamental equations in fluid dynamics as the Navier–Stokes equations, requires special attention: a careful treatment of the pressure and the convective terms is needed, to avoid the onset of spurious oscillations. We advocate the use of multilevel bases (wavelets, or hierarchical finite elements) to build-up consistent (i.e., residual based) stabilization devices. The proposed strategy is particularly well-suited for being applied within an adaptive discretization procedure. Classical SUPG-type methods can be interpreted as special realizations of the multilevel stabilization idea. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Coerciveness; Inf–sup condition; Stabilization; Multilevel bases; Stokes problem; Convection–diffusion problems
1. Introduction Multilevel bases include hierarchical finite elements [5], wavelets and other wavelet-like families [18], as well as the bases associated with sparse grid methods [4] or incremental unknown methods [17]. Although such bases are less simple and easy-to-handle than standard bases such as Lagrange’s, yet they enjoy extra properties, which make their use in numerical simulation very attractive. Significant applications of multilevel finite element and finite difference techniques have
q
This work was partially supported by the European Commission within the TMR project (Training and Mobility for Researchers) Wavelets and Multiscale Methods in Numerical Analysis and Simulation, no. ERB FMRX CT98 0184, and by the Italian funds Murst Cofin’98 Analisi Numerica. E-mail address: [email protected] (C. Canuto). 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 3 - 9
422
C. Canuto / Computers & Fluids 31 (2002) 421–436
been given over the years; the use of wavelet-like methods is currently under investigation, with various promising results [22]. Multilevel bases provide a natural capability of preconditioning symmetric positive-definite problems, as they yield stiffness matrices which are spectrally equivalent to diagonal matrices (see, e.g., Refs. [5,8,23]). The possibility of approximating a function (a signal) in an optimal way by adding and removing details [25] can be exploited in designing efficient reconstruction algorithms, with applications to, e.g., the numerical discretization of conservation laws [1,27]. Details (which form the multiscale representation of an approximate solution, or a residual) are inherently error indicators; sound adaptive strategies can be based on their inspection [19]. Related to this, is the recently proposed and stimulating attempt to separate coherent and incoherent structures in the numerical simulation of a turbulent flow, by means of wavelets [26]. The present contribution focusses on another feature of multilevel bases, namely the possibility of representing inner products and norms of function spaces of any regularity order in an explicit and computable form, i.e., by weighted sums of subsequent details [21,32,33]. This property, which somehow generalized the Parseval identity in Fourier analysis, is very important; indeed, often norms and inner products of fractional and/or negative order are defined in a way that is not at all amenable to an efficient numerical treatment. By means of multilevel bases, these mathematical objects become implementable. Many boundary value problems of interest lack of the good symmetric and positive-definite (coercive) character of classical elliptic problems. This is a major cause of difficulty in their numerical treatment. In fluid dynamics, spurious velocity oscillations may appear in the numerical simulation of flows at high Reynolds numbers; they are allowed by the insufficient control provided by the diffusion terms. Furthermore, pressure instabilities may rise in the simulation of incompressible flows, whenever the discrete trial velocity and pressure spaces fail to satisfy a particular form of coerciveness condition, the well-known LBB or inf–sup condition. Over the years, many different approaches have been proposed to transform a boundary value problem into an equivalent one with enhanced coerciveness; SUPG-type stabilization methods, least-square and Galerkin/least-square methods are among the examples. The general idea of consistent stabilization amounts to augmenting the standard variational formulation of the problem by another variational term, which is in some sense ‘proportional’ to the residual of the equation. The effect is to provide a tighter control on the solution (for instance, in a stronger norm). Classical stabilization is accomplished at the discrete level; however, the problem can also be set in an abstract way on the exact operator equation [3,6]. It is here that multilevel inner products and norms come usefully into play (see, e.g., Refs. [7,11]). In the first part of this paper, we present an abstract framework in which consistent stabilizations of weakly coercive or non-coercive problems can be cast. The setting is particularly geared towards multilevel-based stabilizations. Next, we illustrate two applications to the numerical discretization of model problems of interest in fluid dynamics: the Stokes problem and the convection–diffusion problem. In the first case, the stabilization proposed in Ref. [14] allows a fully independent discretization of velocity and pressure, a feature which is fundamental if the numerical method incorporates an adaptive strategy. Uniform control on the exact L2 -norm of the discrete pressure (up to the physically indetermined constant mode) is guaranteed. For the convection–diffusion problem, following [11–13], we advocate the use of a fractional order norm (precisely, of order 1=2) to get a uniform control on the streamline derivative of the discrete
C. Canuto / Computers & Fluids 31 (2002) 421–436
423
solution. In addition, our functional framework suggests to view classical SUPG methods as truncated versions of the proposed multilevel stabilization. Notation: Throughout the paper, we shall use the notation A K B to indicate A 6 cB for a suitable constant c, independent of the relevant parameters for the particular formula in which the symbol is used (e.g., the diffusion parameter m, the multiscale parameters k, the discretization parameter h). The notation A B means A K B and B K A.
2. An abstract setting Let us suppose that we aim at numerically solving an operator equation of the form Lu ¼ f ;
ð2:1Þ
where L is linear continuous between a Hilbert space V and the dual W 0 of a (possibly different) Hilbert space W. Let h; i denote the duality pairing between W 0 and W. Introducing the continuous bilinear form a : V W ! R;
aðv; wÞ :¼ hLv; wi;
the operator equation can be written variationally as u 2 V : aðu; wÞ ¼ hf ; wi 8 w 2 W :
ð2:2Þ
It is classically known [2] that the well-posedness of this problem is equivalent to the two inf–sup conditions inf sup v2V w2W
aðv; wÞ > 0; kvkV kwkW
inf sup
w2W v2V
aðv; wÞ > 0: kvkV kwkW
ð2:3Þ
Let us introduce a finite dimensional approximation of this problem. To this end, for h 2 ð0; hmax , let Vh V and Wh W be subspaces satisfying dim Vh ¼ dim Wh < þ1, and let us consider the Petrov–Galerkin discretization of Eq. (2.2) uh 2 Vh : aðuh ; wh Þ ¼ hf ; wh i
8 w h 2 Wh :
ð2:4Þ
The solvability of this problem, as well as the behavior of the solution, is related to the (best) continuity constant j > 0 and inf–sup constant ah P 0, which are defined through the inequalities ah kvh kV 6 sup
wh 2Wh
aðvh ; wh Þ 6 jkvh kV kwh kW
8 vh 2 Vh :
ð2:5Þ
Note that j can be chosen independently of h, as a consequence of the continuity of the form a on V W ; in addition, it is not restrictive to assume that j 1, possibly after a proper scaling of the equation. Conversely, the inf–sup conditions (2.3) for the exact problem may imply nothing on ah : neither its independence of h, nor even its being different from 0. If ah ¼ 0, problem (2.4) is not well-posed. This is, e.g., the case of the approximations of the Stokes problem, in which the discrete velocity and pressure spaces do not satisfy a suitable compatibility condition; as a consequence, the discrete pressure solution is defined up to spurious
424
C. Canuto / Computers & Fluids 31 (2002) 421–436
oscillating components (see, e.g., Ref. [9]). On the other hand, if 0 < ah j, problem (2.4) is indeed well-posed, but the classical error estimate j inf ku vh kV ku uh kV 6 1 þ ah vh 2Vh allows the discrete solution to be far apart from the best approximation in Vh of the exact solution. This is the case, e.g., of pure Galerkin approximations to convection–diffusion problems in the convection-dominated regime (see, e.g., Ref. [34]). The previous considerations motivate the development over the years of a great number of strategies, which modify the variational problem (2.4) in a consistent way (i.e., in such a way that the solution of Eq. (2.2) is also the solution of the modified problem). The aim is to improve the constant a~h which appears in the analog of Eq. (2.5) for the modified variational form a~; more precisely, a~h should possibly be bounded away from 0 independently of h, and have the same order of magnitude as the continuity constant j~. We propose an abstract functional framework which leads to stabilization strategies based on the multilevel representation of a Hilbert space. Let us denote by Z 0 a Hilbert space continuously imbedded into W 0 , with inner product ðð; ÞÞZ 0 . Let us set Ve :¼ fv 2 V : Lv 2 Z 0 g equipped with the graph norm kvkV~ :¼ kvkV þ d1=2 kLvkZ 0 : Here, d > 0 is a parameter, whose actual choice will depend on the particular application but whose order of magnitude should be thought to be in the order of unity. From now on, let us suppose that the data f, assigned in Eq. (2.1) as an element of W 0 , be indeed in Z 0 . This immediately implies that the solution u belongs to Ve . e W with continuous imbedding, and a linear Next, let us introduce a second Hilbert space W 0 e operator S 2 Lð W ; Z Þ. By such operator, we define the modified bilinear form e ! R; a~ : Ve W
a~ðv; wÞ :¼ aðv; wÞ þ dððLv; SwÞÞZ 0 ;
as well as the modified linear form e ! R; f~ : W
f~ðwÞ :¼ hf ; wi þ dððf ; SwÞÞZ 0 :
This modification is consistent with Eq. (2.1), in the sense that u also solves the variational problem u 2 Ve : a~ðu; wÞ ¼ f~ðwÞ
e: 8w2W
ð2:6Þ
We now assume that the finite dimensional subspaces Vh V and Wh W previously introduced, e , respectively. Thus, we consider the Petrov–Galerkin approxiare also included into Ve and W mation of problem Eq. (2.1) uh ; wh Þ ¼ f~ðwh Þ u~h 2 Vh : a~ð~
8 wh 2 Wh :
ð2:7Þ
Let us discuss the continuity and inf–sup constants of the form a~. As far as the former one is e concerned, one has for any v 2 Ve and w 2 W
C. Canuto / Computers & Fluids 31 (2002) 421–436
425
j~ aðv; wÞj 6 jkvkV kwkW þ dkLvkZ 0 kSwkZ 0 6 c1 jkvkV kwkW~ þ c2 d1=2 kLvkZ 0 kwkW~ 6 j~kvkV~ kwkW~ ; e W , c2 is the norm of the operator d1=2 S 2 Lð W e ; Z0Þ where c1 is the norm of the imbedding W and e j :¼ maxðc1 j; c2 Þ. Concerning the inf–sup constant, we make the following Assumption 2.1. The space Z 0 and the operator S are chosen in such a way that the inf–sup condition a~kvh kV~ 6 sup
wh 2Wh
a~ðvh ; wh Þ kwh kW~
8 vh 2 Vh ;
ð2:8Þ
is fulfilled with a constant a~ > 0 independent of h, having the same order of magnitude as j~. This yields a constant of order of magnitude 1 in the error estimate ! j~ ku u~h kV~ 6 1 þ inf ku vh kV~ ; a~ vh 2Vh thus expressing the enhanced stabilization property of the modified approximation. In the rest of the paper, we will consider two relevant examples of the above functional setting.
3. The Stokes problem Given a bounded domain X Rn (n P 2) with Lipschitz boundary oX, we consider the Stokes problem with homogeneous boundary conditions: ~ þ rP ¼ ~ DU f in X;
ð3:1Þ
~¼0 rU
ð3:2Þ
in X;
~ ¼ 0 on oX: U ð3:3Þ n 1=2 ~ :¼ ðH 1 ðXÞÞ equipped with the norm k~ vkX~ ¼ ðr~ v; r~ vÞ . Here and in the sequel, we Let us set X 0 indicate by ð; ÞX both the inner product in L2 ðXÞ and the one in ðL2 ðXÞÞn . Let h; i denote the ~. Let M be a fixed subspace of L2 ðXÞ such that ~0 ¼ ðH 1 ðXÞÞn and X duality pairing between X 2 L ðXÞ ¼ M span (1); we denote by kqkM ¼ ðq; qÞ1=2 X the norm of M. 0 ~ ~ We assume that f is given in X , which implies the existence and uniqueness of the solution ~ ; P Þ in X ~ M. Setting ðU w; pÞX þ ðr ~ v; qÞX A½ð~ v; pÞ; ð~ w; qÞ ¼ ðr~ v; r~ wÞX ðr ~ and F ð~ w; qÞ ¼ h~ f;~ wi;
426
C. Canuto / Computers & Fluids 31 (2002) 421–436
~ ; P Þ as the solution of the variational problem we obtain ðU ~; P Þ 2 X ~ M : A½ðU ~ ; P Þ; ð~ ðU w; qÞ ¼ F ð~ w; qÞ
~ M: 8 ð~ w; qÞ 2 X
~ M equipped with Thus, in the abstract setting of the previous section, we choose V ¼ W :¼ X the graph norm; next, we identify L with the Stokes operator, i.e., we set D r L :¼ ; ð3:4Þ r 0 finally, we choose the forms a and f to be A and F, respectively. ~ and Mh M be finite dimensional subspaces. In the adaptive approximation of the ~h X Let X Stokes problem, it may be convenient to allow as much mutual independence as possible in the selection of the active degrees of freedom for velocity and pressure. Consequently, it might happen that the classical inf–sup condition bkqkM 6 sup ~h ~ w2X
ðr ~ w; qÞX k~ wkX~
8 q 2 Mh ;
is not satisfied (i.e., b ¼ 0), or it is satisfied with b > 0 but depending upon the discretization. In this case, the resulting discrete pressure is affected by spurious or weakly-spurious modes. In order to avoid such a drawback, a stabilization procedure of the type (2.6) has been recently proposed in Ref. [14]. Let us briefly recall the idea. Assume that the basis functions in Mh are selected within a Riesz basis Wp ¼ fwpk : k 2 Kp g of M, in the sense that there exists a finite set p p p Kph Kp such that k : k 2 Kh g. We recall that W is a Riesz basis if any q 2 M can be P Mh ¼ spanfw p written as q ¼ k2Kp qk wk with q ¼ ðqk Þk2Kp 2 ‘2 ðKp Þ (the series being convergent in L2 ðXÞ), and additionally kqkM kqk‘2 ðKp Þ
8 q 2 M:
ð3:5Þ
e p of L2 ðXÞ (the subspace of the zero By the Riesz representation theorem, there exists a basis W 0 p mean valued functions), which is biorthogonal to W , in the sense that ðwpk ; w~pk0 ÞX ¼ dk;k0
8 k; k0 2 Kp :
ð3:6Þ
Consequently, any q 2 M can be represented as X qk wpk with qk ¼ ðq; w~pk ÞX : q¼ k2Kp
~ satisfying w k 2 X Let us associate to each k 2 Kp an auxiliary velocity ~ r~ w k ¼ w~pk
in X:
Furthermore, let us assume that X qk ~ w k K kqk‘2 ðKp Þ 8 q ¼ ðqk Þk2Kp 2 ‘2 ðKp Þ: k2Kp
ð3:7Þ
ð3:8Þ
~ X
An example of construction of such auxiliary velocities, starting from a wavelet basis of the pressure space, will be detailed in the next subsection.
C. Canuto / Computers & Fluids 31 (2002) 421–436
427
~ :M!X ~ as Using the auxiliary velocities, we define the operator S X X ~q ¼ q¼ qk wpk 7! S qk ~ w k : k2Kp
ð3:9Þ
k2Kp
Proposition 3.1. There exists a constant c > 0 such that ~ qk~ 6 c kqk kS M X
8 q 2 M;
ð3:10Þ
furthermore, there exists a constant b > 0 such that 2
~ qi P b kqk hrq; S M
8 q 2 M:
ð3:11Þ
Proof. Inequality (3.10) follows from (3.8) and (3.5), whereas (3.11) is a consequence of Eq. (3.7) and again Eq. (3.8). We are now ready to introduce our modified formulation of the Stokes problem, which reads as ~2X ~ and p 2 M such that follows: find U ~ ; r~ ~; ðrU wÞX ðr ~ w; P ÞX ¼ h~ f;~ wi 8 ~ w2X ~ qi ¼ dh~ ~ qi ~ ; qÞ þ dhDU ~ rP ; S ðr U f;S X
ð3:12Þ 8 q 2 M:
ð3:13Þ
The reader should note formal analogies and actual differences of this formulation with the ~h and Galerkin/least-squares finite element method presented in Ref. [28]. Precisely, assume that X Mh are spaces of continuous finite elements on a triangulation Th whose maximum diameter of triangles is h; let hT denote the diameter of the triangle T. That formulation reads as follows: find ~h 2 X ~h and ph 2 Mh such that U ~h ; r~ ~h ; ðrU wh ÞX ðr ~ wh ; Ph ÞX ¼ ð~ f;~ wh ÞX 8 ~ wh 2 X X X 2 ~h rPh ; rqh Þ ¼ d ~h ; qh Þ þ d h ðD U h2T ð~ f ; rqh ÞT ðr U T X T T 2Th
8 qh 2 Mh :
T 2Th
As in our formulation, here the stabilization is achieved through a consistent term (depending on the residual of the momentum equation) added to the continuity equation. On the other hand, our ~ acting on the continuous level is here only mimicked on the discrete level by the operator S operator r on each triangle. Indeed, thinking multilevel, the weights h2T simulate the effect of ~0 , which transforming the elementwise L2 -inner products into a global inner product in the space X ~h rPh þ ~ f and the term rqh . is the correct function space for both the residual DU Going back to our variational problem, we observe that it can be written as e ½ðU ~; P Þ 2 X ~ M: A ~ ; P Þ; ð~ ðU w; qÞ ¼ Fe ð~ w; qÞ
~ M; 8 ð~ w; qÞ 2 X
after defining e ½ð~ ~ qi; A v; rÞ; ð~ w; qÞ :¼ A½ð~ v; rÞ; ð~ w; qÞ þ dhD~ v rr; S ~ qi: Fe ð~ w; qÞ :¼ F ð~ w; qÞ dh~ f;S
ð3:14Þ
428
C. Canuto / Computers & Fluids 31 (2002) 421–436
The proposed stabilization can be cast into the abstract setting of the previous section, in the ~0 M. Denote by J the canonical isomorphism between X ~0 following way. Choose Z 0 :¼ W 0 ¼ X ~ and X , i.e., the inverse of the minus Laplacian operator supplemented by Dirichlet boundary conditions. Then, we can choose the inner product f ; J~ gi þ ðr; qÞX ððð~ f ; rÞ; ð~ g; qÞÞÞX~0 M :¼ h~ ~0 M. The space Ve coincides with V ¼ X ~ M. Furthermore, let us choose W e :¼ W ¼ X ~
in X 0 ~ M!X ~ M as M and let us define S : X ~ q; 0Þ: Sð~ w; qÞ :¼ ðJ1 S ~ M, Then, recalling Eq. (3.4), it is easily seen that, for all ð~ v; rÞ; ð~ w; qÞ 2 X ~ qi: ððLð~ v; rÞ; Sð~ w; qÞÞÞX~0 M ¼ hD~ v rr; S e and f~ :¼ Fe . Assumption (2.1) is fulfilled Thus, Eq. (3.14) coincides with Eq. (2.6) if we set a~ :¼ A provided d is not too large. Indeed, taking into account Eqs. (3.10) and (3.11), the following result can be established (see Ref. [14] for the details of the proof). Theorem 3.2. Let us assume that d satisfy the bound d 6 b =c2 . Then, there exists a constant b > 0 (independent of d) such that inf
sup
ð~ v;rÞ2V ð~ w;qÞ2V
e ½ð~ A v; rÞ; ð~ w; qÞ P b; kð~ v; rÞkX~ M kð~ w; qÞkX~ M
~ M is ~ M or any finite dimensional subspace X ~h Mh ; and the norm in X where V equals either X 2 2 2 vkX~ þ dkrkM . scaled as kð~ v; rÞkX~ M ¼ k~ Using this result, it is then straightforward to obtain optimal a priori and a posteriori error ~ ; P Þ of the Stokes problem ((3.1)–(3.3)) and the solution ð~ estimates between the solution ðU u; pÞ of the Galerkin discretization: e ½ð~ ~h Mh : ~h Mh : A u; pÞ; ð~ w; qÞ ¼ Fe ð~ w; qÞ 8 ð~ w; qÞ 2 X ð~ u; pÞ 2 X Numerical results, obtained by this formulation using the auxiliary stabilizing velocities described below, show good agreement with the theoretical predictions. They can be found in Ref. [14]. 3.1. An example We now indicate how the auxiliary velocities satisfying Eqs. (3.7) and (3.8) can be constructed by means of appropriate wavelet bases. Starting from compactly supported (bi)orthogonal scaling function and wavelet systems on the real line, such as the Daubechies family or any family of B-splines, it is possible to construct biorthogonal scaling function and wavelet systems on the reference interval I ¼ ð0; 1Þ. For the technical details, we refer, e.g., to Ref. [18]. For each level index j P some j0 , there exist families e j ¼ fn~jk ; k 2 Kj g of mutually biorthogonal scaling functions (i.e., Nj ¼ fnjk ; k 2 Kj g and N
C. Canuto / Computers & Fluids 31 (2002) 421–436
R
429
S
I njk n~jk0 ¼ dk;k 0 for all k; k 0 2 Kj ) such that, setting VjI ¼ span Nj , one has VjI Vjþ1 ; 8j, and j VjI e j . Furthermore, there exist is dense in L2 ðIÞ; similar properties hold for the sequence VejI ¼ span N e gjh ; h 2 Hj g of mutually biorthogonal wavelets such families j ¼ fgjh ; h 2 Hj g and j ¼ f~ I I e j , one has V I ¼ V I W I and Ve I ¼ Ve I W e e I , with that, setting Wj ¼ span !j and W j ¼ span ! jþ1 j j jþ1 j j I I I I e e Wj ? V j and W j ? Vj . The construction can be carried on in a way that at least the zeroth-order moment (i.e., the average over I) of each wavelet vanishes; in addition, the dual scaling functions and wavelets belong to H01 ðIÞ. ^ ¼ ð0; 1Þn . Biorthogonal With such a construction at hand, let us first consider the domain X ^ Þ can be defined by taking tensor products of univariate scaling functions and Riesz bases in L2 ðX ^ wavelets. Precisely, for any j P j0 , let WXj denote the set of wavelets w ¼ ni¼1 #i , such that each #i is ~ X^ be either a scaling function n 2 Nj or a wavelet g 2 j , and at least one #i is a wavelet. Let W j ^ ^ Þ, and let UX^ e X be a tensor product basis of ð n Ve I Þ \ L2 ðX defined similarly. Furthermore, let U i¼1 j0 0 j0 j0 be the companion biorthogonal set in ni¼1 VjI0 . Finally, set [ ^ [ ^ ^ ^ e X: e p;X^ ¼ U e X^ [ W WXj ; W Wp;X ¼ UXj0 [ j j0 I
j P j0
j P j0
^ ^ Þ ¼ span W e p;X^ ; more precisely, it is Then, we define M :¼ span Wp;X . One also checks that L20 ðX ^ p p e p;X has the form w~ ¼ n #~i , where #~i ð0Þ ¼ #~i ð1Þ ¼ 0 for all i and there easily seen that any w~k 2 W k R i¼1 exists at least one index, say m ¼ mðkÞ such that I #~m ¼ 0. Then, we associate to w~pk the auxiliary ^ velocity ~ w ;k X defined by ( ^ ^ X ;X ~ where w ; k ¼ ð0; . . . ; wk ; . . . ; 0Þ R x^m ^ ;X ~ xÞ ¼ ð 0 #m ðsÞ dsÞ i6¼m #~i ð^ xi Þ: wk ð^
Condition (3.7) is trivially satisfied, whereas Eq. (3.8) can be checked by a version of the so-called ‘vaguelettes’ Lemma [32]. ^ under a smooth mapping F. Set G ¼ F 1 and let DG denote the Next, let X be the image of X Jacobian matrix of G. Then, we define ^
p;X Wp;X ¼ fwp;X k :¼ wk # Gg
and
e p;X ¼ fw~p;X :¼ jDGjw~p;X^ # Gg; W k k
^ X ~ ;X defined by and we associate to each w~p;X the Piola transform of ~ w ; k , namely, the velocity wk k ^ ;X ; X ~ wk # GÞ. Again, Eqs. (3.7) and (3.8) are satisfied. wk :¼ jDGjDGð~ For more general domains, one can resort to domain decomposition. In this case, pressure are discontinuous between subdomains; this implies that the stabilizing wavelets can be built independently on each subdomain as shown above, and then extended by zero outside the subdomain of definition. This completes the construction of the stabilizing wavelets.
4. The convection–diffusion problem Let X be again a bounded domain contained in Rn (n P 1), with Lipschitz boundary oX. We consider the problem
430
C. Canuto / Computers & Fluids 31 (2002) 421–436
mDu þ a ru þ bu ¼ f ; u ¼ 0;
in X; on oX;
ð4:1Þ
n
where m > 0 is a constant, a 2 ðW 1;1 ðXÞÞ , b 2 L1 ðXÞ and there exists a constant c > 0 (of order of magnitude 1) such that 12r a þ b P c
ð4:2Þ
a:e: in X:
Problem (4.1) is written variationally as in Eq. (2.2), if we assume that f 2 H 1 ðXÞ and we set V ¼ W :¼ H01 ðXÞ, Z Z Z ð4:3Þ aðv; wÞ :¼ m rv rw dx þ a rvw dx þ bvw dx; X
X
X
and hf ; wi :¼ H 1 ðXÞ hf ; wiH 1 ðXÞ : 0
The form a is continuous and coercive on V; precisely, one has for all v 2 V mkrvk2L2 ðXÞ þ ckvk2L2 ðXÞ 6 aðv; vÞ 6 Kðkrvk2L2 ðXÞ þ kvk2L2 ðXÞ Þ
ð4:4Þ
with K proportional to maxðm; kakðL1 ðXÞÞn ; kbkL1 ðXÞ Þ. If m kakðL1 ðXÞÞn , we are precisely in one of the situations considered in Section 2, namely the inf–sup constant (here, the coercivity constant) and the continuity constant are sensibly different from each other. As a consequence, the solution of Problem (4.1) is allowed to exhibit strong gradients, which usually take the form of highly localized transition layers (boundary and interior layers). If a standard Galerkin method is used to approximate that problem, and if functions in Vh do not have the potential to faithfully represent such layers, the corresponding discrete solution exhibits strong gradients as well, but they manifest themselves as oscillations, more or less localized around the layers. These oscillations can be prevented by a stabilization procedure. The key idea underlying the stabilization method proposed in Ref. [11] and developed in Refs. [12,13], is to obtain a better continuity estimate than in Eq. (4.4). Just for the sake of simplicity, let us assume that r a ¼ 0 in X; furthermore, let us set Dv :¼ a rv; and let, as before, ð; ÞX denote the L2 ðXÞ-inner product. The convective term ðDv; wÞX can be bounded as jðDv; wÞX j 6 kDvkL2 ðXÞ kwkL2 ðXÞ
8 v; w 2 H01 ðXÞ;
on the other hand, by partial integration, it can also be bounded as jðDv; wÞX j ¼ j ðv; DwÞX j 6 kvkL2 ðXÞ kDwkL2 ðXÞ : The two bounds and the theory of interpolation of function spaces suggest to look for a more symmetric bound, involving the same norm of both v and w. Such bound should control Dv and Dw in the norm of an intermediate space between L2 ðXÞ and H 1 ðXÞ (remember that if v 2 L2 ðXÞ, then Dv 2 H 1 ðXÞ and kDvkH 1 ðXÞ K kvkL2 ðXÞ ). Taking into account the boundary conditions, a 1=2 natural candidate as the intermediate space is the dual of the space H00 ðXÞ, which is the in-
C. Canuto / Computers & Fluids 31 (2002) 421–436
431
terpolation space of order 1=2 between H01 ðXÞ and L2 ðXÞ (see, e.g. Ref. [31]). For simplicity, let us 1=2 1=2 0 set H1=2 ðXÞ :¼ H00 ðXÞ and H1=2 ðXÞ :¼ ðH00 ðXÞÞ . Furthermore, let ðð; ÞÞH1=2 ðXÞ be an inner product in H1=2 ðXÞ, whose definition will be discussed later on. In the sequel, we shall first present the conceptually simplest way to build-up a stabilized variational formulation of Eq. (2.2), which guarantees an H1=2 ðXÞ-control on the streamline derivative. Next, we shall introduce a more refined and flexible form of the stabilization device. Let us start by considering the space X 1=2 ðXÞ :¼ fv 2 H01 ðXÞ : Lv 2 H1=2 ðXÞg
ð4:5Þ
equipped with the norm kvkX 1=2 ðXÞ :¼ kvkL2 ðXÞ þ m1=2 krvkL2 ðXÞ þ d1=2 kLvkH1=2 ðXÞ ;
ð4:6Þ
where d > 0 is a suitable parameter. Note that this norm controls Lv, rather than Dv. Indeed, we aim at introducing a consistent stabilization of Eq. (2.2), based on the residual f Lv. However, we shall recover the control on kDvkH1=2 ðXÞ from the control on kLvkH1=2 ðXÞ , by some sort of ‘post-processing’ analysis. From now on, we assume that f 2 H1=2 ðXÞ, so that u 2 X 1=2 ðXÞ. At this point, we introduce the stabilized bilinear form A : X 1=2 ðXÞ X 1=2 ðXÞ ! R;
Aðv; wÞ :¼ aðv; wÞ þ dððLv; LwÞÞH1=2 ðXÞ
ð4:7Þ
and the stabilized linear form F : X 1=2 ðXÞ ! R;
F ðwÞ :¼ hf ; wi þ dððf ; LwÞÞH1=2 ðXÞ ;
ð4:8Þ
and we consider the variational problem: find u 2 X 1=2 ðXÞ such that Aðu; wÞ ¼ F ðwÞ
8 w 2 X 1=2 ðXÞ:
ð4:9Þ
e ¼ Thus, with respect to the abstract framework of Section 2, we have Z 0 ¼ H1=2 ðXÞ, Ve ¼ W 1=2 X ðXÞ and S ¼ L. With this setting, we obtain a~ ¼ A and f~ ¼ F . Coming to Assumption 2.1, we note that the form A is uniformly coercive on X 1=2 ðXÞ. Indeed, Aðv; vÞ P mkrvk2L2 ðXÞ þ ckvk2L2 ðXÞ þ dkLvk2H1=2 ðXÞ P akvk2X 1=2 ðXÞ
8 v 2 X 1=2 ðXÞ;
ð4:10Þ
for a constant a > 0 independent of m. On the other hand, as far as continuity is concerned, we have the symmetric bound jAðv; wÞj 6 CkvkX 1=2 ðXÞ kwkX 1=2 ðXÞ
8 v; w 2 X 1=2 ðXÞ;
but the constant C may blow-up with m1 (see Ref. [12, estimate (3.12)]). However, introducing a slightly stronger norm of v yields the uniform (in m) estimate (see Ref. [12, Proposition 4.1]) jAðv; wÞj K ðkvkX 1=2 ðXÞ þ d1=2 kvkH1=2 ðXÞ ÞkwkX 1=2 ðXÞ
8 v; w 2 X 1=2 ðXÞ:
ð4:11Þ
If a Galerkin approximation of Eq. (4.9) is considered, the bounds (4.10) and (4.11) imply a uniform error estimate in the X 1=2 ðXÞ-norm. We shall not detail such analysis here, since a similar result will be obtained later on, for the modified stabilization setting. Let us discuss the definition of the inner product in H1=2 ðXÞ, which so far was not made precise. We propose to use a multilevel decomposition of L2 ðXÞ, generated by two mutually
432
C. Canuto / Computers & Fluids 31 (2002) 421–436
e :¼ fw~k : k 2 Kg. Each k is a couple, k :¼ ðj; kÞ, biorthogonal Riesz bases W :¼ fwk : k 2 Kg and W where j P j0 (for some j0 2 N) is the ‘level’ (or ‘scale’) index and k 2 Kj (for some finite set Kj Z) is the ‘position’ index, associated with the basis function wk . For instance, wk might have (at least locally) the form wk ðxÞ ¼ 2nj=2 wð2j x kÞ, where w is an appropriate ‘generating function’. For notational convenience, we set jkj :¼ j. The elements of both bases are compactly S k similarly, we have supported functions; denoting by Sk the support of wk and defining e S k 2jkj . We recall that biorthogonality means that the relations diam Sk diam e ðwk ; w~k0 ÞX ¼ dk;k0
8 k; k0 2 K
ð4:12Þ
are fulfilled, whereas the Riesz basis condition says that any v 2 L2 ðXÞ can be written as X v¼ ðv; w~k ÞX wk k2K
with kvk2L2 ðXÞ
X
jðv; w~k ÞX j2 :
k2K
In addition, let us suppose that W H01 ðXÞ and that any g 2 H 1 ðXÞ can be written as X g¼ hg; wk iw~k : k2K
Let us introduce the scale of Sobolev spaces Hr ðXÞ, 1 6 r 6 1, defined as 8 r 1=2 < r 6 1; H ðXÞ; > > < 01=2 r ¼ 1=2; H00 ðXÞ; Hr ðXÞ :¼ r > ðXÞ; 0 6 r < 1=2; H > : 0 r ðH ðXÞÞ ; 1 6 r < 0: e allow the characterization of the spaces Hr ðXÞ. Precisely, for We assume that the bases W and W any 0 6 r 6 1, we require that X 22rjkj jðv; w~k ÞX j2 8 v 2 Hr ðXÞ; ð4:13Þ kvk2Hr ðXÞ k2K
and kgk2Hr ðXÞ
X
22rjkj jhg; wk ij2
8 g 2 Hr ðXÞ:
ð4:14Þ
k2K
Families of biorthogonal basis functions in non-trivial domains have been constructed in the last few years; see, e.g. Refs. [15,16,20,22,24]. Inspired by Eq. (4.14), we are led to use the following H1=2 ðXÞ-inner product X ððf ; gÞÞH1=2 ðXÞ :¼ 2jkj hf ; wk ihg; wk i k2K 1
in the stabilized forms (4.7) and (4.8). The parameter d should be taken proportional to kakL1 ðXÞ therein. Obviously, in general infinitely many terms appear on the right-hand side, making the
C. Canuto / Computers & Fluids 31 (2002) 421–436
433
practical implementation of such inner product unfeasible. However, once a finite dimensional trial/test space Vh H01 ðXÞ has been chosen, it is possible ([12, Section 5]) to truncate the inner product to a finite sum (depending on the ‘resolution power’ of Vh ), without affecting the stabilization effect. To be more precise, if Vh contains all basis functions wk with jkj smaller than some jðhÞ, then the sum in Eq. (4.18) can be restricted to those k satisfying jkj P jðhÞ. On the other hand, if we assume that the following (absolutely reasonable) inverse inequality kDvh kL2 ðXÞ K h1=2 kDvh kH1=2 ðXÞ
8 vh 2 Vh ;
holds true, then there exists a level J ðhÞ j log2 hj such that all the details with jkj > J ðhÞ are inessential in the stabilization process. In conclusion, the H1=2 ðXÞ-inner product defined above can be replaced by the truncated inner product X 2jkj hf ; wk ihg; wk i: ððf ; gÞÞh :¼ jðhÞ 6 jkj 6 J ðhÞ
Let us now turn to the more refined version of our multilevel stabilization approach. The key observation is that––thinking multilevel––the stabilization mechanism should provide extra H1=2 ðXÞ-control only on the ‘low-level’ components of Dv, since the ‘high-level’ components are implicitly controlled by the diffusive term mDv. In other words, one should insert information on the local relationship between convection and diffusion. In classical SUPG-stabilizations of finite element methods (see, e.g., Ref. [29]), such information is inserted at the discrete level through a switch in the stabilization parameter, based on a local (in position) Peclet number. In our approach, we are going to act at the continuos level (i.e., before any finite-dimensional discretization) and on the basis of a local (in position but also in ‘scale’) Peclet number. To be precise, for any index k 2 K, let us introduce the multiscale Peclet number
Pek :¼ c
kakL1 ðSk Þ 2jkj m
;
ð4:15Þ
where c is a tuning parameter depending on the multilevel basis W. Next, set ! ( ) 2jkj if Pek > 1 2jkj 22jkj c kakL1 ðS Þ k ; ; ¼ min ; sk :¼ 22jkj c kakL1 ðSk Þ m if Pe 6 1 k m and define the inner product X ððf ; gÞÞs :¼ sk hf ; wk ihg; wk i;
ð4:16Þ
k2K
as well as the associated norm kgks :¼ ððg; gÞÞ1=2 s . Recalling Eq. (4.14), such bilinear form behaves 1=2 ðXÞ-inner product on the ‘low-level’ components of its arguments, and as an like a (scaled) H H 1 ðXÞ-inner product (scaled by m1 ) on the ‘high-level’ components. In particular, Eq. (4.16) is e ¼ X1=2 ðXÞ, where defined on H 1 ðXÞ. Therefore, now we can choose Z 0 ¼ H 1 ðXÞ, so that Ve ¼ W 1=2 1 X ðXÞ :¼ H0 ðXÞ equipped with the norm kvkX1=2 ðXÞ :¼ kvkL2 ðXÞ þ m1=2 krvkL2 ðXÞ þ kLvks : The continuous problem we aim at discretizing reads as follows: find u 2 H01 ðXÞ such that
434
C. Canuto / Computers & Fluids 31 (2002) 421–436
8 w 2 H01 ðXÞ;
Aðu; wÞ ¼ FðwÞ with
A : H01 ðXÞ H01 ðXÞ ! R;
Aðv; wÞ :¼ aðv; wÞ þ ððLv; LwÞÞs
ð4:17Þ
and F : H01 ðXÞ ! R;
FðwÞ :¼ ðf ; wÞ þ ððf ; LwÞÞs :
Note the explicit expression of the stabilizing term appearing in Eq. (4.17): ððLv; LwÞÞs ¼
X
X 22jkj 2jkj aðv; wk Þaðw; wk Þ: aðv; w Þaðw; w Þ þ k k c kakL1 ðSk Þ m Pe 6 1
Pek >1
ð4:18Þ
k
For the new form, one can prove the coerciveness inequality akvk2X1=2 ðXÞ 6 Aðv; vÞ 8 v 2 H01 ðXÞ; for a constant a 1, and the uniform continuity bound jAðv; wÞj K ðkvkX1=2 ðXÞ þ kak1=2 L1 ðXÞ kvkH1=2 ðXÞ ÞkwkX1=2 ðXÞ
8 v; w 2 H01 ðXÞ;
which are similar to the analogous ones (4.10) and (4.11). Next, given any family fVh gh>0 of finite dimensional subspaces of H01 ðXÞ, we consider the stabilized Galerkin approximation: find uh 2 Vh such that Aðuh ; wh Þ ¼ Fðwh Þ
8 wh 2 Vh :
ð4:19Þ
Standard arguments lead to the error bound 1=2
ku uh kX1=2 ðXÞ K inf ðku vh kX1=2 ðXÞ þ kakL1 ðXÞ ku vh kH1=2 ðXÞ Þ: vh 2Vh
ð4:20Þ
Let us assume that the family fVh gh>0 provides classical approximation properties, expressed by the bounds inf kv vh kHr ðXÞ K hsr jvjH s ðXÞ
vh 2Vh
8 v 2 H s ðXÞ \ H01 ðXÞ;
for 0 6 r 6 1 and 1 6 s < some s0 . Then, one can establish [13, Theorem 3.4] the error estimate ku uh kL2 ðXÞ þ m1=2 kru ruh kL2 ðXÞ þ kLu Luh ks K ½Chs þ m1=2 hs1 þ Bhsð1=2Þ jujHs ðXÞ ; ð4:21Þ where B and C are constants depending on the velocity field a (they are surely K 1 if a is locally Lipschitz continuous at the points where it vanishes). As anticipated, from this result we deduce a bound on the H1=2 ðXÞ-norm of Du Duh . Indeed, expressing the negative order norms in the equivalent multilevel form (4.14), one obtains 1=2 1=2 kru ruh kL2 ðXÞ kak1=2 L1 ðXÞ kDu Duh kH1=2 ðXÞ K Cku uh kL2 ðXÞ þ ðkakL1 ðXÞ þ BÞm
þ kLu Luh ks :
C. Canuto / Computers & Fluids 31 (2002) 421–436
435
This and Eq. (4.21) imply that the error on the streamline derivative measured in the H1=2 ðXÞnorm can be bounded in terms of the right-hand side of Eq. (4.21). The stabilized scheme (4.19) is still not implementable, as in the stabilizing term (4.18) infinitely many details appear. However, the sum in Eq. (4.18) can be truncated to a finite number of details, depending on Vh . The truncation can be accomplished precisely in the same way as explained above for the ‘global version’ of our stabilization. 5. Conclusions In our opinion, one of the interests of the proposed multilevel framework is that it sheds new light on classical stabilization strategies such as Galerkin/least squares or SUPG. Usually, these methods add a locally finite number of terms to the pure Galerkin approximation. Such perturbation terms should be interpreted as, or at least related to, a particular truncated version of a Sobolev-type inner product, as they mimic the inner product at the discrete level. For the Stokes problem, the relationship has been discussed in Section 3. For model 1D convection–diffusion problems, the interpretation has been rigorously established in Ref. [12, Section 6]: the classical SUPG stabilization considered in Ref. [29] is shown to coincide with our multilevel-stabilized method (4.19), in which the (truncated) inner product is defined via piecewise linear interpolatory wavelets. The theoretical results briefly described above apply to this situation. The functional framework we have developed indicates that the ‘order 1=2’ regularity setting is appropriate for the correct understanding of numerical discretizations of convection–diffusion problems, in the convection-dominated regime. Other recent and independent results [10,30], which involve order 1=2 norms in the study of stabilized approximations to convection–diffusion problems, confirm our perspective. The extension of our investigation to more complex equations of fluid dynamics, such as the incompressible Navier–Stokes equations, is underway. References [1] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function. II. Nonlinear reconstruction techniques. SIAM J Sci Comput 1999;20:1053–93. [2] Babuska I, Aziz AK. Survey lectures on the mathematical foundation of the finite element method. In: Aziz AK, editor. The mathematical foundation of the finite element method with applications to partial differential equations. New York: Academic Press; 1972. p. 3–359. [3] Baiocchi C, Brezzi F. Stabilization of unstable numerical methods, Atti del Simposio Internazionale Problemi attuali dell’Analisi e della Fisica Matematica. Taormina, 15–17 ottobre 1992 (1993), pp. 59–64. [4] Balder R, Zenger C. The solution of multidimensional real Helmholtz equations on sparse grids. SIAM J Sci Comput 1996;17:631–46. [5] Bank RE. In: Hierarchical bases and the finite element method, Acta Numerica. Cambridge: Cambridge University Press; 1996. p. 1–43. [6] Bertoluzza S. Stabilization by multiscale decomposition. Appl Math Lett 1998;11:129–34. [7] Bramble JH, Lazarov RD, Pasciak JE. A least-squares approach based on a discrete minus one inner product for first order systems. Math Comput 1997;66:935–55. [8] Bramble JH, Pasciak JE, Xu J. Parallel multilevel preconditioners. Math Comput 1990;55:1–22. [9] Brezzi F, Fortin M. Mixed and hybrid finite element methods. New York: Springer; 1991. [10] Brezzi F, Marini D, S€ uli E. Residual-free bubbles for advection–diffusion problems: the general error analysis. Numer Math 2000;85:31–47.
436
C. Canuto / Computers & Fluids 31 (2002) 421–436
[11] Bertoluzza S, Canuto C, Tabacco A. Negative norm stabilization of convection–diffusion problems. Appl Math Lett 2000;13:121–7. [12] Bertoluzza S, Canuto C, Tabacco A. Stable discretizations of convection–diffusion problems via computable negative-order inner products. SIAM J Numer Anal 2000;38:1034–55. [13] Canuto C. Multilevel stabilization of convection–diffusion problems by variable-order inner products. Computing 2001;66:121–38. [14] Canuto C, Masson R. Stabilized wavelet approximations of the Stokes problem. Math Comput 2001;70:1397–416. [15] Canuto C, Tabacco A, Urban K. The wavelet element method. Part I: construction and analysis. Appl Comput Harm Anal 1999;6:1–52. [16] Canuto C, Tabacco A, Urban K. The wavelet element method. Part II: realization and additional features in 2D and 3D. Appl Comput Harm Anal 2000;8:123–65. [17] Chen M, Temam R. Incremental unknowns for solving partial differential equations. Numer Math 1991;59:255–71. [18] Cohen A. Wavelet methods in numerical analysis, In: Ciarlet PG, Lions JL, editors. Handbook of numerical analysis, vol. VII. Amsterdam: Elsevier; 2000. [19] Cohen A, Dahmen W, DeVore R. Adaptive wavelet methods for elliptic operator equations––convergence rates. Math Comp 2001;70:27–75. [20] Cohen A, Masson R. Wavelet adaptive methods for 2nd order elliptic problems. Domain decomposition and boundary conditions. Numer Math 2000;86:193–238. [21] Dahmen W. Stability of multiscale transformations. J Fourier Anal Appl 1996;4:341–62. [22] Dahmen W. In: Wavelet and multiscale methods for operator equations Acta Numerica 6. Cambridge: Cambridge University Press; 1997. p. 55–228. [23] Dahmen W, Kunoth A. Multilevel preconditioning. Numer Math 1992;63:315–44. [24] Dahmen W, Schneider R. Composite wavelet bases for operator equations. Math Comput 1999;68:1533–67. [25] DeVore R. In: Nonlinear approximation Acta Numerica. Cambridge: Cambridge University Press; 1998. p. 53– 150. [26] Farge M, Schneider K, Kevlahan N. Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis. Phys Fluids 1999;11:2187–201. [27] Harten A. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Comm Pure Appl Math 1995;48:1305–42. [28] Hughes TJR, Franca LP, Balestra M. A new finite element formulation of computational fluid dynamics: a stable Petrov–Galerkin formulation of the Stokes problem accomodating equal-order interpolations. Comput Methods Appl Mech Engng 1986;59:85–99. [29] Hughes TJR, Franca LP, Hulbert GM. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput Meth Appl Mech Engng 1989;73: 173–89. [30] Kellogg RB, Stynes M. The intrinsic norm finite element method for convection–diffusion problems, Preprint 1999-1, Department of Mathematics, University College, Cork, Ireland. [31] Lions JL, Magenes E. Problemes aux limites non homogenes. Paris: Dunod; 1968. [32] Meyer Y. Ondelettes et operateurs––tomes 1 et 2. Paris: Hermann; 1990. [33] Oswald P. Multilevel finite element approximation. Theory and applications. Stuttgart: Teubner; 1994. [34] Roos H-G, Stynes M, Tobiska L. Numerical methods for singularly perturbed differential equations. Convection– diffusion and flow problems. Berlin: Springer; 1996.
Multilevel stabilization devices in CFD q Claudio Canuto Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Turin, Italy Received 1 August 2001 Dedicated to Professor Roger Peyret
Abstract The numerical discretization of such fundamental equations in fluid dynamics as the Navier–Stokes equations, requires special attention: a careful treatment of the pressure and the convective terms is needed, to avoid the onset of spurious oscillations. We advocate the use of multilevel bases (wavelets, or hierarchical finite elements) to build-up consistent (i.e., residual based) stabilization devices. The proposed strategy is particularly well-suited for being applied within an adaptive discretization procedure. Classical SUPG-type methods can be interpreted as special realizations of the multilevel stabilization idea. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Coerciveness; Inf–sup condition; Stabilization; Multilevel bases; Stokes problem; Convection–diffusion problems
1. Introduction Multilevel bases include hierarchical finite elements [5], wavelets and other wavelet-like families [18], as well as the bases associated with sparse grid methods [4] or incremental unknown methods [17]. Although such bases are less simple and easy-to-handle than standard bases such as Lagrange’s, yet they enjoy extra properties, which make their use in numerical simulation very attractive. Significant applications of multilevel finite element and finite difference techniques have
q
This work was partially supported by the European Commission within the TMR project (Training and Mobility for Researchers) Wavelets and Multiscale Methods in Numerical Analysis and Simulation, no. ERB FMRX CT98 0184, and by the Italian funds Murst Cofin’98 Analisi Numerica. E-mail address: [email protected] (C. Canuto). 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 3 - 9
422
C. Canuto / Computers & Fluids 31 (2002) 421–436
been given over the years; the use of wavelet-like methods is currently under investigation, with various promising results [22]. Multilevel bases provide a natural capability of preconditioning symmetric positive-definite problems, as they yield stiffness matrices which are spectrally equivalent to diagonal matrices (see, e.g., Refs. [5,8,23]). The possibility of approximating a function (a signal) in an optimal way by adding and removing details [25] can be exploited in designing efficient reconstruction algorithms, with applications to, e.g., the numerical discretization of conservation laws [1,27]. Details (which form the multiscale representation of an approximate solution, or a residual) are inherently error indicators; sound adaptive strategies can be based on their inspection [19]. Related to this, is the recently proposed and stimulating attempt to separate coherent and incoherent structures in the numerical simulation of a turbulent flow, by means of wavelets [26]. The present contribution focusses on another feature of multilevel bases, namely the possibility of representing inner products and norms of function spaces of any regularity order in an explicit and computable form, i.e., by weighted sums of subsequent details [21,32,33]. This property, which somehow generalized the Parseval identity in Fourier analysis, is very important; indeed, often norms and inner products of fractional and/or negative order are defined in a way that is not at all amenable to an efficient numerical treatment. By means of multilevel bases, these mathematical objects become implementable. Many boundary value problems of interest lack of the good symmetric and positive-definite (coercive) character of classical elliptic problems. This is a major cause of difficulty in their numerical treatment. In fluid dynamics, spurious velocity oscillations may appear in the numerical simulation of flows at high Reynolds numbers; they are allowed by the insufficient control provided by the diffusion terms. Furthermore, pressure instabilities may rise in the simulation of incompressible flows, whenever the discrete trial velocity and pressure spaces fail to satisfy a particular form of coerciveness condition, the well-known LBB or inf–sup condition. Over the years, many different approaches have been proposed to transform a boundary value problem into an equivalent one with enhanced coerciveness; SUPG-type stabilization methods, least-square and Galerkin/least-square methods are among the examples. The general idea of consistent stabilization amounts to augmenting the standard variational formulation of the problem by another variational term, which is in some sense ‘proportional’ to the residual of the equation. The effect is to provide a tighter control on the solution (for instance, in a stronger norm). Classical stabilization is accomplished at the discrete level; however, the problem can also be set in an abstract way on the exact operator equation [3,6]. It is here that multilevel inner products and norms come usefully into play (see, e.g., Refs. [7,11]). In the first part of this paper, we present an abstract framework in which consistent stabilizations of weakly coercive or non-coercive problems can be cast. The setting is particularly geared towards multilevel-based stabilizations. Next, we illustrate two applications to the numerical discretization of model problems of interest in fluid dynamics: the Stokes problem and the convection–diffusion problem. In the first case, the stabilization proposed in Ref. [14] allows a fully independent discretization of velocity and pressure, a feature which is fundamental if the numerical method incorporates an adaptive strategy. Uniform control on the exact L2 -norm of the discrete pressure (up to the physically indetermined constant mode) is guaranteed. For the convection–diffusion problem, following [11–13], we advocate the use of a fractional order norm (precisely, of order 1=2) to get a uniform control on the streamline derivative of the discrete
C. Canuto / Computers & Fluids 31 (2002) 421–436
423
solution. In addition, our functional framework suggests to view classical SUPG methods as truncated versions of the proposed multilevel stabilization. Notation: Throughout the paper, we shall use the notation A K B to indicate A 6 cB for a suitable constant c, independent of the relevant parameters for the particular formula in which the symbol is used (e.g., the diffusion parameter m, the multiscale parameters k, the discretization parameter h). The notation A B means A K B and B K A.
2. An abstract setting Let us suppose that we aim at numerically solving an operator equation of the form Lu ¼ f ;
ð2:1Þ
where L is linear continuous between a Hilbert space V and the dual W 0 of a (possibly different) Hilbert space W. Let h; i denote the duality pairing between W 0 and W. Introducing the continuous bilinear form a : V W ! R;
aðv; wÞ :¼ hLv; wi;
the operator equation can be written variationally as u 2 V : aðu; wÞ ¼ hf ; wi 8 w 2 W :
ð2:2Þ
It is classically known [2] that the well-posedness of this problem is equivalent to the two inf–sup conditions inf sup v2V w2W
aðv; wÞ > 0; kvkV kwkW
inf sup
w2W v2V
aðv; wÞ > 0: kvkV kwkW
ð2:3Þ
Let us introduce a finite dimensional approximation of this problem. To this end, for h 2 ð0; hmax , let Vh V and Wh W be subspaces satisfying dim Vh ¼ dim Wh < þ1, and let us consider the Petrov–Galerkin discretization of Eq. (2.2) uh 2 Vh : aðuh ; wh Þ ¼ hf ; wh i
8 w h 2 Wh :
ð2:4Þ
The solvability of this problem, as well as the behavior of the solution, is related to the (best) continuity constant j > 0 and inf–sup constant ah P 0, which are defined through the inequalities ah kvh kV 6 sup
wh 2Wh
aðvh ; wh Þ 6 jkvh kV kwh kW
8 vh 2 Vh :
ð2:5Þ
Note that j can be chosen independently of h, as a consequence of the continuity of the form a on V W ; in addition, it is not restrictive to assume that j 1, possibly after a proper scaling of the equation. Conversely, the inf–sup conditions (2.3) for the exact problem may imply nothing on ah : neither its independence of h, nor even its being different from 0. If ah ¼ 0, problem (2.4) is not well-posed. This is, e.g., the case of the approximations of the Stokes problem, in which the discrete velocity and pressure spaces do not satisfy a suitable compatibility condition; as a consequence, the discrete pressure solution is defined up to spurious
424
C. Canuto / Computers & Fluids 31 (2002) 421–436
oscillating components (see, e.g., Ref. [9]). On the other hand, if 0 < ah j, problem (2.4) is indeed well-posed, but the classical error estimate j inf ku vh kV ku uh kV 6 1 þ ah vh 2Vh allows the discrete solution to be far apart from the best approximation in Vh of the exact solution. This is the case, e.g., of pure Galerkin approximations to convection–diffusion problems in the convection-dominated regime (see, e.g., Ref. [34]). The previous considerations motivate the development over the years of a great number of strategies, which modify the variational problem (2.4) in a consistent way (i.e., in such a way that the solution of Eq. (2.2) is also the solution of the modified problem). The aim is to improve the constant a~h which appears in the analog of Eq. (2.5) for the modified variational form a~; more precisely, a~h should possibly be bounded away from 0 independently of h, and have the same order of magnitude as the continuity constant j~. We propose an abstract functional framework which leads to stabilization strategies based on the multilevel representation of a Hilbert space. Let us denote by Z 0 a Hilbert space continuously imbedded into W 0 , with inner product ðð; ÞÞZ 0 . Let us set Ve :¼ fv 2 V : Lv 2 Z 0 g equipped with the graph norm kvkV~ :¼ kvkV þ d1=2 kLvkZ 0 : Here, d > 0 is a parameter, whose actual choice will depend on the particular application but whose order of magnitude should be thought to be in the order of unity. From now on, let us suppose that the data f, assigned in Eq. (2.1) as an element of W 0 , be indeed in Z 0 . This immediately implies that the solution u belongs to Ve . e W with continuous imbedding, and a linear Next, let us introduce a second Hilbert space W 0 e operator S 2 Lð W ; Z Þ. By such operator, we define the modified bilinear form e ! R; a~ : Ve W
a~ðv; wÞ :¼ aðv; wÞ þ dððLv; SwÞÞZ 0 ;
as well as the modified linear form e ! R; f~ : W
f~ðwÞ :¼ hf ; wi þ dððf ; SwÞÞZ 0 :
This modification is consistent with Eq. (2.1), in the sense that u also solves the variational problem u 2 Ve : a~ðu; wÞ ¼ f~ðwÞ
e: 8w2W
ð2:6Þ
We now assume that the finite dimensional subspaces Vh V and Wh W previously introduced, e , respectively. Thus, we consider the Petrov–Galerkin approxiare also included into Ve and W mation of problem Eq. (2.1) uh ; wh Þ ¼ f~ðwh Þ u~h 2 Vh : a~ð~
8 wh 2 Wh :
ð2:7Þ
Let us discuss the continuity and inf–sup constants of the form a~. As far as the former one is e concerned, one has for any v 2 Ve and w 2 W
C. Canuto / Computers & Fluids 31 (2002) 421–436
425
j~ aðv; wÞj 6 jkvkV kwkW þ dkLvkZ 0 kSwkZ 0 6 c1 jkvkV kwkW~ þ c2 d1=2 kLvkZ 0 kwkW~ 6 j~kvkV~ kwkW~ ; e W , c2 is the norm of the operator d1=2 S 2 Lð W e ; Z0Þ where c1 is the norm of the imbedding W and e j :¼ maxðc1 j; c2 Þ. Concerning the inf–sup constant, we make the following Assumption 2.1. The space Z 0 and the operator S are chosen in such a way that the inf–sup condition a~kvh kV~ 6 sup
wh 2Wh
a~ðvh ; wh Þ kwh kW~
8 vh 2 Vh ;
ð2:8Þ
is fulfilled with a constant a~ > 0 independent of h, having the same order of magnitude as j~. This yields a constant of order of magnitude 1 in the error estimate ! j~ ku u~h kV~ 6 1 þ inf ku vh kV~ ; a~ vh 2Vh thus expressing the enhanced stabilization property of the modified approximation. In the rest of the paper, we will consider two relevant examples of the above functional setting.
3. The Stokes problem Given a bounded domain X Rn (n P 2) with Lipschitz boundary oX, we consider the Stokes problem with homogeneous boundary conditions: ~ þ rP ¼ ~ DU f in X;
ð3:1Þ
~¼0 rU
ð3:2Þ
in X;
~ ¼ 0 on oX: U ð3:3Þ n 1=2 ~ :¼ ðH 1 ðXÞÞ equipped with the norm k~ vkX~ ¼ ðr~ v; r~ vÞ . Here and in the sequel, we Let us set X 0 indicate by ð; ÞX both the inner product in L2 ðXÞ and the one in ðL2 ðXÞÞn . Let h; i denote the ~. Let M be a fixed subspace of L2 ðXÞ such that ~0 ¼ ðH 1 ðXÞÞn and X duality pairing between X 2 L ðXÞ ¼ M span (1); we denote by kqkM ¼ ðq; qÞ1=2 X the norm of M. 0 ~ ~ We assume that f is given in X , which implies the existence and uniqueness of the solution ~ ; P Þ in X ~ M. Setting ðU w; pÞX þ ðr ~ v; qÞX A½ð~ v; pÞ; ð~ w; qÞ ¼ ðr~ v; r~ wÞX ðr ~ and F ð~ w; qÞ ¼ h~ f;~ wi;
426
C. Canuto / Computers & Fluids 31 (2002) 421–436
~ ; P Þ as the solution of the variational problem we obtain ðU ~; P Þ 2 X ~ M : A½ðU ~ ; P Þ; ð~ ðU w; qÞ ¼ F ð~ w; qÞ
~ M: 8 ð~ w; qÞ 2 X
~ M equipped with Thus, in the abstract setting of the previous section, we choose V ¼ W :¼ X the graph norm; next, we identify L with the Stokes operator, i.e., we set D r L :¼ ; ð3:4Þ r 0 finally, we choose the forms a and f to be A and F, respectively. ~ and Mh M be finite dimensional subspaces. In the adaptive approximation of the ~h X Let X Stokes problem, it may be convenient to allow as much mutual independence as possible in the selection of the active degrees of freedom for velocity and pressure. Consequently, it might happen that the classical inf–sup condition bkqkM 6 sup ~h ~ w2X
ðr ~ w; qÞX k~ wkX~
8 q 2 Mh ;
is not satisfied (i.e., b ¼ 0), or it is satisfied with b > 0 but depending upon the discretization. In this case, the resulting discrete pressure is affected by spurious or weakly-spurious modes. In order to avoid such a drawback, a stabilization procedure of the type (2.6) has been recently proposed in Ref. [14]. Let us briefly recall the idea. Assume that the basis functions in Mh are selected within a Riesz basis Wp ¼ fwpk : k 2 Kp g of M, in the sense that there exists a finite set p p p Kph Kp such that k : k 2 Kh g. We recall that W is a Riesz basis if any q 2 M can be P Mh ¼ spanfw p written as q ¼ k2Kp qk wk with q ¼ ðqk Þk2Kp 2 ‘2 ðKp Þ (the series being convergent in L2 ðXÞ), and additionally kqkM kqk‘2 ðKp Þ
8 q 2 M:
ð3:5Þ
e p of L2 ðXÞ (the subspace of the zero By the Riesz representation theorem, there exists a basis W 0 p mean valued functions), which is biorthogonal to W , in the sense that ðwpk ; w~pk0 ÞX ¼ dk;k0
8 k; k0 2 Kp :
ð3:6Þ
Consequently, any q 2 M can be represented as X qk wpk with qk ¼ ðq; w~pk ÞX : q¼ k2Kp
~ satisfying w k 2 X Let us associate to each k 2 Kp an auxiliary velocity ~ r~ w k ¼ w~pk
in X:
Furthermore, let us assume that X qk ~ w k K kqk‘2 ðKp Þ 8 q ¼ ðqk Þk2Kp 2 ‘2 ðKp Þ: k2Kp
ð3:7Þ
ð3:8Þ
~ X
An example of construction of such auxiliary velocities, starting from a wavelet basis of the pressure space, will be detailed in the next subsection.
C. Canuto / Computers & Fluids 31 (2002) 421–436
427
~ :M!X ~ as Using the auxiliary velocities, we define the operator S X X ~q ¼ q¼ qk wpk 7! S qk ~ w k : k2Kp
ð3:9Þ
k2Kp
Proposition 3.1. There exists a constant c > 0 such that ~ qk~ 6 c kqk kS M X
8 q 2 M;
ð3:10Þ
furthermore, there exists a constant b > 0 such that 2
~ qi P b kqk hrq; S M
8 q 2 M:
ð3:11Þ
Proof. Inequality (3.10) follows from (3.8) and (3.5), whereas (3.11) is a consequence of Eq. (3.7) and again Eq. (3.8). We are now ready to introduce our modified formulation of the Stokes problem, which reads as ~2X ~ and p 2 M such that follows: find U ~ ; r~ ~; ðrU wÞX ðr ~ w; P ÞX ¼ h~ f;~ wi 8 ~ w2X ~ qi ¼ dh~ ~ qi ~ ; qÞ þ dhDU ~ rP ; S ðr U f;S X
ð3:12Þ 8 q 2 M:
ð3:13Þ
The reader should note formal analogies and actual differences of this formulation with the ~h and Galerkin/least-squares finite element method presented in Ref. [28]. Precisely, assume that X Mh are spaces of continuous finite elements on a triangulation Th whose maximum diameter of triangles is h; let hT denote the diameter of the triangle T. That formulation reads as follows: find ~h 2 X ~h and ph 2 Mh such that U ~h ; r~ ~h ; ðrU wh ÞX ðr ~ wh ; Ph ÞX ¼ ð~ f;~ wh ÞX 8 ~ wh 2 X X X 2 ~h rPh ; rqh Þ ¼ d ~h ; qh Þ þ d h ðD U h2T ð~ f ; rqh ÞT ðr U T X T T 2Th
8 qh 2 Mh :
T 2Th
As in our formulation, here the stabilization is achieved through a consistent term (depending on the residual of the momentum equation) added to the continuity equation. On the other hand, our ~ acting on the continuous level is here only mimicked on the discrete level by the operator S operator r on each triangle. Indeed, thinking multilevel, the weights h2T simulate the effect of ~0 , which transforming the elementwise L2 -inner products into a global inner product in the space X ~h rPh þ ~ f and the term rqh . is the correct function space for both the residual DU Going back to our variational problem, we observe that it can be written as e ½ðU ~; P Þ 2 X ~ M: A ~ ; P Þ; ð~ ðU w; qÞ ¼ Fe ð~ w; qÞ
~ M; 8 ð~ w; qÞ 2 X
after defining e ½ð~ ~ qi; A v; rÞ; ð~ w; qÞ :¼ A½ð~ v; rÞ; ð~ w; qÞ þ dhD~ v rr; S ~ qi: Fe ð~ w; qÞ :¼ F ð~ w; qÞ dh~ f;S
ð3:14Þ
428
C. Canuto / Computers & Fluids 31 (2002) 421–436
The proposed stabilization can be cast into the abstract setting of the previous section, in the ~0 M. Denote by J the canonical isomorphism between X ~0 following way. Choose Z 0 :¼ W 0 ¼ X ~ and X , i.e., the inverse of the minus Laplacian operator supplemented by Dirichlet boundary conditions. Then, we can choose the inner product f ; J~ gi þ ðr; qÞX ððð~ f ; rÞ; ð~ g; qÞÞÞX~0 M :¼ h~ ~0 M. The space Ve coincides with V ¼ X ~ M. Furthermore, let us choose W e :¼ W ¼ X ~
in X 0 ~ M!X ~ M as M and let us define S : X ~ q; 0Þ: Sð~ w; qÞ :¼ ðJ1 S ~ M, Then, recalling Eq. (3.4), it is easily seen that, for all ð~ v; rÞ; ð~ w; qÞ 2 X ~ qi: ððLð~ v; rÞ; Sð~ w; qÞÞÞX~0 M ¼ hD~ v rr; S e and f~ :¼ Fe . Assumption (2.1) is fulfilled Thus, Eq. (3.14) coincides with Eq. (2.6) if we set a~ :¼ A provided d is not too large. Indeed, taking into account Eqs. (3.10) and (3.11), the following result can be established (see Ref. [14] for the details of the proof). Theorem 3.2. Let us assume that d satisfy the bound d 6 b =c2 . Then, there exists a constant b > 0 (independent of d) such that inf
sup
ð~ v;rÞ2V ð~ w;qÞ2V
e ½ð~ A v; rÞ; ð~ w; qÞ P b; kð~ v; rÞkX~ M kð~ w; qÞkX~ M
~ M is ~ M or any finite dimensional subspace X ~h Mh ; and the norm in X where V equals either X 2 2 2 vkX~ þ dkrkM . scaled as kð~ v; rÞkX~ M ¼ k~ Using this result, it is then straightforward to obtain optimal a priori and a posteriori error ~ ; P Þ of the Stokes problem ((3.1)–(3.3)) and the solution ð~ estimates between the solution ðU u; pÞ of the Galerkin discretization: e ½ð~ ~h Mh : ~h Mh : A u; pÞ; ð~ w; qÞ ¼ Fe ð~ w; qÞ 8 ð~ w; qÞ 2 X ð~ u; pÞ 2 X Numerical results, obtained by this formulation using the auxiliary stabilizing velocities described below, show good agreement with the theoretical predictions. They can be found in Ref. [14]. 3.1. An example We now indicate how the auxiliary velocities satisfying Eqs. (3.7) and (3.8) can be constructed by means of appropriate wavelet bases. Starting from compactly supported (bi)orthogonal scaling function and wavelet systems on the real line, such as the Daubechies family or any family of B-splines, it is possible to construct biorthogonal scaling function and wavelet systems on the reference interval I ¼ ð0; 1Þ. For the technical details, we refer, e.g., to Ref. [18]. For each level index j P some j0 , there exist families e j ¼ fn~jk ; k 2 Kj g of mutually biorthogonal scaling functions (i.e., Nj ¼ fnjk ; k 2 Kj g and N
C. Canuto / Computers & Fluids 31 (2002) 421–436
R
429
S
I njk n~jk0 ¼ dk;k 0 for all k; k 0 2 Kj ) such that, setting VjI ¼ span Nj , one has VjI Vjþ1 ; 8j, and j VjI e j . Furthermore, there exist is dense in L2 ðIÞ; similar properties hold for the sequence VejI ¼ span N e gjh ; h 2 Hj g of mutually biorthogonal wavelets such families j ¼ fgjh ; h 2 Hj g and j ¼ f~ I I e j , one has V I ¼ V I W I and Ve I ¼ Ve I W e e I , with that, setting Wj ¼ span !j and W j ¼ span ! jþ1 j j jþ1 j j I I I I e e Wj ? V j and W j ? Vj . The construction can be carried on in a way that at least the zeroth-order moment (i.e., the average over I) of each wavelet vanishes; in addition, the dual scaling functions and wavelets belong to H01 ðIÞ. ^ ¼ ð0; 1Þn . Biorthogonal With such a construction at hand, let us first consider the domain X ^ Þ can be defined by taking tensor products of univariate scaling functions and Riesz bases in L2 ðX ^ wavelets. Precisely, for any j P j0 , let WXj denote the set of wavelets w ¼ ni¼1 #i , such that each #i is ~ X^ be either a scaling function n 2 Nj or a wavelet g 2 j , and at least one #i is a wavelet. Let W j ^ ^ Þ, and let UX^ e X be a tensor product basis of ð n Ve I Þ \ L2 ðX defined similarly. Furthermore, let U i¼1 j0 0 j0 j0 be the companion biorthogonal set in ni¼1 VjI0 . Finally, set [ ^ [ ^ ^ ^ e X: e p;X^ ¼ U e X^ [ W WXj ; W Wp;X ¼ UXj0 [ j j0 I
j P j0
j P j0
^ ^ Þ ¼ span W e p;X^ ; more precisely, it is Then, we define M :¼ span Wp;X . One also checks that L20 ðX ^ p p e p;X has the form w~ ¼ n #~i , where #~i ð0Þ ¼ #~i ð1Þ ¼ 0 for all i and there easily seen that any w~k 2 W k R i¼1 exists at least one index, say m ¼ mðkÞ such that I #~m ¼ 0. Then, we associate to w~pk the auxiliary ^ velocity ~ w ;k X defined by ( ^ ^ X ;X ~ where w ; k ¼ ð0; . . . ; wk ; . . . ; 0Þ R x^m ^ ;X ~ xÞ ¼ ð 0 #m ðsÞ dsÞ i6¼m #~i ð^ xi Þ: wk ð^
Condition (3.7) is trivially satisfied, whereas Eq. (3.8) can be checked by a version of the so-called ‘vaguelettes’ Lemma [32]. ^ under a smooth mapping F. Set G ¼ F 1 and let DG denote the Next, let X be the image of X Jacobian matrix of G. Then, we define ^
p;X Wp;X ¼ fwp;X k :¼ wk # Gg
and
e p;X ¼ fw~p;X :¼ jDGjw~p;X^ # Gg; W k k
^ X ~ ;X defined by and we associate to each w~p;X the Piola transform of ~ w ; k , namely, the velocity wk k ^ ;X ; X ~ wk # GÞ. Again, Eqs. (3.7) and (3.8) are satisfied. wk :¼ jDGjDGð~ For more general domains, one can resort to domain decomposition. In this case, pressure are discontinuous between subdomains; this implies that the stabilizing wavelets can be built independently on each subdomain as shown above, and then extended by zero outside the subdomain of definition. This completes the construction of the stabilizing wavelets.
4. The convection–diffusion problem Let X be again a bounded domain contained in Rn (n P 1), with Lipschitz boundary oX. We consider the problem
430
C. Canuto / Computers & Fluids 31 (2002) 421–436
mDu þ a ru þ bu ¼ f ; u ¼ 0;
in X; on oX;
ð4:1Þ
n
where m > 0 is a constant, a 2 ðW 1;1 ðXÞÞ , b 2 L1 ðXÞ and there exists a constant c > 0 (of order of magnitude 1) such that 12r a þ b P c
ð4:2Þ
a:e: in X:
Problem (4.1) is written variationally as in Eq. (2.2), if we assume that f 2 H 1 ðXÞ and we set V ¼ W :¼ H01 ðXÞ, Z Z Z ð4:3Þ aðv; wÞ :¼ m rv rw dx þ a rvw dx þ bvw dx; X
X
X
and hf ; wi :¼ H 1 ðXÞ hf ; wiH 1 ðXÞ : 0
The form a is continuous and coercive on V; precisely, one has for all v 2 V mkrvk2L2 ðXÞ þ ckvk2L2 ðXÞ 6 aðv; vÞ 6 Kðkrvk2L2 ðXÞ þ kvk2L2 ðXÞ Þ
ð4:4Þ
with K proportional to maxðm; kakðL1 ðXÞÞn ; kbkL1 ðXÞ Þ. If m kakðL1 ðXÞÞn , we are precisely in one of the situations considered in Section 2, namely the inf–sup constant (here, the coercivity constant) and the continuity constant are sensibly different from each other. As a consequence, the solution of Problem (4.1) is allowed to exhibit strong gradients, which usually take the form of highly localized transition layers (boundary and interior layers). If a standard Galerkin method is used to approximate that problem, and if functions in Vh do not have the potential to faithfully represent such layers, the corresponding discrete solution exhibits strong gradients as well, but they manifest themselves as oscillations, more or less localized around the layers. These oscillations can be prevented by a stabilization procedure. The key idea underlying the stabilization method proposed in Ref. [11] and developed in Refs. [12,13], is to obtain a better continuity estimate than in Eq. (4.4). Just for the sake of simplicity, let us assume that r a ¼ 0 in X; furthermore, let us set Dv :¼ a rv; and let, as before, ð; ÞX denote the L2 ðXÞ-inner product. The convective term ðDv; wÞX can be bounded as jðDv; wÞX j 6 kDvkL2 ðXÞ kwkL2 ðXÞ
8 v; w 2 H01 ðXÞ;
on the other hand, by partial integration, it can also be bounded as jðDv; wÞX j ¼ j ðv; DwÞX j 6 kvkL2 ðXÞ kDwkL2 ðXÞ : The two bounds and the theory of interpolation of function spaces suggest to look for a more symmetric bound, involving the same norm of both v and w. Such bound should control Dv and Dw in the norm of an intermediate space between L2 ðXÞ and H 1 ðXÞ (remember that if v 2 L2 ðXÞ, then Dv 2 H 1 ðXÞ and kDvkH 1 ðXÞ K kvkL2 ðXÞ ). Taking into account the boundary conditions, a 1=2 natural candidate as the intermediate space is the dual of the space H00 ðXÞ, which is the in-
C. Canuto / Computers & Fluids 31 (2002) 421–436
431
terpolation space of order 1=2 between H01 ðXÞ and L2 ðXÞ (see, e.g. Ref. [31]). For simplicity, let us 1=2 1=2 0 set H1=2 ðXÞ :¼ H00 ðXÞ and H1=2 ðXÞ :¼ ðH00 ðXÞÞ . Furthermore, let ðð; ÞÞH1=2 ðXÞ be an inner product in H1=2 ðXÞ, whose definition will be discussed later on. In the sequel, we shall first present the conceptually simplest way to build-up a stabilized variational formulation of Eq. (2.2), which guarantees an H1=2 ðXÞ-control on the streamline derivative. Next, we shall introduce a more refined and flexible form of the stabilization device. Let us start by considering the space X 1=2 ðXÞ :¼ fv 2 H01 ðXÞ : Lv 2 H1=2 ðXÞg
ð4:5Þ
equipped with the norm kvkX 1=2 ðXÞ :¼ kvkL2 ðXÞ þ m1=2 krvkL2 ðXÞ þ d1=2 kLvkH1=2 ðXÞ ;
ð4:6Þ
where d > 0 is a suitable parameter. Note that this norm controls Lv, rather than Dv. Indeed, we aim at introducing a consistent stabilization of Eq. (2.2), based on the residual f Lv. However, we shall recover the control on kDvkH1=2 ðXÞ from the control on kLvkH1=2 ðXÞ , by some sort of ‘post-processing’ analysis. From now on, we assume that f 2 H1=2 ðXÞ, so that u 2 X 1=2 ðXÞ. At this point, we introduce the stabilized bilinear form A : X 1=2 ðXÞ X 1=2 ðXÞ ! R;
Aðv; wÞ :¼ aðv; wÞ þ dððLv; LwÞÞH1=2 ðXÞ
ð4:7Þ
and the stabilized linear form F : X 1=2 ðXÞ ! R;
F ðwÞ :¼ hf ; wi þ dððf ; LwÞÞH1=2 ðXÞ ;
ð4:8Þ
and we consider the variational problem: find u 2 X 1=2 ðXÞ such that Aðu; wÞ ¼ F ðwÞ
8 w 2 X 1=2 ðXÞ:
ð4:9Þ
e ¼ Thus, with respect to the abstract framework of Section 2, we have Z 0 ¼ H1=2 ðXÞ, Ve ¼ W 1=2 X ðXÞ and S ¼ L. With this setting, we obtain a~ ¼ A and f~ ¼ F . Coming to Assumption 2.1, we note that the form A is uniformly coercive on X 1=2 ðXÞ. Indeed, Aðv; vÞ P mkrvk2L2 ðXÞ þ ckvk2L2 ðXÞ þ dkLvk2H1=2 ðXÞ P akvk2X 1=2 ðXÞ
8 v 2 X 1=2 ðXÞ;
ð4:10Þ
for a constant a > 0 independent of m. On the other hand, as far as continuity is concerned, we have the symmetric bound jAðv; wÞj 6 CkvkX 1=2 ðXÞ kwkX 1=2 ðXÞ
8 v; w 2 X 1=2 ðXÞ;
but the constant C may blow-up with m1 (see Ref. [12, estimate (3.12)]). However, introducing a slightly stronger norm of v yields the uniform (in m) estimate (see Ref. [12, Proposition 4.1]) jAðv; wÞj K ðkvkX 1=2 ðXÞ þ d1=2 kvkH1=2 ðXÞ ÞkwkX 1=2 ðXÞ
8 v; w 2 X 1=2 ðXÞ:
ð4:11Þ
If a Galerkin approximation of Eq. (4.9) is considered, the bounds (4.10) and (4.11) imply a uniform error estimate in the X 1=2 ðXÞ-norm. We shall not detail such analysis here, since a similar result will be obtained later on, for the modified stabilization setting. Let us discuss the definition of the inner product in H1=2 ðXÞ, which so far was not made precise. We propose to use a multilevel decomposition of L2 ðXÞ, generated by two mutually
432
C. Canuto / Computers & Fluids 31 (2002) 421–436
e :¼ fw~k : k 2 Kg. Each k is a couple, k :¼ ðj; kÞ, biorthogonal Riesz bases W :¼ fwk : k 2 Kg and W where j P j0 (for some j0 2 N) is the ‘level’ (or ‘scale’) index and k 2 Kj (for some finite set Kj Z) is the ‘position’ index, associated with the basis function wk . For instance, wk might have (at least locally) the form wk ðxÞ ¼ 2nj=2 wð2j x kÞ, where w is an appropriate ‘generating function’. For notational convenience, we set jkj :¼ j. The elements of both bases are compactly S k similarly, we have supported functions; denoting by Sk the support of wk and defining e S k 2jkj . We recall that biorthogonality means that the relations diam Sk diam e ðwk ; w~k0 ÞX ¼ dk;k0
8 k; k0 2 K
ð4:12Þ
are fulfilled, whereas the Riesz basis condition says that any v 2 L2 ðXÞ can be written as X v¼ ðv; w~k ÞX wk k2K
with kvk2L2 ðXÞ
X
jðv; w~k ÞX j2 :
k2K
In addition, let us suppose that W H01 ðXÞ and that any g 2 H 1 ðXÞ can be written as X g¼ hg; wk iw~k : k2K
Let us introduce the scale of Sobolev spaces Hr ðXÞ, 1 6 r 6 1, defined as 8 r 1=2 < r 6 1; H ðXÞ; > > < 01=2 r ¼ 1=2; H00 ðXÞ; Hr ðXÞ :¼ r > ðXÞ; 0 6 r < 1=2; H > : 0 r ðH ðXÞÞ ; 1 6 r < 0: e allow the characterization of the spaces Hr ðXÞ. Precisely, for We assume that the bases W and W any 0 6 r 6 1, we require that X 22rjkj jðv; w~k ÞX j2 8 v 2 Hr ðXÞ; ð4:13Þ kvk2Hr ðXÞ k2K
and kgk2Hr ðXÞ
X
22rjkj jhg; wk ij2
8 g 2 Hr ðXÞ:
ð4:14Þ
k2K
Families of biorthogonal basis functions in non-trivial domains have been constructed in the last few years; see, e.g. Refs. [15,16,20,22,24]. Inspired by Eq. (4.14), we are led to use the following H1=2 ðXÞ-inner product X ððf ; gÞÞH1=2 ðXÞ :¼ 2jkj hf ; wk ihg; wk i k2K 1
in the stabilized forms (4.7) and (4.8). The parameter d should be taken proportional to kakL1 ðXÞ therein. Obviously, in general infinitely many terms appear on the right-hand side, making the
C. Canuto / Computers & Fluids 31 (2002) 421–436
433
practical implementation of such inner product unfeasible. However, once a finite dimensional trial/test space Vh H01 ðXÞ has been chosen, it is possible ([12, Section 5]) to truncate the inner product to a finite sum (depending on the ‘resolution power’ of Vh ), without affecting the stabilization effect. To be more precise, if Vh contains all basis functions wk with jkj smaller than some jðhÞ, then the sum in Eq. (4.18) can be restricted to those k satisfying jkj P jðhÞ. On the other hand, if we assume that the following (absolutely reasonable) inverse inequality kDvh kL2 ðXÞ K h1=2 kDvh kH1=2 ðXÞ
8 vh 2 Vh ;
holds true, then there exists a level J ðhÞ j log2 hj such that all the details with jkj > J ðhÞ are inessential in the stabilization process. In conclusion, the H1=2 ðXÞ-inner product defined above can be replaced by the truncated inner product X 2jkj hf ; wk ihg; wk i: ððf ; gÞÞh :¼ jðhÞ 6 jkj 6 J ðhÞ
Let us now turn to the more refined version of our multilevel stabilization approach. The key observation is that––thinking multilevel––the stabilization mechanism should provide extra H1=2 ðXÞ-control only on the ‘low-level’ components of Dv, since the ‘high-level’ components are implicitly controlled by the diffusive term mDv. In other words, one should insert information on the local relationship between convection and diffusion. In classical SUPG-stabilizations of finite element methods (see, e.g., Ref. [29]), such information is inserted at the discrete level through a switch in the stabilization parameter, based on a local (in position) Peclet number. In our approach, we are going to act at the continuos level (i.e., before any finite-dimensional discretization) and on the basis of a local (in position but also in ‘scale’) Peclet number. To be precise, for any index k 2 K, let us introduce the multiscale Peclet number
Pek :¼ c
kakL1 ðSk Þ 2jkj m
;
ð4:15Þ
where c is a tuning parameter depending on the multilevel basis W. Next, set ! ( ) 2jkj if Pek > 1 2jkj 22jkj c kakL1 ðS Þ k ; ; ¼ min ; sk :¼ 22jkj c kakL1 ðSk Þ m if Pe 6 1 k m and define the inner product X ððf ; gÞÞs :¼ sk hf ; wk ihg; wk i;
ð4:16Þ
k2K
as well as the associated norm kgks :¼ ððg; gÞÞ1=2 s . Recalling Eq. (4.14), such bilinear form behaves 1=2 ðXÞ-inner product on the ‘low-level’ components of its arguments, and as an like a (scaled) H H 1 ðXÞ-inner product (scaled by m1 ) on the ‘high-level’ components. In particular, Eq. (4.16) is e ¼ X1=2 ðXÞ, where defined on H 1 ðXÞ. Therefore, now we can choose Z 0 ¼ H 1 ðXÞ, so that Ve ¼ W 1=2 1 X ðXÞ :¼ H0 ðXÞ equipped with the norm kvkX1=2 ðXÞ :¼ kvkL2 ðXÞ þ m1=2 krvkL2 ðXÞ þ kLvks : The continuous problem we aim at discretizing reads as follows: find u 2 H01 ðXÞ such that
434
C. Canuto / Computers & Fluids 31 (2002) 421–436
8 w 2 H01 ðXÞ;
Aðu; wÞ ¼ FðwÞ with
A : H01 ðXÞ H01 ðXÞ ! R;
Aðv; wÞ :¼ aðv; wÞ þ ððLv; LwÞÞs
ð4:17Þ
and F : H01 ðXÞ ! R;
FðwÞ :¼ ðf ; wÞ þ ððf ; LwÞÞs :
Note the explicit expression of the stabilizing term appearing in Eq. (4.17): ððLv; LwÞÞs ¼
X
X 22jkj 2jkj aðv; wk Þaðw; wk Þ: aðv; w Þaðw; w Þ þ k k c kakL1 ðSk Þ m Pe 6 1
Pek >1
ð4:18Þ
k
For the new form, one can prove the coerciveness inequality akvk2X1=2 ðXÞ 6 Aðv; vÞ 8 v 2 H01 ðXÞ; for a constant a 1, and the uniform continuity bound jAðv; wÞj K ðkvkX1=2 ðXÞ þ kak1=2 L1 ðXÞ kvkH1=2 ðXÞ ÞkwkX1=2 ðXÞ
8 v; w 2 H01 ðXÞ;
which are similar to the analogous ones (4.10) and (4.11). Next, given any family fVh gh>0 of finite dimensional subspaces of H01 ðXÞ, we consider the stabilized Galerkin approximation: find uh 2 Vh such that Aðuh ; wh Þ ¼ Fðwh Þ
8 wh 2 Vh :
ð4:19Þ
Standard arguments lead to the error bound 1=2
ku uh kX1=2 ðXÞ K inf ðku vh kX1=2 ðXÞ þ kakL1 ðXÞ ku vh kH1=2 ðXÞ Þ: vh 2Vh
ð4:20Þ
Let us assume that the family fVh gh>0 provides classical approximation properties, expressed by the bounds inf kv vh kHr ðXÞ K hsr jvjH s ðXÞ
vh 2Vh
8 v 2 H s ðXÞ \ H01 ðXÞ;
for 0 6 r 6 1 and 1 6 s < some s0 . Then, one can establish [13, Theorem 3.4] the error estimate ku uh kL2 ðXÞ þ m1=2 kru ruh kL2 ðXÞ þ kLu Luh ks K ½Chs þ m1=2 hs1 þ Bhsð1=2Þ jujHs ðXÞ ; ð4:21Þ where B and C are constants depending on the velocity field a (they are surely K 1 if a is locally Lipschitz continuous at the points where it vanishes). As anticipated, from this result we deduce a bound on the H1=2 ðXÞ-norm of Du Duh . Indeed, expressing the negative order norms in the equivalent multilevel form (4.14), one obtains 1=2 1=2 kru ruh kL2 ðXÞ kak1=2 L1 ðXÞ kDu Duh kH1=2 ðXÞ K Cku uh kL2 ðXÞ þ ðkakL1 ðXÞ þ BÞm
þ kLu Luh ks :
C. Canuto / Computers & Fluids 31 (2002) 421–436
435
This and Eq. (4.21) imply that the error on the streamline derivative measured in the H1=2 ðXÞnorm can be bounded in terms of the right-hand side of Eq. (4.21). The stabilized scheme (4.19) is still not implementable, as in the stabilizing term (4.18) infinitely many details appear. However, the sum in Eq. (4.18) can be truncated to a finite number of details, depending on Vh . The truncation can be accomplished precisely in the same way as explained above for the ‘global version’ of our stabilization. 5. Conclusions In our opinion, one of the interests of the proposed multilevel framework is that it sheds new light on classical stabilization strategies such as Galerkin/least squares or SUPG. Usually, these methods add a locally finite number of terms to the pure Galerkin approximation. Such perturbation terms should be interpreted as, or at least related to, a particular truncated version of a Sobolev-type inner product, as they mimic the inner product at the discrete level. For the Stokes problem, the relationship has been discussed in Section 3. For model 1D convection–diffusion problems, the interpretation has been rigorously established in Ref. [12, Section 6]: the classical SUPG stabilization considered in Ref. [29] is shown to coincide with our multilevel-stabilized method (4.19), in which the (truncated) inner product is defined via piecewise linear interpolatory wavelets. The theoretical results briefly described above apply to this situation. The functional framework we have developed indicates that the ‘order 1=2’ regularity setting is appropriate for the correct understanding of numerical discretizations of convection–diffusion problems, in the convection-dominated regime. Other recent and independent results [10,30], which involve order 1=2 norms in the study of stabilized approximations to convection–diffusion problems, confirm our perspective. The extension of our investigation to more complex equations of fluid dynamics, such as the incompressible Navier–Stokes equations, is underway. References [1] Arandiga F, Donat R, Harten A. Multiresolution based on weighted averages of the hat function. II. Nonlinear reconstruction techniques. SIAM J Sci Comput 1999;20:1053–93. [2] Babuska I, Aziz AK. Survey lectures on the mathematical foundation of the finite element method. In: Aziz AK, editor. The mathematical foundation of the finite element method with applications to partial differential equations. New York: Academic Press; 1972. p. 3–359. [3] Baiocchi C, Brezzi F. Stabilization of unstable numerical methods, Atti del Simposio Internazionale Problemi attuali dell’Analisi e della Fisica Matematica. Taormina, 15–17 ottobre 1992 (1993), pp. 59–64. [4] Balder R, Zenger C. The solution of multidimensional real Helmholtz equations on sparse grids. SIAM J Sci Comput 1996;17:631–46. [5] Bank RE. In: Hierarchical bases and the finite element method, Acta Numerica. Cambridge: Cambridge University Press; 1996. p. 1–43. [6] Bertoluzza S. Stabilization by multiscale decomposition. Appl Math Lett 1998;11:129–34. [7] Bramble JH, Lazarov RD, Pasciak JE. A least-squares approach based on a discrete minus one inner product for first order systems. Math Comput 1997;66:935–55. [8] Bramble JH, Pasciak JE, Xu J. Parallel multilevel preconditioners. Math Comput 1990;55:1–22. [9] Brezzi F, Fortin M. Mixed and hybrid finite element methods. New York: Springer; 1991. [10] Brezzi F, Marini D, S€ uli E. Residual-free bubbles for advection–diffusion problems: the general error analysis. Numer Math 2000;85:31–47.
436
C. Canuto / Computers & Fluids 31 (2002) 421–436
[11] Bertoluzza S, Canuto C, Tabacco A. Negative norm stabilization of convection–diffusion problems. Appl Math Lett 2000;13:121–7. [12] Bertoluzza S, Canuto C, Tabacco A. Stable discretizations of convection–diffusion problems via computable negative-order inner products. SIAM J Numer Anal 2000;38:1034–55. [13] Canuto C. Multilevel stabilization of convection–diffusion problems by variable-order inner products. Computing 2001;66:121–38. [14] Canuto C, Masson R. Stabilized wavelet approximations of the Stokes problem. Math Comput 2001;70:1397–416. [15] Canuto C, Tabacco A, Urban K. The wavelet element method. Part I: construction and analysis. Appl Comput Harm Anal 1999;6:1–52. [16] Canuto C, Tabacco A, Urban K. The wavelet element method. Part II: realization and additional features in 2D and 3D. Appl Comput Harm Anal 2000;8:123–65. [17] Chen M, Temam R. Incremental unknowns for solving partial differential equations. Numer Math 1991;59:255–71. [18] Cohen A. Wavelet methods in numerical analysis, In: Ciarlet PG, Lions JL, editors. Handbook of numerical analysis, vol. VII. Amsterdam: Elsevier; 2000. [19] Cohen A, Dahmen W, DeVore R. Adaptive wavelet methods for elliptic operator equations––convergence rates. Math Comp 2001;70:27–75. [20] Cohen A, Masson R. Wavelet adaptive methods for 2nd order elliptic problems. Domain decomposition and boundary conditions. Numer Math 2000;86:193–238. [21] Dahmen W. Stability of multiscale transformations. J Fourier Anal Appl 1996;4:341–62. [22] Dahmen W. In: Wavelet and multiscale methods for operator equations Acta Numerica 6. Cambridge: Cambridge University Press; 1997. p. 55–228. [23] Dahmen W, Kunoth A. Multilevel preconditioning. Numer Math 1992;63:315–44. [24] Dahmen W, Schneider R. Composite wavelet bases for operator equations. Math Comput 1999;68:1533–67. [25] DeVore R. In: Nonlinear approximation Acta Numerica. Cambridge: Cambridge University Press; 1998. p. 53– 150. [26] Farge M, Schneider K, Kevlahan N. Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis. Phys Fluids 1999;11:2187–201. [27] Harten A. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Comm Pure Appl Math 1995;48:1305–42. [28] Hughes TJR, Franca LP, Balestra M. A new finite element formulation of computational fluid dynamics: a stable Petrov–Galerkin formulation of the Stokes problem accomodating equal-order interpolations. Comput Methods Appl Mech Engng 1986;59:85–99. [29] Hughes TJR, Franca LP, Hulbert GM. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput Meth Appl Mech Engng 1989;73: 173–89. [30] Kellogg RB, Stynes M. The intrinsic norm finite element method for convection–diffusion problems, Preprint 1999-1, Department of Mathematics, University College, Cork, Ireland. [31] Lions JL, Magenes E. Problemes aux limites non homogenes. Paris: Dunod; 1968. [32] Meyer Y. Ondelettes et operateurs––tomes 1 et 2. Paris: Hermann; 1990. [33] Oswald P. Multilevel finite element approximation. Theory and applications. Stuttgart: Teubner; 1994. [34] Roos H-G, Stynes M, Tobiska L. Numerical methods for singularly perturbed differential equations. Convection– diffusion and flow problems. Berlin: Springer; 1996.