Skip to main content

Certified reduced basis approximation for parametrized partial differential equations and applications

Abstract

Reduction strategies, such as model order reduction (MOR) or reduced basis (RB) methods, in scientific computing may become crucial in applications of increasing complexity. In this paper we review the reduced basis methods (built upon a high-fidelity ‘truth’ finite element approximation) for a rapid and reliable approximation of parametrized partial differential equations, and comment on their potential impact on applications of industrial interest. The essential ingredients of RB methodology are: a Galerkin projection onto a low-dimensional space of basis functions properly selected, an affine parametric dependence enabling to perform a competitive Offline-Online splitting in the computational procedure, and a rigorous a posteriori error estimation used for both the basis selection and the certification of the solution. The combination of these three factors yields substantial computational savings which are at the basis of an efficient model order reduction, ideally suited for real-time simulation and many-query contexts (for example, optimization, control or parameter identification). After a brief excursus on the methodology, we focus on linear elliptic and parabolic problems, discussing some extensions to more general classes of problems and several perspectives of the ongoing research. We present some results from applications dealing with heat and mass transfer, conduction-convection phenomena, and thermal treatments.

1 Introduction and motivation

Although the increasing computer power makes the numerical solution problems of very large dimensions that model complex phenomena essential, a computational reduction is still determinant whenever interested in real-time simulations and/or repeated output evaluations for different values of some inputs of interest. For a general introduction on the development of the reduced basis methods we refer to [13].

In this work we review the reduced basis (RB) approximation and a posteriori error estimation methods for the rapid and reliable evaluation of engineering outputs associated with elliptic and parabolic parametrized partial differential equations (PDEs). In particular, we consider a (say, single) output of interest s(μ)R expressed as a functional of a field variable u(μ) that is the solution of a partial differential equation, parametrized with respect to the input parameter p-vector μ; the input parameter domain - that is, the set of all possible inputs - is a subset D of R p . The input-parameter vector typically characterizes physical properties and material, geometrical configuration, or even boundary conditions and force fields or sources. The outputs of interest are physical quantities or indexes used to measure and assess the behavior of a system, that is, related to fields variables or fluxes, as for example, domain or boundary averages of the field variables, or other quantities such as energies, drag forces, flow rates, and so on. For the sake of simplicity, we consider throughout the paper the case of a linear output of a field variable, that is, s(μ)=l(u(μ)) for a suitable linear operator l(). Finally, the field variablesu(μ) that link the input parameters to the output depend on the selected PDE models and may represent temperature or concentration, displacements, potential functions, distribution functions, velocity or pressure. We thus arrive at an input-output relationship μs(μ), whose evaluation requires the solution of a parametrized PDE.

The reduced basis methodology we recall in this paper is motivated by, and applied within two particular contexts: the real-time context (for example, in-the-field robust parameter-estimation, or nondestructive evaluation); and the many-query context (for example, design or shape optimization, optimal control or multi-model/scale simulation). Both are crucial in view of more widespread application of numerical methods for PDEs in engineering practice and more specific industrial processes. They also feature a remarkable challenge to classical numerical techniques, such as - but not limited to - the finite element (FE) method; in fact, classical FE approximations may require big computational efforts (and also data/memory management) when the dimension N of the discretisation space becomes large. This makes unaffordable both real-time and many-query simulations: hence, looking also for computational efficiency in numerical methods becomes mandatory. The real-time and many-query contexts are often much better served by a model reduction technique such as the reduced basis approximations and associated a posteriori error bound estimation revised in this work. We note, however, that the RB methods do not replace, but rather build upon and are measured - as regards accuracy - relative to, a finite element model: the reduced basis approximates not the exact solution but rather a ‘given’ finite element discretization of (typically) very large dimension N, indicated as a high-fidelity truth approximation. In short, we promote an algorithmic collaboration rather than a computational competition between RB and FE methods.

In this paper we shall focus on the case of linear functional outputs of affinely parametrized linear elliptic and parabolic coercive partial differential equations. This kind of problems - relatively simple, yet relevant to many important applications in transport (for example, steady/unsteady conduction, convection-diffusion), mass transfer, and more generally in continuum mechanics - proves a convenient expository vehicle for the methodology, with the aim of stressing on the potential impact on possible industrial applications, dealing with optimization for devices and/or processes, diagnosis, control.

We provide here a short table of contents for the remainder of this review paper. For a wider framework on the position occupied by reduced basis method compared with other reduced order modelling (ROM) techniques and their current developments and trends, see [1]. After a brief historical excursus, we present in Section 2 the state of the art of the reduced basis method, presenting the essential components of this approach. We describe the affine linear elliptic and parabolic coercive settings in Section 3, discussing briefly admissible classes of piecewise-affine geometry and coefficients. In Sections 4 and 5 we present the essential components of the reduced basis method: RB Galerkin projection and optimality; greedy sampling procedures; an Offline-Online computational stratagem. In Section 6 we recall rigorous and relatively sharp a posteriori error bounds for RB approximations of field variables and outputs of interest. In Section 7 we briefly discuss several extensions of the methodology to more general and difficult classes of problems and applications, while in Section 8 we introduce three ‘working examples’ which shall serve to illustrate the RB formulation and its potential. In the last Section 9 we provide some future perspectives.

Although this paper focuses only on the affine linear elliptic and parabolic coercive cases - in order to allow to catch all the main ingredients - the reduced basis approximation and associated a posteriori error estimation methodology is much more general; nevertheless, many problems can successfully be faced in the even simplest affine case.

2 State of the art of the methodology

In this section we briefly review the current landscape starting from a brief historical excursus, introduce the essential RB ingredients and provide several references for further inquiry.

2.1 Computational opportunities and collaborations

The development of the reduced basis methodology can be viewed as a response to the issues described before, to address a significative computational reduction and improvement in computational performances. However, the parametric real-time and many-query contexts represent also computational opportunities, since an important role in the RB paradigm and computational stratagem is played by the parametric setting. In particular:

(i) Our attention is restricted to a typically smooth and rather low-dimensional parametrically induced manifold M, spanned by the set of fields engendered as the input varies over the parameter domain: for example, in the elliptic case

M={u(μ)X:μD},

where X is a suitable functional space. Clearly, generic approximation spaces are unnecessarily rich and hence unnecessarily expensive within the parametric framework. Our approach is premised upon a classical finite element method ‘truth approximation’ space X N X of (typically very large) dimension N; the RB method consists in a low-order approximation of the ‘truth’ manifold M N (see Figure 1) given by

M N ={ u N (μ) X N :μD}.
(1)

Several classical RB proposals focus on the truth manifold M N ; much of what we present shall be relevant to any of these reduced basis spaces/approximations.

Figure 1
figure 1

In the case of a single parameter, the parametrically induced manifold M N X N is a one-dimensional filament; the bullets represent the FE solutions used as basis functions. Indeed, the red dotted line denotes all the possible RB solutions, obtained as combinations of the basis functions.

(ii) Under suitable assumptions, the parametric setting enables to decouple the computational effort in two stages: a very extensive (parameter independent) pre-processing performed Offline once that prepares the way for subsequent very inexpensive calculations performed Online for each new input-output evaluation required. In the real-time or many-query contexts, where the goal is to achieve a very low marginal cost per input-output evaluation, we can accept an increased ‘Offline’ cost - not tolerable for a single or few evaluations - in exchange for greatly decreased ‘Online’ cost for each new/additional input-output evaluation.

2.2 A brief historical path

Reduced Basis discretization is, in brief, a Galerkin projection on an N-dimensional approximation space that focuses on the parametrically induced manifold M N . We restrict the attention to the Lagrange reduced basis spaces, which are based on the use of ‘snapshot’ FE solutions of the PDEs, corresponding to certain (properly selected) parameter values, as global approximation basis functions previously computed and stored; other possible approaches, such as Taylor [4] or Hermite spaces [5], take into account also partial derivatives of these basis solutions.

Initial ideas grew out of two related research topics dealing with linear/nonlinear structural analysis in the late 70’s: the need for more effective many-query design evaluation and more efficient parameter continuation methods [68]. The first work presented in these early somewhat domain-specific contexts were soon extended to (i) general finite-dimensional systems as well as certain classes of ODEs/PDEs [912], and (ii) a variety of different reduced basis approximation spaces - in particular Taylor and Lagrange and more recently Hermite expansions. The next decade saw further expansion into different applications and classes of equations, such as fluid dynamics and, more specifically, the incompressible Navier-Stokes equations [1316].

However, in these early methods, the approximation spaces tended to be rather local and typically low-dimensional in parameter (often a single physical parameter), due also to the absence of a posteriori error estimators and effective sampling procedures. It is clear that in higher-dimensional parameter domains the ad hoc reduced basis predictions ‘far’ from any sample points can not necessarily be trusted, and hence a posteriori error estimators combined with efficient parametric space exploration techniques are crucial to guarantee reliability, accuracy and efficiency.

Much current effort in the last ten years in the RB framework has thus been devoted to the development of (i) a posteriori error estimation procedures - and in particular rigorous error bounds for outputs of interest - and (ii) effective sampling strategies, in particular for higher dimensional parameter domains [17, 18]. The a posteriori error bounds are of course mandatory for rigorous certification of any particular RB Online output prediction. Not only, an a priori theory for RB approximations is also available, dealing with a class of single parameter coercive problems [19] and more recently extended also to the multi-parameter case [20].

However, the error estimators also play an important role in effective (greedy) sampling procedures [1, 18]: they allow us to explore efficiently the parameter domain in search of most representative ‘snapshots’, and to determine when we have just enough basis functions. We note here that greedy sampling methods are similar in objective to, but very different in approach from, more well-known Proper Orthogonal Decomposition (POD) methods [21]; the former are usually applied in the (multi-dimensional) parameter domain, while the latter are most often applied in the (one-dimensional) temporal domain. An efficient combination of the two techniques greedy-POD in parameter-time has been proposed [22, 23] and is currently used for the treatment of parabolic problems [24]; see Section 5.2.

Concerning instead computational reduction and decoupling stratagems, early work on the RB method certainly exploited - but not fully - the Offline-Online procedure. In particular, early RB approaches did not fully decouple the underlying FE approximation - with space of very high dimension N - from the subsequent reduced basis projection and evaluation - of very low dimension N. Consequently, the computational savings provided by RB treatment (relative to classical FE evaluation) were typically rather modest [4, 7, 10]despite the very small size of the RB linear systems. Much work has thus been devoted to full decoupling of the FE and RB spaces through Offline-Online procedures, above all concerning the efficient a posteriori error estimation: the complexity of the Offline stage depends on N; the complexity of the Online stage - solution and/or output evaluation for a new value of μ - depends only on N and Q (used to measure the parametric complexity of the operator and data, as defined below). In this way, in the Online stage we can reach the accuracy of a high-fidelity FE model but at the very low cost of a reduced-order model.

In the context of affine parameter dependence, in which the operator is expressible as the sum of Q products of parameter-dependent functions and parameter-independent operators (see Section 3), the Offline-Online idea is quite self-apparent and has been naturally exploited [16, 25] and extended more recently in order to obtain efficient a posteriori error estimation. In the case of nonaffine parameter dependence the development of Offline-Online strategies is even more challenging and only in the last few years effective procedures have been studied and applied [26] to allow more complex parametrizations; clearly, Offline-Online procedures are an important element both in the real-time and the many-query contexts. We recall that also historically [9] RB methods have been built upon, and measured (as regards accuracy) relative to, underlying finite element discretizations. However, spectral element approaches [27, 28], finite volume [22], and other traditional discretization methods may be considered too.

2.3 Essential RB components

The essential components of the reduced basis method, which will be analyzed in detail along the next sections, can be summarized as below.

  1. (i)

    Rapidly convergent global reduced basis (RB) approximations - (Galerkin) projection onto a (Lagrange) space X N N spanned by solution of the governing partial differential equation at N (optimally) selected points S N in the parameter set D. Typically, N will be small, as we focus attention on the (smooth) low-dimensional parametrically-induced manifold of interest. The RB approximations to the field variable and output will be denoted u N (μ) and s N (μ), respectively.

  2. (ii)

    Rigorous a posteriori error estimation procedures that provide inexpensive yet sharp bounds for the error in the RB field-variable approximation, u N (μ), and output(s) approximation, s N (μ). Our error indicators are rigorous upper bounds for the error (relative to the FE truth field u N (μ) and output s N (μ)=l( u N (μ)) approximation, respectively) for all μD and for all N. Error estimators are also employed during the greedy procedure [1] to construct optimal RB samples/spaces ensuring an efficient and well-conditioned RB approximation.

  3. (iii)

    Offline/Online computational procedures - decomposition stratagems which decouple the generation and projection stages of the RB approximation: very extensive (μ-independent) pre-processing performed Offline once that prepares the way for subsequent inexpensive calculations performed Online for each new input-output evaluation required.

3 Elliptic & parabolic parametric PDEs

We introduce the formulation of affinely parametrized linear elliptic/parabolic coercive problems; the methodology addressed in this work is intended for heat and mass convection/conduction problems. For the sake of simplicity, we consider only compliant outputs, referring to Section 7 for the treatment of general (non-compliant) outputs and the extensions to other classes of equations.

3.1 Elliptic coercive parametric PDEs

We consider the following problem: Given μD R p , evaluate the output of interest

s(μ)=(u(μ)),
(2)

where u(μ)X(Ω) satisfies

a(u(μ),v;μ)=f(v),vX(Ω).
(3)

Ω is a suitably regular bounded spatial domain in R d (for d=2 or 3), X=X(Ω) is a suitable Hilbert space; a(,;μ) and f(;μ) are the bilinear and linear forms, respectively, associated with the PDE. We shall exclusively consider second-order PDEs, and hence ( H 0 1 ( Ω ) ) ν X(Ω) ( H 1 ( Ω ) ) ν , where ν=1 (respectively, ν=d) for a scalar (respectively, vector) field; here L 2 (Ω) is the space of square integrable functions over Ω, H 1 (Ω)={v|v L 2 (Ω),v ( L 2 ( Ω ) ) d }, H 0 1 (Ω)={v H 1 (Ω):v | Ω =0}. We denote by ( , ) X the inner product associated with the Hilbert space X, whose induced norm X = ( , ) X is equivalent to the usual ( H 1 ( Ω ) ) ν norm. Similarly, (,) and denote the L 2 (Ω) inner product and induced norm, respectively.

We shall assume that the bilinear form a(,;μ):X×XR is continuous and coercive over X for all μ in D, that is,

γ ( μ ) : = sup w X sup v X a ( w , v ; μ ) w X v X < + , μ D ,
(4)
α 0 > 0 : α ( μ ) : = inf w X a ( w , w ; μ ) w X 2 α 0 , μ D .
(5)

Finally, f() and () are linear continuous functionals over X; we assume - solely for simplicity of exposition - that f and are independent of μ. Under these standard hypotheses on a and f, (3) admits a unique solution. For the sake of simplicity,1 we shall further presume for most of this paper that we are in ‘compliance’ case [1]. In particular, we assume that (i) a is symmetric - a(w,v;μ)=a(v,w;μ)w,vXμD - and furthermore (ii) =f. We shall make one last assumption, crucial to Offline-Online procedures, by assuming that the parametric bilinear form a is ‘affine’ in the parameter μ: for some finite Q a a(,;μ) can be expressed as

a(w,v;μ)= q = 1 Q Θ a q (μ) a q (w,v),
(6)

for given smooth μ-dependent functions Θ a q 1q Q a , and continuous μ-independent bilinear forms a q 1q Q a (in the compliant case the a q are additionally symmetric). Under this assumption, M N defined by (1) lies on a smooth p-dimensional manifold in X N . In actual practice, f may also depend affinely on the parameter: in this case, f(v;μ) may be expressed as a sum of Q f products of μ-dependent functions and μ-independent X-bounded linear forms. As we shall see in the following, the assumption of affine parameter dependence is broadly relevant to many instances of both property and geometry parametric variation. Nevertheless, this assumption may be relaxed [26], as detailed in Section 7.

3.2 Parabolic coercive parametric PDEs

We also consider the following parabolic model problem: Given μD R p , evaluate the output of interest

s(t;μ)=(u(t;μ)),tI=[0, t f ],
(7)

where u(μ) C 0 (I; L 2 (Ω)) L 2 (I;X) is such that

m( u t (t;μ),v;μ)+a(u(t;μ),v;μ)=g(t)f(v),vX,tI,
(8)

subject to initial condition u(0;μ)= u 0 L 2 (Ω); g(t) L 2 (I) is called control function. In addition to the previous assumptions (4)-(6), we shall assume that a(,;μ) - which represents convection and diffusion - is time-invariant; moreover, m(,;μ) - which represents ‘mass’ or inertia - is assumed to be time-invariant, symmetric, and continuous and coercive over L 2 (Ω), with coercivity constant

σ 0 :σ(μ):= inf w X m ( w , w ; μ ) w X 2 σ 0 ,μD.
(9)

Finally, we assume that also m(,;μ) is ‘affine in parameter’, that is, it can be expressed as

m(w,v;μ)= q = 1 Q m Θ m q (μ) m q (w,v),
(10)

for given smooth parameter-dependent functions Θ m q , 1 q Q m , and continuous parameter-independent bilinear forms m q , 1 q Q m , for suitable integer Q m .

3.3 Parametrized formulation

We now describe a general class - through not the most general one - of elliptic and parabolic problems which honors the hypotheses previously introduced; for simplicity we consider a scalar field (ν=1) in two space dimension (d=2). We shall first define an ‘original’ problem (subscript o), posed over the parameter-dependent domain Ω o = Ω o (μ); we denote X o (μ) a suitable Hilbert space defined on Ω o (μ). In the elliptic case, the original problem reads as follows: Given μD, evaluate

s o (μ)= l o ( u o (μ)),

where u o (μ) X o (μ) satisfies

a o ( u o (μ),v;μ)= f o (v),v X o (μ).

In the same way, for the parabolic case we have: Given μD, evaluate

s o (t;μ)= l o ( u o (t;μ)),

being u o (μ) C 0 (I; L 2 (Ω)) L 2 (I; X o (μ)) such that

m o ( u o t (t;μ),v;μ)+ a o ( u o (t;μ),v;μ)=g(t) f o (v),v X o (μ),tI.

The RB framework requires a reference (μ-independent) domain Ω in order to compare, and combine, FE solutions that would be otherwise computed on different domains and grids. For this reason, we need to map Ω o (μ) to a reference domain Ω= Ω o ( μ ref ) μ ref D, in order to get the ‘transformed’ problem (2)-(3) or (7)-(8) - which is the point of departure of RB approach - for elliptic and parabolic case, respectively. The reference domain Ω is thus related to the original domain Ω o (μ) through a parametric mapping T(;μ), such that Ω o (μ)=T(Ω;μ). It remains to place some restrictions on both the geometry (that is, on Ω o (μ)) and the operators (that is, a o m o f o l o ) such that (upon mapping) the transformed problem satisfies the hypotheses introduced above - in particular, the affinity assumption (6), (10). To this aim, a domain decomposition is useful [1].

We first consider the class of admissible geometries. In order to build a parametric mapping related to geometrical properties, we introduce a conforming domain decomposition of Ω o (μ),

Ω o (μ)= l = 1 L dom Ω o l (μ),
(11)

consisting of mutually nonoverlapping open subdomains Ω o l (μ), s.t. Ω o l (μ) Ω o l (μ)=, 1l< l L dom . If related to geometrical properties used as input parameters (for example, lengths, thicknesses, diameters or angles) the definition of parametric mappings can be done in a quite intuitive fashion.2 In the following we will identify Ω l = Ω o l ( μ ref ), 1l L dom , and denote (11) the ‘RB triangulation’; it will play an important role in the generation of the affine representation (6), (10). Hence, original and reference subdomains must be linked via a mapping T(;μ): Ω l Ω o l (μ), 1l L dom , such that

Ω o l (μ)= T l ( Ω l ;μ),1l L dom ;
(12)

these maps must be individually bijective, collectively continuous, and such that T l (x,μ)= T l (x;μ), x Ω l Ω l , for 1l< l L dom .

Here we consider the affine case, where the transformation is given, for any μD and x Ω l , by

T i l (x,μ)= C i l (μ)+ j = 1 d G i j l (μ) x j ,1id,
(13)

for given translation vectors C l :D R d and linear transformation matrices G l :D R d × d . The linear transformation matrices can effect rotation, scaling and/or shear and have to be invertible. The associated Jacobians can be defined as J l (μ)=|det( G l (μ))|, 1l L dom .

We next introduce the class of admissible operators. We may consider the associated bilinear forms

a o (w,v;μ)= l = 1 L dom Ω o l ( μ ) [ w x o 1 | w x o 2 | w ] K o , l (μ) [ v x o 1 v x o 2 v ] ,
(14)

where K o , l :D R 3 × 3 , 1l L dom , are prescribed coefficients.3 In the parabolic case, we also may consider

m o (w,v;μ)= l = 1 L dom Ω o l ( μ ) w M o , l (μ)v,
(15)

where M o , l :DR represents the identity operator. Similarly, we require that f o () and l o () are written as

f o (v)= l = 1 L dom Ω o l ( μ ) F o , l (μ)v, l o (v)= l = 1 L dom Ω o l ( μ ) L o , l (μ)v,

where F o , l :DR and L o , l :DR, for 1l L dom , are prescribed coefficients. By identifying u(μ)= u o (μ)T(;μ) in the elliptic case (resp. u(t;μ)= u o (t;μ)T(;μ)t>0 in the parabolic case), and tracing (14) back on the reference domain Ω by the mapping T(;μ), it follows that the transformed bilinear form a(,;μ) can be expressed as

a(w,v;μ)= l = 1 L dom Ω l [ w x 1 | w x 2 | w ] K l (μ) [ v x 1 v x 2 v ] ,
(16)

where K l :D R 3 × 3 , 1l L dom , is a parametrized tensor given by

K l (μ)= J l (μ) G l (μ) K o , l (μ) ( G l ( μ ) ) T

and G l :D R 3 × 3 is given by

G l (μ)=( ( G l ( μ ) ) 1 0 0 1 ),1l L dom .

In the same way, the transformed bilinear form m(,;μ) can be expressed as

m(w,v;μ)= l = 1 L dom Ω l w M l (μ)v,
(17)

where M l :DR, 1l L dom , M l (μ)= J l (μ) M o , l (μ). The transformed linear forms can be expressed similarly as

f(v)= l = 1 L dom Ω F l (μ)v,l(v)= l = 1 L dom Ω L l (μ)v,

where F l :DR and L l :DR are given by F l (μ)= J l (μ) M o , l (μ), L l = J l (μ) L o , l (μ), for 1l L dom . Hence, the original problem has been reformulated on a reference configuration, resulting in a parametrized problem where the effect of geometry variations is traced back onto its parametrized transformation tensors. The affine formulation (6) (resp. (6) and (10)) can then be derived by simply expanding the expression (16) (and (17)) in terms of the subdomains Ω l and the different entries of K i j l . This results, for example, in

a(w,v;μ)= K 1 11 (μ) Ω 1 w x 1 v x 1 + K 1 12 (μ) Ω 1 w x 1 v x 2 +.

The affine representation is now clear: for each term in (18) the (parameter-independent) integral represents a q (w,v), while the (parameter-dependent) prefactor represents Θ q (μ); the bilinear form m admits a similar treatment. The process by which we map this original problem to the transformed problem can be largely automated [1]. There are many ways in which we can relax the given assumptions and thus treat an even broader class of problems; for example, we may consider ‘elliptical’ or ‘curvy’ triangular subdomains [1]; we may consider non-time-invariant bilinear forms a and m; we may consider coefficient functions K M which are polynomial in the spatial coordinate (or more generally approximated by the Empirical Interpolation Method [26]). Some generalizations will be addressed in Section 7 and can be pursued by modification of the method presented in Section 4: in general, increased complexity in geometry and operator will result in more terms in affine expansions - larger - with a corresponding increase in the reduced basis (Online) computational costs.

4 The reduced basis method

We discuss in this section all the details related to the construction of the reduced basis approximation in both the elliptic and the parabolic case, for rapid and reliable prediction of engineering outputs associated with parametrized PDEs.

4.1 Elliptic case

We assume that we are given a FE approximation space X N of (typically very large) dimension N. Hence, the FE discretization of problem (2)-(3) [29, 30] is as follows: given μD, evaluate

s N (μ)=( u N (μ)),
(18)

where u N (μ) X N satisfies

a( u N (μ),v;μ)=f(v),v X N .
(19)

We then introduce, given a positive integer N max , an associated sequence of (what shall ultimately be reduced basis) approximation spaces: for N=1,, N max X N N is a N-dimensional subspace of X N ; we further suppose that they are nested (or hierarchical), that is, X 1 N X 2 N X N max N X N ; this condition is fundamental in ensuring (memory) efficiency of the resulting RB approximation. We recall from Section 2 that there are several classical RB proposals - Taylor, Lagrange, and Hermite spaces - as well as many different approaches, such as POD spaces. Even if we focus on Lagrange RB spaces, much of what is presented in this paper - in particular, concerning the discrete formulation, Offline-Online procedures and a posteriori error estimation - shall be relevant to any of these RB spaces/approximations, even if they are not of immediate application in industrial problems (where we want to preserve the Offline-Online procedure and hierarchical spaces).

In order to define a (hierarchical) sequence of Lagrange spaces X N N , 1N N max , we first introduce a ‘master set’ of properly selected parameter points μ n D, 1n N max . We then define, for given N{1,, N max }, the Lagrange parameter samples

S N ={ μ 1 ,, μ N },
(20)

and associated Lagrange RB spaces

X N N =span{ u N ( μ n ),1nN};
(21)

the u N ( μ n ), 1n N max , are often referred to as ‘(retained) snapshots’ of the parametric manifold M N and are obtained by solving the FE problem (19) for μ n , 1n N max . It is clear that, if indeed the manifold is low-dimensional and smooth, then we would expect to well approximate any member of the manifold - any solution u N (μ) for some μ in D - in terms of relatively few retained snapshots. However, we must ensure that we can choose a good combination of the available retained snapshots; represent the retained snapshots in a stable RB basis, efficiently obtain the associated RB basis coefficients; and finally choose the retained snapshots (that is, the sample S N max ) in an optimal way. The sampling strategy used to build the set S N will be discussed in Section 5.

4.1.1 Galerkin projection

For our particular class of equations, Galerkin projection is arguably the best approach. Given μD, evaluate (recalling the compliance assumption)

s N N (μ)=f( u N N (μ)),
(22)

where u N N (μ) X N N X N (or more precisely, u X N N N (μ) X N N ) satisfies

a( u N N (μ),v;μ)=f(v),v X N N .
(23)

We immediately obtain the classical optimality result in the energy norm:4

(24)

in the energy norm, the Galerkin procedure automatically selects the best combination of snapshots; moreover, we have that

(25)

that is, the output converges as the ‘square’ of the energy error. Although this latter result depends critically on the compliance assumption, extension via adjoint approximations to the non-compliant case is possible; we discuss this further in Section 7.

We now consider the discrete equations associated with the Galerkin approximation (23). First of all, we apply the Gram-Schmidt process with respect to the ( , ) X inner product to snapshots u N ( μ n ), 1n N max , to obtain mutually ( , ) X -orthonormal basis functions ζ n N , 1n N max . Then, the RB solution can be expressed as:

u N N (μ)= m = 1 N u N m N (μ) ζ m N ;
(26)

by taking v= ζ n N , 1nN, into (23) and using (26), we obtain the RB ‘stiffness’ equations

m = 1 N a( ζ m N , ζ n N ;μ) u N m N (μ)=f( ζ n N ),
(27)

for the RB coefficients u N m N (μ), 1m,nN; we can subsequently evaluate the RB output as

s N N (μ)= m = 1 N u N m N (μ)f( ζ m N ).
(28)

4.1.2 Offline-Online procedure

The system (27) is nominally of small size: a set of N linear algebraic equations in N unknowns. However, the formation of the stiffness matrix, and indeed the load vector, involves entities ζ n N , 1nN, associated with our N-dimensional FE approximation space. Fortunately, we can appeal to affine parameter dependence to construct very efficient Offline-Online procedures. In particular, system (27) can be expressed, thanks to (6), as

m = 1 N ( q = 1 Q Θ q (μ) a q ( ζ m N , ζ n N )) u N m N (μ)=f( ζ n N ),

for 1nN. The equivalent matrix form is

( q = 1 Q a Θ a q (μ) A N q ) u N (μ)= f N ,
(29)

where ( u N ( μ ) ) m = u N m N (μ) and

( A N q ) m n = a q ( ζ m N , ζ n N ), ( f N ) n =f( ζ n N ),

for 1m,n N max . Since each basis function ζ n N belongs to the FE space X N , they can be written as

ζ n N = i = 1 N ζ n i N ϕ i ,1n N max ,

that is, as a linear combination of the FE basis functions { ϕ i } i = 1 N ; therefore, the RB ‘stiffness’ matrix can be assembled once the corresponding the RB ‘stiffness’ matrix can be assembled. Then, by denoting

Z=[ ζ 1 || ζ N ] R N × N ,1N N max ,

we have that

A N q = Z T A N Z q , f N = Z T F N ,

being

( A N q ) i j = a q ( ϕ j , ϕ i ), ( F N ) i =f( ϕ i )

the structures given by the FE discretization. In this way, computation entails an expensive μ-independent Offline stage performed only once and an Online stage for any chosen parameter value μD. During the former the FE structures { A N q } q = 1 Q a and F N , as well as the snapshots { u N ( μ n ) } n = 1 N max and the corresponding orthonormal basis { ζ n N } n = 1 N max , are computed and stored. In the latter, for any given μ, all the Θ q a (μ) coefficients are evaluated, and the N×N linear system (29) is assembled and solved, in order to get the RB approximation u N N (μ). Then, the RB output approximation is obtained through the simple scalar product (37). Although being dense (rather than sparse as in the FE case), the system matrix is very small, with a size independent of the FE space dimension N.

The Online operation count is O(Q N 2 ) to get and O( N 3 ) to invert the matrix in (29), and finally O(N) to effect the inner product (37). The Online storage is - thanks to the hierarchy assumption - only O(Q N max 2 )+O( N max ): for any given N, we may extract the necessary RB N×N matrices (respectively, N-vectors) as principal submatrices (respectively, principal subvectors) of the corresponding N max × N max (respectively, N max ) quantities. The Online (marginal) cost (operation count and storage) to evaluate μ s N N (μ) is thus independent of N (see Figure 2).

Figure 2
figure 2

Comparison between the finite element and the reduced basis approximation frameworks: δ τ N is the marginal computational time for a single FE field/output approximation, while δ τ N is the marginal computational time for a single RB field/output (Online) evaluation, provided the data structures assembled and stored during the Offline stage (courtesy A. Patera, http://augustine.mit.edu).

4.2 Parabolic case

We next introduce the finite difference in time and finite element (FE) in space discretization [29, 30] of the parabolic problem (8). We first divide the time interval I into K subintervals of equal length Δt= t f /K and define t k =kΔt0kK, and define the FE approximation space X N . Hence, given μD, we look for u N k (μ)X0kK, such that

1 Δ t m ( u N k ( μ ) u N k 1 ( μ ) , v ; μ ) + a ( u N k ( μ ) , v ; μ ) = g ( t k ) f ( v ) , v X N , 1 k K ,
(30)

subject to initial condition ( u N 0 ,v)=( u 0 ,v)v X N . We then evaluate the output (recalling the compliance assumption): for 0kK

s N k (μ)=f( u N k (μ)).
(31)

We shall sometimes denote u N k (μ) as u N ( t k ;μ) and s N k (μ) as s N ( t k ;μ) to more clearly identify the discrete time levels. Under the coercivity assumption (9) of the bilinear form a(,;μ) and the smoothness assumption of Θ a , m q (μ) coefficients,

M N K ={ u N k (μ):1kK,μD},
(32)

the analogous entity of (32) in the parabolic case, lies on a smooth (p+1)-dimensional manifold in X N .

Equation (30) - Backward Euler-Galerkin discretization of (8) - shall be our point of departure: we shall presume that Δt is sufficiently small and N is sufficiently large such that u N ( t k ;μ) and s N ( t k ;μ) are effectively indistinguishable from u( t k ;μ) and s( t k ;μ), respectively. The development readily extends to Crank-Nicholson or higher order discretization; for purposes of exposition, we consider the simple Backward Euler approach.

The RB approximation in this case [24, 31] is based on RB spaces X N N 1N N max , generated by a sampling procedure which combines spatial snapshots in time and parameter - u N k (μ) - in an optimal fashion (see Section 5). Given μD, we now look for u N k (μ) X N N 0kK, such that

1 Δ t m ( u N k ( μ ) u N k 1 ( μ ) , v ; μ ) + a ( u N k ( μ ) , v ; μ ) = g ( t k ) f ( v ) , v X N N , 1 k K ,
(33)

subject to ( u N 0 (μ),v)=( u N 0 ,v)v X N N . We then evaluate the associated output: for 0kK

s N k (μ)=f( u N k (μ)).
(34)

We shall sometimes denote u N k (μ) as u N ( t k ;μ) and s N k (μ) as s N ( t k ;μ) to more clearly identify the discrete time levels. (Note that all the RB quantities should bear a N - X N N u N N k (μ) s N N k (μ) - since the RB approximation is defined in terms of the truth discretization; however, for clarity of exposition, we shall typically suppress this superscript.)

We now develop the algebraic equations associated with (33)-(34). First of all, the RB approximation u N k (μ) X N N shall be expressed as

u N k (μ)= m = 1 N u N m k (μ) ζ m N ,
(35)

given a set of mutually ( , ) X orthogonal basis functions ζ n N X N , 1n N max , and corresponding (hierarchical) RB spaces

X N =span{ ξ n ,1nN},1N N max .

By taking v= ζ n N , 1nN, into (33) and using (35), we obtain:

1 Δ t m = 1 N m ( ζ m N , ζ n N ; μ ) u N m k ( μ ) + m = 1 N a ( ζ m N , ζ n N ; μ ) u N m k ( μ ) = f ( ζ n N ) + 1 Δ t m = 1 N m ( ζ m N , ζ n N ; μ ) u N m k 1 ( μ ) ,
(36)

for the RB coefficients u N m N (μ), 1m,nN; we can subsequently evaluate the RB output as

s N k (μ)= m = 1 N u N m k (μ)f( ζ m N ).
(37)

The equivalent matrix form is

( q = 1 Q a Θ a q ( μ ) A N q + 1 Δ t q = 1 Q m Θ m q ( μ ) M N q ) u N ( μ ) = f N + 1 Δ t q = 1 Q m Θ m q ( μ ) M N q ,
(38)

where ( u N k ( μ ) ) m = u N m k (μ) and

( M N q ) m n = m q ( ζ m N , ζ n N ),1m,n N max ;

other terms are the same as in the elliptic case (see Sections 4.1.1-4.1.2). Moreover, also the RB mass terms can be computed from the FE mass terms as

M N q = Z T M N Z q ,where ( M N q ) i j = m q ( ϕ j , ϕ i ),

being { ϕ i } i = 1 N the basis of the FE space X N .

The Offline-Online procedure is now straightforward; in particular, the unsteady case is very similar to the steady case discussed before. There are a few new twists: as regards storage, we must now append to the elliptic Offline dataset an affine development for the mass matrix M N q , 1q Q m , associated with the unsteady term; as regards computational complexity, we must multiply the elliptic operation counts by K to arrive at O(K N 3 ) (in fact, O(K N 2 ) for a linear time-invariant system) for the Online operation count, where K is the number of time steps (recall that in actual practice the ‘truth’ is discrete in time). Thus, the Online evaluation of s N (μ) remains independent of N even in the unsteady case.

5 Sampling strategies

We now review two sampling strategies used for the construction of RB spaces: a greedy procedure for the elliptic case and a combined POD-greedy procedure for the parabolic case. Let us denote by Ξ a finite sample of points in D, which shall serve as surrogates for D in the calculation of errors (and error bounds) over the parameter domain.

5.1 Elliptic case

We denote the particular samples which shall serve to select the RB space - or ‘train’ the RB approximation - by Ξ train . The cardinality of Ξ train will be denoted | Ξ train |= n train . We note that although the ‘test’ samples Ξ serve primarily to understand and assess the quality of the RB approximation and a posteriori error estimators, the ‘train’ samples Ξ train serve to generate the RB approximation. The choice of n train and Ξ train thus have important Offline and Online computational implications. Moreover, let us denote ε tol a chosen tolerance for the stopping criterium of the greedy algorithm.

The greedy sampling strategy can be implemented as follows:

As we shall describe in detail in Section 6, Δ N (μ) is a sharp, (asymptotically) inexpensive a posteriori error bound for u N (μ) u X N N N (μ) X .

Roughly, at iteration N the greedy algorithm appends to the retained snapshots that particular candidate snapshot - over all candidate snapshots u N (μ)μ Ξ train - which is (predicted5 by the a posteriori error bound to be the) least well approximated by (the RB prediction associated to) X N 1 N . We refer to [32] for a general analysis of the greedy algorithm and related convergence rates.

5.2 Parabolic case

The temporal evolution case is quite different: the greedy approach [31] can encounter difficulties best treated by incorporating elements of the POD selection process [22]. Our sampling method thus combine the POD in t k - to capture the causality associated with the evolution equation - with the greedy procedure in μ [1, 18, 31] - to treat efficiently the higher dimensions and more extensive ranges of parameter variation.

To begin, we summarize the basic POD optimality property: given J elements w j X N , 1jJ, POD({ w 1 ,, w J },M) returns M<J ( , ) X -orthonormal functions { χ m ,1mM} such that the space P M =span{ χ m ,1mM} is optimal, that is,

P M = arginf Y M span { w j , 1 j J } ( 1 J j = 1 J inf v Y M w j v X 2 ) 1 / 2 ,

where Y M denotes an M-dimensional linear space.

To initiate the POD-greedy sampling procedure we must specify Ξ train , an initial sample S ={ μ 0 } and a tolerance ε tol . The algorithm depends on two suitable integers M 1 and M 2 (the criterium behind their setting is addressed later) and reads as follows:

As we shall describe in detail in Section 6, Δ N ( t k ;μ) provides a sharp inexpensive a posteriori error bound for u N ( t k ;μ) u N N ( t k ;μ) X . In practice, we exit the POD-greedy sampling procedure at N= N max N max , 0 for which a prescribed error tolerance is satisfied: to wit, we define

ε N , max = max μ Ξ train Δ N ( t K ;μ),

and terminate when ε N , max ε tol . Note, by virtue of the final re–definition, the POD-greedy generates hierarchical spaces X N , 1N N max , which is computationally very advantageous.

We choose M 1 to satisfy an internal POD error criterion based on the usual sum of eigenvalues and ε tol ; we choose M 2 M 1 to minimize duplication in the RB space. It is important to note that the POD-greedy method readily accommodates a repeat μ in successive greedy cycles - new information will always be available and old information rejected; in contrast, a pure greedy approach in both t and μ[31], though often generating good spaces, can ‘stall’. Furthermore, since the POD is conducted in only one (time) dimension - with the greedy addressing the remaining (parameter) dimensions - the procedure remains computationally feasible even for large parameter domains and very extensive parameter train samples (and in particular in higher parameter dimensions).

Concerning the computational aspects, the crucial point is that the operation count for the POD-greedy algorithm is additive and not multiplicative in n train and N; in contrast, in a pure POD approach, we would need to evaluate the FE ‘truth’ solution at the n train candidate parameter values. As a result, in the POD-greedy approach we can take n train relatively large: we can thus anticipate RB spaces and approximations that provide rapid convergence uniformly over the parameter domain.

6 A posteriori error estimation

Effective a posteriori error bounds for field variables and outputs of interest are crucial for both the efficiency and the reliability of RB approximations. As regards efficiency, a posteriori error estimation permits us to (inexpensively) control the error, as well as to minimize the computational effort by controlling the dimension of the RB space. Not only, in the greedy algorithm the application of error bounds (as surrogates for the actual error) allows significantly larger training samples Ξ train D and a better parameter space exploration at greatly reduced Offline computational cost. Concerning reliability, a posteriori error bounds allows a confident exploitation of the rapid predictive power of the RB approximation. By means of an efficient a posteriori error bound, we can make up for an error quantification for each new parameter value μ in the online stage and thus can make sure that feasibility (and safety/failure) conditions are verified.

The motivations for error estimation in turn place requirements on the error bounds. First, the error bounds must be rigorous - valid for all N and for all parameter values in the parameter domain D: non-rigorous error ‘indicators’ may suffice for adaptivity during basis assembling, but not for reliability. Second, the bounds must be reasonably sharp: an overly conservative error bound can yield inefficient approximations (N too large) or even dangerous suboptimal engineering results (unnecessary safety margins). And third, the bounds must be very efficient: the Online operation count and storage to compute the RB error bounds - the marginal average cost - must be independent of N (and commensurate with the cost associated with the RB output prediction).

6.1 Elliptic case

Let us now consider a posteriori error bounds for the field variable u N N (μ) and the output s N N (μ) in the elliptic case (22)-(23). We introduce two basic ingredients of our error bounds: the error residual relationship and coercivity lower bounds.

6.1.1 Basic ingredients

The central equation in a posteriori theory is the error residual relationship. In particular, it follows from the problem statements for u N (μ), (19), and u N N (μ), (23), that the error e(μ):= u N (μ) u N N (μ) X N satisfies

a(e(μ),v;μ)=r(v;μ),v X N .
(39)

Here r(v;μ) ( X N ) (the dual space to X N ) is the residual,

r(v;μ):=f(v;μ)a( u N N (μ),v;μ),v X N .
(40)

Indeed, (39) directly follows from the definition (40), f(v;μ)=a( u N (μ),v;μ), v X N , bilinearity of a, and the definition of e(μ). It shall prove convenient to introduce the Riesz representation of r(v;μ): e ˆ (μ) X N satisfies

( e ˆ ( μ ) , v ) X =r(v;μ),v X N .
(41)

This allows us to write the error residual equation (39) as

a(e(μ),v;μ)= ( e ˆ ( μ ) , v ) X ,v X N ,
(42)

and it follows that the dual norm of the residual can be evaluated through the Riesz representation:

r(;μ) ( X N ) := sup v X N r ( v ; μ ) v X = e ˆ (μ) X ;
(43)

this shall prove to be important for the Offline-Online stratagem developed in Section 6.1.3 below.

As a second ingredient, we need a positive, parametric lower bound function α LB N (μ) for α N (μ), the FE coercivity constant6 defined as

α N (μ)= inf w X N a ( w , w ; μ ) w X 2 ;
(44)

hence, we introduce

0< α LB N (μ) α N (μ)μD,
(45)

where the online computational time to evaluate μ α LB N (μ) has to be independent of N in order to fulfill the efficiency requirements on the error bounds articulated before. An efficient algorithm for the computation of α LB N (μ) is given by the so-called Successive Constraint Method (SCM), widely analyzed in [1, 33, 34]. Moreover, the SCM algorithm - which is based on the successive solution of suitable linear optimization problems - has been developed for the special requirements of the RB method; it thus features an efficient Offline-Online strategy, making the Online calculation complexity independent of N - a fundamental requisite.

6.1.2 Error bounds

We define error estimators for the solution in the energy norm and for the output as

Δ N (μ):= e ˆ (μ) X / ( α LB N ( μ ) ) 1 / 2 ,
(46)

and

Δ N s (μ):= e ˆ (μ) X 2 Δ N 2 (μ)/ α LB N (μ),
(47)

respectively. We next introduce the effectivities associated with these error estimators as

and

η N s (μ):= Δ N s (μ)/( s N (μ) s N N (μ)),

respectively. Clearly, the effectivities are a measure of the quality of the proposed estimator: for rigor, we shall insist upon effectivities ≥1; for sharpness, we desire effectivities as close to unity as possible. We can prove7[1] that for any N=1,, N max , the effectivities satisfy

1 η N ( μ ) γ ( μ ) α LB N ( μ ) , μ D ,
(48)
1 η N s ( μ ) γ ( μ ) α LB N ( μ ) , μ D ,
(49)
γ(μ)

being defined in (4). It is important to observe that the effectivity upper bounds, (48) and (49), are independent of N, and hence stable with respect to RB refinement.

6.1.3 Offline-Online for e ˆ (μ) X computation

The error bounds of the previous section are of no utility without an accompanying Offline-Online computational approach.

The computationally crucial component of all the error bounds of the previous section is e ˆ (μ) X , the dual norm of the residual. To develop an Offline-Online procedure we first expand the residual (40) according to (26) and (6):

r ( v ; μ ) = f ( v ) a ( n = 1 N u N n N ( μ ) ζ n N , v ; μ ) = f ( v ) n = 1 N u N n N ( μ ) a ( ζ n N , v ; μ ) = f ( v ) n = 1 N u N n N ( μ ) q = 1 Q Θ q ( μ ) a q ( ζ n N , v ) .
(50)

If we insert (50) in (41) and apply linear superposition, we obtain

( e ˆ ( μ ) , v ) X =f(v) q = 1 Q n = 1 N Θ q (μ) u N n N (μ) a q ( ζ n N ,v),

or

e ˆ (μ)=C+ q = 1 Q n = 1 N Θ q (μ) u N n N (μ) L n q ,

where ( C , v ) X =f(v), v X N , that is, C is the Riesz representation of f, and ( L n q , v ) X = a q ( ζ n N ,v), v X N , 1nN, 1qQ, that is, L n q is the Riesz representation of A n q ( X N ) defined as A n q (v)= a q ( ζ n N ,v), v X N . We denote the C, L n q , 1nN, 1qQ, as FE ‘pseudo’-solutions, that is, solutions of ‘associated’ FE Poisson problems. We thus obtain

(51)

from which we can directly calculate the requisite dual norm of the residual through (43).

The Offline-Online decomposition is now clear. In the Offline stage we form the μ-independent quantities. In particular, we compute the FE ‘pseudo’-solutions C, L n q , 1n N max , 1qQ, and store ( C , C ) X , ( C , L n q ) X , ( L n q , L n q ) X , 1n, n N max , 1q, q Q. The Offline operation count depends on N max , Q, andN.

In the Online stage, given any ‘new’ value of μ - and Θ q (μ), 1qQ, u N n N (μ), 1nN - we simply retrieve the stored quantities ( C , C ) X , ( C , L n q ) X , ( L n q , L n q ) X , 1n, n N, 1q, q Q, and then evaluate the sum (51). The Online operation count, and hence also the marginal cost, is O( Q 2 N 2 ) - and independent ofN.8

6.2 Parabolic case

In this section we deal with a posteriori error estimation in the reduced basis context for affinely parametrized parabolic coercive PDEs. As for the elliptic case, to construct the a posteriori error bounds we need two ingredients. The first ingredient is the dual norm of the residual

ε N ( t k ;μ)= sup v X N r N ( v ; t k ; μ ) v X ,1kK,
(52)

where r N (v; t k ;μ) is the residual associated with the RB approximation (33)

r N ( v ; t k ; μ ) = g ( t k ) f ( v ) 1 Δ t m ( u N k ( μ ) u N k 1 ( μ ) , v ; μ ) a ( u N k ( μ ) , v ; μ ) , v X N , 1 k K .
(53)

The second ingredient is a lower bound for the coercivity constant α N (μ), 0< α LB N (μ) α N (μ), μD.

We can now define our error bounds in terms of these two ingredients; in fact, it can readily be proven [22, 31] that for all μD and all N

(54)
| s N k ( μ ) s N k ( μ ) | Δ N s k ( μ ) , 1 k K ,
(55)

where Δ N k (μ) Δ N ( t k ;μ) and Δ N s k (μ) Δ N s ( t k ;μ) are given by

Δ N k ( μ ) = ( Δ t α LB N ( μ ) k = 1 k ε N 2 ( t k ; μ ) ) 1 / 2 ,
(56)
Δ N s k ( μ ) = ( Δ N k ( μ ) ) 2 .
(57)

(We assume for simplicity that u N 0 X N ; otherwise there will be an additional contribution to Δ N k (μ).)

Even if based on the same components as in the elliptic case, now the Construction-Evaluation procedure for the error bound is a bit more involved. The necessary computations for the Offline and Online stages - by construction rather similar to the elliptic case - are discussed in details, for example, in [24]. We consider here only the decomposition for the dual norm of the residual [31]. We first invoke duality, our RB expansion, the affine parametric dependence of a and m, and linear superposition to express

ε N 2 ( t k ; μ ) = Q N f f + n = 1 N ( q = 1 Q a Θ a q ( μ ) u N n k ( μ ) Q N n q f a + 1 Δ t q = 1 Q m Θ m q ( μ ) ϕ N n k ( μ ) Q N n q f m ) + n , n = 1 N , N ( q , q = 1 Q a , Q a Θ a q ( μ ) Θ a q ( μ ) u N n k ( μ ) u N n k ( μ ) Q N n n q q a a + n , n = 1 N , N ( + 1 ( Δ t ) 2 q , q = 1 Q m , Q m Θ m q ( μ ) Θ m q ( μ ) ϕ N n k ( μ ) ϕ N n k ( μ ) Q N n n q q m m + n , n = 1 N , N ( + 1 Δ t q , q = 1 Q a , Q m Θ a q ( μ ) Θ m q ( μ ) u N n k ( μ ) ϕ N n k ( μ ) Q N n n q q a m ) ,
(58)

for 1kK, where ϕ N n k (μ):= u N n k (μ) u N n k 1 (μ) and Q N f f = ( z f , z f ) X Q N n q f a =2 ( z n q a , z f ) X 1q Q a 1nN Q N n q f m =2 ( z n q m , z f ) X 1q Q m 1nN Q N n n q q a a = ( z n q a , z n q a ) X 1q, q Q a 1n, n N Q N n n q q a m =2 ( z n q a , z n q m ) X 1q Q a 1 q Q m 1n, n N, and Q N n n q q m m = ( z n q m , z n q m ) X 1q, q Q m 1n, n N. Here the z f z n q a z n q m are solutions to time-independent and μ-independent ‘Poisson’ problems: ( z f , v ) X =f(v)v X N ( z n q a , v ) X = a q ( ξ n ,v)v X N 1nN1q Q a , and ( z n q m , v ) X = m q ( ξ n ,v)v X N 1nN1 q Q m .

The Construction-Evaluation decomposition is now clear. In the μ-independent construction stage we find z f , z a , z m , and the inner products Q N max f f , Q N max f a , Q N max f m , Q N max a a , Q N max m m , and Q N max a m at (considerable) computational cost O( Q a Q m N max N ). In the μ-dependent Evaluation stage - performed many times - we simply perform the sum (58) from the stored inner products in O( ( 1 + Q m N + Q a N ) 2 ) operations per time step and hence O( ( 1 + Q m N + Q a N ) 2 K) operations in total. The crucial point, again, is that the cost and storage in the Evaluation phase - the marginal cost for each new value of μ - is independent of N: thus we can not only evaluate our output prediction but also our rigorous output error bound very rapidly in the parametrically interesting contexts of real-time or many-query investigation.

7 Extensions to more general problems

We now briefly discuss some extensions of the reduced basis methodology presented in Section 4 to address more general classes of problems, also to face industrial problems of a certain degree of complexity.

7.1 Non-compliant problems

For the sake of simplicity, we addressed in Section 4 the RB approximation of affinely parametrized coercive problems in the compliant case. We now consider the elliptic case and the more general non-compliant problem: given μD, find

s(μ)=(u(μ)),
(59)

where u(μ)X satisfies

a(u(μ),v;μ)=f(v),vX.
(60)

We assume that a is coercive and continuous (and affine, (6)) but not necessarily symmetric. We further assume that both and f are bounded functionals but we no longer require =f.9 Following the methodology (and the notation) addressed in Section 4, we can readily develop an a posteriori error bound for s N (μ): by standard arguments [1, 2]

| s N (μ) s N N (μ)| ( X N ) Δ N (μ),

where and Δ N (μ) is given by (46). We denote the method already illustrated as ‘primal-only’. Although for many outputs primal-only is perhaps the best approach (each additional output, and associated error bound, is a simple ‘add-on’), this approach has two deficiencies:

  1. (i)

    we loose the ‘quadratic convergence’ effect (25) for outputs (unless =f and a is symmetric);

  2. (ii)

    the effectivities Δ N s (μ)/|s(μ) s N (μ)| may be unbounded: if =f then we know, from (25), that |s(μ) s N (μ)| e ˆ (μ) X 2 and hence Δ s (μ)/|s(μ) s N (μ)|1/ e ˆ (μ) X as N, that is, the effectivity of the output error bound (47) tends to infinity as (N and) u N pr N (μ) u N (μ). We may expect similar behavior for any ‘close’ to f: the failing is that (47) does not reflect the contribution of the test space to the convergence of the output.

The introduction of RB primal-dual approximation will take care of the previous issue - and ensure a stable limit N. We thus introduce the dual problem associated to , that reads as follows: find ψ(μ)X such that

a(v,ψ(μ);μ)=(v),vX;

ψ is denoted the ‘adjoint’ or ‘dual’ field. Let us define the RB spaces for the primal and the dual problem, respectively:

X N pr N , pr = span { u N ( μ k , pr ) ζ k N , 1 k N pr } , X N du N , du = span { Ψ N ( μ k , du ) , 1 k N du } ;

for 1 N pr N pr , max 1 N du N du , max . For our purposes a single FE space suffices for both the primal and dual, even if in actual practice the FE primal and dual spaces may be different. The resulting RB approximation u N pr N X N pr N , pr Ψ N du X N du du solve

a ( u N pr N ( μ ) , v ; μ ) = f ( v ) , v X N pr N , pr , a ( v , Ψ N du N ( μ ) ; μ ) = ( v ) , v X N du N , du ;

then, the RB output can be evaluated as [35]

s N pr , N du N (μ)=( u N pr N ) r pr ( Ψ N du N ;μ),

where

r pr (v;μ)=f(v)a( u N pr N ,v;μ), r du (v;μ)=(v)a(v, Ψ N du N ;μ)

are the primal and the dual residual. In particular, in the non-compliant case, the output error bound takes the form

Δ N s (μ) r pr ( ; μ ) ( X N ) ( α LB N ( μ ) ) 1 / 2 r du ( ; μ ) ( X N ) ( α LB N ( μ ) ) 1 / 2 .
(61)

We thus recover the ‘quadratic’ output effect; note that the Offline-Online procedure is very similar to the ‘primal-only’ case, but now we need to do everything both for primal and dual; moreover, we need to evaluate both a primal and a dual residual for the a posteriori error bounds, but at a reasonable computational cost and by reusing the same computational framework built and set for the ‘primal-only’ approach. Error bounds related to the gradient of computed quantities, such as velocity and pressure in potential flows problems, have been addressed in [36]. For parabolic problems, the treatment of non-compliant outputs follows the same strategy; we only remark that the dual problem in this case shall evolve backward in time [31].

7.2 Nonaffine and noncoercive problems

In this section we introduce wider classes of problems to be treated with the reduced basis method: nonaffine problems and noncoercive problems in order to provide a general framework for the methodology. The reader interested in numerical applications may go directly to Section 8 without affecting the understanding of the subsequent sections.

7.2.1 Nonaffine problems

The assumption of affine parametric dependence - expressed by conditions (6) and (10) - is of fundamental importance in order to exploit the Offline-Online stratagem and then minimize the marginal cost associated with each input-output evaluation. However, also nonaffine problems, that is, problems in which conditions (6) and (10) are not still valid, can be efficiently treated in the reduced basis framework. In this case, we rely on the Empirical Interpolation Method (EIM) [26, 37, 38], which is an interpolation method for parametric functions based on adaptively chosen interpolation points and global shape functions.

In practice, if the problem is not affinely parametrized (for example, when the geometric transformation (12) has a more general expression than in (13), or the physical coefficients appearing in the tensor K o , l are nonaffine functions of x and μ), the parametrized tensors in (16) and (17) depend both on the parameter μ and the spatial coordinate x. In this case, the operators can not be expressed as in (18) - and ultimately as (6) and (10). Hence, we need an additional pre-processing, before the FE assembling stage, in order to recover the affinity assumption. According to EIM, each component K i j l (x,μ) is approximated by an affine expression given by

K ˜ i j l (x,μ)= k = 1 K i j l a β k i j l (μ) η k i j l (x)+ ε a i j l (x,μ);
(62)

the same approximation is set up for the components of the M i j l (x,μ) tensor in the parabolic case:

M ˜ i j l (x,μ)= k = 1 K i j l m γ k i j l (μ) ϕ k i j l (x)+ ε m i j l (x,μ).
(63)

All the coefficients β k i j l ’s, γ k i j l ’s, η k i j l ’s and ϕ k i j l ’s are efficiently computable scalar functions and the error terms are guaranteed to be under some tolerance,

ε a i j l (;μ) ε tol EIM , ε m i j l ε tol EIM ,μD.

In this way, we can identify the μ-dependent coefficients in the developments (62), (63) as the coefficients Θ a q (μ) (resp. Θ m q (μ)) in (6) and (10), that is, Θ a q (μ)= β k i j l (μ), Θ m q (μ)= γ k i j l (μ), being q a condensed index for (i,j,k,l), while the μ-independent functions will be treated as pre-factors in the integrals which give the μ-independent bilinear forms a q (w,v) (resp. m q (w,v)).

We refer the reader to [26] and [39] for details on EIM procedures for nonaffine problems. The nonaffine treatment is really important since many problems involving more complex geometrical parametrizations and/or more complex physical instances (that is, non-homogeneous or non-isotropic properties in materials) are hold by nonaffine parametric dependence.

7.2.2 Noncoercive problems

The reduced basis framework can be effectively applied also to problems involving operators which do not satisfy the (quite strict) coercivity assumption [18]; this is the case, for example, of the (Navier)-Stokes problem, where stability is in fact fulfilled in the more general sense of the inf-sup condition [29]. For the sake of simplicity, we restrict our considerations to the elliptic (scalar) case (2)-(3). We assume that the (parametrized) bilinear form a(,;μ): X 1 × X 2 R is continue and satisfies the more general inf-sup condition:

β 0 >0:β(μ):= inf w X 1 sup v X 2 a ( w , v ; μ ) w X 1 v X 2 β 0 ,μD.
(64)

In this case the finite element (and thus the subsequent reduced basis) approximation is based on a more general Petrov-Galerkin approach. Given two FE spaces X 1 , N X 1 X 2 , N X 2 , the FE approximation u N (μ) X 1 , N satisfies

a( u N (μ),v,μ)=f(v),v X 2 , N ,

and the output can be evaluated as10

s N (μ)=l( u N (μ)).

In order to have a stable FE approximation, we require that exists β 0 0 such that

β N (μ)= inf w X 1 , N sup v X 2 , N a ( w , v ; μ ) w X 1 v X 2 β 0 μD.
(65)

This condition can be reformulated in terms of the so-called inner supremizer operator T μ : X 1 , N X 2 , N

( T μ w , v ) X 2 =a(w,v;μ),w X 1 , N ,v X 2 , N ;

by Cauchy-Schwarz inequality and taking v= T μ w, we have that for any w X 1

a(w, T μ w;μ)β(μ)w X 1 T μ w X 2 .

The reduced basis approximation inherits the same Petrov-Galerkin structure; in order to guarantee its stability, we need to introduce two different spaces (note that the second is μ-dependent):

X N 1 =span{u( μ n ),1nN}, X N 2 , μ =span{ T μ u( μ n ),1nN},

for 1N N max ; then u N N (μ) X N 1 satisfies

a( u N N (μ),v;μ)=f(v),v X N 2 , μ ,

and

s N (μ)=l( u N (μ)).

If we define

β N (μ) inf w X N 1 sup v X N 2 , μ a ( w , v ; μ ) w X 1 v X 2 ,
(66)

we obtain

u N (μ) u N N (μ) X (1+ γ β N ( μ ) ) inf w N X N 1 u N (μ) w N X 1 ,

which is the analogue of (24) for noncoercive problems. In this case we can show that β N (μ) β N (μ)μD; this property, which yields the stability of the RB approximation, is not automatically satisfied by a (simple) Galerkin formulation; hence, we need to enforce this property through the introduction of a Petrov-Galerkin framework. Observe that approximation is provided by X N 1 and stability (through β N ) by X N 2 , μ .

The Offline-Online computational strategem, as well as the a posteriori error estimation, are based on the same arguments described in Section 6 for the coercive case; we remark that also the inner supremizer operator can be written in the affine form under the affinity assumption (6) on a(,;μ). In particular, from (66), we can easily prove that

Δ N (μ) e ˆ ( μ ) X β N L B ( μ ) ,

where β N L B (μ) is a lower bound of inf-sup constant (65) and can be computed by means of the same SCM procedure used for the lower bound of coercivity constants [34, 40].

An interesting case of noncoercive problems is given by Stokes problems where approximation stability is guaranteed by the fullfillment of an equivalent inf-sup stability condition on the pressure term with RB approximation spaces properly enriched [41, 42]. Error bounds can be developed in the general noncoercive framework [40] or with a penalty setting [43].

8 Working examples

Reduced basis methods have already been and may be applied in many problems of industrial interest: material sciences and linear elasticity [17, 4446], heat and mass transfer [4750], acoustics [51], potential flows [36], micro-fluid dynamics [40], electro-magnetism [52]; for examples of implementation of some worked problems in the mentioned fields, see [53, 54] for a versatile setting.

In many of these problems there are physical or engineering parameters which characterize the problem but also geometrical parameters holding a Cartesian geometrical setting; this configuration is quite typical for industrial devices, and plants and related constructions and products. More complex geometrical parametrizations will be briefly considered in Section 9, involving, for example, biomedical devices and/or aerodynamic shapes.

We discuss in this section11 three working examples of industrial interest, dealing with different heat or mass transfer problems. The first example deals with forced steady heat conduction/convection; the second application deals with a transient heat treatment, while the third one is an example of a (simple) coupled problem, dealing with the transient evolution of the concentration field near the surface of a body immersed into a fluid flowing across a channel. All numerical details concerning the construction of RB spaces and computational costs are reported in Table 1.

Table 1 Numerical details for the test cases presented.

8.1 A ‘Couette-Graetz’ conduction-convection problem

This problem deals with forced steady heat convection combined with heat conduction in a straight duct, whose walls can be kept at fixed temperature or insulated or characterized by heat exchange. The flow has an imposed temperature at the inlet and a known convection field (a Couette flow, that is, a given linear velocity profile [55]). From the engineering point of view, this example describes a class of heat transfer problems in fluidic devices with a versatile configuration. In particular, Péclet number as a measure of axial transport velocity field (modeling the physics of the problem) and the length of the non-insulated portion of the duct are only two of the possible parameters to be varied in order to extract average temperatures. Also discontinuities in Neumann boundary conditions (different heat fluxes) and thermal boundary layers are interesting phenomena to be studied.

We consider the physical domain Ω o (μ) shown in Figure 3; all lengths are non-dimensionalized with respect to a unity length h ˜ (dimensional channel width); moreover, let us denote k ˜ the dimensional (thermal) conductivity coefficient for the air flowing in the duct, ρ ˜ its density and c ˜ p the specific heat capacity under constant pressure. We introduce the (thermal) diffusion coefficient D ˜ = k ˜ / ρ ˜ c ˜ p , as well as the Péclet number, given by the ratio Pe= U ˜ h ˜ / D ˜ , being U ˜ the reference dimensional velocity for the convective field. We consider here P=2 parameters: μ 1 is the length of the non-insulated bottom portion of the duct (unity heat flux), while μ 2 represents the Péclet number; the parameter domain is given by D=[1,10]×[0.1,100].

Figure 3
figure 3

‘Couette-Graetz’ conduction-convection problem: parametrized geometry and domain boundaries.

The solution u(μ), defined as the non-dimensional temperature u(μ)=(τ τ i n )/ τ i n (where τ is the dimensional temperature, τ i n is the dimensional temperature of the air at the inflow and in the first portion of the duct) satisfies the following steady advection-diffusion equation:

{ 1 μ 2 Δ u ( μ ) + x 2 x 1 u ( μ ) = 0 in Ω o ( μ ) , 1 μ 2 u n ( μ ) = 0 on Γ 1 Γ 3 , 1 μ 2 u n ( μ ) = 1 on Γ 2 , u ( μ ) = 0 on Γ 4 Γ 5 Γ 6 ,

with summation (i,j=1,2) over repeated indices; hence, we impose the temperature at the top walls and in the ‘inflow’ zone of the duct (Γ6), while we consider an insulated wall (zero heat flux on Γ1 and Γ3) or heat exchange at a fixed rate (that is, unity on Γ2) on other boundaries. We note that the forced convection field is given by a linear velocity profile x 2 U ˜ (Couette type flow). The output of interest is the average temperature of the fluid on the non-insulated portion of the bottom wall of the duct, given by

s(μ):= T av (μ)= 1 μ 1 Γ 2 u(μ).

This problem is then mapped to the fixed reference domain Ω and discretized by piecewise linear finite elements; the dimension of the corresponding space is N=5,433. Since we are in a non-compliant case, a further dual problem has to be solved in order to obtain better output evaluations and related error bounds, see Section 7.1. In particular, we show in Figure 4 the lower bound of the coercivity constant of the bilinear form associated to our problem.

Figure 4
figure 4

‘Couette-Graetz’ conduction-convection problem: lower bound of the coercivity constant α LB N (μ) as a function of μ.

We plot in Figure 5 the convergence of the greedy algorithm for the primal and the dual problem, respectively; with a fixed tolerance ε tol = 10 2 , N pr , max =21 and N du , max =30 basis have been selected, respectively. In Figure 6 the selected parameter values S N pr for the primal and S N du for the dual problems, respectively, are shown; in each case Ξ train is a uniform random sample of size n train =1,000. Moreover, in Figure 7 some representative solutions (computed for N= N max ) for selected values of parameters are reported.

Figure 5
figure 5

‘Couette-Graetz’ conduction-convection problem: relative errors max μ Ξ train ( Δ N pr (μ)/ u N pr N (μ) X ) and max μ Ξ train ( Δ N du (μ)/ ψ N du N (μ) X ) as a function of N pr and N du for the RB approximations computed during the greedy procedure, for the primal (left) and the dual (right) problem, respectively. Here Ξ train is a uniform random sample of size n train =1,000 and the RB tolerance is ϵ tol = 10 2 .

Figure 6
figure 6

‘Couette-Graetz’ conduction-convection problem: selected parameter values S N pr for the primal (left) and S N du for the dual (right) in the parameter space.

Figure 7
figure 7

‘Couette-Graetz’ conduction-convection problem: representative solutions for μ=(1,0.1), μ=(1,100) (top), μ=(10,0.1), μ=(10,100) (bottom).

The thermal boundary layer looks very different in the four cases. In particular, higher variations of temperature, as well as large gradients along the lower wall - are remarkable for higher Péclet number, when forced convection dominates steady conduction; moreover, the standard behavior of boundary layer width - usually given by O(1/Pe) - is captured correctly. In Figure 8 the RB evaluation (for N= N max ) of the output of interest is reported as a function of the parameters, as well as the related error bound. As we can see, for low values of μ 2 (Péclet number) the dependence of the output on μ 1 (geometrical aspect) is rather modest; for high values of μ 2 , instead, the output shows a larger variations wih respect to μ 1 . In the same way, for longer/shorter channels the dependence on the Péclet number is higher/lower.

Figure 8
figure 8

‘Couette-Graetz’ conduction-convection problem: computed RB output (left) and related error bound (right) as functions of μ in the parameter space.

8.2 A transient thermal treatment problem

This problem considers a transient thermal treatment on a sectional slice of a railroad rail. Heat treatment is a method used to alter the physical, and sometimes chemical, properties of a material, which involves the use of heating or chilling, normally to extreme temperatures, to achieve a desired result such as hardening or softening of a material. Heat treatment techniques include annealing, case hardening, precipitation strengthening, tempering, and quenching. Although the most common application is metallurgical, heat treatments are also used in the manufacturing of many other materials.

We consider here P=2 parameters: μ 2 is a geometrical parameter representing the thickness of the web connecting the top and the bottom of the railroad rail slice (see Figure 9), while μ 1 denotes the non-dimensional Biot number, given by Bi h ˜ c d ˜ / k ˜ . We assume that the railroad rail slice has thermal conductivity k ˜ and we characterize the heat transfer coefficient between the railroad section and the fluid surrounding the railroad rail slice itself by a heat transfer coefficient h ˜ c ; moreover, d ˜ denotes the height of the slice of the railroad rail. The parameter domain is given by D=[0.01,10]×[0.02,0.2].

Figure 9
figure 9

Heat treatment problem: parametrized geometry and domain boundaries.

The (non-dimensional) temperature distribution is denoted u(μ) (the dependence of time is omitted for sake of simplicity) and is defined in terms of dimensional temperature as u(μ)=(τ τ init )/( τ env τ init ) where τ is the dimensional temperature, τ init the initial dimensional temperature (at t=0) and τ env is the dimensional temperature of the fluid surrounding the railroad slice (at every time) and the (asymptotic) temperature at the end of the treatment.

The governing equation for u(μ,t) is the following time-dependent linear PDE: for t[0,T]

{ u ( μ ) t Δ u ( μ ) = 0 in Ω o ( μ ) , u ( μ , t = 0 ) = 0 in Ω o ( μ ) , u n + μ 1 u ( μ ) = μ 1 g ( t ) on Ω o ( μ ) .

The inhomogeneous Robin conditions correspond to the heat exchange between the railroad rail slice section and the fluid used for the thermal treatment. Here the control input g(t) is a function of time t; the problem considers any square-integrable function for g(t). In practice, the PDE is replaced by a discrete-time (backward Euler [30]) approximation with time-steps of size Δt=0.005. Note that the final time is T=0.75 and that the number of time-steps is n t =150; the spatial discretization is made by piecewise linear finite elements, whose corresponding space dimension is N=16,737. Our output of interest is the average temperature all over the piece of railroad rail slice, given by

s(μ)= 0 T (h(t) Ω o ( μ ) u(μ))dt,

where h(t) is a function of time t; the problem considers any function (including Dirac delta) for h(t).

In Figure 10 we plot the lower bound of the coercivity constant of the bilinear form associated to the problem. As in the previous case, a further dual problem has to be solved in order to obtain better output evaluations and related error bounds. We show in Figure 11 the convergence of the greedy algorithm for the primal and the dual problem, respectively; with a fixed tolerance ε tol = 10 2 , N pr , max =22 and N du , max =6 basis have been selected, respectively.

Figure 10
figure 10

Heat treatment problem: lower bound of the coercivity constant α LB N (μ) as a function of μ.

Figure 11
figure 11

Heat treatment problem: relative errors max μ Ξ train ( Δ N pr (μ)/ u N pr N (μ) X ) and max μ Ξ train ( Δ N du (μ)/ ψ N du N (μ) X ) as a function of N pr and N du for the RB approximations computed during the greedy procedure, for the primal (left) and the dual (right) problem, respectively. Here Ξ train is a uniform random sample of size n train =1,000 and the RB tolerance is ε tol = 10 2 .

In Figures 12 and 13 some representative solutions for selected values of parameters are reported, for both t=Δt and t=T. In particular, two different heat treatments have been investigated: heating and cooling process. In the first case, we have imposed a thermal flux g(t)=10t, while in the second case g(t)=10t. We can remark more sensible variations of temperature all over the body for larger values of μ 1 (Biot number); moreover, the behavior of the temperature changes strongly between narrower and larger configurations.

Figure 12
figure 12

Heat treatment problem (heating): representative solutions for μ=(10,0.02) and μ=(10,0.2), at time t=Δt (top) and t=T (bottom).

Figure 13
figure 13

Heat treatment problem (cooling): representative solutions for μ=(10,0.02) and μ=(10,0.2), at time t=Δt (top) and t=T (bottom).

Concerning the output (67), two cases have been taken into account: a distributed (in time) output - corresponding to h(t)=1 - given by the integral of the temperature in time and space, and a concentrated (in time) output - corresponding to h(t)=δ(t) - given by the spatial integral of temperature at each timestep. In Figures 14 and 15 the RB evaluation (for N= N max ) of these two outputs of interest are reported, as well as the related error bounds. Higher values of the output are obtained with larger values of the two parameters; moreover, keeping the geometry fixed, variations w.r.t. Biot number in output values are of about one order of magnitude.

Figure 14
figure 14

Heat treatment problem: RB distributed output (left) and related error bound (right) as functions of μ in the parameter space.

Figure 15
figure 15

Heat treatment problem: RB concentrated outputs and related error bounds as functions of time t, for μ=(0.01,0.02), μ=(0.01,0.2) (left), μ=(10,0.02), μ=(10,0.2) (right).

8.3 A transient (coupled) diffusion-transport problem around a cylinder

The problem represents the transient evolution of a concentration field near the surface of a body (a two-dimensional cylinder section) immersed into a fluid flowing into a channel. The mass (for example, of oxygen or drug) can be released or absorbed through the body surface within the surrounding fluid. This is a well-known mass transfer problem in the design and sizing of substances diffusers used for many industrial, civil and, more recently, biomedical applications (drug and/or oxygen release, stent design); in the same way, it can be seen as an heat transfer problem through an heat exchanger [56].

The problem is described by the coupling of an unsteady mass (or heat) transfer phenomenon (or substance release) by diffusion (or conduction) into a body and by transport (or convection) phenomena inside the field where the fluid is flowing; the transport field is given, for example, by a potential solution (see, for example, [55]).

We consider the physical domain Ω o (μ) shown in Figure 16, non-dimensionalized with respect to R ˜ , the unit radius of the cylinder immersed in the fluid. Moreover, we denote D ˜ the dimensional mass diffusion coefficient, U ˜ a reference dimensional velocity for transport field, and we introduce the Péclet number as Pe= U ˜ R ˜ / D ˜ , while time is non-dimensionalized by the quantity R ˜ 2 / D ˜ .

Figure 16
figure 16

Diffusion-transport problem around a cylinder: parametrized geometry and domain boundaries.

In this problem the boundary segments Γ1, Γ7 are curved (all other boundary segments are straight lines) and they represent the semi-circular section of the cylinder immersed in the flow (thanks to symmetry the problem can be simplified by considering just ‘half’ configuration). The segments Γ1, Γ7 are given by the parametrization

[ x 1 x 2 ] = [ 0 0 ] + [ 1 0 0 1 ] [ 1 0 0 1 ] [ cos ( t ) sin ( t ) ] ,

where for Γ1, , for Γ7, .

We consider here only one parameter μ 1 , the Péclet number, which is given by the ratio between the transport and diffusion terms; the parameter domain is given by D=[0.1,100]. The solution is characterized by the (adimensional) concentration u(μ,t)=(c c init )/ c inlet , being c the dimensional concentration, c init the initial dimensional concentration (at t=0), and inlet the dimensional concentration imposed at the inflow (at every time step). The governing equation for u(μ,t) is the following time-dependent linear PDE: for t[0,T]

u ( μ ) t 1 μ 1 Δ u ( μ ) + ( v r sin ( θ ) v θ cos ( θ ) ) x 1 u ( μ ) + ( v r sin ( θ ) + v θ cos ( θ ) ) x 2 u ( μ ) = 0 in Ω o ( μ ) ,
u ( μ , t = 0 ) = 0 in Ω o ( μ ) , u ( μ ) = 0 on Γ 3 Γ 4 Γ 5 , 1 μ 1 u n ( μ ) = 0 on Γ 2 Γ 6 , 1 μ 1 u n ( μ ) = g ( t ) on Γ 1 Γ 7 ;

the control input g(t) is a (square-integrable) function of time t. The potential velocity field (ideal inviscid fluid) is given in polar coordinates by ( v r , v θ ), being [55]

v r = ( 1 r 0 2 r 2 ) cos ( θ ) , v θ = ( 1 + r 0 2 r 2 ) sin ( θ ) ,

where r= x 1 2 + x 2 2 r 0 = R ˜ =1 and θ=arcsin( x 2 / x 1 2 + x 2 2 ).

In practice, the PDE is replaced by a discrete-time (backward Euler) approximation with time steps of size Δt=0.01; note that the final time is T=1 and that the number of time steps is n t =100. The spatial discretization is made by piecewise linear finite elements, whose corresponding space dimension is N=13,976.

Our output of interest is the average concentration on the cylinder surface, given by

(67)

where h(t) may be a function of time t. As for the two previous cases, we deal with a non-compliant problem, for which the dual problem has to be introduced and solved. In Figure 17 we plot the lower bound of the coercivity constant of the bilinear form associated to the problem. We show in Figure 18 the convergence of the greedy algorithm for the primal and the dual problem, respectively; with a fixed tolerance ε tol = 10 2 , N pr , max =17 and N du , max =17 basis have been selected, respectively.

Figure 17
figure 17

Diffusion-transport problem around a cylinder: lower bound of the coercivity constant α LB N (μ) as a function of μ 1 .

Figure 18
figure 18

Diffusion-transport problem around a cylinder: relative errors max μ Ξ train ( Δ N pr (μ)/ u N pr N (μ) X ) and max μ Ξ train ( Δ N du (μ)/ ψ N du N (μ) X ) as a function of N pr and N du for the RB approximations computed during the greedy procedure, for the primal (left) and the dual (right) problem, respectively. Here Ξ train is a uniform random sample of size n train =1,000 and the RB tolerance is ε tol = 10 2 .

In Figures 19 and 20 some representative solutions at time t=T, for selected values of the parameter, show the physical convective phenomena at different Péclet numbers, to underline the different nature of the problem: from a diffusion dominated (lower Péclet number) to a transport dominated (higher Péclet number) problem. Two different cases have been analyzed, concerning the mass transfer through the cylindrical body: in the first case, we have imposed a mass flux g(t)=10t (substance release by the cylinder), while in the second case g(t)=10t (substance absorption through the cylinder). In any case, higher values of concentration and higher gradients are obtained for larger Peclet numbers: absorption or release are more effective when transport dominates over diffusion.

Figure 19
figure 19

Diffusion-transport problem around a cylinder: representative solutions for μ 1 =0.1 and μ 1 =100 at time t=T, g(t)=10t.

Figure 20
figure 20

Diffusion-transport problem around a cylinder: representative solutions for μ 1 =0.1 and μ 1 =100 at time t=T, g(t)=10t.

In the following Figure 21 the behavior of the (RB evaluation of) output (67) is shown, as well as the related error bounds (magnified by a factor 10), in the case of heat emission (Figure 19); we have considered a concentrated (in time) output (corresponding to h(t)=δ(t)), given by the (spatial) average of the concentration on the cylinder at each timestep. According to the behavior of solutions, we obtain higher values of the output when μ 1 increases.

Figure 21
figure 21

Diffusion-transport problem around a cylinder: RB concentrated outputs and related error bounds as functions of time t, for μ 1 =0.1,1,100.

8.4 Computational aspects

We conclude this section by discussing some computational aspects related to the three numerical examples presented above, and showing how reduced basis techniques allow a substantial reduction of computational work. We recall that, in order to obtain a rapid and reliable procedure, we are interested in (i) the minimization of the (marginal) cost associated with each input-output evaluation as well as in (ii) the possibility to provide a certification of each reduced approximation, both with respect to the corresponding finite element approximation.

All the details are reported in Table 1. Compared to the corresponding FE approximation, RB Online evaluations of field variables and outputs enable a computational speedup, defined as S= t FE / t RB online , of about two orders of magnitude. In particular, the average time over 2,500 Online output evaluations is of 0.107 for the first Couette-Graetz problem, of 0.198 for the second heat treatment problem, as well as of 0.158 for the third diffusion-transport problem. Note that the times related to the RB Online evaluation take into account also the a posteriori error estimation for solution and output. This great computational advantage is due, basically, to the reduction in linear system dimensions, and finally in the huge dimensional reduction - N vsN - between RB spaces and corresponding FE spaces. For the three cases considered, this ratio goes from 260 (first case) to 820 (third case). Thanks to the Gram-Schmidt orthonormalization, the condition number of the RB matrices is limited to O( 10 2 ), while without this procedure it will go to O( 10 14 ).

In the end, we take into account also the time spent for the Offline construction and storage; this allows to determine the break-even point, given by Q BE = t RB offline / t FE . In particular, we obtain a break-even point of O( 10 2 ) in the three cases, which can be considered acceptable whenever interested either in the real-time context, or in the limit of many queries. The performances described in Table 1 are valid even if we consider a higher number of parameters (for example, with P between 10 and 25, see [57]).

9 Perspectives and ongoing research

We end this review paper dedicated to applications of reduced basis method in an industrial framework by putting current methodology development in perspective.

9.1 Extension to complex problems

Growing research areas are devoted to the following kind of problems.

(i) Nonlinear problems: the reduced basis framework and related model–reduction approaches are well developed for linear parametrized partial differential equations. They can be effectively applied also to nonlinear problems [37, 58, 59], even if this in turn introduces both numerical and theoretical complications, and many open research issues are still to be faced. Classical problems arising in applied sciences are, for example, Navier-Stokes/Boussinesq and Burgers’ equations in fluid mechanics [1618, 47, 48, 60, 61] and nonlinear elasticity in solid mechanics.

First of all, computational complexity is increasing at both the Offline and the Online stage: we need to solve nonlinear problems of big dimension O(N) during the RB space generation, as well as nonlinear problems of reduced dimension O(N) for each Online evaluation; in both the cases, classical iterative procedure - such as fixed point or Newton-type algorithms - can be used. A posteriori error bounds introduced for linear problems can be effectively extended to steady nonlinear problems (see for example, [62] for steady incompressible Navier-Stokes equations). However, the most important challenge deals with the reliability and/or the certification of the methodology in the unsteady - parabolic - problems [23, 63]: in these cases exponential instability seriously compromises a priori and a posteriori error estimates, yielding to bounds which are limited to modest (final) times and modest Reynolds numbers. More precisely, stability considerations limit the product of the final time and the Reynolds number [64].

(ii) Problems dealing with (homogeneous or even heterogeneous) couplings in a multiphysics setting and based on domain decomposition techniques: a domain decomposition approach [29, 65] combined with reduced basis method has been successfully applied in [27, 28, 66] and further extensions are foreseen [67]. A coupled multiphysics setting has been proposed for simple fluid-structure interaction problems [68, 69].

(iii) Optimal control [7073], shape optimization, inverse and design problems [74, 75] as many-query applications have been and are subject to extensive research, which is of interest also in an industrial context. One of the main goals of this field is the study of efficient techniques to deal with geometrical parameters, in order to keep the number of parameters reasonable but also to guarantee versatility in the parametrization in order to treat and represent complex shapes. Recent works [7679] deal with free-form deformation techniques combined with empirical interpolation in bio-medical and aerodynamic problems.

(iv) Another growing field is related with the development and application of the reduced basis methodology to the quantification of uncertainty [24, 80, 81].

9.2 Efficiency improvement in RB methodology

The efforts are also aimed at improving the computational performance in three-dimensional settings to have a more efficient implementation of the Offine ‘construction stage’ (for example, on high-performance parallel supercomputers) and more and more attractive real-time applications such as the ones currently available on smartphones [82].

Improvements in the efficiency of parameters space exploration are also crucial; see, for example, modified greedy algorithms and combined adaptive techniques [83], such as ‘hp’ RB method [84, 85]. At the same time, (i) improvements in the a posteriori error bounds for nonaffine problems [38]; (ii) reduction of the complexity of the parametrized operators and more efficient estimation of lower bounds of stability factors (that is, coercivity or inf-sup constants) for complex nonaffine problems [86]; or (iii) more specialized RB spaces [87] are under investigation.

Footnotes

1This assumption will greatly simplify the presentation while still exercising most of the important RB concepts; furthermore, many important engineering problems are in fact ‘compliant’.

2These regions can represent different material properties, but they can also be used for algorithmic purposes to ensure well-behaved mappings.

3Here, for 1l L dom , K o , l :D R 3 × 3 is a given SPD matrix (which in turn ensures coercivity of the bilinear form): the upper 2×2 principal submatrix of K o , l is the usual tensor conductivity/diffusivity; the (3,3) element of K o , l represents the identity operator (‘mass matrix’) and is equal to M o , l ; and the (3,1), (3,2) (and (1,3), (2,3)) elements of K o , l - which we can choose here as zero thanks to the current restriction to symmetric operators - permit first derivative terms to take into consideration transport/convective terms.

4Under the coercivity and the symmetry assumptions, the bilinear form a(,;μ) defines a (energy) scalar product given by ( ( w , v ) ) μ :=a(w,v;μ)w,vX; the induced energy norm is given by .

5Clearly the accuracy and cost of the a posteriori error estimator Δ N (μ) are crucial to the success of the greedy algorithm.

6As we assumed that the bilinear form is coercive and the FE approximation spaces are conforming, it follows that α N (μ)α(μ) α 0 >0, μD.

7Similar results can be obtained for the a posteriori error bounds in the X norm.

8It thus follows that the a posteriori error estimation contribution to the cost of the greedy algorithm of Section 5 is O(Q N max N )+O( Q 2 N max 2 N)+O( n train Q 2 N max 3 ): we may thus choose N and n train independently (and large).

9Typical output fuctionals correspond to the ‘integral’ of the field u(μ) over an area or line (in particular, boundary segment) in Ω ¯ . However, by appropriate lifting techniques, ‘integrals’ of the flux over boundary segments can also be considered.

10We pursue here just a primal approximation, however we can readily extend the approach to a primal-dual formulation as described for coercive problems in Section 7.1.

11All over the section, Ω o (μ) denotes the original (physical) domain, whose generic point is indicated as x=( x 1 , x 2 ); for the sake of simplicity, we formulate all the problems in the original domain, but remove all the subscripts o . Moreover, a tilde ~ denotes dimensional quantities, while the absence of a tilde signals a non-dimensional quantity.

References

  1. Rozza G, Huynh P, Patera A: Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch. Comput. Methods Eng. 2008, 15: 229–275. 10.1007/s11831-008-9019-9

    Article  MATH  MathSciNet  Google Scholar 

  2. Patera, A., Rozza, G.: Reduced Basis Approximation and a posteriori Error Estimation for Parametrized Partial Differential Equations. Version 1.0, MIT. http://augustine.mit.edu (2006)

    Google Scholar 

  3. Prud’homme C, Rovas D, Veroy K, Maday Y, Patera A, Turinici G: Reliable real-time solution of parametrized partial differential equations: reduced-basis output bounds methods. J. Fluids Eng. 2002, 124: 70–80. 10.1115/1.1448332

    Article  Google Scholar 

  4. Porsching TA: Estimation of the error in the reduced basis method solution of nonlinear equations. Math. Comput. 1985,45(172):487–496. 10.1090/S0025-5718-1985-0804937-0

    Article  MATH  MathSciNet  Google Scholar 

  5. Ito K, Ravindran S: A reduced-order method for simulation and control of fluid flow. J. Comput. Phys. 1998,143(2):403–425. 10.1006/jcph.1998.5943

    Article  MATH  MathSciNet  Google Scholar 

  6. Almroth BO, Stern P, Brogan FA: Automatic choice of global shape functions in structural analysis. AIAA J. 1978, 16: 525–528. 10.2514/3.7539

    Article  Google Scholar 

  7. Noor A: Recent advances in reduction methods for nonlinear problems. Comput. Struct. 1981, 13: 31–44. 10.1016/0045-7949(81)90106-1

    Article  MATH  Google Scholar 

  8. Noor A: On making large nonlinear problems small. Comput. Methods Appl. Mech. Eng. 1982, 34: 955–985. 10.1016/0045-7825(82)90096-2

    Article  MATH  Google Scholar 

  9. Fink JP, Rheinboldt WC: On the error behavior of the reduced basis technique for nonlinear finite element approximations. Z. Angew. Math. Mech. 1983, 63: 21–28. 10.1002/zamm.19830630105

    Article  MATH  MathSciNet  Google Scholar 

  10. Porsching TA, Lee MYL: The reduced-basis method for initial value problems. SIAM J. Numer. Anal. 1987, 24: 1277–1287. 10.1137/0724083

    Article  MATH  MathSciNet  Google Scholar 

  11. Barrett A, Reddien G: On the reduced basis method. Z. Angew. Math. Mech. 1995,75(7):543–549. 10.1002/zamm.19950750709

    Article  MATH  MathSciNet  Google Scholar 

  12. Rheinboldt WC: On the theory and error estimation of the reduced basis method for multi-parameter problems. Nonlinear Anal. 1993,21(11):849–858. 10.1016/0362-546X(93)90050-3

    Article  MATH  MathSciNet  Google Scholar 

  13. Gunzburger, M.D.: Finite Element Methods for Viscous Incompressible Flows. Academic Press, (1989)

    MATH  Google Scholar 

  14. Ito K, Ravindran S: A reduced basis method for control problems governed by PDEs. In Control and Estimation of Distributed Parameter System Edited by: Desch W., Kappel F., Kunisch K.. 1998, 153–168.

    Chapter  Google Scholar 

  15. Ito K, Ravindran S: Reduced basis method for optimal control of unsteady viscous flows. International Journal of Computational Fluid Dynamics 2001,15(2):97–113. 10.1080/10618560108970021

    Article  MATH  MathSciNet  Google Scholar 

  16. Peterson J: The reduced basis method for incompressible viscous flow calculations. SIAM J. Sci. Stat. Comput. 1989,10(4):777–786. 10.1137/0910047

    Article  MATH  Google Scholar 

  17. Nguyen, N.C., Veroy, K., Patera, A.T.: Certified real-time solution of parametrized partial differential equations. In: Yip, S. (ed.) Handbook of Materials Modeling, pp. 1523–1558. Springer (2005)

    Google Scholar 

  18. Veroy, K., Prud’homme, C., Rovas, D.V., Patera, A.: A posteriori error bounds for reduced basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations. In: Proceedings of the 16th AIAA Computational Fluid Dynamics Conference (2003). Paper 2003–3847

    Book  Google Scholar 

  19. Maday Y, Patera A, Turinici G: A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations. J. Sci. Comput. 2002,17(1–4):437–446.

    Article  MATH  MathSciNet  Google Scholar 

  20. Buffa, A., Maday, Y., Patera, A., Prud’homme, C., Turinici, G.: A priori convergence of the greedy algorithm for the parametrized reduced basis. M2AN Math. Model. Numer. Anal. (2009), submitted

    Google Scholar 

  21. Holmes P, Lumley J, Berkooz G: Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge University Press, UK; 1996.

    Book  MATH  Google Scholar 

  22. Haasdonk B, Ohlberger M: Reduced basis method for finite volume approximations of parametrized linear evolution equations. M2AN Math. Model. Numer. Anal. 2008,42(2):277–302. 10.1051/m2an:2008001

    Article  MATH  MathSciNet  Google Scholar 

  23. Nguyen N, Rozza G, Patera A: Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers’ equation. Calcolo 2009,46(3):157–185. 10.1007/s10092-009-0005-x

    Article  MATH  MathSciNet  Google Scholar 

  24. Nguyen N, Rozza G, Huynh P, Patera A: Reduced basis approximation and a posteriori error estimation for parametrized parabolic PDEs; application to real-time Bayesian parameter estimation. In Large-Scale Inverse Problems and Quantification of Uncertainty. Edited by: Biegler L., Biros G., Ghattas O., Heinkenschloss M., Keyes D., Mallick B., Marzouk Y., Tenorio L., van Bloemen Waanders B., Willcox K.. John Wiley & Sons, Ltd, UK; 2010:151–178. Chap. 8 Chap. 8

    Chapter  Google Scholar 

  25. Balmes E: Parametric families of reduced finite element models: theory and applications. Mech. Syst. Signal Process. 1996,10(4):381–394. 10.1006/mssp.1996.0027

    Article  Google Scholar 

  26. Barrault M, Maday Y, Nguyen N, Patera A: An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. C.R. Math. Acad. Sci. Paris, Series I 2004,339(9):667–672. 10.1016/j.crma.2004.08.006

    Article  MATH  MathSciNet  Google Scholar 

  27. Løvgren AE, Maday Y, Rønquist EM: A reduced basis element method for the steady Stokes problem. M2AN Math. Model. Numer. Anal. 2006,40(3):529–552. 10.1051/m2an:2006021

    Article  MathSciNet  Google Scholar 

  28. Løvgren AE, Maday Y, Rønquist EM: The Reduced Basis Element Method for Fluid Flows. In Analysis and Simulation of Fluid Dynamics. Advances in Mathematical Fluid Dynamics. Birkhäuser, Boston; 2007:129–154.

    Chapter  Google Scholar 

  29. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations. Springer-Verlag, (1994)

    MATH  Google Scholar 

  30. Quarteroni, A.: Numerical Models for Differential Problems, Series MS&A, vol. 2 Springer (2009)

    Book  Google Scholar 

  31. Grepl M, Patera AT: A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations. M2AN Math. Model. Numer. Anal. 2005, 39: 157–181. 10.1051/m2an:2005006

    Article  MATH  MathSciNet  Google Scholar 

  32. Binev, P., Cohen, A., Dahmen, W., Devore, R., Petrova, G., Wojtaszczyk, P.: Convergence rates for greedy algorithms in reduced basis methods. (2010), in preparation

    Google Scholar 

  33. Huynh P, Rozza G, Sen S, Patera A: A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability costants. C. R. Acad. Sci. Paris, Series I 2005, 345: 473–478.

    Article  MathSciNet  Google Scholar 

  34. Huynh P, Knezevic D, Chen Y, Hesthaven J, Patera A: A natural-norm successive constraint method for inf-sup lower bounds. Comput. Methods Appl. Mech. Eng. 2010,199(29–32):1963–1975. 10.1016/j.cma.2010.02.011

    Article  MATH  MathSciNet  Google Scholar 

  35. Pierce N, Giles M: Adjoint recovery of superconvergent functionals from PDE approximations. SIAM Rev. 2000,42(2):247–264. 10.1137/S0036144598349423

    Article  MathSciNet  Google Scholar 

  36. Rozza G: Reduced basis approximation and error bounds for potential flows in parametrized geometries. Commun. Comput. Phys. 2011, 9: 1–48.

    MATH  MathSciNet  Google Scholar 

  37. Grepl M, Maday Y, Nguyen N, Patera A: Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations. ESAIM Math. Modelling Numer. Anal. 2007,41(3):575–605. 10.1051/m2an:2007031

    Article  MATH  MathSciNet  Google Scholar 

  38. Eftang J, Grepl M, Patera A: A posteriori error bounds for the empirical interpolation method. C.R. Math. Acad. Sci. Paris, Series I 2010,348(9–10):575–579. 10.1016/j.crma.2010.03.004

    Article  MATH  MathSciNet  Google Scholar 

  39. Nguyen N: A posteriori error estimation and basis adaptivity for reduced-basis approximation of nonaffine-parametrized linear elliptic partial differential equations. J. Comput. Phys. 2007, 227: 983–1006. 10.1016/j.jcp.2007.08.031

    Article  MATH  MathSciNet  Google Scholar 

  40. Rozza, G., Huynh, P., Manzoni, A.: Reduced basis approximation and error bounds for Stokes flows in parametrized geometries: roles of the inf-sup stability constants. (2010), submitted

    Google Scholar 

  41. Rozza G, Veroy K: On the stability of the reduced basis method for Stokes equations in parametrized domains. Comput. Methods Appl. Mech. Eng. 2007,196(7):1244–1260. 10.1016/j.cma.2006.09.005

    Article  MATH  MathSciNet  Google Scholar 

  42. Rozza G: Reduced basis methods for Stokes equations in domains with non-affine parameter dependence. Comput. Vis. Sci. 2009,12(1):23–35. 10.1007/s00791-006-0044-7

    Article  MathSciNet  Google Scholar 

  43. Gerner, A., Veroy, K.: Reduced basis a posteriori error bounds for the Stokes equations in parametrized domains: a penalty approach. Math. Mod. and Meth. in Appl. Sc., (2011), in press

    Google Scholar 

  44. Milani R, Quarteroni A, Rozza G: Reduced basis method for linear elasticity problems with many parameters. Comput. Methods Appl. Mech. Eng. 2008, 197: 4812–4829. 10.1016/j.cma.2008.07.002

    Article  MATH  MathSciNet  Google Scholar 

  45. Huynh P, Patera A: Reduced basis approximation and a posteriori error estimation for stress intensity factors. Int. J. Numer. Methods Eng. 2007,72(10):1219–1259. 10.1002/nme.2090

    Article  MATH  Google Scholar 

  46. Huynh, P., Rozza, G.: Reduced basis method and a posteriori error estimation: application to linear elasticity problems. (2011), submitted

    Google Scholar 

  47. Deparis S, Rozza G: Reduced basis method for multi-parameter-dependent steady Navier-Stokes equations: applications to natural convection in a cavity. J. Comput. Phys. 2009,228(12):4359–4378. 10.1016/j.jcp.2009.03.008

    Article  MATH  MathSciNet  Google Scholar 

  48. Rozza G, Nguyen C, Patera A, Deparis S: Reduced basis methods and a posteriori error estimators for heat transfer problems. Proceedings of HT2009, 2009 ASME Summer Heat Transfer Conference, S. Francisco, USA 2009. Paper HT 2009–88211 Paper HT 2009–88211

    Google Scholar 

  49. Rozza G, Nguyen C, Huynh P, Patera A: Real-time reliable simulation of heat transfer phenomena. Proceedings of HT2009, 2009 ASME Summer Heat Transfer Conference, S. Francisco, USA 2009. Paper HT 2009–88212 Paper HT 2009–88212

    Google Scholar 

  50. Gelsomino F, Rozza G: Comparison and combination of reduced order modelling techniques in 3D parametrized heat transfer problems. Math. Comput. Model. Dyn. Syst. 2011, 17: 371–394. 10.1080/13873954.2011.547672

    Article  MATH  MathSciNet  Google Scholar 

  51. Sen S, Veroy K, Huynh P, Deparis S, Nguyen N, Patera A: ‘Natural norm’ a posteriori error estimators for reduced basis approximations. J. Comput. Phys. 2006, 217: 37–62. 10.1016/j.jcp.2006.02.012

    Article  MATH  MathSciNet  Google Scholar 

  52. Chen Y, Hesthaven J, Maday Y, Rodríguez J: Certified reduced basis methods and output bounds for the harmonic Maxwell’s equations. SIAM J. Sci. Comput. 2010,32(2):970–996. 10.1137/09075250X

    Article  MATH  MathSciNet  Google Scholar 

  53. rbMIT Software: MIT http://augustine.mit.edu/methodology/methodology_rbMIT_System.htm (2007–2011)

  54. Reduced Basis at MIT: MIT http://augustine.mit.edu/methodology.htm (2007–2011)

  55. Panton, R.L.: Incompressible Flow, 3rd edn. John Wiley & Sons, Inc. (2005)

    Google Scholar 

  56. Incropera, F.P., DeWitt, D.P.: Fundamentals of Heat and Mass Transfer. John Wiley & Sons (1990)

    Google Scholar 

  57. Sen S: Reduced basis approximation and a posteriori error estimation for many-parameter heat conduction problems. Numer. Heat Transf., Part B, Fundam. 2008,54(5):369–389. 10.1080/10407790802424204

    Article  Google Scholar 

  58. Canuto C, Tonn T, Urban K: A posteriori error analysis of the reduced basis method for non-affine parameterized nonlinear PDEs. SIAM J. Numer. Anal. 2009,47(3):2001–2022. 10.1137/080724812

    Article  MATH  MathSciNet  Google Scholar 

  59. Jung N, Haasdonk B, Kröner D: Reduced basis method for quadratically nonlinear transport equations. International Journal of Computing Science and Mathematics 2009,2(4):334–353. 10.1504/IJCSM.2009.030912

    Article  MATH  MathSciNet  Google Scholar 

  60. Veroy K, Prud’homme C, Patera AT: Reduced-basis approximation of the viscous Burgers equation: rigorous a posteriori error bounds. C. R. Acad. Sci. Paris, Série I 2003,337(9):619–624.

    Article  MATH  MathSciNet  Google Scholar 

  61. Quarteroni A, Rozza G: Numerical solution of parametrized Navier-Stokes equations by reduced basis methods. Numer. Methods Partial Differ. Equ. 2007,23(4):923–948. 10.1002/num.20249

    Article  MATH  MathSciNet  Google Scholar 

  62. Veroy K, Patera A: Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds. Int. J. Numer. Methods Fluids 2005, 47: 773–788. 10.1002/fld.867

    Article  MATH  MathSciNet  Google Scholar 

  63. Knezevic D, Nguyen N, Patera A: Reduced basis approximation and a posteriori error estimation for the parametrized unsteady Boussinesq equations. Math. Models Methods Appl. Sci. 2011, 21: 1415–1442. 10.1142/S0218202511005441

    Article  MATH  MathSciNet  Google Scholar 

  64. Johnson C, Rannacher R, Boman M: Numerics and hydrodynamic stability: toward error control in computational fluid dynamics. SIAM J. Numer. Anal. 1995,32(4):1058–1079. 10.1137/0732048

    Article  MATH  MathSciNet  Google Scholar 

  65. Quarteroni, A., Valli, A.: Domain Decomposition Methods for Partial Differential Equations. Oxford University Press, (1999)

    MATH  Google Scholar 

  66. Løvgren, A.E., Maday, Y., Rønquist, EM: The reduced basis element method: offline-online decomposition in the nonconforming, nonaffine case. In: Hesthaven, J.S., Rønquist, E.M. (ed.) Spectral and High Order Methods for Partial Differential Equations. Selected papers from the ICOSAHOM ’09 Conference, June 22–26, Trondheim, Norway. Lecture Notes in Computational Science and Engineering, vol. 76, pp. 247–254. Springer (2011)

    Chapter  Google Scholar 

  67. Iapichino, L., Quarteroni, A., Rozza, G.: A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks. (2011), submitted

    Google Scholar 

  68. Lassila T, Rozza G: Model reduction of steady fluid-structure interaction problems with free-form deformations and reduced basis method. Proceedings of 10th Finnish Mechanics Days, Jyvaskyla, Finland 2009, 454–465.

    Google Scholar 

  69. Lassila, T., Quarteroni, A., Rozza, G.: A reduced model with parametric coupling for fluid-structure interaction problems. (2011), submitted

    Google Scholar 

  70. Quarteroni A, Rozza G, Quaini A: Reduced basis methods for optimal control of advection-diffusion problem. In Advances in Numerical Mathematics Edited by: Fitzgibbon W., Hoppe R., Periaux J., Pironneau O., Vassilevski Y.. 2007, 193–216.

    Google Scholar 

  71. Tonn T, Urban K, Volkwein S: Optimal control of parameter-dependent convection-diffusion problems around rigid bodies. SIAM J. Sci. Comput. 2010,32(3):1237–1260. 10.1137/08074194X

    Article  MATH  MathSciNet  Google Scholar 

  72. Dedè L: Reduced basis method and a posteriori error estimation for parametrized linear-quadratic optimal control problems. SIAM J. Sci. Comput. 2010,32(2):997–1019. 10.1137/090760453

    Article  MATH  MathSciNet  Google Scholar 

  73. Tonn T, Urban K, Volkwein S: Comparison of the reduced basis and POD a posteriori error estimators for an elliptic linear-quadratic optimal control problem. Math. Comput. Model. Dyn. Syst. 2011, 17: 355–369. 10.1080/13873954.2011.547678

    Article  MATH  MathSciNet  Google Scholar 

  74. Rozza, G.: Shape design by optimal flow control and reduced basis techniques: applications to bypass configurations in haemodynamics. Ph.D. thesis, N. 3400, École Polytechnique Fédérale de Lausanne (2005)

    Google Scholar 

  75. Rozza G: On optimization, control and shape design of an arterial bypass. Int. J. Numer. Methods Fluids 2005,47(10–11):1411–1419. 10.1002/fld.888

    Article  MATH  MathSciNet  Google Scholar 

  76. Lassila T, Rozza G: Parametric free-form shape design with PDE models and reduced basis method. Comput. Methods Appl. Mech. Eng. 2010, 199: 1583–1592. 10.1016/j.cma.2010.01.007

    Article  MATH  MathSciNet  Google Scholar 

  77. Rozza, G., Lassila, T., Manzoni, A.: Reduced basis approximation for shape optimization in thermal flows with a parametrized polynomial geometric map. In: Hesthaven, J.S., Rønquist, E.M. (ed.) Spectral and High Order Methods for Partial Differential Equations. Selected papers from the ICOSAHOM ’09 Conference, June 22–26, Trondheim, Norway. Lecture Notes in Computational Science and Engineering, vol. 76, pp. 307–315. Springer (2011)

    Chapter  Google Scholar 

  78. Manzoni, A., Quarteroni, A., Rozza, G.: Shape optimization for viscous flows by reduced basis methods and free-form deformation. (2011), submitted

    Google Scholar 

  79. Rozza G, Manzoni A: Model order reduction by geometrical parametrization for shape optimization in computational fluid dynamics. In Proceedings of ECCOMAS CFD 2010, V European Conference on Computational Fluid Dynamics, Lisbon, Portugal Edited by: Pereira J.C.F., Sequeira A.. 2010.

    Google Scholar 

  80. Boyaval S, Le Bris C, Maday Y, Nguyen N, Patera A: A reduced basis approach for variational problems with stochastic parameters: Application to heat conduction with variable Robin coefficient. Comput. Methods Appl. Mech. Eng. 2009,198(41–44):3187–3206. 10.1016/j.cma.2009.05.019

    Article  MATH  MathSciNet  Google Scholar 

  81. Huynh, P., Knezevic, D., Patera, A.: Certified reduced basis model characterization: a frequentistic uncertainty framework. Comput. Methods Appl. Mech. Engrg (2011), submitted

    Google Scholar 

  82. Huynh P, Knezevic D, Peterson J, Patera A: High-fidelity real-time simulation on deployed platforms. Comput. Fluids 2011, 43: 74–81. 10.1016/j.compfluid.2010.07.007

    Article  MATH  MathSciNet  Google Scholar 

  83. Haasdonk, B., Dihlmann, M., Ohlberger, M.: A training set and multiple bases generation approach for parametrized model reduction based on adaptive grids in parameter space. (2010), submitted

    Google Scholar 

  84. Eftang J, Patera A, Rønquist E: An “hp” certified reduced basis method for parametrized elliptic partial differential equations. SIAM J. Sci. Comput. 2010,32(6):3170–3200. 10.1137/090780122

    Article  MATH  MathSciNet  Google Scholar 

  85. Eftang J, Knezevic D, Patera A: An “hp” certified reduced basis method for parametrized parabolic partial differential equations. Math. Comput. Model. Dyn. Syst. 2011, 17: 395–422. 10.1080/13873954.2011.547670

    Article  MATH  MathSciNet  Google Scholar 

  86. Lassila T, Rozza G: Model reduction of semiaffinely parametrized partial differential equations by two-level affine approximation. C.R. Math. Acad. Sci. Paris, Series I 2011,349(1–2):61–66. 10.1016/j.crma.2010.11.016

    Article  MATH  MathSciNet  Google Scholar 

  87. Eftang, J., Huynh, P., Knezevic, D., Patera, A.: A two-step certified reduced basis method. J. Sci. Comput. (2011), submitted

    Google Scholar 

Download references

Acknowledgements

A special thank to Prof. Anthony Patera (MIT) for deep insights and several discussions and ideas. We acknowledge all the people who have contributed to the rbMIT package (beta version) used for RB computations presented in this work. In particular, we thank Dr. N.C. Nguyen (MIT) and Dr. D.B.P. Huynh (MIT) for their insights and contributions, as well as Dr. T. Lassila (EPFL) for his feedbacks and suggestions.

This work has been supported in part by the Swiss National Science Foundation (Project 200021-122136) and by the Progetto Roberto Rocca (MIT-Politecnico di Milano).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Gianluigi Rozza.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Authors’ original file for figure 22

Authors’ original file for figure 23

Authors’ original file for figure 24

Authors’ original file for figure 25

Authors’ original file for figure 26

Authors’ original file for figure 27

Authors’ original file for figure 28

Authors’ original file for figure 29

Authors’ original file for figure 30

Authors’ original file for figure 31

Authors’ original file for figure 32

Authors’ original file for figure 33

Authors’ original file for figure 34

Authors’ original file for figure 35

Authors’ original file for figure 36

Authors’ original file for figure 37

Authors’ original file for figure 38

Authors’ original file for figure 39

Authors’ original file for figure 40

Authors’ original file for figure 41

Authors’ original file for figure 42

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Quarteroni, A., Rozza, G. & Manzoni, A. Certified reduced basis approximation for parametrized partial differential equations and applications. J.Math.Industry 1, 3 (2011). https://doi.org/10.1186/2190-5983-1-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-5983-1-3

Keywords