On the verification of model reduction methods based on the proper generalized decomposition

Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

On the verification of model reduction methods based on the proper generalized decomposition Pierre Ladevèze a,b,⇑, Ludovic Chamoin a a b

LMT-Cachan (ENS Cachan/CNRS/Paris 6 Univ./PRES UniverSud Paris) 61 Avenue du Président Wilson, 94235 CACHAN Cedex, France EADS Foundation Chair ‘‘Advanced Computational Structural Mechanics’’, France

a r t i c l e

i n f o

Article history: Received 30 September 2010 Received in revised form 21 February 2011 Accepted 28 February 2011 Available online 4 March 2011 Keywords: Model reduction Verification Error estimation Proper generalized decomposition Separated representation Constitutive relation error

a b s t r a c t In this work, we introduce a consistent error estimator for numerical simulations performed by means of the proper generalized decomposition (PGD) approximation. This estimator, which is based on the constitutive relation error, enables to capture all error sources (i.e. those coming from space and time numerical discretizations, from the truncation of the PGD decomposition, etc.) and leads to guaranteed bounds on the exact error. The specificity of the associated method is a double approach, i.e. a kinematic approach and a unusual static approach, for solving the parameterized problem by means of PGD. This last approach makes straightforward the computation of a statically admissible solution, which is necessary for robust error estimation. An attractive feature of the error estimator we set up is that it is obtained by means of classical procedures available in finite element codes; it thus represents a practical and relevant tool for driving algorithms carried out in PGD, being possibly used as a stopping criterion or as an adaptation indicator. Numerical experiments on transient thermal problems illustrate the performances of the proposed method for global error estimation. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Nowadays, numerical simulation constitutes a common tool in science and engineering activities. It is especially used for prediction and decision making, or simply for a better understanding of physical phenomena. However, in order to give an accurate representation of the real world, a large set of parameters may need to be introduced in the mathematical models involved in the simulation, which leads to important (and often overwhelming) computational efforts. This is for example the case when dealing with models that aim at taking uncertainties into account, or in optimization problems. A drawback of such complex multi-parameter models is the fact that they usually lead to a huge number of degrees of freedom (due to the so-called curse of dimensionality) so that they cannot be tackled with classical brute force methods. Therefore, alternative numerical approaches are necessary. For that purpose, model reduction methods exploit the fact that the response of complex models can often be approximated with a reasonable precision by the response of a surrogate model, seen as the projection of the initial model on a low-dimensional reduced basis [19,4,2,15,40,12,39,16,3,41]. Model reduction methods distinguish ⇑ Corresponding author at: LMT-Cachan (ENS Cachan/CNRS/Paris 6 Univ./PRES UniverSud Paris) 61 Avenue du Président Wilson, 94235 CACHAN Cedex, France. E-mail addresses: [email protected] (P. Ladevèze), chamoin@lmt. ens-cachan.fr (L. Chamoin). 0045-7825/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.02.019

themselves by the way of defining and constructing the reduced functional basis. During the last few years, an appealing model reduction technique, based on separation of variables within a spectral resolution approach, has received a growing interest and has been the object of intensive research works. Even though it is not a new approach, and is not restricted to Computational Mechanics (see [5] for an application in quantum chemistry), it is particularly well-suited for the resolution of highly multidimensional problems currently considered. Indeed, the associated complexity scales linearly with the number of dimensions, instead of the exponential growing characteristic of mesh-based discretization strategies. This model reduction strategy implies the solution of a few tractable problems in order to get a correct approximation of the exact solution. The proper orthogonal decomposition (POD) [20,12] is a wellknown technique for separation of variables which has been widely used in various scientific domains. It is related to other decomposition techniques such as the Karhunen–Loeve expansion [18] or the singular value decomposition (SVD) [14]. All these POD-like techniques are referred as a posteriori as they require some knowledge (at least partial) on the solution of the problem. Even though POD was preliminary used in the context of data analysis (the objective being a low-dimensional approximate description of a high-dimensional process, see [6]) as well as in image processing applications (the objective being to extract the most relevant part of an image in order to compress the data), it is now used as a model reduction

2033

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

technique [19,5,40,2]. The associated procedure implies the resolution, in what is called a learning phase or snapshot, of a problem which is less costly than the original problem, and whose solution is decomposed using POD. Then, the initial model is reduced through its projection over the resulting basis, leading to an inexpensive resolution. However, the quality of the obtained solution is highly dependent on the choices made during the learning phase, which is one of the main drawbacks of the method. A different technique for separation of variables has been introduced as radial loading approximation in [21]. It is nowadays named proper generalized decomposition (PGD) [10] as it can be seen as a POD extension. The PGD approximation does not require any knowledge on the solution (it is thus referred as a priori) and does not use any orthogonality property. It operates in an iterative strategy in which a set of simple problems, that can be seen as pseudo eigenvalue problems, need to be solved. PGD has been developed during years for solving time-dependent nonlinear problems in Structural Mechanics [21,22,35,31,37] as a keypoint of the socalled LATIN method. Extension to stochastic problems, initially under the name Generalized Spectral Decomposition, has been done in [32,33]. The PGD approach has also lead to major breakthroughs in multidimensional or inverse problems [9,11]. The PGD approximation drastically reduces computation and storage cost compared to the brute force approximation. For instance, considering a standard transient model defined in a 3D physical space, and discretizing the time domain with N time steps (N can be millions), usual incremental strategies require the solution of N (usually nonlinear) 3D space problems. However, if the space–time PGD (or radial approximation) is considered, about k  m 3D problems should be solved for computing the space functions and k  m 1D problems for computing the time functions (m being the number of modes considered in the PGD approximation, and k the number of iterations). As k  m  100 in practical applications, the computing time savings can reach many orders of magnitude. Even though PGD is usually very effective, a major drawback is that it does not include, until now, any robust error estimator that could give an idea of the quality of the approximate solution. Basic results on a priori error estimation for representations using separation of variables can be found in [22]. A first attempt to compute a posteriori error estimates was also presented in [1] when solving multidimensional models with PGD representations; an error estimator, given with respect to a specific output of interest, was computed at each iteration to control both the quality of added terms in the decomposition and the convergence of the approximate solution. However, the method is not robust enough as it merely defines a goal-oriented error indicator without providing for strict bounds. Furthermore, error on the adjoint problem is not clearly analyzed as a very fine solution of the adjoint problem is taken, which is usually not available in practical applications. In this paper, the goal is not to construct a new PGD approximation, but to verify in a robust manner computations performed with PGD. For that, we introduce an error estimation procedure based on the constitutive relation error [27] which leads to a guaranteed and relevant error evaluation. What is specific here is the way to construct the required admissible fields, and the procedure used to split contributions of the various error sources; the methodology lies on a double standard kinematic PGD computation, associated with kinematic and static formulations, which can be easily handled using commercial finite element codes. The obtained data, as well as the prescribed data on displacements, tractions, and body forces define the starting point of the error estimation method. The resulting error estimator takes all error sources into account, particularly that related to classical discretizations (use of the finite element method (FEM), time integration schemes) and that related to the truncation of the decomposition sum. Therefore, this error estimator and inherent error indicators

provide useful information for stopping criteria and adaptive procedures, which leads to an optimal approach without ever doing the full computation. The methodology for error estimation is applied here to the case of linear problems with radial loading approximation, taking a transient diffusion problem as an example. Let us finally note that the approach we propose can be generalized to a broader set of problems. First, it is possible to consider other mechanical models, provided operators defining constitutive relations present sufficient convexity; it thus includes (visco-)elasticity, (visco-)plasticity, and transient visco-dynamics models for which guaranteed error bounds (based on admissible fields) can be computed [27,23]. Moreover, it is not limited to PGD but applies to the class of model reduction methods that can provide an equilibrated (in the finite element sense) solution. The paper is organized as follows: after this introduction, the reference diffusion problem is presented in Section 2; Section 3 gives a general review on decomposition of the solution using PGD; the global error estimation method we propose is introduced in Section 4; numerical results on 1D and 2D applications are shown in Section 5 before drawing concluding remarks in Section 6. 2. Reference problem and notations We consider a transient diffusion problem defined on an open bounded domain X  Rd (d = 1, 2, 3), with boundary @ X, over a time interval I ¼ ½0; T. We assume that a prescribed zero temperature is applied on a part @ uX – ; of the boundary, and that the domain is subjected to a time-dependent thermal loading that consists of: (i) a given thermal flux qd(x, t) on @ qX  @ X, with @ uX \ @ qX = ; and @ u X [ @ q X ¼ @ X; (ii) a source term fd(x, t) in X (see Fig. 1). The material that composes X is assumed to be isotropic. Furthermore, for the sake of simplicity, we consider that initial conditions are zero. The problem thus consists of finding the temperature-flux pair (u(x, t), q(x, t)), with ðx; tÞ 2 X  I , that verifies:  the thermal constraints:

u ¼ 0 in @ u X  I ;

ð1Þ

 the equilibrium equations:

@u ¼ $  q þ fd @t

in X  I ;

q  n ¼ qd

in @ q X  I;

ð2Þ

 the constitutive relation:

q ¼ $u in X  I ;

ð3Þ

 the initial conditions:

uðx; 0þ Þ ¼ 0 8x 2 X;

ð4Þ

n denoting the outgoing normal to X. In the following, in order to be consistent with other linear problems encountered in Mechanics (linear elasticity for instance), we carry out the change of variable q ? q which leads, in particular, to the new constitutive relation q ¼ $u.

qd

Ω fd

Fig. 1. Representation of the reference diffusion problem.

2034

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

Remark 1. The methodology we develop in this paper would apply to cases where material parameters (diffusion coefficient, density, specific heat capacity) are piecewise constant in X. However, for what we aim at showing in this paper, and for the sake of simplicity, we consider here a parameter-free problem. Defining V ¼ H10 ðXÞ ¼ fv 2 H1 ðXÞ; v j@u X ¼ 0g, the weak formulation in space of the diffusion problem consists of finding u(x, t), with uð; tÞ 2 V for all t 2 I , such that:

bðu; v Þ ¼ lðv Þ 8v 2 V;

8t 2 I ;

ð5Þ

ujt¼0þ ¼ 0; where bilinear form b(, ) and linear form l() are defined as:

bðu; v Þ ¼

Z  X

 @u v þ $u  $v dX; @t

lðv Þ ¼

Z

fd v dX 

Z

X

qd v dS:

@q X

ð6Þ As regards the space–time weak formulation, we introduce the functional spaces T ¼ L2 ðI Þ and L2 ðI ; VÞ ¼ V T . We therefore search the solution u 2 L2 ðI ; VÞ such that u_ 2 L2 ðI; L2 ðXÞÞ and

Bðu; v Þ ¼ Lðv Þ 8v 2 L2 ðI; VÞ with

Bðu; v Þ ¼

ð7Þ

Bðum ; v m Þ ¼ Lðv m Þ 8v m 2 V m T :

T

bðu; v Þdt þ

Z

uðx; 0þ Þv ðx; 0þ ÞdX;

Lðv Þ ¼

Z

T

lðv Þdt:

 Z m  Z X k_ i wi wj dX þ ki $wi  $wj dX ¼ lðwj Þ 8t 2 I ;

ð8Þ Problem (7) is classically solved using the FEM in space associated with a particular time integration scheme, or a (discontinuous) Galerkin approximation in time. The exact solution of (7), which is usually out of reach, is denoted (uex, qex). To close this section, we introduce the following L2-norms and inner products, for functions u and v in V T :

Z

hu; v iX ¼

Z

hhu; v ii ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kv kX ¼ hv ; v iX ;

uv dX;

X T

uv dt; 0

Z

0

T

Z

kv kI ¼

uv dXdt;

X

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hv ; v iI ;

ð9Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jjjv jjj ¼ hhv ; v ii:

ð12Þ

ki ð0þ Þhwi ; wj iX ¼ 0:

Let us notice that initial conditions in (7) are considered in a weak sense only, contrary to initial conditions in (12). In practice, we use the formulation of (12) to take initial conditions into account. Another point of view consists of first constructing a lowdimensional subspace T m ¼ spanfkj gm j¼1  T . The Galerkin approximation um 2 V T m is then defined as the one which verifies:

Bðum ; v m Þ ¼ Lðv m Þ 8v m 2 V T m :

m  X

aij

i¼1

ð13Þ

Z

wi w dX þ bij

X

Z X



$wi  $w dX ¼ cj lðw Þ 8w 2 V

ð14Þ

with aij ¼ hk_ i ; kj iI ; bij ¼ hki ; kj iI , and cj ¼ h1; kj iI . In practice, this last approach is not used for large scale applications due to the prohibitive computational cost required for solving (14).

3. Decomposition of the solution with PGD

3.2. Particular case: the POD

3.1. Separated representation Solutions of parameterized problems describe a set which, without being known, often has a small complexity in the sense that a correct representation of this set can be performed from linear combinations of a few relevant elements in the set (even though this is not the case with hyperbolic problems such as transient dynamics problems); this is the main motivation for model reduction methods. We restrict here to the case of model reduction based on a separated representation of the solution; it is presented for the diffusion problem presented in Section 2, the associated solution depending on a space parameter x and a time parameter t. The basic idea in separated representation is to construct the solution u(x, t), over the space–time domain, as a sum of m products of space and time functions (m 2 N is the order of the representation): m X

X

This leads to the mapping Wm = S(Km) which can be interpreted as a system of m coupled PDEs in space (for j = 1, 2, . . . , m):

X

uðx; tÞ  um ðx; tÞ ¼

m X i¼1

0

X

ð11Þ

This leads to the mapping Km = T(Wm) in which the time functions basis Km is obtained solving the system of m ODEs in time (for j = 1, 2, . . . , m):

i¼1

Z 0

hu; v iI ¼

Remark 2. The separated representation proposed in (10) is not the only one that can be considered. First, when domain X presents sufficient symmetries, it would be valuable to decompose the space dependency into products of one-dimensional space functions Xi(x)Yi(y)Zi(z). Second, in cases material parameters can vary (and are thus real parameters of the problem), functions of these material parameters could also be considered in the decomposition. The methodology we propose here for error estimation can be applied to these last two separated representations; however, a specific treatment is necessary to obtain an admissible solution, and will be the topic of a forthcoming paper. Classical model reduction techniques start with the construction of the low-dimensional subspace V m ¼ spanfwj gm j¼1  V and define the Galerkin approximation um 2 V m T as the one which verifies:

wi ðxÞki ðtÞ ¼ Wm ðxÞ  Km ðtÞ;

ð10Þ

i T

where Wm(x) = [w1(x), w2(x), . . . , wm(x)] (resp. Km(t) = [k1(t), k2(t), . . . , km(t)]T) is a low-dimensional reduced basis of space functions (resp. time functions).

A widely used separated representation is the Proper Orthogonal Decomposition (POD), also known as Karhunen–Loeve Decomposition [18,20], Principal Component Decomposition [17] or, in the finitedimensional case, Singular Value Decomposition (SVD) [14]. When u is known, this optimal separated representation (tensor product approximation) can be defined as the one which minimizes the distance to the solution with respect to a given metric. Taking the L2-norm jjjjjj on V T defined in (9), we obtain the following definition of POD:

ðWm ; Km Þ ¼ argminW m ;K m jjju  W m  K m jjj2 :

ð15Þ

Noticing that inner product hh, ii has the separation property hhw1 k1 ; w2 k2 ii ¼ hw1 ; w2 iX hk1 ; k2 iI , problem (15) leads to an eigenvalue problem of the form:

GðwÞ ¼ rw

ð16Þ

that defines the space basis functions wi. Indeed, the minimization problem:

min jjju  wkjjj2

wðxÞ;kðtÞ

ð17Þ

2035

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

yields after differentiation with respect to w and k:

hu; wiX kðtÞ ¼ ; hw; wiX

hu; kiI wðxÞ ¼ : hk; kiI

3.3. The standard PGD approximation

ð18Þ

Using a substitution method, problem (17) comes down to solving:

    hu; wiX hu; wiX hu; wiX u; ¼ ; w; hw; wiX I hw; wiX hw; wiX I

ð19Þ

which can be put under the form of the classical eigenvalue problem (16). Remark 3. A similar eigenvalue problem can be obtained for k. G is called the spatial correlation operator of u; it is a symmetric positive linear operator. Under classical regularity assumptions on u, it is also a compact operator, so that classical spectral theory applies: the eigenvalue problem has a countable sequence of eigensolutions with positive eigenvalues ri and orthogonal eigenfunctions {wi}. The eigenpairs (wi, ri) being sorted by decreasing eigenvalues, an optimal separated representation of u can be obtained defining time functions ki as:

ki ðtÞ ¼

hwi ; uð; tÞiX 2

kwi kX

" lim jjju  um jjj2 ¼ lim kuk2  m!1

Bðum1 þ wk; wk þ w kÞ ¼ Lðwk þ w kÞ 8k 2 T ;

m X

8w 2 V: ð22Þ

We can then define two applications:  Sm : T ! V is the application which maps a time function k into a space function w = Sm(k) defined by

Bðum1 þ wk; w kÞ ¼ Lðw kÞ 8w 2 V:

ð20Þ

:

Consequently, jjjwi ki jjj ¼ kwi kX  kki kI ¼ ror verifies the convergence property:

m!1

We now introduce the recently called Proper Generalized Decomposition (PGD) technique [10,28,11,31,34] which constitutes an a priori construction of a separated representation of the solution u. An attractive feature is that it does not require any knowledge on u; neither functions ki(t) nor functions wi(x) are initially given, and both families are computed on the fly. In this section, we give a classical definition of PGD, called progressive Galerkin-based PGD. We assume that a separated variables decomposition um1 = Wm1  Km1 of order (m  1) is known. A new couple ðw; kÞ 2 V  T is defined for the order m decomposition as the pair which verifies the following double Galerkin orthogonality criterion:

ð23Þ

It is associated with a time-independent space problem;

pffiffiffiffiffi

ri , and the truncation er T m : V ! T is the application which maps a space function w into a time function k = Tm(w) defined by

#

ri ¼ 0:

ð21Þ

Bðum1 þ wk; wk Þ ¼ Lðwk Þ 8k 2 T :

i¼1

Functions wi and ki being orthogonal with respect to inner products h, iX and h; iT , respectively, the decomposition is called a biorthogonal decomposition of u. Until now, POD has been mainly considered as an a posteriori reduction method, extracting information from data. But it can also be used as an a priori model reduction method, such as in optimization [41], fluid dynamics [12,16], or thermal conduction [2]. In this framework, as information on the solution u needs to be known in order to perform POD, a partial, total or even approximate resolution of the reference problem in a preliminary stage called learning phase or snapshot is carried out to construct the POD basis. Equations of the problem are then projected on the extracted POD basis which leads to the solution of reduced order problems similar to (11) or (13). Remark 4. There are two ways to realize the snapshots: the first one consists in using the full model and to perform a fine calculation over a limited time interval (more generally, with a limited number of parameters); the second one consists in realizing a set of calculations over the full time interval but with a coarse mesh. Remark 5. In the finite dimensional case where we possess discretized values of u at Nx space points and Nt time points, the solution of the POD minimization problem is given by the Singular Value Decomposition (SVD) of the so-called snapshots matrix A defined as Aij = u(xi, tj), 1 6 i 6 Nx, 1 6 j 6 Nt. The SVD consists of decomposing the real matrix A (of dimensions Nx  Nt) under the form A ¼ URVT where U and V are orthogonal matrices, respectively of dimension Nx  Nx and Nt  Nt, and R is a diagonal matrix of dimension Nx  Nt that contains positive values r1, . . . , rp called singular values of A. A priori model reduction with POD can be useful in many multiparameter problems, such as optimization where we usually solve many problems on quasi-similar configurations. However, a drawback is the choice of the snapshots as the initially constructed basis may not be optimal and has strong consequences on the effectivity of the model reduction.

ð24Þ

It is associated with a time problem (scalar ODE). A couple (w, k) verifies (22) if and only if w = Sm(k) and k = Tm(w), which is a nonlinear problem. It is shown in [34] that the progressive Galerkin-based PGD can be interpreted in terms of pseudo eigenvalue problems, functions w and k being respectively the dominant eigenfunctions of mappings Gm :¼ Sm  T m and e m :¼ T m  Sm . This interpretation is fruitful in the sense that it alG lows to propose dedicated algorithms, inspired from classical algorithms for eigenvalue problems, for the construction of the decomposition. One possible algorithm, denoted power iterations algorithm, is described below and will be used in the numerical experiments. It is based on power-type iterations; starting from an initial space function w(0) = Sm(k(0)) (in practice, k(0) is generated randomly), the sequence w(k+1) = Gm(w(k)) is computed leading to the algorithm: ? ? ? ? ? ? ? ? ? ? ?

for m = 1 to mmax do define k(0) (initialization) for k = 1 to kmax do compute w(k) = Sm(k(k1)) normalize w(k) (kw(k)kX = 1) compute k(k) = Tm(w(k)) check convergence of (w(k)k(k)) end for set wm = w(k) and km = k(k) set um = um1 + wmkm and check convergence end for

We therefore remark that um is built by means of the solution of a few space and time problems, without any knowledge on u: the first step w(k) = Sm(k(k1)) of the iterative strategy, which is the most costly, consists of solving a space problem with fixed time function k(k1); the second step k(k) = Tm(w(k)) consists of solving a time problem over the whole time domain I , with fixed space function w(k). In practice, these two problems are solved using classical discretization techniques i.e. the FEM and a given time integration scheme.

2036

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

Remark 6. In the previous algorithm, there is no need for normalization of functions w and k. However, we choose here to normalize space functions w for practical reasons. Remark 7. At a given order m, space function wm may be orthogonalized with respect to the existing space basis {w1, w2, . . . , wm1} (using the Gram-Schmidt method for instance). This procedure usually improves the results as it avoids numerical pollutions due to bad conditioning (possible dependency between functions wi). The power iterations algorithm does not usually lead to an optimal decomposition um = Wm  Km in the sense of the Galerkin projection. An optimal Galerkin PGD could be defined by imposing the residual to be simultaneously orthogonal to reduced spaces V m ¼ spanðWm Þ  V and T m ¼ spanðKm Þ  T , i.e. verifying the two Galerkin orthogonality criteria (11) and (13). The associated algorithm would then be: ? for m = 1 to mmax do ? define Kð0Þ m (initialization) ? for k = 1 to kmax do

ðk1Þ ? compute WðkÞ m ¼ S Km

ðkÞ ? compute KðkÞ m ¼ T Wm

? ? ? ?

check convergence end for ðkÞ set um ¼ WðkÞ m  Km and check convergence end for

This algorithm leads to optimal sets of functions in the sense of Galerkin projection and leads to a decomposition um which can be significantly better than the decomposition obtained with a progressive definition of the PGD. Nevertheless, each iteration requires the application of mappings T and S; the application of mapping T is relatively cheap since it corresponds to a set of ODEs (formulation of the initial problem on a reduced basis of space functions). However, the application of mapping S is generally very costly for large scale applications (2D, 3D) since it involves the solution of a set of coupled PDEs. A simple modification of the initial power iterations algorithm is nevertheless possible to build a better decomposition than the progressive PGD without applying mapping S [28,34,37]. It consists in introducing the application of mapping T in order to update the whole set of time functions Km after each construction of a new couple (wm, km). This leads to the following algorithm, which corresponds to the power iterations algorithm with an additional updating of time functions. ? ? ? ? ? ?

for m = 1 to mmax do perform steps 2 to 8 of the power iterations algorithm set wm = w(k) compute Km = T(Wm) set um = Wm  Km and check convergence end for

This new decomposition, which is better quality than the one given by the initial progressive PGD, will be used in the following. Remark 8. Two other algorithms have been recently proposed for PGD approximations: - the first one, introduced in [33] for stochastic PDEs, is inspired from the Arnoldi algorithm for solving classical eigenproblems. This algorithm allows the capture of an approximation of the dominant eigenspace of G without applying the mapping S. For many applications, power iterations algorithm with update

and Arnoldi’s algorithm lead to very similar convergence properties. Arnoldi’s algorithm is usually more efficient since it requires the solution of only m uncoupled PDEs in order to build an order m decomposition; - the second one, introduced in [34], comes from a new definition of PGD which allows to improve the convergence properties of Galerkin-based PGD with respect to a given metric. This new definition can be interpreted as a PGD based on Petrov–Galerkin criteria, where orthogonality of the residual is imposed with respect to another set of space and time functions, which are solutions of an adjoint problem. It is shown that this so-called dual-based PGD leads to a decomposition which is closer to the classical POD. Concurrently to the Galerkin-based PGD, another definition of the PGD based on a minimal residual criterion has been developed in [28]. The residual RðuÞ 2 L2 ðI ; VÞ of (7) is defined as

hhv ; RðuÞii :¼ Lðv Þ  Bðu; v Þ :¼ hhv ; w  BðuÞii 8v 2 L2 ðI ; VÞ; ð25Þ 2

where w 2 L ðI ; VÞ and operator B are defined by using Riesz representations in Hilbert space L2 ðI ; VÞ. In the natural definition of the Minimal Residual PGD, an optimal couple (wm, km) is then the one which minimizes the residual norm:

ðwm ; km Þ 2 argminðw;kÞ jjjRðum1 þ wkÞjjj2 :

ð26Þ

An advantage of the minimal residual formulation is that the convergence with m of the decomposition um is monotonic, if convergence is evaluated with the residual norm. It is thus a robust construction of a separated representation which can be used in cases where Galerkin-based PGD fails. However, it has several drawbacks: - the resulting decomposition may present very poor convergence properties with respect to usual norms. Taking the natural norm of L2(X) L2([0, T]) usually leads to bad convergence rates; - the formulation requires much more computational cost than the Galerkin PGD. Consequently, the PGD based on a minimal residual formulation should be avoided in practical applications for which Galerkinbased PGDs work. 4. Procedure for global error upper bound in the PGD approximation In this part, we develop an error estimation method based on the constitutive relation error and applied to simulations carried out with PGD. Considering the diffusion problem introduced in Section 2, we first present two PGD approaches that provide for kinematically (or thermally) and statically admissible solutions of the problem; these approaches are both based on classical kinematic Galerkin formulations. We then introduce the error estimator and analyze its properties. Let us notice that denotations kinematic and static refer, in a more general setting, to primal and dual variables, respectively. 4.1. Kinematic PGD approach We consider the space–time weak formulation (7). In order to compute an approximation of u, we use the progressive PGD approach defined in Section 3.3, associated with the power iterations algorithm and updating of the time functions. Knowing um1, this approach first consists of finding the pair (w, k) such that:

2037

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

Bðwk; wk þ w kÞ ¼ Rm1 ðwk þ w kÞ 8k 2 T ;

8w 2 V;

ð27Þ

where Rm1() = L()  B(um1, ) is the residual at order m  1; this is solved using the power iterations algorithm. After convergence of (w(k), k(k)), we take wm = w(k) and update all the time functions using the mapping Km = T(Wm). We then obtain the order m PGD approximation um of u:

um ¼ Wm  Km :

Dt



qðuh; m Þ  $v dX ¼ lðv Þ 

X

Z X

h;Dt @um v dX 8v 2 V; @t

8t 2 I;

nor Z Z h;Dt @um Dt



qðuh; v dX 8v 2 V h ; m Þ  $v dX ¼ lðv Þ  @t X X

8t 2 I :

ð37Þ

ð38Þ

ð28Þ

In practice, the power iterations algorithm requires the solution of two kinds of problems:  mapping w = Sm(k) leads to a space problem that consists of finding w 2 V such that

Z

Z

faS ww þ bS $w  $w gdX

X

¼ lðw Þ 

m1 X i¼1

Z





aS;i wi w þ bS;i $wi  $w dX

Remark 9. We will see that in cases where the KA solution h;Dt h;Dt ðum ; qðum ÞÞ also verifies the equilibrium equation in the FE sense, the procedure used to define a guaranteed error estimator is very simple. In the following, we define by Q (resp. Qh ) the space of fluxes which verify (37) (resp. (38)). We introduce in the next section a static PGD approach that provides for a flux that belongs to Qh .

ð29Þ 4.2. Static PGD approach

X

with

_ ki ; b ¼ kkk2 ; lðÞ ¼ hlðÞ; ki ; aS ¼ hk; S I I I aS;i ¼ hk_ i ; kiI ; bS;i ¼ hki ; kiI :

We define an order m PGD approximation of the flux as:

ð30Þ

qm ðx; tÞ ¼ qh0 ðx; tÞ þ

ðaS M þ bS KÞX ¼ F

ð31Þ

with

wh ¼ NT X; M ¼ NT NdX; X Z NT DT DNdX ðD : space derivative operatorÞ; K¼

where Um(x) = [/1(x), /2(x), . . . , /m(x)]T and Nm(t) = [n1(t), n2(t), . . . , nm(t)]T are the reduced bases of space and time functions, respectively.

qh0 ðx; tÞ is a particular flux which belongs to Qh i.e. which verifies (38), for each time t 2 I , and space functions /i(x) are statically admissible to zero in the FE sense i.e. they verify:

X m1 X

ð39Þ

Remark 10. The flux field could be represented using a PGD order mq different from the PGD order mu used to represent the temperature field. However, we consider mq = mu = m in the whole paper.

Z

F ¼ Ffd  Fqd 

/i ðxÞni ðtÞ ¼ qh0 ðx; tÞ þ Um ðxÞ  Nm ðtÞ;

i¼1

h

An approximate solution w of w is computed using the FEM with associated FE space V h  V; it leads to a discretized system of the form:

m X

ðaS;i M þ bS;i KÞXi ;

i¼1

ð32Þ

Z

/i  $v dX ¼ 0 8v 2 V h ;

8t 2 ½0; T:

ð40Þ

X

 mapping k = Tm(w) leads to a time problem that consists of finding k 2 T such that

kð0Þ ¼ 0;

aT k_ þ bT k ¼ dT

ð33Þ

fd ðx; tÞ ¼ F d ðxÞjf ðtÞ;

with

aT ¼

Z

w2 dX;

bT ¼

Z

X

$w  $wdX;

X

Dt

An approximate solution k integration scheme.

ð34Þ of k is computed using a given time

M ij

¼

Z X

M K_ m þ K Km ¼ F

wi wj dX;

K ij

¼

Z X

ð35Þ

Z

$zhf  $v dX ¼

Z

F d v dX 8v 2 V h ; Z Z $zhq  $v dX ¼  Q d v dS 8v 2 V h ; X @q X Z Z $shi  $v dX ¼  wi v dX 8v 2 V h ði ¼ 1; . . . ; mÞ; X

X

$wi  $wj dX;

¼ lðwj Þ:

ð36Þ

An approximate solution KDmt of this system is easily obtained using a time integration scheme. We eventually obtain an order m PGD approximation of the P h Dt Dt form uh; ¼ m i¼1 wi ðxÞki ðtÞ. It is kinematically admissible (KA) in m the sense that it verifies Eqs. (1) and (4) of the initial problem. HowDt h;Dt h;Dt h;Dt ever, the associated pair ðuh; m ; qðum ÞÞ, with qðum Þ ¼ $um , is usually not statically admissible, even in a FE sense. Indeed, it does verify neither:

ð42Þ

X

qh0 ðx; tÞ can be defined as

"

F j

ð41Þ

the computation of can be easily performed with a classical FE procedure. Indeed, computing functions ðzhf ; zhq ; sh1 ; sh2 ; . . . shm Þ 2 ðV h Þmþ2 such that: X

As regards the updating of the time functions at order m, it merely leads to an order m system of the form:

Kmjt¼0 ¼ 0;

qd ðx; tÞ ¼ Q d ðxÞjq ðtÞ;

qh0 ðx; tÞ

dT ¼ lðwÞ  bðum1 ; wÞ:

with

The associated space is denoted Qh0 . On the one hand, assuming that the external loading (fd, qd) can be written with the radial approximation:

qh0 ðx; tÞ

¼ $

zhf ðxÞ

jf ðtÞ þ $

zhq ðxÞ

jq ðtÞ þ

m X

#

$

shi k_ i ðtÞ

:

ð43Þ

i¼1

On the other hand, noticing that the dual Galerkin formulation of the diffusion problem is to find q 2 Q  T such that:

Z 0

T

Z

q  q dXdt ¼ 0 8q 2 Q0

ð44Þ

X

P and assuming that qm1 ¼ qh0 þ m1 i¼1 /i ni is known, the progressive PGD at order m (with updating of time functions) first consists of searching a couple (/, n) that verifies:

2038

Z

T

0

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

Z X

ðqm1 þ /nÞ  ð/ n þ /n ÞdXdt ¼ 0 8/ 2 Qh0 ;

8 n 2 T

ð45Þ

Similarly to the kinematic approach, the associated power iterations algorithm requires the solution of two kinds of problems:  a space problem of the form / ¼ e S m ðnÞ such that:

Z

X

fa/ þ bg  / dX ¼ 0 8/ 2 Qh0

ð46Þ

with a ¼ knk2I and b ¼ hn; qm1 iI . Using duality arguments, we deduce that there is a function wh ðxÞ 2 V h such that:

a/ þ b ¼ $wh :

ð47Þ

Therefore, as / should belong to Q h0 ; wh should verify:

Z

X

1

a



$wh  b  $w dX ¼ 0 8w 2 V h :

ð48Þ

An important point is that this last problem can be solved with the standard FEM, and / is eventually deducted using (47). e m ð/Þ such that:  a time problem of the form n ¼ T

Z

ðqm1 þ /nÞ  /dX ¼ 0 8t 2 I :

ð49Þ

^ h;Dt Þ verifying (2) is said statically admissible (SA); ^ h;Dt ; q - a pair ðu ^ h;Dt is the most technical point. It consists the construction of q of using an energetic relation, denoted prolongation condition, Dt involving the field qh; coming from the static PGD approach. m Dt Due to the fact that qh; 2 Qh , the prolongation condition enables m to compute, on finite element edges, some flux densities which are equilibrated with the external loading (fd(t),qd(t)) (provided that this external loading can be exactly represented in the FE ^ h;Dt is constructed space). In a second stage, an admissible flux q from the densities by solving local problems at the element level. Full details on the technique used to construct the equilibrated ^ h;Dt can be found in [27,25,26,29,38]; here, we spedensities and q cifically use the new technique introduced in [29,38] which defines local problems on patches of elements, and represents a nice compromise between robustness, computational cost, and implementation facilities. We also refer to [36,13] for other approaches enabling to construct a statically admissible field. Dt ^ h; Remark 11. In practice, q is constructed as m

Dt ^ h0 þ ^ h; ¼q q m

m X

^ h nDt ; / i i

i¼1

X

This provides for the explicit expression of n(t) (after normalization of /);

nðtÞ ¼ hqm1 ð; tÞ; /iX :

ð50Þ

In a second step, the vector Nm(t) of time functions is updated e ðUm Þ. The power iterations algousing the global mapping Nm ¼ T rithm which is employed to compute the static PGD version of q thus reads: ? ? ? ? ? ? ? ? ? ? ? ? ?

compute qh0 for m = 1 to mmax do define n(0) (initialization) for k = 1 to kmax do compute /ðkÞ ¼ e S m ðnðk1Þ Þ normalize /(k) (k/(k)kX = 1) e m ð/ðkÞ Þ compute nðkÞ ¼ T check convergence of (/(k)n(k)) end for set /m = /(k) and nm = n(k) e ðUm Þ (updating of time functions) compute Nm ¼ T set qm ¼ qh0 þ Um  Nm and check convergence end for

Of course, space and time problems associated with mappings e e m are again approximately solved, so that the computed S m and T P h Dt Dt decomposition is noted qh; ¼ qh0 þ m m i¼1 /i ni . Dt h;Dt Eventually, the pair ðuh; ; q Þ is SA in the FE sense, i.e. it verifies m m the FE equilibrium equations. 4.3. Definition of an admissible pair The global error estimator we set up in this section is based on the constitutive relation error (CRE) [24,27]. The use of the CRE requires the computation of an admissible temperature-flux pair, de^ h;Dt Þ in the following, associated to the computed PGD ^ h;Dt ; q noted ðu ^ h;Dt Þ is admissible in the sense that it ^ h;Dt ; q solution. The pair ðu should verify (1), (2) and (4). Its construction is performed from Dt h;Dt the approximate solution ðuh; m ; qm Þ at hand: ^ h;Dt verifying (1) and (4) is said kinemati- a temperature field u cally admissible (KA); a simple choice consists in taking the temDt perature field uh; coming from the kinematic PGD approach; m

^ ðzhf Þjf ðtÞ þ q ^ ðzhq Þjq ðtÞ þ ^ h0 ¼ q q

m X

^ ðshi Þk_ Di t ðtÞ; q

ð51Þ

i¼1

^ h are SA and SA0 fields verifying the prolongation ^ h0 and / where q i condition with qh0 and /hi , respectively. The associated computational cost is (m + 2)  (CvNv + CeNe), where Nv (resp. Ne) is the number of vertices (resp. elements) in the mesh, and Cv (resp. Ce) is the computational cost associated to a problem on a patch of elements (resp. to a problem at the element level). 4.4. Calculation of the constitutive relation error and associated error upper bound ^ Þ (the previous computations ^; q From any admissible solution ðu Dt ^ h;Dt ^ Þ ¼ ðu ^; q ^ h; lead to take ðu ; q Þ), the CRE is defined as a measure m m ECRE of the non-verification of (3) by this admissible pair. It reads:

^Þ ^; q E2CRE ðu

Z

T 0

Z X

^  $u ^  $u ^  $u ^ Þ  ðq ^ÞdXdt ¼ jjjq ^ jjj2q : ðq

ð52Þ

A fundamental property of the CRE is its relation with the exact solution (uex, qex):

^ jjj2q þ Gðuex  u ^Þ ^jjj2u þ jjjqex  q ^ Þ ¼ E2CRE ðu ^; q jjjuex  u with jjj  jjj2u ¼ P 0.

RT R 0

X

^Þ ¼ $ðÞ  $ðÞdXdt and Gðuex  u

ð53Þ R

X ðu

ex

^ Þ2jt¼T dX u

Proof.

^  $u ^  qex þ $uex  $u ^ k2q ¼ kq ^ k2q kq ^ k2q þ kuex  u ^ k2u  2 ¼ kqex  q ^ k2q þ kuex  u ^ k2u þ 2 ¼ kqex  q ¼ kqex 

^ k2q q

þ kuex 

^ k2u u

Z

^ Þ  $ðuex  u ^ ÞdX ðqex  q

X

Z

d þ dt

X

Z

^ Þ ex @ðuex  u ^ ÞdX ðu  u @t

X

^Þ2jt dX ðuex  u

and the result yields integrating (54) over [0, T].

ð54Þ

h

Eq. (53), which can be seen as the Prager–Synge equality for ^ jjju 6 ECRE time-dependent problems, directly shows that jjjuex  u so that the measure ECRE is an upper bound of the global error (given in the energy norm). It takes into account all error sources,

2039

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

we introduce error indicators which enable to split contributions of the different error sources, and thus to drive adaptive procedures effectively.

qd(t)

4.5. Error indicators for the PGD approximation

L

In order to optimize the adaptive process which aims at decreasDt ^ h; ing the error jjjuex  u m jjju , one should be able to estimate the contributions brought by the various error sources. More specifically, we wish to split the contribution coming from the time/space discretization (use of a space mesh in the FE method, use of a time integration scheme) and that coming from the PGD representation (truncation of the sum, assuming that error due to iteration stopping is neglectable); this is carried out defining appropriate error indicators.

Fig. 2. The 1D diffusion problem.

i.e. space and time discretizations as well as truncation of the sum Dt ^ h;Dt ^ Þ ¼ ðu ^; q ^ h; in the PGD representation when choosing ðu m ; qm Þ. The computational cost required to compute ECRE is proportional to Ng  (Np + 1) where Ng (resp. Np) is the number of Gauss points in the mesh (resp. the number of time steps). In the next section,

4

0.6

3.5

0.5

3

0.4

u(L/2,t)

u(x,T)

2.5 2 1.5

0.3 0.2 0.1

1

0

0.5 0 0

1

2

3

4

5 x

6

7

8

9

−0.1

10

0

1

2

3

4

5 t

6

7

8

9

10

0.9

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

ψ

0.9

0.8

ψ

0.9

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0

0

1

2

3

4

5

x

6

7

8

9

10

0

1

2

3

4

5

x

6

7

8

9

0

10

5

5

5

4.5

4.5

4.5

4

4

4

3.5

3.5

3.5

3

3

2.5

2.5

λ

3 2.5

λ

λ

ψ

Fig. 3. Exact solution of the problem: space evolution at t = T (left), and time evolution at x = L/2 (right).

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

1

2

3

4

5

t

6

7

8

9

0

10

0

1

2

3

4

5

t

6

7

8

9

0

10

0

1

2

3

4

5

6

7

8

9

10

0

1

2

3

4

5

6

7

8

9

10

x

t

ðkÞ

ðkÞ

Fig. 4. Evolution of modes with the number k of iterations, for k = 1, 2, 3 (from left to right): space function w1 (upper row) and time function k1 (lower row).

Table 1 Thermic energy of the first PGD mode with respect to the iteration number k. k

1

2

3

4

5

6

7

8

9

10

Energy

14.4185

11.4975

12.1713

12.0278

12.0589

12.0522

12.0536

12.0533

12.0534

12.0533

2040

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047 0.9

4.5

0.8

4

0.7

3.5 3

0.5

2.5 λ

ψ

0.6

0.4

2

0.3

1.5

0.2

1

0.1

0.5 0 0

1

2

3

4

5 x

6

7

8

9

0

10

0

0.4

−0.1

0.2

−0.2

0

−0.3

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

λ

ψ

0.6

0

−0.2

−0.4

−0.4

−0.5

−0.6

−0.6

−0.8

0

1

2

3

4

5 x

6

7

8

9

−0.7

10

0.5

0.2 0.15 0.1 0.05

0 λ

ψ

0 −0.05 −0.1 −0.5

−0.15 −0.2 −0.25

−1

0

1

2

3

4

5 x

6

7

8

9

−0.3

10

0.8

0.06

0.6

0.04

0.4

0.02

0.2

0 −0.02 λ

ψ

0 −0.2

−0.04

−0.4

−0.06

−0.6

−0.08

−0.8

−0.1

−1

−0.12

0

1

2

3

4

5 x

6

7

8

9

10

0.8

0.07 0.06

0.6

0.05

0.4

0.04 0.03

0

λ

ψ

0.2

0.02 0.01

−0.2

0

−0.4

−0.01

−0.6

−0.02

−0.8

−0.03 0

1

2

3

4

5 x

6

7

8

9

10

Fig. 5. Space functions wm (left) and time functions km (right) with respect to order m, for m = 1, 2, 3, 4, 5 (from top to bottom).

The idea to define such error indicators, which is not restricted to the PGD framework and has been applied in other works (see [27] for instance), is quite simple: we first define the solution uh,Dt

of the discrete problem obtained discretizing the initial diffusion problem (5) in space and time (with the same discretizations as Dt those used to compute uh; m ). Using for instance a classical FE

2041

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047 2

0

3

1.8 1.6

2.5

−0.5

1.4

2

1.2

−1

1

1.5

0.8

−1.5 1

0.6 0.4

−2

0.5

0.2 0 9

9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 x

−2.5 0

1

2

3

4

5 x

6

7

8

9

10

Fig. 6. Representation of fluxes required to compute qh0 . From left to right:

dzhq dx

1

;

dsh1 , dx

2

3

and

dsh2 . dx

4

5 x

6

7

8

9

10

Remark 12. Using a similar procedure, we could also split contribution of the error coming from space discretization and that coming from time discretization. In order to do so, we should introduce a reference problem obtained by discretizing (5) in space only. The associated solution uh should then verify a system of the form:

0.4 0.35 0.3

relative error

0 0

0.25

Uhjt¼0 ¼ 0;

_ h þ KUh ¼ Fh MU

8t 2 T :

ð60Þ

Comparing values of ECRE,dis and ECRE,PGD enables to assess the critical error source and to drive the adaptive process effectively.

0.2 0.15

5. Numerical results

0.1 0.05 0 1

2

3

4

5

6

7

8

9

10

m Fig. 7. Evolution of ECRE with respect to the PGD order m.

method and a forward Euler time scheme, u a system of the form:

U1h ¼ 0;

M

Upþ1 h

 Dt

Uph

þ KUph ¼ Fph

h,Dt

In all the following examples, the power iterations algorithm is initialized with k(0)(t) = t (n(0)(t) = t). We use the FEM with linear elements to solve problems in space, and a forward Euler scheme to solve problems in time. The number of iterations performed in the power iterations algorithm is fixed to kmax = 4. 5.1. 1D problem

is obtained solving

8p 2 ½1; P  1;

ð55Þ

where P is the number of time steps. We can then write:

  h;Dt Dt uex  um ¼ uex  uh;Dt þ uh;Dt  uh; ; m

ð56Þ

We consider a 1D diffusion problem over the space–time domain ]0, L[  [0, T]. A homogeneous Dirichlet boundary condition is imposed at x = 0 whereas a given time-dependent flux qd(t) is applied at x = L (see Fig. 2). For the numerical study, we take: L = 10, T = 10, and a constant loading qd(t) = 1. The (quasi-) exact solution of the problem, obtained from an overkill solution with 100 elements and 10000 time steps, is given in Fig. 3.

which leads, using the Galerkin orthogonality, to the result: Dt 2 ex h;Dt 2 Dt 2 jjjuex  uh; jjju þ jjjuh;Dt  uh; m jjju ¼ jjju  u m jjju ; ex

ð57Þ

h,Dt

where jjju  u jjju is the contribution of the error coming from h;Dt space and time discretizations alone, and jjjuh;Dt  um jjju is the contribution of the error coming from the PGD representation alone. The second contribution can be easily estimated, without any additional computation, defining problem (55) as the reference problem and using the CRE. In that framework, an associated admis^ h;Dt Þ (in the sense of (55)) should be defined; it is di^ h;Dt ; q sible pair ðu rectly obtained at each time step as a direct post-processing Dt h;Dt (interpolation) of the solution ðuh; m ; qm Þ at hand (the latter solution is admissible in a discrete sense). We can thus compute the estimate:

^ h;Dt  $u ^ h;Dt jjjq ECRE;PGD ¼ jjjq h;Dt

ð58Þ

h;Dt um jjju

of the error jjju  coming from the PGD representation. In a second time, (57) enables to define the estimate ECRE,dis of the error jjjuex  uh,Dtjjju coming from space and time discretizations; it reads:

ECRE;dis ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E2CRE  E2CRE;PGD :

ð59Þ

5.1.1. Convergence of the PGD procedure We first analyze the convergence of the power iterations algorithm. For the approximate solution of space and time problems involved in this algorithm, the structure is divided into Ne = 10 elements and the time domain into Np = 100 time steps. For the ðkÞ first order (m = 1), we plot in Fig. 4 the space function w1 and time ðkÞ function k1 with respect to k (=1, 2, 3). ðkÞ ðkÞ We give in Table 1 the associated thermal energy of w1 k1 with respect to k 2 [1, 10]. We observe that the convergence is reached quite rapidly after a few iterations (k  4). In the following computations, we thus take kmax = 4. We now analyze the convergence of the PGD solution with respect to the number m of modes included in the PGD representation. We plot in Fig. 5 the modes (wm, km) obtained after applying the power iterations algorithm, with updating of the time functions at mmax = 5. As regards the static approach, we plot in Fig. 6 some of the dzh dsh fields dxq and dxi ði ¼ 1; . . . ; mÞ which enable to compute qh0 2 Qh . We easily show that terms /i are zero (which is a particular case

2042

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047 0.9

5

0.8

4.5 4

0.7

3.5 3

0.5 λ

ψ

0.6

0.4

2

0.3

1.5

0.2

1

0.1

0.5

0

0

1

2

3

4

5 x

6

7

8

9

0

10

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0 −0.1

0.35

−0.2 −0.3 −0.4

0.25 λ

ψ

1

0.4

0.2

−0.7

0.1

−0.8

0.05 0

−0.5 −0.6

0.15

−0.9

0

1

2

3

4

5 x

6

7

8

9

−1

10

0

0 −0.05

−0.1

−0.1

−0.2

−0.15 −0.2

−0.3 λ

ψ

0

0.45

0.3

−0.25

−0.4

−0.3

−0.5

−0.35 −0.4

−0.6

−0.45

−0.7 0

1

2

3

4

5 x

6

7

8

9

−0.5 0

10

0.45

2

3

4

5 t

6

7

8

9

10

0.05

0.4

0

0.35

−0.05

0.3

−0.1

0.25

−0.15

0.2

−0.2

0.15

−0.25

0.1

−0.3

0.05

−0.35

0

1

λ

ψ

2.5

0

1

2

3

4

5 x

6

7

8

9

−0.4

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0.5

0

0.4

−0.2

0.3

−0.4

0.2 0.1

ψ

λ

−0.6

0

−0.8

−0.1

−1

−0.2 −0.3

−1.2 −1.4 0

−0.4 1

2

3

4

5 x

6

7

8

9

10

−0.5

Fig. 8. Space functions wm (left) and time functions km (right) with respect to order m, for m = 1, 2, 3, 4, 5 (from top to bottom), without neither orthogonalization nor updating.

2043

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047 0.9

5

0.8

4.5 4

0.7

3.5 3

0.5 λ

ψ

0.6

0.4

2.5 2

0.3

1.5

0.2

1

0.1

0.5

0

0

1

2

3

4

5 x

6

7

8

9

0

10

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

0

1

2

3

4

5 t

6

7

8

9

10

−0.08 0

1

2

3

4

5 t

6

7

8

9

10

0.6

0.6

0.4

0.4

0.2

0.2 0 λ

ψ

0

0

−0.2

−0.2

−0.4

−0.4

−0.6

−0.6

−0.8

0

1

2

3

4

5 x

6

7

8

9

−0.8

10

0.5

0.15 0.1 0.05 0

0 λ

ψ

−0.05 −0.1

−0.15 −0.5

−0.2 −0.25 −0.3

−1

0

1

2

3

4

5 x

6

7

8

9

10

−0.35

0.8

0.06

0.6

0.04

0.4

0.02 0 −0.02

0

−0.04

λ

ψ

0.2

−0.2

−0.06

−0.4

−0.08

−0.6

−0.1

−0.8

−0.12

−1

0

1

2

3

4

5 x

6

7

8

9

10

−0.14

0.8

0.04

0.6

0.02

0.4 0

0

−0.02

λ

ψ

0.2

−0.2

−0.04

−0.4 −0.06

−0.6 −0.8 0

1

2

3

4

5 x

6

7

8

9

10

Fig. 9. Space functions wm (left) and time functions km (right) with respect to order m, for m = 1, 2, 3, 4, 5 (from top to bottom), with orthogonalization of space modes but no updating of time modes.

2044

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

of this 1D problem). Futhermore the only admissible flux field reads:

0.4 basic ortho opti

0.35

^m ðx; tÞ ¼ qd ðtÞ  q

relative error

0.3

Z

k_ i ðtÞ

L

wi ðsÞds:

ð61Þ

x

i¼1

We plot in Fig. 7 the value of the relative error estimate ^jjjq , given by the CRE, with respect to the PGD order ECRE ¼ ECRE =jjjq m which is used (m 2 [1, 10]). This error estimator takes into account all error sources (space and time discretizations, truncation of the PGD sum). When order grows, we observe that ECRE tends to an asymptotic value which is an estimate of the discretization error (due to space and time discretizations only); we will see in a following example how indicators enable to assess this last error.

0.25 0.2 0.15 0.1 0.05 0

m X

1

2

3

4

5

6

7

8

9

5.1.2. Influence of orthogonalization of space modes and updating of time functions We now observe how orthogonalization of space modes wm(x) and updating of time functions km(t) improve the quality of the PGD approximation. We represent in Fig. 8 (resp. Fig. 9) the five

10

m Fig. 10. Evolution of ECRE with respect to the PGD order m, for the basic PGD, the PGD with orthogonalization only, and the optimized PGD. 0

2.5

25

−0.1 −0.2

2

20

−0.3 u(x,T)

qd

−0.5 −0.6

u(L/2,t)

1.5

−0.4

15 10

0.5

−0.7 −0.8

1

5

0

−0.9 −1

0

0 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t

1

2

3

4

5 x

6

7

8

9

−0.5

10

0

1

2

3

4

5 t

6

7

8

9

10

9

10

Fig. 11. Linear time evolution of qd(t) (left), and associated exact solution of the problem: in space at t = T (center), in time at x = L/2 (right). 1

0.2

0.8

0.1

0.6

0.01

u(x,T)

0 −0.2 −0.4

u(L/2,t)

−0.1

0.2

−0.2 −0.3

0.005 0

−0.4

−0.6

−0.005

−0.5

−0.8 0

1

2

3

4

5 t

6

7

8

9

−0.6

10

0

1

2

3

4

5 x

6

7

8

9

−0.01 0

10

1

2

3

4

5 t

6

7

8

Fig. 12. Cyclic time evolution of qd(t) (left), and associated exact solution of the problem: in space at t = T (center), in time at x = L/2 (right). 0.04

1

0.035

0.9

0.03

0.8 0.7

0.025

relative error

relative error

qd

0.015

0

0.4

−1

0.02

0.02 0.015

0.6 0.5 0.4 0.3

0.01

0.2 0.005 0

0.1 1

2

3

4

5

m

6

7

8

9

10

0

1

2

3

4

5

m

6

7

8

9

10

Fig. 13. Evolution of ECRE with respect to the PGD order m: for a linear evolution of qd(t) (left), for a cyclic evolution of qd(t) (right).

2045

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

5.1.3. Influence of the loading evolution We consider in this section some loading evolutions which are different from the previous case qd(t) = 1. More precisely, we consider a linear evolution and a cyclic evolution. A representation of the time evolution of qd(t), as well as the quasi-exact solution of the reference diffusion problem, are given in Fig. 11 (resp. Fig. 12) for the linear evolution (resp. cyclic evolution) of the loading. We give in Fig. 13 the evolution of the relative error estimate ECRE with respect to the PGD order m 2 [1, 10], for the linear and cyclic evolution of qd(t). In the two cases, we observe that the error due to PGD reaches convergence for m  4, that means that few modes are actually necessary to get a correct approximation of the solution. However, we verify that discretization errors are larger when considering a cyclic loading than when considering a linear loading (due to time discretization mainly).

space function ψ

0.14

3

0.12

2.5

0.1 2

λ

1.5

0.08 0.06

1

0.04

0.5

0.02 0 0

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

t

space function ψ

0 −0.005

1

−0.01

0

−0.015

−2

−0.02

λ

−1

−0.025

−3

−0.03

−4

−0.035 −0.04

−5

−0.045 0

1

2

3

4

5

t

space function ψ

0.015

2 1

0.01

0

0.005

−1 −2

0

λ

first modes obtained without orthogonalization of the wm(x) nor updating of the km(t) (resp. with orthogonalization of the wm(x) but no updating of the km(t)). We compare in Fig. 10 the evolution with m of the relative error estimate ECRE in the three various approaches we considered to obtained the modes, i.e. the basic PGD (no orthogonalization, no updating), the PGD with orthogonalization only, and the optimized PGD (orthogonalization and updating). On the one hand, we clearly observe that orthogonalization of space modes wm(x) enables to avoid numerical effects which may be due to depency between modes (see Fig. 8) and the fact that these modes are not computed exactly. On the other hand, we also observe that updating of time functions provides for a better quality PGD representation of the solution. In the following, we use the optimal PGD procedure, with orthogonalization of space modes and updating of time functions.

−0.005

−3

−0.01

−4

−0.015

−5 −6

−0.02

−7

−0.025

space function ψ

0

1

2

3

4

5

6

7

8

9

10

−0.02 0

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

t

0.015 3

0.01

2

0.005

1

0

0

λ

5.2. 2D problem

−1

5.2.1. Convergence of the PGD procedure We represent in Fig. 15 the modes (wm, km) obtained for mmax = 5. We plot in Fig. 16 the value of the relative error estimate

−2

−0.01

−3 −4

−0.015

t

space function ψ

0.035 3

0.03

2

0.025 0.02

1 0 −1

0.015

λ

Let us now consider the 2D structure of Fig. 14, which presents two rectangular holes in which a fluid circulates. Using symmetries, we keep a quarter of the 2D domain (here, we model the upper right quarter of the plate alone) that we denote X. It is subjected over [0, T] to a given time-dependent flux qd(t) on the hole boundary, a zeroflux on the symmetry planes, and a given temperature ud = 0 on the other boundaries. A time-independent source term of the form fd(x, y) = 200xy is also applied in X. For the numerical study, we take T = 10 and consider a constant loading of the form qd = 1. The space mesh used to discretize X consists of Ne = 50 quadrangular elements, and the time interval [0, T] is divided into Np = 1000 time steps.

−0.005

0.01 0.005

−2

0

−3

−0.005

−4

−0.01 −0.015 0

1

2

3

4

5

t

Fig. 15. Space functions wm (left) and time functions km (right) with respect to order m, for m = 1, 2, 3, 4, 5 (from top to bottom).

y x

Fig. 14. The 2D domain with rectangular holes (left), and associated FE mesh with applied loading (right).

2046

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047

6. Conclusions

0.4 0.35

We presented a method that provides for a guaranteed and robust error estimator in numerical simulations performed using PGD. The estimator, based on the CRE, lies on a double PGD approach of the problem which enables to compute admissible fields by means of classical tools incorporated in FE codes. Indicators on the various error sources have also been derived in order to lead adaptive approaches in an optimized manner; the adaptation issue for PGD computations will be addressed in a forthcoming paper. The error estimation method is quite general: it could be applied to other parameterized problems such as stochastic or optimization problems. Furthermore, following [15,23,7,8,30], the obtained error estimates could be used to compute guaranteed error bounds on specific quantities of interest; goal-oriented error estimation in the PGD framework is a work in progress.

relative error

0.3 0.25 0.2 0.15 0.1 0.05 0

1

2

3

4

5

6

7

8

9

10

m

References

Fig. 16. Evolution of ECRE with respect to the PGD order m.

0.4 E2CRE E2PGD

0.35

2

2 −EPGD ECRE

0.3 0.25 0.2 0.15 0.1 0.05 0

1

2

3

4

5

6

7

8

9

10

m

Fig. 17. Evolution of estimates E2CRE ; E2CRE;PGD , and E2CRE;dis with respect to the PGD order m.

^ jjjq with respect to the PGD order m which is used ECRE ¼ ECRE =jjjq (m 2 [1, 10]). Again, this error estimator (which takes into account all error sources) reaches convergence at m = 4. The asymptotic value is an estimate of the discretization error alone; in the next section, we compute error indicators that enable to assess this last error at an early stage of the computation. 5.2.2. Computation of error indicators Using the procedure defined in Section 4.5, we introduce a reference problem obtained from a time and space discretization of the initial diffusion problem. After computing associated admissible fields, we obtain the estimate ECRE,PGD of the error coming from the PGD representation alone. From ECRE and ECRE,PGD, we then deduce an estimate ECRE,dis of the time/space discretization error alone. The three estimates are given in Fig. 17 with respect to the PGD order m. We observe that indicator E2CRE;dis provides for a relevant assessment of the discretization error, even for small values of m. Moreover, we can see that indicator ECRE,dis becomes larger than ECRE,PGD for m P 3; this piece of information could be used in an optimized adaptive process in which we would first increase the number of modes, up to m = 3, before refining the time and space discretizations to get a better quality solution.

[1] A. Ammar, F. Chinesta, P. Diez, A. Huerta, An error estimator for separated representations of highly multidimensional models, Computer Methods in Applied Mechanics and Engineering 199 (2010) 1872–1880. [2] R.A. Bialecki, A.J. Kassab, A. Fic, Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis, International Journal for Numerical Methods in Engineering 62 (6) (2005) 774–797. [3] S. Boyaval, C. Le Bris, Y. Maday, N.C. Nguyen, A.T. Patera, A reduced basis approach for variational problems with stochastic parameters: application to heat conduction with variable robin coefficient, Computer Methods in Applied Mechanics and Engineering 198 (41–44) (2009) 3187–3206. [4] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 1–123. [5] E. Cances, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational quantum chemistry: a primer Handbook of numerical analysis, in: P. Ciarlet, C. Le Bris (Eds.), Handbook of Numerical Analysis, vol. 10, 2003, pp. 3–270. [6] A. Chatterjee, An introduction to the proper orthogonal decomposition, Current Science 78 (7) (2000) 808–817. [7] L. Chamoin, P. Ladevèze, Bounds on history-dependent or independent local quantities in viscoelasticity problems solved by approximate methods, International Journal for Numerical Methods in Engineering 71 (12) (2007) 1387–1411. [8] L. Chamoin, P. Ladevèze, A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems, Computer Methods in Applied Mechanics and Engineering 197 (9–12) (2008) 994–1014. [9] F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, F. Boust, Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization, Computer Methods in Applied Mechanics and Engineering 197 (2008) 400–413. [10] F. Chinesta, P. Ladevèze, A. Ammar, E. Cueto, A. Nouy, Proper generalized decomposition in extreme simulations: towards a change of paradigm in computational mechanics?, IACM Expressions 26/09 (2009) 2–7 [11] F. Chinesta, A. Ammar, E. Cueto, Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models, Archives of Computational Methods in Engineering 17 (4) (2010) 327–350. [12] L. Cordier, M. Bergmann, Réduction dynamique par décomposition orthogonale aux valeurs propres (POD). OCET, 2006. [13] R. Cottereau, P. Diez, A. Huerta, Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method, Computational Mechanics 44 (4) (2009) 533–547. [14] G.H. Golub, C.F.V. Loan, Matrix Computations, Third ed., Johns Hopkins University Press, Baltimore, 1996. [15] M.A. Grepl, A.T. Patera, A posteriori error bounds for reduced-basis approximation of parametrized parabolic partial differential equations, ESAIM-Mathematical Modelling and Numerical Analysis (M2AN) 39 (1) (2005) 157–181. [16] M.D. Gunzburger, J.S. Peterson, J.N. Shadid, Reduced-order modeling of timedependent PDEs with multiple parameters in the boundary data, Computer Methods in Applied Mechanics and Engineering 196 (4–6) (2007) 1030–1047. [17] I. Jolliffe, Principal Component Analysis, Springer Verlag, New York, 1986. [18] K. Karhunen, Zur spektraltheorie stochastisher prozesse, Annales Academiae Scientiarum Fennicae Series (1946) 34. [19] P. Krysl, S. Lall, J.E. Marsden, Dimensional model reduction in nonlinear finite element dynamics of solids and structures, International Journal for Numerical Methods in Engineering 51 (2001) 479–504. [20] K. Kunish, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numerische Mathematik 90 (1) (2001) 148–177. [21] P. Ladevèze, The large time increment method for the analysis of structures with nonlinear constitutive relation described by internal variables, Comptes Rendus Académie des Sciences, Paris 309 (II) (1989) 1095–1099 (in french). [22] P. Ladevèze, Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation, Springer, 1998.

P. Ladevèze, L. Chamoin / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2032–2047 [23] P. Ladevèze, Strict upper error bounds for calculated outputs of interest in computational structural mechanics, Computational Mechanics 42 (2) (2008) 271–286. [24] P. Ladevèze, D. Leguillon, Error estimate procedure in the finite element method and application, SIAM Journal of Numerical Analysis 20 (3) (1983) 485–509. [25] P. Ladevèze, P. Marin, J.P. Pelle, J.L. Gastine, Accuracy and optimal meshes in finite element computation for nearly incompressible materials, Computer Methods in Applied Mechanics and Engineering 94 (3) (1992) 303–314. [26] P. Ladevèze, E.A.W. Maunder, A general method for recovering equilibrating element tractions, Computer Methods in Applied Mechanics and Engineering 137 (1996) 111–151. [27] P. Ladevèze, J.-P. Pelle, Mastering Calculations in Linear and Nonlinear Mechanics, Springer, NY, 2004. [28] P. Ladevèze, J.C. Passieux, D. Néron, The LATIN multiscale computational method and the proper generalized decomposition, Computer Methods in Applied Mechanics and Engineering 199 (21) (2009) 1287–1296. [29] P. Ladevèze, L. Chamoin, E. Florentin, A new non-intrusive technique for the construction of admissible stress fields in model verification, Computer Methods in Applied Mechanics and Engineering 199 (9–12) (2010) 766–777. [30] P. Ladevèze, L. Chamoin, Calculation of strict error bounds for finite element approximations of non-linear pointwise quantities of interest, International Journal for Numerical Methods in Engineering 84 (2010) 1638–1664. [31] D. Néron, P. Ladevèze, Proper generalized decomposition for multiscale and multiphysics problems, Archives of Computational Methods in Engineering 17 (4) (2010) 351–372. [32] A. Nouy, A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations, Computer Methods in Applied Mechanics and Engineering 196 (45–48) (2007) 4521–4537.

2047

[33] A. Nouy, Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms, Computer Methods in Applied Mechanics and Engineering 197 (2008) 4718–4736. [34] A. Nouy, A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations, Computer Methods in Applied Mechanics and Engineering 199 (23–24) (2010) 1603– 1626. [35] A. Nouy, P. Ladevèze, Multiscale computational strategy with time and space homogenization: a radial-type approximation technique for solving microproblems, International Journal for Multiscale Computational Engineering 2 (4) (2004) 557–574. [36] N. Pares, P. Diez, A. Huerta, Subdomain-based flux-free a posteriori error estimators, Computer Methods in Applied Mechanics and Engineering 195 (4– 6) (2006) 297–323. [37] J.C. Passieux, P. Ladevèze, D. Néron, A scalable time–space multiscale domain decomposition method: adaptive time scales separation, Computational Mechanics 46 (4) (2010) 621–633. [38] F. Pled, L. Chamoin, P. Ladevèze, On the techniques for constructing admissible stress fields in model verification: performances on engineering examples, International Journal for Numerical Methods in Engineering (2010). online. [39] D.V. Rovas, L. Machiels, Y. Maday, Reduced-basis output bound methods for parabolic problems, IMA Journal of Numerical Analysis 26 (3) (2006) 423– 445. [40] D. Ryckelynck, A priori hyperreduction method: an adaptive approach, Journal of Computational Physics 202 (2005) 346–366. [41] F. Schmidt, N. Pirc, M. Mongeau, F. Chinesta, Efficient mold danscooling optimization by using model reduction, Technical Report, University of Toulouse, 2010.