On the stability of the reduced basis method for Stokes equations in parametrized domains

Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260 www.elsevier.com/locate/cma

On the stability of the reduced basis method for Stokes equations in parametrized domains Gianluigi Rozza b

a,*

, Karen Veroy

b

a CMCS – Modeling and Scientific Computing, EPFL – Ecole Polytechnique Fe´de´rale de Lausanne, Station 8, CH1015 Lausanne, Switzerland Mechanical Engineering Department, MIT – Massachusetts Institute of Technology, 77 Mass Avenue, Room 3-264, Cambridge, MA 02139, USA

Received 30 March 2006; received in revised form 6 September 2006; accepted 7 September 2006

Abstract We present an application of reduced basis method for Stokes equations in domains with affine parametric dependence. The essential components of the method are (i) the rapid convergence of global reduced basis approximations – Galerkin projection onto a space WN spanned by solutions of the governing partial differential equation at N selected points in parameter space; (ii) the off-line/on-line computational procedures decoupling the generation and projection stages of the approximation process. The operation count for the on-line stage – in which, given a new parameter value, we calculate an output of interest – depends only on N (typically very small) and the parametric complexity of the problem; the method is thus ideally suited for the repeated and rapid evaluations required in the context of parameter estimation, design, optimization, and real-time control. Particular attention is given (i) to the pressure treatment of incompressible Stokes problem; (ii) to find an equivalent inf–sup condition that guarantees stability of reduced basis solutions by enriching the reduced basis velocity approximation space with the solutions of a supremizer problem; (iii) to provide algebraic stability of the problem by reducing the condition number of reduced basis matrices using an orthonormalization procedure applied to basis functions; (iv) to reduce computational costs in order to allow real-time solution of parametrized problem.  2006 Elsevier B.V. All rights reserved. Keywords: Parametrized Stokes equations; Reduced basis methods; Approximation stability; Inf–sup condition; Supremizer; Galerkin approximation; Algebraic stability; Gram–Schmidt basis orthogonalization

1. Introduction The optimization, control, design and characterization of an engineering component or system requires the prediction of certain ‘‘quantities of interest,’’ or performance metrics, which we denote outputs – for example velocity field, maximum stresses, maximum temperatures, heat transfer, flow rates, vorticity, or lifts and drags. These outputs are typically expressed as functionals of field variables associated with parametrized partial differential equations (P2DE) which describe the physical behavior of the system. The parameters, which we shall denote inputs, serve to *

Corresponding author. Tel.: +41 21 6932733; fax: +41 21 6934303. E-mail addresses: gianluigi.rozza@epfl.ch (G. Rozza), veroy@alum. mit.edu (K. Veroy). 0045-7825/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.09.005

identify a particular ‘‘configuration’’ of the components: these inputs may represent design or decision variables, such as geometry, or characterization variables, such as physical properties – for example in inverse design problems. We thus get an implicit input–output relationship whose evaluation demands the solution of the underlying partial differential equations. The development of computational methods is permitting rapid and reliable evaluation of this input–output relationship in the design, optimization and control contexts. The approach used is based on the reduced basis method, first introduced in the late 1970s for nonlinear structural analysis, and later developed more broadly in the 1980s and 1990s (see [1,3] for general framework, [24] for extension to multi-parameter problems, [4,16] for application in nonlinear problems). Recent applications of the method

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

1245

are contained in [14,18–20,25,33]. For the general a priori convergence properties of reduced basis method see, for example, [11,12]. In the application we use global approximation spaces and we exploit off-line/on-line computational decompositions (see [1] for application of this strategy within the reduced basis context). These steps allow us – for the class of ‘‘parameter-affine’’ problems – to reliably decouple the generation and projection stages of reduced basis approximation, thereby effecting computational economies of several orders of magnitude. In this paper we deal with the application and investigation of reduced basis techniques for Stokes equations in parametrized domains, focusing our attention to pressure treatment, to approximation stability (by the introduction of an equivalent inf–sup condition) and to algebraic stability (by achieving condition number reduction and basis orthonormalization). Preliminary applications of reduced basis techniques to incompressible viscous flow problems are given in [15,10] and also [9] for optimal control problems; however in these contributions the pressure approximation is not considered and stability of solutions is not discussed. More recent works dealing with divergence-free velocity spaces, physical parameters and focused on a posteriori error estimates are the ones by Patera et al. [31,14,32]. After this introduction, in Section 2 we formulate the problem for parametrized Stokes equations. In Section 3 we introduce the Stokes reduced basis formulation and analyze the stability of the approximation. At the end of this section we introduce an application of the method to a ‘‘T-bypass’’ configuration. In Section 4 we study a procedure to control N more tightly and to apply an adaptive procedure for the basis construction using a suitable error projection. In Section 5 we analyze the inf– sup condition which is necessary for the stability of reduced basis approximation. Then in Section 6 we present some numerical results for three different test cases based on the same geometrical configuration. In Section 7 we deal with algebraic stability problem and we present possible solutions to achieve it by orthonormalization. In Section 8 we mention possible developments for viscous flow and shape optimization problem in haemodynamics within the reduced basis framework.

bw [ C b in . After introbD ¼ C such that the Dirichlet portion C ^ ducing a lift function L^gin (extension of non-homogeneous boundary conditions to the interior of the domain) such ^^g and ^uj ¼ 0, the weak formulation of (1) that ^u ¼ ~u  L CD in b  H 1^ ð XÞ, b (where H 1^ ð XÞ b ¼ reads: find ^u 2 Yb ¼ H 1C^ D ð XÞ CD CD b : u ¼ 0 on C b such that b D g), p^ 2 M b ¼ L2 ð XÞ fu 2 H 1 ð XÞ 8 Z Z > > ^ dX w dX  ^pr  w m r^u  r^ > > > ^ ^ X X > > > Z < ^ i; 8^ ^ dX þ h Fb 0 ; w w 2 Yb ; ¼ ^f  w ð2Þ > ^ X > > > Z > > > > b 0 ; ^qi; 8^q 2 M b; : ^qr  ^u dX ¼ h G

2. The parametrized Stokes problem

We want to build a system of parametrized partial differential equations (P2DEs) depending on a set of geometrical parameters l ¼ fl1 ; . . . ; lP g 2 Dl  RP , that we call parameters. Then (7) is traced back to a reference domain b r into Xr. For by an affine mapping of the sub-domains X r r b ^ 2 X , r = 1, . . . , R, its image x 2 X is given by any x

We consider the following steady Stokes problem [5] in a b  R2 with boundary C: b domain X 8 ^ b mM~ u þ r^p ¼ f in X; > > > > > b > > > u > : m o~ b out ; ^ p^ n ¼ 0 on C o^ n where ~ u is the velocity, ^p the pressure, ^f a force field, m a kinematic viscosity and ^ n is the normal unit vector. The b is split into three components, C b w; C b in ; C b out , boundary o X

^ X

b 0 are terms due to non-homogeneous Dirichlet Fb 0 ; G b in , on C b out we have a boundary condition (~u ¼ ^gin ) on C stress-free Neumann condition. We assume that the dob ¼ SR X b r , so that main is made up of R sub-domains: X r¼1 the bilinear and linear forms of the weak formulation read for 1 6 i, j 6 2 and ^mi;j ¼ mdi;j as ^i ¼ hc A^u; w

R Z X

^mij

^r X

r¼1

o^u o^ w b  d X; o^xi o^xj

ð3Þ

where summation over i and j is understood, b p; w ^i ¼  h B^

R Z X r¼1

^r X

b ^ d X; ^pr  w

ð4Þ

^ i ¼ h Fb s ; w ^ i þ h Fb 0 ; w ^i h Fb ; w

ð5Þ

and ^i ¼ h Fb s ; w

R Z X r¼1

^r X

b h Fb 0 ; w ^i ^ d X; f^ w ð6Þ

^^g ; w ^ i; ¼ h c AL in 0 b ; ^qi ¼ h B^ b q; L ^^g i: hG in

Then we write ( b p; w ^ i þ h B^ ^ i ¼ h Fb ; w ^ i; hc A^u; w b q; ^ui ¼ h G b 0 ; ^qi; h B^

8^ w 2 Yb ;

b: 8^q 2 M

^Þ ¼ Gr ðlÞ^ x ¼ Gr ðl; x x þ gr ;

1 6 r 6 R;

ð7Þ

ð8Þ

we thus write o oxj o o ¼ ¼ Grji ðlÞ o^xi o^xi oxj oxj

ð9Þ

1246

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

and we get in the reference domain X: hAu; wi ¼ R Z  X ou  r 1 ow Gii0 ðlÞ^mi0 j0 Grjj0 ðlÞ detðGr ðlÞÞ dX; 8w 2 Y ; r oxi oxj r¼1 X ð10Þ hBp; wi ¼ R Z   ow X j p Grij ðlÞ detðGr ðlÞÞ1 dX; 8w 2 Y ;  r ox i r¼1 X ð11Þ hF ; wi ¼ hF s ; wi þ hF 0 ; wi; ð12Þ where hF s ; wi ¼

R Z   X 1 f^ r detðGr ðlÞÞ w dX; r¼1

Xr

^^g ; wi; hF ; wi ¼ hAL in

ð13Þ

^^g i: hG ; qi ¼ hBq; L in

ties and finite dimensional sub-spaces like Vh  V and Mh  M for velocity (uh) and pressure (ph), respectively. A necessary condition for the well posedness of this problem is the so-called inf–sup condition (LBB) [23]: 9b0 > 0:

ðT l q; wÞY ¼ Bðl; q; wÞ; so that T l q ¼ arg sup w2Y h

b2 ðlÞ ¼ inf

q2M h

ð14Þ

1 6 i; j 6 2; r ¼ 1; . . . ; R: ð15Þ

Furthermore, we may define H

qði;j;rÞ

ðlÞ ¼

Z

ou ow  dX; X oxi oxj Z owi hBsði;j;rÞ p; wi ¼  p dX r oxj X

mrij ðlÞ;

hA

qði;j;rÞ

u; wi ¼

ð16Þ

r

Usði;j;rÞ ðlÞ ¼ vrij ðlÞ;

ð17Þ

for 1 6 r 6 R, 1 6 i, j 6 2 (q and s are condensed indexes of i, j, r quantities) and we apply affine decomposition: Aðl; u; wÞ ¼

Qa X

q

q

H ðlÞAðu; wÞ ;

ð18Þ

Us ðlÞBðp; wÞs

ð19Þ

q¼1

Bðl; p; wÞ ¼

Qb X s¼1

for X  Rd¼2 , Qa may be as large as d · d · d · R and Qb as large as d · d · R. This splitting of the operators into a part which is parameter-dependent and into a part parameter-independent (defined and computed once in the reference domain) is crucial for the computational efficiency of the method. The Stokes problem rewritten on the reference domain X reads: find (u(l),p(l)) 2 Y · M such that  Aðl; uðlÞ; wÞ þ Bðl; pðlÞ; wÞ ¼ hF ; wi; 8w 2 Y ; Bðl; q; uðlÞÞ ¼ hG0 ; qi;

8w 2 Y h ;

Bðl; q; wÞ : kwkY

ð22Þ

ð23Þ

ðT l q; T l qÞY kqk2M

ð24Þ

;

by definition (22) it follows that kT l qk2Y ¼ Bðl; q; T l qÞ;

The tensors for pressure and divergence forms are 1

8l 2 Dl :

To verify it let us introduce a supremizer operator Tl: Mh ! Yh defined as follows:

mrij ðlÞ ¼ Grii0 ðlÞ^mi0 j0 Grjj0 ðlÞ detðGr ðlÞÞ1 ; 1 6 i; i0 ; j; j0 6 2; r ¼ 1; . . . ; R:

Bðl; q; wÞ P b0 ; kwkY kqkM

It is possible to demonstrate that

The transformation tensors for the bilinear viscous terms are defined as follows:

vrij ðlÞ ¼ Grij detðGr ðlÞÞ ;

q2M h w2Y h

ð21Þ

0

0

bðlÞ ¼ inf sup

8q 2 M: ð20Þ

In our numerical approximation the Stokes problem has been solved by the Galerkin-Finite Element Method, see [6–8]. With the subscript h we indicate discretized quanti-

and by applying (23) and (21) we have  2 ðT l q; T l qÞY kT l qk2Y ðT l q; wÞY ¼ inf ¼ inf sup inf 2 q2M h q2M h kqk2 q2M h w2Y kwk kqk kqkM h Y M M  2 Bðl; q; wÞ ¼ inf sup ¼ b2 ðlÞ; q2M h w2Y kwk kqk h Y M this proves (24). 3. The reduced basis formulation of the Stokes equations In the reduced basis approximation we take some ‘‘l’’ samples S lN ¼ fl1 ; . . . ; lN g, where ln 2 Dl , n = 1, . . . , N and we solve N times the Stokes problem (20) using the Galerkin-Finite Element Method to obtain uh(l) and ph(l). The reduced basis pressure space is MN = span{nn, n = 1, . . . , N}, where nn = ph(ln). We can build the reduced basis velocity space as Y lN ¼ spanffn ; n ¼ 1; . . . ; N ; T l nn ; n ¼ 1; . . . ; N g, where fn = uh(ln). The reduced basis approximation problem reads: find ðuN ðlÞ; pN ðlÞÞ 2 Y lN  M N such that (

Aðl; uN ðlÞ; wÞ þ Bðl; pN ðlÞ; wÞ ¼ hF ; wi; Bðl; q; uN ðlÞÞ ¼ hG; qi;

8w 2 Y lN ;

8q 2 M N : ð25Þ

Problem (25) is subject to an equivalent reduced basis inf–sup condition. Lemma 3.1. We define bN ðlÞ ¼ inf sup q2M N

l w2Y N

Bðl; q; wÞ : kwkY kqkM

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

8 2N N X X > > l > A u ðlÞ þ Blil pNl ðlÞ ¼ F li ; > Nj > < j¼1 ij l¼1

Then the following inequalities hold bN ðlÞ P bðlÞ P b0 > 0;

1247

8l 2 Dl :

2N > X > > > Bljl uNj ðlÞ ¼ Gll ; > :

Proof

1 6 i 6 2N ; ð27Þ

1 6 l 6 N;

j¼1

bðlÞ ¼ inf sup

q2M h w2Y h

Bðl; q; wÞ Bðl; q; wÞ 6 inf sup kwkY kqkM q2M N w2Y h kwkY kqkM

where the sub-matrices A and B are given by

Bðl; q; T l qÞ Bðl; q; wÞ ¼ inf 6 inf sup l l q2M N kT qk kqk q2M N w2Y kwkY kqkM Y M

Alij

¼

k¼1 b

where Qb ¼ Qb þ 1; UQ ¼ 1. For n = 1, . . . , N: rkn ¼ 0; for k ¼ 1; . . . ; Qb ; rQ b n ¼ fn ¼ uh ðln Þ: k

ðrkn ; wÞY ¼ BðnnN ; wÞ ;

8w 2 Y h ; for k ¼ 1; . . . ; Qb ; ð26Þ

rQ b n ¼ 0: In this way we have compacted the reduced basis velocity space made up of velocity (fn) and supremizer solutions (Tlnn). For a new ‘‘l’’ we want a solution given by a combination of previously computed stored solutions as basis functions: ! Qb 2N X X k uN ðlÞ ¼ uNj ðlÞ U ðlÞrkj ;

pN ðlÞ ¼

k

Blil

¼

Qb X Qb X

0

k

Uk ðlÞUk ðlÞBðrk0 i ; nl Þ ;

k¼1 k 0 ¼1

1 6 i 6 2N ; 1 6 l 6 N ; and F li

¼

Qb X k 0 ¼1 0

0

Uk ðlÞhF ; rk0 i i;

Gll ¼ hG ; nl i;

1 6 i 6 2N ;

1 6 l 6 N:

Finally, problem (27) can be written in compact form as       F uN A B ¼ : ð28Þ  T G 0 pN B This linear system, whose unknowns are the coefficients of the linear combination of previously computed off-line solutions, has the same structure of a finite element Stokes problem. Using reduced basis we deal with a matrix of considerably smaller dimension (of order of N) and with full matrices (instead of sparse ones). 3.1. Reduced basis on-line complexity

For n = N + 1, . . . , 2N:

N X

00

1 6 i; j 6 2N ;



To demonstrate the lemma above we have applied the fact that MN  Mh, the definition of the supremizer and the fact that the velocity space is enriched by supremizers, respectively. For further considerations on supremizer operator and the reduced basis framework, see [25] for the general non-coercive problem and [34] especially for Helmholtz and Burgers equations. We rewrite for computational convenience Y lN using the affine dependence of Bðl; q; wÞ on the parameter and the P b linearity of Tl: T l n ¼ Qq¼1 Uq ðlÞT q n for any n and l, which allows us to write ( b ) Q X l Y N ¼ span Uk ðlÞrkn ; n ¼ 1; . . . ; 2N ;

j¼1

0

Hk ðlÞUk ðlÞUk ðlÞAðrk0 i ; rk00 j Þ ;

k¼1 k 0 ¼1 k 00 ¼1

N

¼ bN ðlÞ:

Qa X Qb X Qb X

We have the following computational costs to build (online) reduced basis matrices, given also the supremizer com2 ponent in the velocity space: OðQa ðQb Þ 4N 2 Þ for sub-matrix 2 b b b A, OðQ Q 2N Þ for B, OðQ 2N Þ for F and O(27N3) for the inversion of the full reduced basis matrix (28). Note that the approach presented to build reduced basis velocity space is only one of the possible solutions and the space Y lN , in this case, is depending on the value of the on-line l parameter. Other options that avoid ill conditioning problems characterized by different computational costs are available and will be presented in Section 7. 3.2. An application to bypass configuration

k¼1

pNl ðlÞnl ;

l¼1

whose weights uNj and pNl are given by the following reduced basis linear system (in this case summation over i and j is no more understood to better explain how to build the on-line system):

After the general presentation of the reduced basis method for Stokes equations we introduce an application of interest used to carry out some preliminary numerical results described in the following sections and to introduce outputs of interest. In Fig. 1 we represent a ‘‘T-bypass’’ configuration, used, for example, to extract preliminary results about a possible optimized configuration for an arterial bypass and to study

1248

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

2

1.5

1.0

0.5

0

–0.5

–1 –1.5

–1

–0.5

0

0.5

1.0

1.5

b the four sub-domains and the problem parameters; C bw ¼ C b n ðC b in [ C b out Þ (left). Reference domain X (right). Fig. 1. Physical domain X:

the sensitivity of some geometrical parameters, see also [27]. We introduce a vector of parameters l ¼ ft; D; L; S; H ; hg 2 Dl  RP , P = 6, where Dl is the Cartesian product: [tmin, tmax] · [Dmin, Dmax] · [Lmin, Lmax] · [Smin, Smax] · [Hmin, Hmax] · [hmin, hmax], representing the range of variation of our parameters of interest during bypass optib and the reference one X mization. The physical domain X have been divided into four sub-domains (R = 4). The transformation tensors for the bilinear viscous terms (14): 1

mrij ðlÞ ¼ Grii0 ðlÞ^mi0 j0 Grjj0 ðlÞ detðGr ðlÞÞ ; 1 6 i; i0 ; j; j0 6 2; r ¼ 1; . . . ; R; assume this formulation: 2 3 2S t  tan h 6 H 7 6D m1 ¼ m 4 5 ; m2 ¼ m4 1 þ tan2 h 0  tan h H t 2t 3 2L 3 0 0 6D 7 6 7 m3 ¼ m4 ; m4 ¼ m4 D 5 5: D D 0 0 t L The tensors for pressure and divergence forms vrij ðlÞ

¼

Grij

are given  t 1 v ¼ 0  t v3 ¼ 0

1

r

detðG ðlÞÞ ;

by  H tan h H

 0 ; D

v4 ¼



v ¼ 

L 0

D S

3 7 5;

ð29Þ

S

0

0

D

 0 : D

As outputs of interest we may consider, for example, the mean value of the velocity components: Z u1 ðlÞ dX R X XrZ s1 ðlÞ ¼ ; r¼1 dX r Z X ð33Þ u ðlÞ dX R 2 X XrZ s2 ðlÞ ¼ : r¼1 dX Xr

ð30Þ

ð31Þ

or wall shear stress: Z ouðlÞ ^  t dC; ð35Þ ss ðlÞ ¼ m o^n Cw where ^t is the tangential unit vector. In the last two examples we may introduce a dual residual correction based on an adjoint problem to improve output accuracy (see, for example [29]).

ð32Þ

4. Off-line optimized basis assembling: adaptive procedure

(15):

 ;

3.3. Outputs of interest

Other outputs are derived from velocity, such as vorticity:  R Z  X ou2 ðlÞ ou1 ðlÞ sv ðlÞ ¼  dX; ð34Þ ox1 ox2 Xr r¼1

1 6 i; j 6 2; r ¼ 1; . . . ; R

2

;

0

nents we have Qa = 20 and Qb = 9. Please note that the tensors mr should be considered for each component (u1, u2)T of the velocity u.

In this case Qa and Qb of (18) and (19) are as large as 32 (d · d · d · R) and 16 (d · d · R), respectively, but due to the fact that the tensors mr and vr have some zero compo-

An adaptive procedure based on H1 maximum relative projected error for velocity EH 1 has been developed to optimize basis assembling, optionally we can also consider and

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

combine L2 maximum relative projected error for pressure EL2 . We underline that, given the higher powers of N that appear in our cost computing estimation, it is crucial (both for on-line and off-line effort) to control N more tightly. To this end, we may gainfully apply an off-line assembling procedure adaptively, see also [34,17]. We first construct, offline, an approximation that, over most of the domain, exhibits an error  (EH 1 or EL2 or both) less than prior : we d begin with a first point l1 (S N 0 ¼1 ¼ fl1 g); we next evaluate error N 0 ¼1 ðlÞ over a large test sample of parameter points Pprior in Dl , denoted with ; we then choose for l2 (and 1 2 0 0 hence Pprior S N ¼2 ¼ fl ; l g) the maximizer of N ¼1 ðlÞ over . We repeat P this process until the maximum of N 0 ¼N prior ðlÞ over prior is lower than prior . Then, on-line, d given a new value of the parameter, l, and an error tolerance post d ðlÞ, we essentially repeat this adaptive process – but now our sample points are drawn from S N prior , and the test sample is a singleton – l. Typically we choose prior  post d d ðlÞ since our test is not exhaustive; and therefore, typically, Npost(l)  Nprior. With the adaptive process we get higher accuracy at lower N: modest reductions in N can translate into measurable performance improvements. This procedure is very important not only to get a computationally cheaper and faster procedure but also to avoid ill conditioning in matrix assembling procedures. Error projection procedure, described below, has permitted us to have an off-line adaptive (and optimized) assembling procedure which is very fast and inexpensive, without reduced basis matrix assembling procedures, but based only on the solution of a linear system.

R where ð X ni nj dXÞ is the pressure mass matrix R P P b pij ¼ NP NP pi ð Ul Uk dXÞpj , 1 6 i, j 6 N, i.e., M l¼1 k¼1 l X k pTi M p pj . 4.2. H1 velocity error projection Given a new l and hence a new approximated velocity solution uh(l) we solve the following linear system  Z 2N Z X fi fj dX þ rfi  rfj dX cj X j ZX Z ¼ fi uh ðlÞ dX þ rfi  ruh ðlÞ dX; 1 6 i 6 2N ; X

X

ð37Þ where fi are 2N velocity solutions used as basis, made up of finite element approximated velocities and supremizer soluP tions, written as fi ¼ N u k¼1 ik Wk on N velocity nodes by Wk finite element velocity local functions. Expanding the system (37) by finite element approximation (on N nodes) we get, for 1 6 i 6 2N: " Z  2N X N X N X uil Wl Wk dX ujk j¼1

þ

l¼1

X

k¼1

2N X

N X

N X

j¼1

l¼1

k¼1

¼

N X

N X

m¼1 n¼1 N X N X m¼1 n¼1

X

X

wherePni are N pressure solutions used as basis given by N ni ¼ k¼1P pik Uk , computed solving a finite element problem on NP pressure nodes of the mesh; Uk are finite element local functions for the pressure. Expanding (36) by finite element approximation we can write Z  NP X NP N X X p il Ul Uk dX pjk cj l¼1

j¼1

¼

X

k¼1

NP X NP X m¼1 n¼1

Z pim

X

 Um Un dX phn ðlÞ;

Z  N X N X 2 ¼ kph ðlÞkL2  ci ni nj dX cj ; j¼1

X

 Wm Wn dX uhn ðlÞ X

Z uim

 rWm  rWn dX uhn ðlÞ;

X

ðuTi ðM þ KÞuj Þc ¼ uTi ðM þ KÞuh ðlÞ; b u c ¼ Fu : M Once we have c we write the H1 velocity error EH 1 as Z  2N X 2N X 2 E2H 1 ¼ kuh ðlÞkH 1  ci fi fj dX þ rfi  rfj dX cj ; i¼1

where ð PN PN

j¼1

X

R

fi fj dXÞ is the velocity mass matrix X R R T u l¼1 k¼1 il ð X Wl Wk dXÞujk , i.e., uj M ui and ð X rfi  R PN PN rfj dXÞ ¼ l¼1 k¼1 uil ð X rWl  rWk dXÞujk the velocity stiffness matrix, i.e., uTj K ui . 5. The computation of the constant b of the inf–sup condition

Once we have c we write the L2 pressure error EL2 as

i¼1

X

Z

1 6 i 6 N;

and using compact notation: ðpTi M p pj Þc ¼ pTi M p ph ðlÞ; b p c ¼ Fp : M

E2L2

uil

and in compact notation:

Given a new l and hence a new approximated ph(l) pressure solution we solve the following linear system: Z N Z X ni nj cj dX ¼ ni ph ðlÞ dX; 1 6 i 6 N ; ð36Þ j

 # rWl  rWk dX ujk cj

Z

uim

þ 4.1. L2 pressure error projection

1249

The calculation of the constant of the inf–sup condition for finite element (b) and reduced basis method (bN) has been carried out as a test to guarantee approximation stability: • b is related to the generalized eigenvalue problem (see [7,13,23]): Sz ¼ kM p z

ð38Þ

1250

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

Table 1 The constants b and bN for several values of the parameter D D

b

bN, N = 1

bN, N = 2

bN, N = 3

0.63 0.68 0.73 0.78 0.83 0.88 0.93 0.98 1.03 1.08

2.0249 2.0286 2.0274 2.0257 2.0237 2.0216 2.0192 2.0167 2.0139 2.011

3.9096 3.9316 3.952 3.971 3.9888 4.0054 4.0208 4.0353 4.0487 4.0612

2.7895 2.8708 2.9456 3.0145 3.0781 3.137 3.1918 3.243 3.291 3.3361

2.531 2.6625 2.7992 2.9386 3.0729 3.0729 3.186 3.264 3.3466 3.3759

6.1. First test: homogeneous Dirichlet boundary conditions, forced flow, three varying parameters

Other parameters are frozen.

where S = GTK1RG, K is the finite element (velocity) stiffness matrix ð X rWl  rWk dXÞ,R G is the velocitydivergence finite element matrix ð X Ul r RWk dXÞ, Mp is the finite element pressure mass matrix ð X Ul Uk dXÞ. Precisely, b is given by pffiffiffiffiffiffiffiffi b ¼ kmin ; ð39Þ S is symmetric and positive definite and its eigenvalues are real and positive. • To calculate bN we have a generalized eigenvalues problem on reduced basis (full) matrices (previously defined): b p zN S N z N ¼ kN M T

basis systems have been used: both iterative (such as Bi-CGSTAB, Conjugate Gradient Method, GMRES) and direct, see [22]. The reduced basis solutions have been compared directly with the true finite element solutions by computing H1 relative error for velocity and L2 relative error for pressure.

ð40Þ 1

with SN = B A B symmetric and positive definite, whose eigenvalues are real and positive. This time pffiffiffiffiffiffiffiffiffiffi bN ¼ kN min : ð41Þ • We have tested that the inequalities: bN ðlÞ P bðlÞ P b0 > 0; 8l 2 Dl are satisfied. Table 1 shows a test driven to calculate the constant b of the inf–sup condition studying a one-parameter varying configuration (the parameter D, for example) for different N.

6. Some preliminary numerical results Several numerical tests were carried out to develop all the Stokes reduced basis toolbox and to study different related aspects: affine mapping transformation in reference domains, off-line optimization for basis assembling procedure, approximation and algebraic stability (inf–sup condition fulfillment and condition number control). Different meshes have been used (from coarse to fine). Using finer mesh the CPU time (on a Pentium IV, 2 GHz) increases a lot during computational procedures such as off-line calculations, i.e., reduced basis approximation spaces and matrix assembling, error estimation and, eventually, orthonormalization, but it is maintained reasonable during on-line calculation. A mesh adaptation procedure in some zones has been implemented. Taylor–Hood P2  P1 elements (for velocity and pressure respectively) have been used [23]. Different solver for finite element and reduced

The first test we introduce deals with a forced Stokes flow in the domain considered with all zero Dirichlet boundary conditions (flow in a cavity), varying only three parameters available (depending by one non-dimensional quantity denoted with s). Data and relationships for geometrical parameters used to solve the problem follow (see also Fig. 1). (I) The varying parameters are t = 3  2s, S = s and L = s. (II) The viscosity is m = 0.04 m2 s1, the force field is ^f ¼ ð0; 10x1 ÞT m s2 in the physical domain X. b (III) The non-dimensional parameter s is ranging in [0.1, 1.45]. All other parameters are frozen: H = 1, D = 1 and the angle h = 0. Fig. 2 shows examples of flow solutions: velocity (absolute value) for two configurations with smin and smax; Fig. 3 shows the error reduction using the adaptive optimized procedure of Section 4 during basis assembling procedure, based on H1 relative (projected) error minimization, on the left, and the parameters distribution during basis assembling procedure, on the right. Fig. 4 shows true relative error reduction (H1 for velocity and L2 for pressure) considering a great number of different geometrical configurations. Table 2 shows error reduction (and its magnitude); Table 3 shows the proofs of approximation stability of reduced basis formulation computing bN, the inf–sup equivalent constant and its comparison with b from Galerkin (FEM) approximation. Fig. 5 shows the error R on a possible output of interest, for example sðlÞ ¼ X f  uðlÞ dX. 6.2. Second test: mixed Dirichlet/Neumann boundary conditions, five varying parameters The second test deals with a flow in the considered domain imposing Neumann homogeneous boundary conditions on inflow Cin and outflow Cout and homogeneous Dirichlet condition on the wall Cw. We deal with five varying parameters (bound by three relationships depending on two non-dimensional quantities s and q), only h is fixed to zero. Data and values used follow: • The parameters are t = 3  2s, S = s, L = s, D = q and H = 2  q.

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

1251

Fig. 2. Velocity solution (absolute value) for s = 0.1 (left) and s = 1.45 (right), m s1 · 102.

102

3

100

2.5

10–2

2

t,L,S

H1 error

S L t

10–4

1.5

10–6

1

10–8

0.5

10–10

0 2

4

6

8 N

10

12

14

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

N

Fig. 3. Left: maximum H1 relative projected error (see Section 4.2) on velocity at each iteration during adaptive basis assembling. Right: parameters distribution during off-line reduced basis assembling procedure (D and H are fixed).

• The viscosity is 0.04 m2 s1, while the force field is b f = (0, 10)T m s2 in true domain X. • The parameters range is s 2 [0.1, 1.45] and q 2 [0.1, 1.9]. Fig. 6 shows an example of flow solution (velocity absolute value) for a certain parameters combination (q = s = 1.0); Fig. 7 shows on the left the off-line parameters distribution during basis assembling procedures based on the reduction of the H1 relative projected error using the adaptive procedure; on the right we have maximum and mean output error Ds over a large number ofR different configurations, considering, for example, sðlÞ ¼ X f uðlÞ dX . Fig. 8 shows the true relative error reduction (H1 for velocity and L2 for pressure) considering a great number of different geometrical configurations. Table 4 shows error reduction (and its magnitude) and Table 5 the errors over the output of interest.

6.3. Third test: homogeneous Dirichlet boundary conditions, forced flow, six varying parameters The third test we carried out considered a forced Stokes flow in the domain imposing all zero Dirichlet boundary conditions, varying all the six parameters available (depending by two non-dimensional quantities: s, q), including the angle h. Data and relationship used in this test are reported below: • The parameters are t = 3  2s, S = s, L = s, D = q, H = 2  q and the angle is h. • The viscosity is 0.04 m2 s1, the force field is ^f ¼ ð0; 10x1 ÞT m s2 in the true domain X. b • The parameters range are s 2 [0.1, 1.45], q 2 [0.1, 1.9] and for the angle h 2 [0, p/3].

1252

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260 100

100

max mean

10–1

max mean

10–2

10–2 10–3

10–4 L2 error

H1 error

10–4 10–5 10–6

10–6

10–8 10–7 10–8 10

10–10

–9

10–10 0 10

101 N

102

10–12 100

101 N

102

Fig. 4. Left: H1 relative errors on velocity, maximum and mean values over a large test sampling (50 configurations). Right: relative L2 pressure errors, maximum and mean values over a large test sampling (50 configurations); errors in log–log scale.

Table 2 Table of H1 and L2 relative errors on velocity and pressure, respectively, 50 test configurations, three parameters varying N

Relative error H1 max

Relative error H1 mean

Relative error L2 max

Relative error L2 mean

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

9.8900e001 2.1583e001 2.3301e002 3.7543e003 3.2670e003 1.4271e004 8.4314e005 4.3142e005 2.4504e006 1.6485e006 1.8195e007 5.3852e008 2.4053e008 1.2923e009 1.0702e009

5.9138e001 7.8312e002 1.0174e002 8.7714e004 5.0383e004 3.6446e005 1.8559e005 7.4860e006 5.2102e007 2.5155e007 1.7635e008 6.9109e009 2.9963e009 2.5982e010 1.6946e010

2.5790e001 2.0387e002 5.9742e003 4.4999e004 3.2752e004 2.8591e005 2.7009e005 7.0025e007 2.2660e007 1.3400e007 8.8023e008 8.6392e010 1.5813e010 1.0129e010 8.7277e011

1.0878e001 1.1826e002 1.4583e003 5.5973e005 3.3635e005 1.7688e006 1.5313e006 1.6084e007 4.1354e008 2.8537e008 4.5140e009 1.5917e010 4.1002e011 1.4957e011 1.1840e011

Table 3 bN equivalent inf–sup constant N

bN

N

bN

N

bN

1 2 3 4 5

6.5012e+000 5.8811e+000 5.7978e+000 5.7876e+000 5.7344e+000

6 7 8 9 10

5.7051e+000 5.7048e+000 5.7047e+000 5.7056e+000 5.6309e+000

11 12 13 14 15

5.6263e+000 5.6245e+000 5.6087e+000 5.5935e+000 5.4651e+000

b = 5.0164 (s = 0.76).

Fig. 9 shows the projected error reduction (as shown in Section 4.2) using the adaptive procedure in basis assembling, based on H1 error minimization and parameters distribution during basis assembling procedure. Fig. 10 shows true error reduction (H1 for velocity and L2 for pressure) considering a great number of different geometrical configurations. The plateau in pressure error plot is due to the fact that we are optimizing reduced basis velocity approximation space with adaptive procedures, but nothing is

done with pressure. Table 6 shows error reduction (and its magnitude) Rand Table 7 the error on an output of interest: sðlÞ ¼ X f  uðlÞ dX. A second option is carried out using and adaptive off-line procedure based on total projected error (velocity and pressure). Fig. 11 shows projected error reduction using the adaptive procedure during basis assembling. Fig. 12 shows error reduction (H1 for velocity and L2 for pressure) considering different geometrical configurations. In this case the plateau in the pressure error behavior disappears because we are optimizing reduced basis velocity and pressure approximation spaces with adaptive procedures minimizing off-line error.

6.4. On computational costs In Table 8 we report some on-line computational costs (CPU time) versus N. The finite element solution is obtained with a mean CPU time of 112.23, adopting a

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260 101

1025.818

100

1025.816

10–1

1253

1025.814

10–2 S–SN

1025.812 10–3 1025.81 10–4 1025.808 10–5 1025.806

10–6

1025.804

10–7 4

6

8

10

12

6

14

7

8

9

10

11

12

13

14

15

N

Fig. 5. Left: error on functional output Ds = (s  sN); s ¼

R X

f  u dX. Right: convergence of sN to s (*) versus N (s = 0.356).

mesh with 1278 elements, see Fig. 1. Note the computational savings in the on-line step once we have assembled reduced basis matrix. In the same table we have inserted H1 error indication to compare the computational costs and precision reached. The example refers to the case of a single varying parameter D.

7. On algebraic and approximation stability To control the condition number of reduced basis matrix we have adopted an orthonormalization procedure applied to velocity and pressure basis functions. After orthonormalization (to achieve algebraic stability) we have to satisfy the approximation stability condition on bN, the equivalent reduced basis inf–sup constant. But if we apply the orthonormalization algorithm to reduced basis approximation spaces, assembled as proposed in Section 3, we may lose the validity of Lemma 3.1 to guarantee the

Fig. 6. Example of flow solution (velocity absolute value, m s1 · 102) in the parametrized domain.

101

3

max mean

100

2.5

2 S,L,t,H,D

S–SN

10–1

10–2

1.5

10–3

1

10–4

0.5

10–5 100

D H L S t

101 N

102

0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 N

R Fig. 7. Left: functional output error Ds = (s  sN); s ¼ X f  u dX: maximum and mean error over large test sample configurations. Right: parameters distribution during off-line reduced basis optimized assembling procedure.

1254

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260 101

max mean

100

100

10–1

10–1 L2 error

H1 error

101

10–2

max mean

10–2

10–3

10–3

10–4

10– 4

10–5

100

101

10–5

102

100

101

102

N

N

Fig. 8. Left: relative H1 velocity errors: maximum and mean values over a large test sampling (50 configurations). Right: relative L2 pressure errors: maximum and mean values over a large test sampling (50 configurations); errors in log–log scale.

Table 4 Table of H1 relative errors on velocity and L2 relative errors on pressure, 50 test configurations, five varying parameters, N 6 20 N

Relative error H1 max

Relative error H1 mean

Relative error L2 max

Relative error L2 mean

5 10 15 20

2.93e001 1.33e002 2.98e003 4.61e004

7.86e002 4.47e003 6.65e004 8.01e005

3.09e002 6.79e003 8.16e004 6.24e004

6.14e003 6.96e004 9.021e005 1.60e005

Table 5 Ds = s  sN max and mean error over output s for 50 configurations

7.1. Orthonormalization: Gram–Schmidt (GS) algorithm

N

Ds max

N

Ds mean

5 10 15 20

9.62e002 2.49e002 2.09e003 7.68e004

5 10 15 20

4.73e002 5.27e003 5.63e004 5.35e005

We recall very briefly the main step of GS orthonormalization:

stability of the approximation. For this reason we are going to propose other options in building the reduced basis velocity space.

• Given: zj, j = 1, . . . , N a family vector of functions; P z • we obtain lj ¼ kP jj zjj k, where Pj is an orthogonal projector onto the orthogonal complement of hl1, l2, . . . , lj1i (i.e., the space generated by lv for v = 1, . . . , j  1). • Each lj is orthogonal to l1, l2, . . . , lj1 and lies in hz1, z2, . . . , zji.

101

3

D L S t H θ

2.5 100

Parameters

H1 max error

2

10–1

1.5

1 10

–2

0.5

10–3 100

101 N

102

0

0 1 2 3 4 5 6 7 8 9 10111213141516171819202122232425 N

Fig. 9. Left: H1 velocity (projected) error reduction during basis assembling. Right: parameters distribution during off-line reduced basis assembling.

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

1255

101

101

mean max

mean max

100 100

L2 error

H1 error

10–1 10–1

10–2

10–2 10–3

10–3

10–4

10–4

100

101 N

102

10–5 100

101

102

N

Fig. 10. Left: relative H1 velocity errors: maximum and mean values over a large test sampling (90 configurations.). Right: relative L2 pressure errors: maximum and mean error over a large test sampling (90 configurations); errors in log–log scale.

Table 6 Table of velocity relative errors H1 and pressure L2, 90 test configurations, 6(3) parameters, N 6 25 N

Relative error H1 max

Relative error H1 mean

Relative error L2 max

Relative error L2 mean

5 10 15 20 25

1.60e001 1.01e001 1.348e002 7.82e003 4.79e003

5.65e002 1.60e002 4.14e003 1.49e003 9.37e004

9.55e001 9.48e001 7.02e003 5.64e003 3.63e003

7.08e002 5.74e002 2.17e004 1.20e004 6.66e005

Table 7 Maximum and mean error on functional output Ds = (s  sN): R s ¼ X f  u dX for 90 configurations N

Ds max

N

Ds mean

5 10 15 20 25

1.34e001 1.26e001 4.064e003 2.50e003 1.38e003

5 10 15 20 25

3.26e002 6.91e003 6.87e004 2.92e004 2.31e004

ðrn ; vÞH 1 ¼

Qb X

q

U ðlÞðrqn ; vÞH 1 ¼

q¼1

Qb X

q

Uq ðlÞBðnnN ; vÞ ;

q¼1 2

8v 2 ðH 1CD ðXÞÞ ; in fact ðrn ; vÞH 1 ¼ Bðl; nnN ; vÞ;

ð42Þ

but we recall that for each q component it holds • • • •

Each kljk = 1. P j ¼ I  Lj1 LTj1 ; Lj1 = {l1, . . . , lj1}; P j zj ¼ zj  ðlT1 lj Þl1  ðlT2 zj Þl2  . . .  ðlTj1 zj1 Þlj1 ; {l1, . . . , lN} is an orthogonal basis of hz1, z2, . . . , zNi.

The norm kÆk used here is the Y = (H1(X))2 for velocity (and supremizers) and L2(X) for pressure. The scalar product lTi zj is the one induced by the functional space and the norm we use. The orthonormalization procedure has been applied to reduced basis functions. For velocity and pressure the procedure is standard. For the supremizer we have to introduce some considerations. Referring to previous supremizer formulation of (26) and to the compact notation already introduced in the reference domain, thanks to the linearity of supremizer operator and to the affine composition property we have for n = N + 1, . . . , 2N:

ðrqn ; vÞH 1 ¼ BðnnN ; vÞq

ð43Þ

and by the affine composition/decomposition of the operator we write Bðl; n; vÞ ¼

Qb X

q

Uq ðlÞBðn; vÞ :

q¼1

At this point we have two possibilities (referring to nth supremizer rn, n = N + 1, . . . , 2N) in applying orthonormalization: (i) an orthonormalization (GS) directly on rn done on-line (being rn dependent on l) to obtain r? n as new element (basis function) to enrich the reduced basis velocity space:

1256

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260 101

3

D theta H L S t

2.5 100

Parameters

Total Error

2

10–1

1.5

1 10–2 0.5

10–3 100

101

0

102

0

2

4

6

N

8

10 N

12

14

16

18

20

Fig. 11. Left: H1 velocity and L2 pressure (projected) errors reduction during basis assembling. Right: corresponding parameters distribution during offline reduced basis assembling.

101

101

max mean

max mean

100

100

10–1

L2 error

H1 error

10–1 10–2 10–3

10–2

10–4 10–3 10–5 10–4 100

101

102

10–6 0 10

101 N

N

102

1

Fig. 12. Left: relative H velocity errors: maximum and mean error over a large test sampling (90 configurations). Right: relative L2 pressure errors: maximum and mean error over a large test sampling (90 configurations); errors in log–log scale.

Table 8 On-line computational costs (CPU time) varying N (and H1 velocity error) and comparison with computational cost of FEM solution (%) N

CPU time

H1 error

%

N

CPU time

H1 error

%

1 3 5 7 9 11 13 15 17 19

1.597 3.225 5.908 6.68 8.442 10.615 14.38 20.601 24.625 24.025

6.24e1 6.58e3 1.92e4 5.27e5 6.83e7 2.66e8 1.25e8 1.38e12 2.20e13 2.19e14

1.4 2.9 5.2 5.9 7.52 9.45 12.81 18.4 21.9 21.4

2 4 6 8 10 12 14 16 18 20

2.445 4.937 5.908 7.55 9.314 14.11 17.895 24.646 23.614 25.257

2.16e1 5.30e4 1.006e4 3.36e5 1.11e7 1.98e8 2.91e12 4.05e13 2.74e14 1.39e15

2.2 4.4 5.2 6.7 8.3 12.57 15.9 21.9 21.1 22.5

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

P? n rn ; kP ? n rn k X b  Q P? U ðlÞr q qn n q¼1 X b  ; r? n ¼ Q ? kP n Uq ðlÞrqn k q¼1

mizer on its component rkn (before summation) to preserve Lemma 3.1. To achieve this goal we may introduce two further different options in assembling the supremizer solutions for stabilization procedure by building in a different way the reduced basis velocity approximation space so that we may guarantee both approximation and algebraic stability.

r? n ¼

T P? n ¼ I  Ln1 Ln1 ;

? Li ¼ fr? 1 ; . . . ; ri g;

(ii) an orthonormalization (GS) on components rqn made off-line (rqn are not depending on l) to get r? qn : r? qn

¼

1257

7.2.1. First option We have considered the following reduced basis spaces: Y N ¼ spanfrn ; n ¼ 1; . . . ; N Qb g;

P? qn rqn kP ? qn rqn k

;

where Qb ¼ Qb þ 1. For n = 1, . . . , N:

T P? qn ¼ I  Lqðn1Þ Lqðn1Þ ;

? Lqi ¼ fr? q1 ; . . . ; rqi g;

rn ¼ fn ¼ uðln Þ:

Fig. 13 shows the reduction of the condition number of the matrix of the reduced basis linear system by orthonormalizing reduced basis functions (using method (i) for supremizer, including also orthonormalization for velocity and pressure basis functions). A very interesting property (visible also in Fig. 13) is that the condition number is limited and bounded after orthonormalization and not depending on N. 7.2. Approximation stability: other supremizer options By orthonormalizing the supremizer solutions according to the approach (i) we could lose approximation stability (guaranteed by Lemma 3.1) in the attempt of preserving algebraic stability by reducing the condition number. We can orthonormalize just using method (i) the pressure n and the velocity f basis functions and not the supremizer rn and use the approach (ii) to orthogonalize the supre-

1018

For n ¼ N þ 1; . . . ; N Qb , condensing index m and k in n, we have ðrn ; wÞY ¼ ð~ rmk ; wÞY , where ð~ rmk ; wÞY ¼ Bðnm ; wÞk ; uN ðlÞ ¼

N Qb X

8w 2 Y ;

k ¼ 1; . . . ; Qb ; m ¼ 1; . . . ; N ;

uNj ðlÞrj ;

j¼1

pN ðlÞ ¼

N X

pNl ðlÞnl ;

l¼1

the reduced basis system becomes 8 N Qb N X X > > > Alij uNj ðlÞ þ Blil pNl ðlÞ ¼ F i ; > > < j¼1 l¼1 > N Qb > X > > > Bljl uNj ðlÞ ¼ Gl ; :

1 6 i 6 N Qb ;

1 6 l 6 N;

j¼1

ð44Þ

2

x 1019

no orthog. GS orthog.

no ortho. GS ortho

1.8

16

10

1.6 1014 1.4

cond.num.

cond.num.

1012 10

10

108

1.2 1 0.8 0.6

106 0.4 104

0.2

102

0 4

6

8

10 N

12

14

12

13

14

15

16

17

18

19

20

N

Fig. 13. Orthonormalization: bounded condition number of reduced basis Stokes linear system matrix with complete orthonormalized basis (left) and with a partial orthonormalization only on velocity and pressure (right), but not on supremizer. In each picture there is a comparison with condition number without any orthonormalization.

1258

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

where Alij

¼

For n = N + 1, . . . ,2N: Qa X

k

k

k

H ðlÞAðri ; rj Þ ;

1 6 i; j 6 N Qb ;

rQ b n ¼ 0:

k¼1

Blil ¼

Qb X

k

Uk ðlÞBðri ; nl Þ ;

1 6 i 6 N Qb ; 1 6 l 6 N ;

k¼1

F i ¼ hF ; ri i;

1 6 i 6 N Qb ; Gl ¼ hG0 ; nl i; 1 6 l 6 N :

In this case the basis is no longer l (on-line) dependent, the reduced basis velocity space, enriched by supremizers, has a bigger dimension (N Qb > 2N ) than previously. The compu2 tational costs are as follows: OðQa ðQb Þ N 2 Þ for sub-matrix b 2 A, OðQb Q N Þ for B, OðQb N Þ for F, but the cost for inversion of the full reduced basis matrix (28) increases now to 3 OððQb þ 1Þ N 3 Þ. This approach has the big advantage to preserve Lemma 3.1, to let us apply orthonormalization (method (ii)) and to preserve the mentioned lemma also after orthonormalization. In conclusion this approach is the best if we want to be sure to preserve both algebraic and approximation stability after orthonormalization, if necessary. Fig. 14 shows errors behavior by testing this first new supremizer option over a large test sampling (zero Dirichlet conditions and three geometrical parameters). 7.2.2. Second option Another approach is based on the idea that supremizers are built upon summation using the same lj values used to store velocity fj(lj) and pressure solutions nj(lj) (also in this case the basis for velocity is not dependent on the on-line value of l and it is completely assembled off-line): ( ) Qb X k n Y N ¼ span rn ¼ U ðl Þrkn ; n ¼ 1; . . . ; 2N ;

The reduced basis solution is given by ! Qb 2N X X uNj ðlÞ Uk ðlj Þrkj ; uN ðlÞ ¼ k¼1

j¼1

pN ðlÞ ¼

N X

pNl ðlÞnl ;

l¼1

by solving the system: 8 2N N X X > > > Alij uNj ðlÞ þ Blil pNl ðlÞ ¼ F i ; > > < j¼1 l¼1 2N > X > > > Bljl uNj ðlÞ ¼ Gl ; > :

Alij

¼

Qa X

1 6 l 6 N;

k

Hk ðlÞAðri ; rj Þ

k¼1

¼

Qa X Qb X Qb X

0

00

Hk ðlÞUk ðli ÞUk ðlj ÞAðrk0 i ; rk00 j Þk ;

k¼1 k 0 ¼1 k 00 ¼1

1 6 i; j 6 2N ; Blil ¼

Qb X

Uk ðlÞBðri ; nl Þ

k

k¼1

¼

Qb X Qb X

0

k

Uk ðlÞUk ðli ÞBðrk0 i ; nl Þ ;

0

k¼1 k ¼1

1 6 i 6 2N ; 1 6 l 6 N ; F i ¼ hF ; ri i ¼

for k ¼ 1; . . . ; Qb ; rQ b n ¼ fn ¼ uðln Þ:

Gl ¼ hG0 ; nl i;

Qb X

0

Uk ðli ÞhF ; rk0 i i;

1 6 i 6 2N ;

k 0 ¼1

101

1 6 i; l 6 N :

100 mean max

mean max

100

10–1

10–1

10–2 L2 error

H1 Error

rkn ¼ 0;

ð45Þ

where

¼ 1. For n = 1, . . . , N:

where Qb ¼ Qb þ 1; U

1 6 i 6 2N ;

j¼1

k¼1 Qb

8w 2 Y ; for k ¼ 1; . . . ; Qb ;

ðrkn ; wÞY ¼ BðnnN ; wÞ ;

10–2

10–3

10–3

10–4

10–4

10–5

10–5 100

101

10–6 100

N

101 N

1

2

Fig. 14. Supremizer: first option, H relative error (for velocity) and L relative error (for pressure) using 50 configurations.

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

This option is also competitive as regards the computational costs dealing with 3N · 3N reduced basis matrices (28) instead of ðQb þ 1ÞN  ðQb þ 1ÞN matrix (usually ðQb þ 1Þ 3). We have the following computational costs to build reduced basis matrices, given also the supremizer components in the velocity space: O(Qa4N2) for sub-matrix A, O(Qb2N2) for B, O(N) for F and O(27N3) for the inversion of the full reduced basis matrix (28). Using this option we cannot demonstrate that Lemma 3.1 is still valid (even without orthonormalization). We have tested numerically this option and we can argue that also this approximation is reasonably stable. Numerical results are shown in Fig. 15 where we have reported errors behavior always testing a large number of configurations (Dirichlet conditions and three varying parameters as in Section 6.3). Convergence is very fast in this case. Fig. 16 shows the condition number reduction of reduced basis linear system and a comparison of condition number of the matrix without orthonormalization.

8. Perspectives Development guidelines are devoted to the application of reduced basis to Navier–Stokes equations in parametrized domains and in problems involving non-affine mapping dependence (i.e., shape design problem). See [31] for Navier–Stokes reduced basis formulation and [30,2] for non-affine parametric dependence in elliptic problems. The latter reference especially introduces an efficient reduced basis discretization procedure replacing non-affine coefficient functions with a collateral reduced basis expansion which permits an off-line/on-line computational decomposition by a stable and inexpensive interpolation procedure. This approach allows us the introduction of curved elements in parametrized geometries. The interest is to expand and apply efficient, accurate and real-time reduced basis techniques on bio-mechanics problems (i.e., bio-medical devices optimization such as bypass [21,26,27]). The most important steps under investigation

100

100 mean max

10–2

10–2

10–4

10–4 L2 error

H1 Error

mean max

10–6

10–6

10–8

10–8

10–10

10–10

10–12 0 10

101 N

1259

10–12 100

102

101 N

102

Fig. 15. Supremizer: second option, H1 relative error (velocity) and L2 relative error (pressure) over 50 configurations.

1018

1200

mean max

mean max

1016 1000 1014

cond. numb.

cond.numb.

800

600

1012 1010 108

400 106 200

0 0

104 10 5

10

2

0

5

10

15

15

N

Fig. 16. Condition number (max and mean) for the Stokes reduced basis system, with the new supremizer second option. On the left results with orthonormalization, on the right without it.

1260

G. Rozza, K. Veroy / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1244–1260

are (i) the use of a great number of geometrical parameters, (ii) the non-affine formulation for Stokes reduced basis problem ([28]) and (iii) the study of a posteriori error bounds for Stokes equations. Acknowledgements We acknowledge Prof. A.T. Patera of MIT, Prof. A. Quarteroni of EPF-Lausanne and MOX-Politecnico di Milano and Prof. Y. Maday of University Paris VI, Dr. D. Rovas (University of Illinois), Dr. C. Prud’homme (EPFL) and Dr. N.C. Nguyen (National University of Singapore) for suggestions, very helpful discussion, comments and insights, Prof. F. Saleri (MOX-Politecnico di Milano) for providing the original finite element Stokes solver. G. Rozza acknowledges the support provided through the European Community’s Human Potential Programme under contract HPRN-CT-2002-00270 HaeMOdel. This work was supported also by Swiss National Science Foundation (PBEL2-111646), by DARPA and ASFOR under grant F4920-03-1-0356, by DARPA and GEAE under grant F49620-03-1-0439 and by Singapore–MIT Alliance (SMA). References [1] E. Balmes, Parametric families of reduced finite element models: theory and applications, Mech. Syst. Signal Process. 10 (4) (1996) 381–394. [2] M. Barrault, Y. Maday, N.C. Nguyen, A.T. Patera, An ‘‘empirical interpolation’’ method: application to efficient reduced-basis discretization of partial differential equations, C.R. Acad. Sci. Paris, Ser. I 339 (2004) 667–672. [3] A. Barrett, G. Reddien, On the reduced basis method, Z. Angew. Math. Mech. 75 (7) (1995) 543–549. [4] J.P. Fink, W.C. Rheinboldt, On the error behavior of the reduced basis technique for nonlinear finite element approximations, Z. Angew. Math. Mech. 63 (1) (1983) 21–28. [5] G.P. Galdi, An introduction to the mathematical theory of the Navier–Stokes equations, Linearized Steady Problem, vol. I, Springer-Verlag, New York, 1994. [6] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations, Springer-Verlag, Berlin, 1986. [7] P.M. Gresho, R.L. Sani, Incompressible Flow and the Finite Elements Method, Wiley, New York, 2000. [8] M.D. Gunzburger, Finite Element Method for Viscous Incompressible Flows: A Guide to Theory, Practice, and Algorithms, Academic Press, Boston, 1989. [9] K. Ito, S.S. Ravindran, A reduced basis method for control problems governed by PDEs, in: W. Desch, F. Kappel, K. Kunish (Eds.), Control and Estimation of Distributed Parameter System, Birka·user, Basel, 1998, pp. 153–168. [10] K. Ito, S.S. Ravindran, A reduced-order method for simulation and control of fluid flow, J. Comput. Phys. 143 (2) (1998) 403– 425. [11] Y. Maday, A.T. Patera, G. Turinici, Global a priori convergence theory for reduced-basis approximations of single-parameter symmetric coercive elliptic partial differential equations, C.R. Acad. Sci. Paris, Se´r. I 335 (2002) 1–6. [12] Y. Maday, A.T. Patera, G. Turinici, A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations, J. Sci. Comput. 17 (1) (2002) 437–446.

[13] D.S. Malkus, Eigenproblems associated with the discrete LBB condition for incompressible finite elements, Int. J. Engrg. Sci. 19 (1981) 1299–1310. [14] N.C. Nguyen, K. Veroy, A.T. Patera, Certified real-time solution of parametrized partial differential equations, in: R. Catlow, H. Shercliff, S. Yip (Eds.), Handbook of Materials Modeling, Kluwer Academic Publishing (Springer), Dordrecht, 2005. [15] J.S. Peterson, The reduced basis method for incompressible viscous flow calculations, SIAM J. Sci. Stat. Comput. 10 (4) (1989) 777–786. [16] T.A. Porsching, Estimation of the error in the reduced basis method solution of nonlinear equations, Math. Comput. 45 (172) (1985) 487–496. [17] C. Prud’homme, Adaptive reduced basis space generation and approximation, submitted for publication. [18] C. Prud’homme, D. Rovas, K. Veroy, Y. Maday, A.T. Patera, G. Turinici, Reliable real-time solution of parametrized partial differential equations: reduced-basis output bound methods, J. Fluids Engrg. 172 (2002) 70–80. [19] C. Prud’homme, A.T. Patera, Reduced-basis output bounds for approximately parametrized elliptic coercive partial differential equations, Comput. Visual. Sci. 6 (2–3) (2004) 147–162. [20] C. Prud’homme, D. Rovas, K. Veroy, A.T. Patera, Mathematical and computational framework for reliable real-time solution of parametrized partial differential equations, M2AN 36 (5) (2002) 747–771. [21] A. Quarteroni, G. Rozza, Optimal control and shape optimization of aorto-coronaric bypass anastomoses, M3AS 13 (12) (2003) 1801– 1823. [22] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, Springer, New York, 2000. [23] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, Berlin, 1994. [24] W.C. Rheinboldt, On the theory and error estimation of the reduced basis method for multi-parameter problems, Nonlin. Anal. Theory, Methods Appl. 21 (11) (1993) 849–858. [25] D. Rovas, Reduced-basis output bound methods for parametrized partial differential equations, PhD Thesis, MIT – Massachusetts Institute of Technology, February 2003. [26] G. Rozza, On optimization, control and shape design of an arterial bypass, Int. J. Numer. Methods Fluids 47 (10–11) (2005) 1411–1419. [27] G. Rozza, Real time reduced basis techniques for arterial bypass geometries, in: K.J. Bathe (Ed.), Computational Fluid and Solid Mechanics, Elsevier, Amsterdam, 2005, pp. 1283–1287. [28] G. Rozza, Reduced basis methods for Stokes equations in domains with non-affine parameter dependence, EPFL-IACS report 06.2005, Comput. Visual. Sci., submitted for publication. [29] G. Rozza, Shape design by optimal flow control and reduced basis techniques: applications to bypass configurations in haemodynamics, PhD Thesis No. 3400, EPFL, Ecole Polytechnique Fe´de´rale de Lausanne, 2005. [30] Y. Solodukhov, Reduced-basis methods applied to locally non-affine problems, PhD thesis, MIT, Massachusetts Institute of Technology, 2004. [31] K. Veroy, A.T. Patera, Certified real-time solution of the parametrized steady incompressible Navier–Stokes equations: rigorous reduced-basis a posteriori error bounds, Int. J. Numer. Methods Fluids 47 (8–9) (2005) 773–788. [32] K. Veroy, A.T. Patera, Reduced-basis approximation of the viscosityparametrized incompressible Navier–Stokes equations: rigorous a posteriori error bounds, in: Proceedings Singapore–MIT Alliance Symposium, January 2004. [33] K. Veroy, Reduced-basis methods applied to problems in elasticity, PhD Thesis, MIT – Massachusetts Institute of Technology, June 2003. [34] K. Veroy, C. Prud’homme, D. Rovas, A.T. Patera, A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations, AIAA American Institute of Aeronautics and Astronautics, Paper 2003–3847, 2003.