Reduced basis approximation and a

Immanuel Maiera · Gianluigi Rozzab · Bernard
Haasdonka
Reduced basis approximation and a-posteriori error
estimation for the coupled Stokes-Darcy system
Stuttgart, February 2014
a
Institute for Applied Analysis and Numerical Simulation, University of Stuttgart,
Pfaffenwaldring 57, 70569 Stuttgart, Germany
Tel.: +49-711-685-65509, Fax: +49-711-685-65507
E-mail: {immanuel.maier,haasdonk}@mathematik.uni-stuttgart.de
b
SISSA MathLab,
International School for Advanced Studies,
Office A-435, Via Bonomea 265, 34136 Trieste, Italy
Tel.: +39 040 3787 451, Fax: +39 040 3787 528
E-mail: [email protected]
Abstract The coupling of a free flow with a flow through porous media has many potential applications in several fields related with computational science and engineering, such as blood flows, environmental problems or food technologies. We present a reduced basis method for such coupled problems.
The reduced basis method is a model order reduction method applied in the context of parametrized
systems. Our approach is based on a heterogeneous domain decomposition formulation, namely the
Stokes-Darcy problem. Thanks to an offline/online-decomposition, computational times can be drastically reduced. At the same time the induced error can be bounded by fast evaluable a-posteriori error
bounds. In the offline-stage the proposed algorithms make use of the decomposed problem structure.
Rigorous a-posteriori error bounds are developed, indicating the accuracy of certain lifting operators
used in the offline-stage as well as the accuracy of the reduced coupled system. Numerical experiments
dealing with groundwater flow scenarios demonstrate the efficiency of the approach.
Keywords Reduced basis method · Stokes flow · Porous medium equation · Domain decomposition ·
Non-coercive problem · Error estimation
Stuttgart Research Centre for Simulation Technology (SRC SimTech)
SimTech – Cluster of Excellence
Pfaffenwaldring 5a
70569 Stuttgart
[email protected]
www.simtech.uni-stuttgart.de
Immanuel Maiera et al.
2
1 Introduction
The fundamental idea of reduced basis (RB) methods [1, 23–25, 29] is approximating the solution manifold of a parametric partial differential equation (PDE) by a low-dimensional linear space, commonly
spanned by solution instantiations (snapshots) for carefully selected parameters. This RB space is built
up on a high-dimensional approximation space, from which the solutions are taken, for example a finite
element (FE) space. The efficiency of the method highly relies on the affine decomposition property of
the weak form. It ensures an efficient decomposition of the computations in an offline- and online-stage.
In the offline-stage the RB space is built up once for all RB simulations. This is a costly process, but
independent of the ultimate simulation requests. In the online-stage RB solutions can be computed
for many parameters, but the computational load is independent of the high dimension. The quality
of the approximation can be certified by efficiently computable a-posteriori error estimates.
The advantages of the method are worth while for several settings, where many parameter requests
are posed (many-query context) or where fast solution evaluations are needed (real-time context). For
example RB methods were applied to optimization problems [5, 22, 27]. The range of RB methods
includes time dependent problems [10, 11] as well as noncoercive and nonlinear problems [25, 33]. More
specifically, much investigation was done in solving saddle-point problems evolving from fluid simulations [9, 14, 15, 17, 27, 28]. The idea to include domain decomposition methods [26, 31] in the computational framework of the RB method was introduced in [18, 19], following a Mortar-like approach for
a coercive and linear PDE and extended in [14, 15, 17] to the Stokes equations. Other approaches are
based on the Schur-complement in domain decomposition [8, 12, 32] and on the Dirichlet-Neumann
solution scheme [20]. A method suited for the simulation of fluidic networks was developed in [14].
In this paper a RB method for the heterogeneous coupling of the Stokes equations and the porous
medium equation [6, 7, 16] is presented. The known procedures for homogeneous couplings of noncoercive PDEs [14] are extended to a more general setting. A subdomain-lifting framework is introduced
and integrated into the Greedy-algorithm [29] for the basis generation, yielding independent processes
on the subdomains. The algorithm takes account of the fact that the complexity of the problem with
respect to the parameter is expected to be smaller on the interface between the subdomains than on the
whole domain. Regarding a-posteriori error estimation of the global error, the standard non-coercive
error estimate [25] can be very pessimistic due to a great difference of the stability constants of the
subproblems. We present a modified error bound, which is more suited in that case and refines the
error bounds of [9].
The outline of the paper is as follows. We introduce the coupled system of PDEs and the weak
formulation in Section 2. In Section 3 the considered parameter dependency is explained and the
resulting parametrized problem is detailed. The analysis of the problem is summarized in Section
4. Then we introduce the reduced basis formulation together with the basis generation procedure in
Section 5 and give the definitions and proofs of the error estimates in Section 6. An investigation
by numerical tests dealing with groundwater flows is portrayed in Section 7. We conclude with some
remarks on the results and future work.
2 Problem formulation
We consider the coupled Stokes-Darcy system [6, 7, 16] in a two-dimensional domain Ω ⊂ R2 , decomposed into two subdomains Ωi , i = 1, 2 with Ω1 ∩ Ω2 = ∅ and Ω 1 ∪ Ω 2 = Ω. We denote by
Γ = ∂Ω1 ∩ ∂Ω2 the common boundary or interface, depicted in Figure 1, with associated normal
unit vector n pointing from Ω1 to Ω2 and tangent unit vector τ . All boundaries are assumed to be
Lipschitz-boundaries.
The differential multidomain formulation is as follows. On Ω1 we consider a Stokes flow: Find
velocity u and pressure p such that
−ν∆u + ∇p = f
∇·u= 0
u = uin
u=0

in Ω1 , 


in Ω1 , 
on Γ1in ,



on Γ1 ,
(1)
RB approximation for the coupled Stokes-Darcy system
3
Fig. 1 Illustration of the computational domains and boundaries.
and on Ω2 the porous medium equation: Find the pressure head ϕ such that

−∇ · (K∇ϕ) = 0
in Ω2 , 

−K∇ϕ · n2 = 0
on Γ2 ,


ϕ = ϕb on Γ2b .
(2)
Here ν is the kinematic viscosity, K is the hydraulic conductivity tensor, f a source term, uin the
inflow data, ϕb the Darcy pressure on the bottom and ni the outward normal unit vector on ∂Ωi ,
i = 1, 2. The Darcy velocity is given by
1
uD = − K∇ϕ ,
n
where n is the volumetric porosity constant. The systems are coupled through the following interface
conditions on Γ :

1

u · n = − K∇ϕ · n on Γ ,


n



∂u
−νn ·
+ p = gϕ
on Γ ,
(3)
∂n


√



k
∂u

τ·
on Γ .
u·τ =−
αBJ
∂n
The first two conditions represent the balance of mass and forces. The last one was proposed based
on experiments and is a boundary condition rather than a transmission condition. Here g denotes
the gravity acceleration, αBJ the Beavers-Joseph constant depending on the permeable material and
k = ν1 τ · K · τ . We introduce the weak formulation with the definition of the functional spaces:

X1,u = {v ∈ (H 1 (Ω1 ))2 : (v · n1 )|Γ1 = 0 ∧ (v · n1 )|Γ1in = 0} ,




in

X1,u
= {v ∈ (H 1 (Ω1 ))2 : (v · n1 )|Γ1 = 0 ∧ v|Γ = 0} ,


2
(4)
X1,p = L (Ω1 ) ,


1


X2 = {ψ ∈ H (Ω2 ) : ψ|Γ2b = 0} ,




b
1
X2 = {ψ ∈ H (Ω2 ) : ψ|Γ = 0} .
The inflow data is lifted into the interior of Ω1 by means of a continuous extension operator Ein :
1/2
1/2
in
(H00 (Γ1in ))2 → X1,u
with (Ein v)|Γ1in = v for v ∈ (H00 (Γ1in ))2 . We also introduce a continuous
Immanuel Maiera et al.
4
extension operator Eb : H 1/2 (Γ2b ) → X2b , such that (Eb ψ)|Γ2b = ψ for ψ ∈ H 1/2 (Γ2b ). The weak solution
of (1)–(3) is given by (u + Ein uin , p, ϕ + Eb ϕb ) with (u, p, ϕ) ∈ X1,u × X1,p × X2 such that

A1 (u, v) + B1 (v, p) + C1 (ϕ, v) = F1,u (v), ∀v ∈ X1,u ,

B1 (u, q) = F1,p (q), ∀q ∈ X1,p ,
(5)


A2 (ϕ, ψ) + C2 (u, ψ) = F2 (ψ),
∀ψ ∈ X2 ,
where
A1 (v, w) =
Z
Ω1
ν∇v : ∇w dx +
Z
αBJ
ν √ (v · τ )(w · τ ) ds,
k
Γ
Z
gψ v · n ds,
C1 (ψ, v) =
Z
q ∇ · v dx,
B1 (v, q) = −
Ω1
Z
f · v dx − A1 (Ein uin , v),
F1,u (v) =
Ω1
Γ
F1,p (q) = −B1 (Ein uin , q)
and
A2 (ψ, χ) =
Z
Ω2
C2 (v, ψ) = −
Z
∇ψ · K∇χ dx,
Γ
nψ v · n ds,
F2 (ψ) = −A2 (Eb ϕb , ψ) .
3 The parametrized problem
We introduce a parameter vector µ ∈ P, P ⊂ RP , on which the material properties and also geometric
quantities depend. Here we follow the methodology for RB approximation in parametrized geometries
developed earlier [27, 28]. We suppose that Ω = Ω(µ) is obtained from a reference domain Ωref
through a piecewise affine geometry transformation. More precisely, let reference subdomains Ω ref =
SRi
r
Ω 1,ref ∪ Ω 2,ref , Ω1,ref ∩ Ω2,ref = ∅ with decompositions Ω i,ref = r=1
Ω i,ref for i = 1, 2 be given. Such
decompositions originate from coarse triangulations, as is depicted in Figure 2. On each subdomain
r
r
Ωi,ref
, r = 1, . . . , Ri , i = 1, 2 acts a bijective affine transformation Tir (·; µ) : Ωi,ref
→ Ω(µ) given
through
Tir (x; µ) = Gri (µ)x + tri (µ),
with a transformation matrix Gri : P → R2×2 and a translation vector tri : P → R2 . The transformation
mappings have to fullfill the following continuity condition:
r
′
Tir (x; µ) = Tir′ (x; µ),
r′
∀x ∈ Ω i,ref ∩ Ω i′ ,ref ,
for all r = 1, . . . , Ri , r′ = 1, . . . , Ri′ and i, i′ = 1, 2. Then our assumption reads
Ω i = Ωi (µ) =
Ri
[
Ωir (µ),
i = 1, 2
r=1
r
with Ωir (µ) = Tir (Ωi,ref
; µ) for all r = 1, . . . , Ri , i = 1, 2. Subsequently we define Γref = ∂Ω1,ref ∩∂Ω2,ref
r
r
and Γref = Γref ∩ Ω1,ref and have
Γ = Γ (µ) =
R1
[
r ; µ).
T1r (Γref
r=1
Regarding the material constants, we assume ν = ν(µ), n = n(µ), k = k(µ) and αBJ = αBJ (µ) to be
constant in space. Let µ = (µgeo , µmat ), where µgeo bundles the parameters shaping the geometrical
domain and µmat bundles the parameters representing material properties. Then ν(µ) = ν(µmat ),
RB approximation for the coupled Stokes-Darcy system
5
Fig. 2 Computational domain with parametrized properties: Here Ω1 is above
Γ and Ω2 below Γ . The widths of areas A and B are given by 0.5 + µA , respectively 0.5 + µB , where µ = (µA , µB ) ∈ [0, 0.5]2 . This is achieved with R1 = 3
affine transformations in the upper domain and R2 = 7 affine transformations
in the lower domain.
n(µ) = n(µmat ), k(µ) = k(µmat ) and αBJ (µ) = αBJ (µmat ). Conversely, we allow heterogeneity in
K = K(µ) away from the interface Γ (µ) (thus k(µ) remains constant) and so
−1
x ∈ Ωir (µ),
K(x; µ) = Kref ((Tir (·; µgeo )) (x); µmat ,
for all r = 1, . . . , Ri , i = 1, 2.
Problem (5) is solved by tracing it back to Ωref . As a consequence, in addition to the dependency on
the material parameters, we get dependency of (5) on the geometry parameters due to transformation
formulas.
Remark 1 For the sake of readability we do not introduce new notations for the function spaces on
Ωref , but understand the definitions (4) to be done on Ωref .
The resulting problem on Ωref then reads: For µ ∈ P find (u(µ), p(µ), ϕ(µ)) ∈ X1,u × X1,p × X2 such
that
A1 (u(µ), v; µ) + B1 (v, p(µ); µ) + C1 (ϕ(µ), v; µ) = F1,u (v; µ),
B1 (u(µ), q; µ) = F1,p (q; µ),
A2 (ϕ(µ), ψ; µ) + C2 (u(µ), ψ; µ) = F2 (ψ; µ),
∀v ∈ X1,u ,
∀q ∈ X1,p ,
∀ψ ∈ X2 ,
with
A1 (v, w; µ) =
B1 (v, q; µ) =
C1 (ψ, v; µ) =
F1,u (v; µ) =
R1 Z
X
r=1
r
(µ)∇wj dx +
(∇vj )T JA
1
r
Ω1,ref
R1 Z
X
r=1
r
Ω1,ref
R1 Z
X
r=1
r
Γref
R1 Z
X
r=1
R1 Z
X
r=1
r
Γref
r
(µ)(v · τ )(w · τ ) ds,
KA
1
q(JBr 1 (µ)∇) · v dx,
KCr1 (µ)ψv · n ds,
r
Ω1,ref
r
(µ)f · v dx − A1 (Ein uin , v; µ),
KF
1
F1,p (q; µ) = −B1 (Ein uin , q; µ),
(6a)
(6b)
(6c)
Immanuel Maiera et al.
6
where we omitted summation over j (j = 1, 2) and
r
(µ) = ν(µ)| det(Gr1 (µ))|Gr1 (µ)−1 Gr1 (µ)−T ,
JA
1
αBJ (µ) r
r
|G1 (µ)τ |,
JBr 1 (µ) = | det(Gr1 (µ))|Gr1 (µ)−T ,
(µ) = ν(µ) p
KA
1
k(µ)
r
(µ) = | det(Gr1 (µ))|
KF
1
KCr1 (µ) = g|Gr1 (µ)τ |,
and on the second subdomain
A2 (ψ, χ; µ) =
R2 Z
X
r
Ω2,ref
r=1
C2 (v, ψ; µ) = −
Z
r
Γref
r
(µ)∇χ dx,
(∇ψ)T JA
2
KCr2 (µ)ψ v · n ds,
F2 (ψ; µ) = −A2 (Eb ϕb , ψ; µ) ,
where
r
(µ) = | det(Gr2 (µ))|Gr2 (µ)−1 Kref (µ)Gr2 (µ)−T ,
JA
2
KCr2 (µ) = n(µ)|Gr2 (µ)τ |.
Here, τ is the tangential unit vector on Γref . We suppose that this procedure yields an affinely
parametrized system [23, 25, 29]. This means that decompositions of the following type can be applied:
QA1
A1 (v, w; µ) =
X
q
ΘA
(µ)Aq1 (v, w),
1
q=1
q
and parameter-in-dependent bilinear
with rapidly evaluable parameter-dependent coefficients ΘA
1
q
forms A1 . The size QA1 of the decomposition is preferably small. Analogous representations must
hold for A2 , B1 and Ci , i = 1, 2 and for the linear forms F1,u , F1,p and F2 .
If all parameter-dependent constants allow an affine decomposition, writing down that very decompositions for the linear and bilinear forms is straightforward, except for the terms involving |Gri (µ)τ |,
i = 1, 2, r = 1, . . . , Ri . For our purposes it is sufficient to assume that the tangential vector τ is
r
constant on Γref
, for all r = 1, . . . , R1 . For more general geometrical configurations one has to consider
the Empirical Interpolation Method (EIM) [2] to get approximative affine decompositions for those
functions.
4 Analysis of the problem
The wellposedness of system (5) has been studied in [6]. The proof makes use of the Brezzi stability
theory for saddle-point problems [3, 4]. We introduce here the properties of the involved linear and bilinear forms of the parametrized system, which ensure wellposedness. The function spaces are equipped
with the following norms: k · kX1,u = k · k(H01 (Ω1,ref ))2 , k · kX1,p = k · kL2 (Ω1,ref ) and k · kX2 = k · kH01 (Ω2,ref ) .
For all µ ∈ P the symmetric bilinear form A1 (·, ·; µ) is continuous and coercive on X1,u , that is
γA1 (µ) = sup
sup
v∈X1,u w∈X1,u
αA1 (µ) =
inf
v∈X1,u
A1 (v, w; µ)
< ∞,
kvkX1,u kwkX1,u
A1 (v, v; µ)
> 0,
kvk2X1,µ
∀µ ∈ P,
∀µ ∈ P.
(7)
Furthermore, supµ∈P γA1 (µ) < ∞ and inf µ∈P αA1 (µ) > 0 is supposed. The same properties are
assumed to be valid for A2 (·, ·; µ) on X2 × X2 with corresponding constants γA2 (µ) and αA2 (µ). The
bilinear form B1 (·, ·; µ) is continuous and inf-sup stable on X1,u × X1,p , that is
β(µ) = inf
sup
q∈X1,p v∈X1,u
B1 (v, q; µ)
> 0,
kvkX1,u kqkX1,p
∀µ ∈ P
(8)
RB approximation for the coupled Stokes-Darcy system
7
and we suppose inf µ∈P β(µ) > 0. The bilinear forms Ci (·, ·; µ), i = 1, 2 are continuous on X1,u × X2
respectively X2 × X1,u with corresponding constants γCi (µ) and supµ∈P γCi (µ) < ∞ for i = 1, 2.
Beyond those standard assumptions, we additionally require the following condition: For all µ ∈ P
there exist constants αCi (µ) > 0, i = 1, 2, such that
αC1 (µ)C1 (ψ, v; µ) + αC2 (µ)C2 (v, ψ; µ) = 0,
∀v ∈ X1,u , ψ ∈ X2 .
(9)
Last but not least the linear forms F1,u (·; µ), F1,p (·; µ) and F2 (·; µ) are continuous on X1,u , X1,p and
X2 , respectively. We report now the a-priori result [6] adapted to our notation.
Proposition 1 For every µ ∈ P there exists a unique solution (u(µ), p(µ), ϕ(µ)) ∈ X1,u × X1,p × X2
of the system (6). The following a-priori estimates are valid:
√ 1
1
γu (µ) + αu (µ) 2
CF ,
k(u(µ), ϕ(µ))kX1,u ×X2 ≤
2CF (µ) +
αu (µ)
β(µ)
√
1
γu (µ)
kp(µ)kX1,p ≤
2 1+
CF ,u (µ)
αC1 (µ)β(µ)
αu (µ)
γu (µ)(γu (µ) + αu (µ))
CF ,P (µ) ,
+
αu (µ)β(µ)
with
αu (µ) = min αCi (µ)αAi (µ),
i=1,2
γu (µ) = 2 max max αCi (µ)γAi (µ), max αCi (µ)γCi (µ) ,
i=1,2
i=1,2
n
o
′
CF ,u (µ) = max αC1 (µ)kF1,u (·; µ)kX1,u
, αC2 (µ)kF2 (·; µ)kX2′ ,
(10)
(11)
′
CF ,P (µ) = kF1,p (·; µ)kX1,p
.
5 Reduced basis approximation
We maintain all the previous notations, despite considering the spaces X1,u , X1,p and X2 to be Finite
Element (FE) approximation spaces now. We introduce N1,u = dim(X1,u ), N1,p = dim(X1,p ) and N2 =
dim(X2 ). Let us assume that we are given low-dimensional RB approximation spaces XN,1,u ⊂ X1,u ,
XN,1,p ⊂ X1,p and XN,2 ⊂ X2 with dimensions N1,u = dim(XN,1,u ) ≪ N1,u , N1,p = dim(XN,1,p ) ≪
N1,p and N2 = dim(XN,2 ) ≪ N2 . For µ ∈ P the solution of (6) is approximated by the solution
(uN (µ), pN (µ), ϕN (µ)) ∈ XN,1,u × XN,1,p × XN,2 of the following problem:

A1 (uN (µ), v; µ) + B1 (v, pN (µ); µ)



+C1 (ϕN (µ), v; µ) = F1,u (v; µ), ∀v ∈ XN,1,u, 
(12)
B1 (uN (µ), q; µ) = F1,p (q; µ), ∀q ∈ XN,1,p , 



A2 (ϕ(µ)N , ψ; µ) + C2 (uN (µ), ψ; µ) = F2 (ψ; µ),
∀ψ ∈ XN,2 .
The wellposedness of problem (6) is almost ensured by the assumptions of Section 4. The missing
assumption is the inf-sup stability of B1 (·, ·; µ), because it is not naturally inherited by the subspaces.
We want to ensure that
βN (µ) =
inf
sup
q∈XN,1,p v∈XN,1,µ
B1 (v, q; µ)
> 0,
kvkXN,1,u kqkXN,1,p
∀µ ∈ P
plus inf µ∈P βN (µ) > 0. This stability can be achieved by the so-called supremizer enrichment [9, 27,
28, 30]. For µ ∈ P the supremizer operator is given through
T (µ) : X1,p → X1,u : q 7→ T (µ)q :
(T (µ)q, v) = B1 (v, q; µ),
∀v ∈ X1,u .
(13)
Immanuel Maiera et al.
8
There are multiple options concerning Offline/Online-decomposition and orthonormalization of the
resulting bases, for which stability has been proven theoretically or at least has been observed in practice
during computation. The procedure we choose will be clearly explained after the basis generation
method is presented in the ongoing section. In order to formulate our basis generation method, some
further definitions are needed. We define the trace operators
1/2
γ1 : X1,u → H00 (Γref ) : v 7→ v|Γref · n,
1/2
γ2 : X2 → H00 (Γref ) : ψ 7→ ψ|Γref
and the function spaces
X1,Γ = {γ1 v : v ∈ X1,u },
X2,Γ = {γ2 ϕ : ϕ ∈ X2 }.
From the definition of the coupling functions Ci (·, ·; µ) it is clear that they are determined by above
defined traces of their arguments. Thus Ci (·, ·; µ) can be evaluated at every ν ∈ X1,Γ and η ∈ X2,Γ and
we use the abbrevation C1 (η, ·; µ) for C1 (ψ, ·; µ), where ψ ∈ X2 is an arbitrary function with γ2 ψ = η
and C2 (ν, ·; µ) for C2 (v, ·; µ), where v ∈ X1,u , γ1 v = ν, respectively. If the solution of (6) is known at
the interface, it can be reconstructed by the following problem specific lifting operators:
E1 (µ) : X2,Γ → X1,u × X1,p : η 7→ (E1,u (µ)η, E1,p (µ)η) :
A1 (E1,u (µ)η, v; µ) + B1 (v, E1,p (µ)η; µ) = F1,u (v; µ)
− C1 (η, v; µ),
B1 (E1,u η(µ), q; µ) = F1,p (q; µ),



∀v ∈ X1,u ,

∀q ∈ X1,p . 
(14)
∀ψ ∈ X2 .
(15)
E2 (µ) : X1,Γ → X2 : ν 7→ E2 (µ)ν :
A2 (E2 (µ)ν, ψ; µ) = F2 (ψ; µ) − C2 (ν, ψ; µ),
So it follows (u(µ), p(µ), ϕ(µ)) = (E1,u (µ)γ2 ϕ(µ), E1,p (µ)γ2 ϕ(µ), E2 (µ)γ1 u(µ)). Now, if the RB spaces
yield good approximations to subproblems (14) and (15) for all traces of solutions {(u(µ), p(µ), ϕ(µ)) :
µ ∈ P}, the approximation quality of (12) is expected to be fairly good, too. To develop further this
analysis, let us assume that we are given sets of interface modes ΞN,i,Γ ⊂ Xi,Γ , i = 1, 2:
k
ΞN,i,Γ = ξN,i,Γ
, k = 1, . . . , N1,Γ ,
i = 1, 2.
Those interface modes can be gained from global snapshots: ΞN,1,Γ = {γ1 u(µ) : µ ∈ S1,Γ } and
ΞN,2,Γ = {γ2 ϕ(µ) : µ ∈ S2,Γ } with Si,Γ ⊂ P, |Si,Γ | = Ni,Γ for i = 1, 2. Other options one can think
of are solutions to a generalized eigenproblem on Γref (see [8, 12]) or Fourier interface functions (see
Section 7.3 and [14]).
Remark 2 The motivation of this construction is that the variation of the solutions to (6) on the interface Γ is potentially smaller than on the whole domain Ω. While using the standard RB approach the
snapshots have to represent the solution on the whole domain accurate enough, here we are only aiming
at the representation of the solutions values on the interface. Once the interface modes are given, the
RB spaces are generated via the lifting operators. In this step, which is explained in Algorithm 1 below, only subproblems have to be solved, in contrast to global coupled problems in a “straightforward”
RB procedure. This procedure yields a certain flexibility in choosing between different interface mode
options. If it is known a-priori, that a specific parameter has hardly any influence on the solutions
interface values, the amount of snapshots computed for the interface modes can be reduced drastically.
RB approximation for the coupled Stokes-Darcy system
9
k
′
′
Now we define RB approximations to Eik (µ) = Ei (µ)ξN,i
′ ,Γ , k = 1, . . . , Ni′ ,Γ , 1 ≤ i, i ≤ 2, i 6= i on
XN,1,u × XN,1,p and XN,2 :
k
k
k
EN,1
(µ) = (EN,1,u
(µ), EN,1,p
(µ)) ∈ XN,1,u × XN,1,p :
k
k
A1 (EN,1,u
(µ), v; µ) + B1 (v, EN,1,p
(µ); µ) = F1,u (v; µ)
k
− C1 (ξN,2,Γ
, v; µ),
∀v ∈ XN,1,u,
k
B1 (EN,1,u
(µ), q; µ) = F1,p (q; µ),
∀q ∈ XN,1,p .
k
EN,2
(µ) ∈ XN,2 :
k
k
A2 (EN,2
(µ), ψ; µ) = F2 (ψ; µ) − C2 (ξN,1,Γ
, ψ; µ),
∀ψ ∈ XN,2 .
We now can formulate the basis generation algorithm for the RB problem (12) using the notation
X1 = X1,u × X1,p and XN,1 = XN,1,u × XN,1,p . It is based on the Greedy-algorithm for RB approximations [29] and especially on the Greedy-procedures realized in [14] in the context of parametric
boundary conditions.
Algorithm 1 Let ΞN,i,Γ , i = 1, 2 be given. In addition, we assume that initial RB spaces XN,1,u,
XN,1,p and XN,2 are given. We further assume that error indicators Iik (µ) for the lifting errors kEik (µ)−
k
EN,i
(µ)kXi , k = 1, . . . , Ni′ ,Γ , i = 1, 2 and an error indicator I(µ) for the error k(u(µ), p(µ), ϕ(µ)) −
(uN (µ), pN (µ), ϕN (µ))kX1 ×X2 are given. These can be the errors themselves or appropriate error
estimates. The following Greedy-procedure enlarges the RB spaces step by step. We use a finite training
set of parameters Ptrain ⊂ P, tolerances ǫi > 0 for the lifting errors and a tolerance ǫ > 0 for the error
of the full system.
parfor i = 1, 2
while max(µ,k)∈Ptrain ×{1,...,Ni,Γ } Iik (µ) > ǫi
if maxµ∈Ptrain I(µ) < ǫ
return
end
(µ∗i , ki∗ ) = arg max(µ,k)∈Ptrain ×{1,...,Ni,Γ } Iik (µ)
k1∗
k1∗
(i = 1): XN,1,u ← XN,1,u ⊕ span{EN,1,u
(µ∗1 ), T (µ∗1 )EN,1,p
(µ∗1 )}
k∗
1
XN,1,p ← XN,1,p ⊕ span{EN,1,p
(µ∗1 )}
k2∗
(i = 2): XN,2 ← XN,2 ⊕ span{EN,2(µ∗2 )}
end
end
Here, parfor indicates a possible parallel processing of the loop and T denotes the supremizer
operator defined in (13). The initialization of the algorithm requires generating ΞN,i,Γ , i = 1, 2 by
computing global solutions of (6) and initializing the RB spaces by performing the enrichment steps
of the algorithm on the empty spaces with random parameters (µ∗i , ki∗ ), i = 1, 2.
6 A-posteriori error estimation
As error indicators Iik (µ) and I(µ) in Algorithm 1 one could consider the errors itself, but this would
make an exhaustive search unrealistic, since Ptrain is a big set covering P well. Thus a-posteriori error
bounds are used, which are fastly evaluable. We can go back to known estimates [25] for kE1k (µ) −
k
k
(µ)kX2 (coercive case). We do not go into detail on
EN,1
(µ)kX1 (non-coercive case) and kE2k (µ) − EN,2
that, but introduce the notation
k
k
πN,i
(µ) := Eik (µ) − EN,i
(µ),
(16)
Immanuel Maiera et al.
10
and assume that estimates
k
k
ΠN,i
(µ) ≥ kπN,i
(µ)kXi ,
(17)
for k = 1, . . . , Ni,Γ , i = 1, 2, µ ∈ P are given. In the following we discuss a-posteriori error estimation
for the full solution (uN (µ), pN (µ), ϕN (µ)).
For convenience we first introduce some notations that will occur often in this section. We denote
by z LB a computable lower bound and by z UB a computable upper bound of any respective constant
quantity z. For any vector space Z, we denote by Z ′ the dual space of Z.
Let us introduce the notation X = X1,u × X1,p × X2 . Now we define the bilinear form A : X × X ×
P → R through
A(V, W; µ) := αC1 (µ)(A1 (v, w; µ) + B1 (w, q; µ) + B1 (v, r; µ) + C1 (ψ, w; µ))
+ αC2 (µ)(A2 (ψ, χ; µ) + C2 (v, χ; µ))
and the linear form F : X × P → R through
F (V; µ) := αC1 (µ)(F1,u (v; µ) + F1,p (q; µ)) + αC2 (µ)F2 (ψ; µ),
where V = (v, q, ψ) and W = (w, r, χ). Then problem (6) can be written as
ˆ
A(U(µ),
V; µ) = F (V; µ),
∀V ∈ X,
ˆ
where U(µ)
= (u(µ), p(µ), ϕ(µ)) ∈ X for µ ∈ P. We introduce the Babuˇska inf-sup condition, which
is known to be fullfilled under our assumptions [34]:
βA (µ) := inf sup
V∈X W∈X
A(V, W; µ)
0
≥ βA
> 0,
kVkX kWkX
∀µ ∈ P
(18)
and define the residual R : X × P → R through
ˆ N (µ), V; µ),
R(V, µ) := F (V; µ) − A(U
ˆ N (µ) = (uN (µ), pN (µ), ϕN (µ)) ∈ X is the reduced basis solution of problem (12). Then the
where U
following result holds:
ˆ
ˆ N (µ) can be bounded in the X-norm by the
Proposition 2 [25, 28] The error eN (µ) = U(µ)
−U
a-posteriori error estimate
ˆ
ˆ N (µ)kX ≤ ∆N (µ) := kR(·; µ)kX ′ .
kU(µ)
−U
LB (µ)
βA
Proof From (18) we get
ˆ
ˆ
ˆ
ˆ N (µ)kX ≤ sup A(U(µ) − UN (µ), V; µ)
βA (µ)kU(µ)
−U
kVkX
V∈X
ˆ
ˆ N (µ), V; µ)
A(U(µ),
V; µ) − A(U
= sup
kVkX
V∈X
ˆ N (µ), V; µ)
F (V; µ) − A(U
= sup
kVkX
V∈X
R(V; µ)
= sup
= kR(·; µ)kX ′
V∈X kVkX
⊓
⊔
RB approximation for the coupled Stokes-Darcy system
11
The drawback of the estimator ∆N (µ) is that it does not take account of the possibility that the
LB
stability constants of the subproblems have different magnitudes. The resulting problem is, that βA
can be very small. Better would be the weighting of residual norms corresponding to the errors on the
subdomains with constants related to the subproblems. In the following an error estimate is developed
that partly solves this problem.
Let us indroduce the notation XU = X1,u × X2 and XP = X1,p , which corresponds to “gathering”
the elliptic parts of the subproblems. Now we define the bilinear form AU : XU × XU × P → R through
AU (V, W, µ) = αC1 (µ) (A1 (v, w; µ) + C1 (ψ, w; µ))
+ αC2 (µ) (A2 (ψ, χ; µ) + C2 (v, χ; µ)) ,
and the linear form FU : XU × P → R through
FU (V, µ) = αC1 (µ)F1,u (v; µ) + αC2 (µ)F2 (ψ; µ),
where V = (v, ψ) and W = (w, χ). Additionaly we define B : XU × XP × P → R through
B(V, Q; µ) = αC1 (µ)B1 (v, q; µ)
and FP : XP × P → R through
FP (Q; µ) = αC1 (µ)F1,p (q; µ),
where Q = q. Then problem (6) can be written as
AU (U(µ), V; µ) + B(V, P (µ); µ) = FU (V; µ),
B(U(µ), Q; µ) = FP (Q; µ),
∀V ∈ XU ,
∀Q ∈ XP ,
(19)
where U(µ) = (u(µ), ϕ(µ)) ∈ XU and P (µ) = p(µ) ∈ XP for µ ∈ P. This reformulation as a saddlepoint problem is the foundation of the a-priori results reported in Section 4. There, constants αU (µ),
(10) and γU (µ), (11) were introduced, which in fact are the coercivity and continuity constants of
AU (·, ·; µ).
We define the residual RU : XU × P → R through
RU (V; µ) = FU (V; µ) − AU (UN (µ), V; µ) − B(V, PN (µ); µ),
where UN (µ) = (uN (µ), ϕN (µ)) and PN (µ) = pN (µ) is the reduced basis solution of problem (12).
Additionaly we define residuals R1,u : X1,u × P → R and R2 : X2 × P → R and R1,p : X1,p × P → R
through
R1,u (v; µ) = F1,u (v; µ) − A1 (uN (µ), v; µ) − B1 (v, pN (µ); µ) − C1 (ϕN , v; µ),
R1,p (q; µ) = F1,p (q; µ) − B1 (uN (µ), q; µ),
R2 (ψ; µ) = F2 (ψ; µ) − A2 (ϕN (µ), ψ; µ) − C2 (uN (µ), ψ; µ).
It holds RU (V; µ) = αC1 (µ)R1,u (v; µ)+αC2 (µ)R2 (ψ; µ) for V = (v, ψ) ∈ XU . With all the definitions
we now introduce new error bounds, which extend the theory of error approximation for parametrized
saddle-point problems of [9]. This means that the structure of the bounds is maintained, but the single
components are refined due to our additional subproblem structure.
Lemma 1 The error eN,U (µ) = U(µ) − UN (µ) can be bounded in the XU -norm by the a-posteriori
error estimate kU(µ) − UN (µ)kXU ≤ ∆N,U (µ) with
!
s
√
′
2
2˜
c1 (µ) kR1,p (·; µ)kX1,p
max
∆
(µ)
+
1
+
,
(20)
∆N,U (µ) :=
N,U,i
i=1,2
β LB (µ)
αLB
αLB
U (µ)
U (µ)
with
∆N,U,1 (µ) :=
s
αC1 (µ)
′
kR1,u (·; µ)kX1,u
,
αLB
A1 (µ)
∆N,U,2 (µ) :=
s
αC2 (µ)
kR2 (·; µ)kX2′
αLB
A2 (µ)
Immanuel Maiera et al.
12
and
UB
(µ) .
(µ), αC2 (µ)γCUB
c˜1 (µ) = max αC1 (µ)γA
2
1
The error eN,P (µ) = P (µ) − PN (µ) can be bounded in the XP -norm by the a-posteriori error estimate
kP (µ) − PN (µ)kXP ≤ ∆N,P (µ) :=
′
kR1,u (·; µ)kX1,u
β LB (µ)
+
√
∆N,U (µ)
,
2˜
c2 (µ) LB
β (µ)
(21)
with
UB
(µ)
(µ), γCUB
c˜2 (µ) = max γA
1
1
The quantities αAi (µ), αC1 (µ) and β(µ) have been introduced in Section 4, equations (7), (9) and (8).
Proof We define the linear operator B(µ) : X1,u → X1,p through
(B(µ)v, q) = B1 (v, q; µ),
∀v ∈ X1,u , q ∈ X1,p .
The error eN,U (µ) ∈ XU can be written as eN,U (µ) = (eN,u (µ), eN,ϕ(µ)). Further, we exploit the
0
⊥
unique decomposition eN,u (µ) = e0N,u (µ) + e⊥
N,u (µ), where eN,u (µ) ∈ ker(B(µ)) and eN,u (µ) ∈
⊥
ker(B(µ))⊥ and define e0N,U (µ) = (e0N,u (µ), eN,ϕ (µ)) and e⊥
N,U (µ) = (eN,u (µ), 0). For the sake of
readability we omit the µ-dependence in the following. It holds
0
AU (e0N,U , e0N,U ) = AU (eN,U , e0N,U ) − AU (e⊥
N,U , eN,U )
= FU (e0N,U ) − AU (UN , e0N,U )
0
− B(e0N,U , P ) − AU (e⊥
N,U , eN,U ).
Here we used (19). Now we add 0 = B(e0N,U , P − PN ) to the left hand side and get
0
AU (e0N,U , e0N,U ) = RU (e0N,U ) − AU (e⊥
N,U , eN,U )
0
= αC1 R1,u (e0N,u ) + αC2 R2 (eN,ϕ ) − AU (e⊥
N,U , eN,U ).
(22)
The residual terms can be estimated by
αC1 R1,u (e0N,u ) + αC2 R2 (eN,ϕ)
′
ke0N,u kX1,u + αC2 kR2 (·)kX2′ keN,ϕkX2
≤ αC1 kR1,u (·)kX1,u
s
q
αC1
′
≤
kR
(·)k
αC1 A1 (e0N,u , e0N,u )
1,u
X1,u
αLB
A1
s
q
αC2
′
kR
(·)k
+
αC2 A2 (eN,ϕ , eN,ϕ)
2
X
2
αLB
A2
q
q
0
0
≤ max ∆N,U,i
αC1 A1 (eN,u , eN,u ) + αC2 A2 (eN,ϕ, eN,ϕ )
i=1,2
q
√
≤ 2 max ∆N,U,i AU (e0N,U , e0N,U ).
i=1,2
√ √
The last line follows from a + b ≤ 2 a2 + b2 for a, b > 0, the definitions of AU , e0N,U and (9). The
rear term in (22) can be estimated via
⊥
⊥
0
0
−AU (e⊥
N,U , eN,U ) = − αC1 A1 (eN,u , eN,u ) − αC2 C2 (eN,u , eN,ϕ )
0
UB
ke⊥
≤ αC1 γA
N,u kX1,u keN,u kX1,u
1
ke⊥
+ αC2 γCUB
N,u kX1,u keN,ϕ kX2
2
√
⊥
ckeN,U kXU ke0N,U kXU
≤ 2˜
s
q
2
c˜ke⊥
AU (e0N,U , e0N,U ).
≤
N,U kXU
LB
αU
RB approximation for the coupled Stokes-Darcy system
Combining those two results and dividing by
13
q
AU (e0N,U , e0N,U ) yields
q
√
γ UB
AU (e0N,U , e0N,U ) ≤ 2 max ∆N,U,i + qU ke⊥
N,U kXU .
i=1,2
αLB
U
With the coercivity of AU we have
ke0N,U kXU
≤
s
UB
2
γU
max
∆
+
ke⊥
N,U,i
N,U kXU .
i=1,2
αLB
αLB
U
U
(23)
Now we estimate the error e⊥
N,U . Equation (8) implies [3, 4]
β(µ) =
inf
sup
v∈ker(B(µ))⊥ q∈X1,p
B1 (v, q; µ)
.
kvkX1,u kqkX1,p
This yields
⊥
ke⊥
N,U kXU = keN,u kX1,u ≤
=
=
1
β LB
sup
q∈X1,p
B1 (e⊥
N,u , q)
kqkX1,p
1
B1 (eN,u , q)
sup
β LB q∈X1,p kqkX1,p
1
F1,p (q) − B1 (uN , q)
1
′
sup
,
= LB kR1,p (·)kX1,p
β LB q∈X1,p
kqkX1,p
β
using e0N,u ∈ ker B(µ) and (6b). Combining the last estimate with (23) via the triangle-inequality
keN,U kXU ≤ ke0N,U kXU + ke⊥
N,U kXU yields the desired result (20).
Secondly, we define eN,p = p − pN and get from (8):
keN,P kXP = keN,p kX1,p ≤
B1 (v, eN,p )
1
sup
β LB v∈X1,u kvkX1,u
(24)
Using (6a) we have
B1 (v, eN,p) = F1,u (v) − A1 (u, v) − C1 (ϕ, v) − B1 (v, pN )
= R1,u (v) − A1 (u − uN , v) − C1 (ϕ − ϕN , v)
UB
UB
′
kvkX1,u
kϕ
−
ϕ
k
≤ kR1,u (·)kX1,u
+
γ
ku
−
u
k
+ γA
N
X
N
X
1,u
1,u
C1
1
√
UB
′
}keN,U kXU kvkX1,u .
, γCUB
≤ kR1,u (·)kX1,u
+ 2 max{γA
1
1
Inserting this estimate into (24) together with (20) yields the desired estimate (21).
⊓
⊔
LB
Remark 3 The quantities β LB (µ) and αLB
Ai (µ), i = 1, 2, respectively βA (µ) can be computed with the
LB
Successive constraint method [13]. A recipe for the computation of βA
(µ) is worked out in [21], in
µ∗
particular a lower bound for the surrogate inf-sup constant βA (µ), defined through
∗
µ
βA
(µ) = inf sup
V∈X W∈X
A(V, W; µ)
,
ˆ
kT (µ∗ )VkX kWkX
is computed, where µ∗ is a fixed parameter and Tˆ (µ∗ ) : X → X is the global supremizer operator
defined through
(Tˆ (µ∗ )V, W)X = A(V, W, µ∗ ),
∀W ∈ X.
Since the error eN,P (µ) is traced back to eN,U (µ), the estimate ∆N,U (µ) is sufficient to indicate
the global approximation quality. But of course the global error eN (µ) itself can be bounded with the
above estimates via eN (µ)2 = eN,U (µ)2 + eN,P (µ)2 .
Immanuel Maiera et al.
14
Fig. 3 Illustration of geometry for Model 1. (FF = Free
Flow, PM = Porous Medium)
7 Numerical results
7.1 A model with parametrized material properties
The first model we consider has a fixed computational domain without geometry parameters. The
model configurations we choose are Ω = [0, 4] × [0, 2], Ω1 = [0, 4] × [1, 2], Ω2 = [0, 4] × [0, 1], Γ1in =
{0} × [1, 2] and Γ2b = [0, 4] × {0}. We consider a heterogeneous porous medium and therefore introduce
a decomposition Ω 2 = Ω 2,U ∪ Ω 2,L depicted in Figure 3. The parameter vector lives in two dimensions:
µ = (µ1 , µ2 ) ∈ P ⊂ R2 with P = [1, 10]2. We introduce the kinematic viscosity of the fluid:
ν(x; µ) = 2µ1 · 10−3 ,
x ∈ Ω, µ ∈ P,
the Beavers-Joseph constant:
αBJ (x; µ) =
5 · 10−2
,
√
µ1
x ∈ Γ, µ ∈ P
and the hydraulic conductivity tensor of the porous medium equation:
5
χΩ2,L (x) Id2 ,
K(x; µ) = 1 · 10−5 χΩ2,U (x) +
µ2
x ∈ Ω2 , µ ∈ P,
where χΩ2,U/L (·) are the characteristic functions and Id2 is the unity matrix in R2×2 . All other data
functions are parameter independent: g = 9.807, n = 0.4, ϕb = 100 and
x ∈ Γ1in ,
uin (x) = 1 · 10−2 1 − 0.26(2 − x2 )2 ,
T
x ∈ Ω1 .
f (x) = 0, −1 · 10−4 x1 ,
At the corner of Γ1in , Γ the inflow conditions and the transmission conditions have to fullfill a compability condition, which is ensured by the data for every µ ∈ P. Two representative solutions of the
model are shown in Figures 4 and 5. Note that the plots are separated due to different value ranges.
7.2 A model with parametrized geometry
We also consider a model with parametrized geometry to investigate our algorithm. The computational
domain was already depicted in Figure 2. We have Ω1 = [0, 2] × [3.5, 5]\B and Ω2 = [0, 2] × [0, 3.5]\B,
where B = B(µ) = [1.5 − µB , 2] × [3, 4]. We assign a lower permeability to area A = A(µ) =
[0, µA ] × [1, 2]:
x ∈ Ω2 , µ ∈ P.
K(x; µ) = 1 · 10−5 χΩ2 /A (x) + 5 · 10−2 χA (x) Id2 ,
RB approximation for the coupled Stokes-Darcy system
15
Fig. 4 Sample solutions of Model 1: µ = (1, 1), left: pressures, right: velocities.
Fig. 5 Sample solutions of Model 1: µ = (10, 10), left: pressures, right: velocities.
Fig. 6 Sample solutions of Model 2: µ = (0, 0), left: pressures, right: velocities.
The parameter vector µ = (µA , µB ) lives in P = [0, 0.5]2 . A representative solution of Model 2 is
shown in Figure 6.
This model is inspired by biological applications, where flows through bloodvessels are investigated.
We performed numerical tests with both models and obtained similar results. The results presented in
the following were obtained with Model 1.
Immanuel Maiera et al.
16
|S | = |P
i,G
| = 100
|S | = |P
train
i,G
|S | = 36
i,G
2
10
i,G
|S | = 25
|S | = 25
2
i,G
i,G
10
|S | = 16
|S | = 16
i,G
|S | = 9
i,G
i,G
max errors
max errors
i,G
|S | = 9
0
10
−2
10
| = 100
train
|S | = 36
0
10
−2
−4
10
10
−6
10
−4
10
5
10
15
20
25
30
Number of iterations in basis generation algorithm
2
4
6
8
10
12
14
Number of iterations in basis generation algorithm
Fig. 7 Progress of the training error E in the Greedy-algorithm for different number of interface modes. Left:
k
k
. Here SI,G stands for S1,Γ = S2,Γ .
Iik = kπN,i
kXi , right: Iik = ΠN,i
Table 1 Final basis size N = N1,u + N1,p + N2 and maximal training error E for different number of interface
modes and different error indicators.
|Si,Γ |
Iik
k
kπN,i
kXi :
=
100
36
25
16
9
N
E
59
57
54
47
40
3.11 · 10−7
6.58 · 10−6
1.76 · 10−4
1.00 · 10−2
3.22 · 10−1
Iik
k
ΠN,i
:
=
N
E
43
40
32
35
25
9.40 · 10−4
3.00 · 10−3
1.10 · 10−2
2.20 · 10−2
3.30 · 10−1
7.3 Basis generation process
To examine the basis generation with Algorithm 1 we choose a training set Ptrain of 100 uniformly
distributed parameters in P. To simplify the notation we introduce
E := max keN (µ)kX ,
Ei :=
∆ := max ∆N (µ),
∆i :=
µ∈Ptrain
µ∈Ptrain
max
k
kπN,i
(µ)kXi ,
max
k
ΠN,i
(µ),
(µ,k)∈Ptrain ×{1,...,Ni,Γ }
(µ,k)∈Ptrain ×{1,...,Ni,Γ }
k
k
where πN,i
and ΠN,i
are defined in (16) and (17), respectively. We expect to reach the best possible
accuracy with respect to the error E on Ptrain when taking all snapshots Si,Γ = Ptrain for i = 1, 2.
Figure 7 shows the progress of the training error E in Algorithm 1 for different sized sets Si,Γ using the
k
k
errors kπN,i
(µ)kXi as well as the estimates ΠN,i
(µ) as error indicators. The tolerance ǫ of Algorithm 1
is chosen, such that it is never reached and so the algorithm continues until the ǫi criteria are violated.
The results are compared in Table 1.
When the number of interface modes is reduced, the error of the full system stagnates, whereas
the algorithm continues because the lifting approximations still improve. One could terminate earlier
in this case to get lower dimensions of the RB spaces. Figure 8 illustrates this observation in the case
of 25 interface modes, where
Erel := max
µ∈Ptrain
keN (µ)kX
,
ˆ
kU(µ)k
X
Ei,rel :=
max
(µ,k)∈Ptrain ×{1,...,Ni,Γ }
k
kπN,i
(µ)kXi
kEik (µ)kXi
.
An inaccuracy of the solver used to compute the residual components at around 1 · 10−12 reflects itself
RB approximation for the coupled Stokes-Darcy system
17
8
10
E
Erel
2
10
6
E
10
1
0
∆
E
10
4
E
10
max errors
10
2,rel
−4
10
max errors
2
E
1
1
1,rel
−2
E
∆
E
E
∆
2
2
2
10
0
−6
10
−8
10
10
−2
10
−10
−4
10
10
5
10
15
20
25
30
Number of iterations in basis generation algorithm
5
10
15
20
Number of iterations in basis generation algorithm
Fig. 8 Progress in Greedy-algorithm using |Si,Γ | = 25. Left: Lifting errors compared to errors of the global
k
system using Iik = kπN,i
kXi . Right: Lifting errors and estimates compared to errors and error estimates of the
k
global system using Iik = ΠN,i
.
in an inaccuracy of the error estimates at around 1 · 10−4 . One can see that the error of the full system
is dominated by the worse of the two lifting errors. Here, the bad effectivity of the estimate ∆N (µ)
clearly sticks out.
We investigate the time performance of our method in generating a basis of size N = 30 with an
accuracy of 1 · 10−4 on the training parameter set. The time needed for the offline stage is 2699s. The
ˆ N = (uN (µ), pN (µ), ϕN (µ) of RB problem
online stage includes the computation of the solution U
(12) and the computation of the error estimate ∆N (µ) defined in Proposition 2. Out of 10 executions
for random parameters, the average time needed for the online stage is 0.15s, whereas the average time
ˆ N = (u(µ), p(µ), ϕ(µ)), fulfilling equation (6), is 7.14s.
needed for computing the detailed solution U
As a result, the average speed up factor is 47.43. More significance has the break even quantity, which
also takes account of the offline computational time. Using the RB method pays off after 386 solution
evaluations. This value is also important for a comparison with a “straightforward” RB procedure. It
will further decrease, if the offline computational time decreases.
7.3.1 Basis enrichment with Fourier interface functions
The accuracy of a basis generated with 9 or 16 snapshots can be enriched without including further
snapshot information but using Fourier functions as interface modes:
k
ξN,i,Γ
(x) = cos π(k − 1)IΓ−1 (x) ,
k = 1, . . . , Ni,Γ , i = 1, 2
with a bijective mapping IΓ : [−0.5, 0.5] → Γref . Here ΞN,1,Γ = ΞN,2,Γ . This option is inspired by
[14]. Using this option, we are able to further decrease the RB approximation error. We introduce
the notations Φ9 , Φ16 , Φ25 and Φ36 for the bases generated in section 7.3 with |Si,Γ | = 9, 16, 25, 36
k
and Iik = ΠN,i
. Now we perform Algorithm 1 with initial basis Φ9 , respectively Φ16 and the Fourier
interface modes defined above. Figure 9 shows that the accuracies gained from 25 and 36 snapshots
can also be reached with less snapshots plus an additional Fourier mode basis enrichment algorithm.
Although needing to compute fewer global snapshots is very advantageous, the possibly larger basis
size resulting from a basis enrichment could affect the simulation times in the online-stage. The decision
between the two possibilities depends on the requirements to the method in different practical cases.
7.4 Comparison of error estimators
In Section 6 we already noted that a very small stability factor could be a problem with the estimate ∆N (µ). In the right plot of Figure 8 we observed the resulting badness of effectivities. Now
Immanuel Maiera et al.
18
0
10
φ9 enrichment
φ16 enrichment
φ25
−1
φ36
max errors
10
−2
10
−3
10
−4
10
0
5
10
15
20
25
Number of iterations in basis enrichment algorithm
30
Fig. 9 Progress of the training error E
in the basis enrichment algorithm of the
bases Φ9 and Φ16 with Fourier interface
modes. For comparison the final training
error E of Φ25 and Φ36 is displayed.
6
10
4
10
2
||eN(µk)||X
10
∆ (µk)
N
∆impr(µk)
0
10
N
−2
10
−4
10
0
20
40
60
k: test parameter index
80
100
Fig.
10 Evaluation
of
the
error
keN (µk )kX
and
the
estimates ∆N (µk ) and ∆impr
(µk )
=
N
(∆N,U (µk )2 + ∆N,P (µk )2 )1/2 at 100
randomly chosen test parameters µk ,
k = 1, . . . , 100.
we investigate the improved bounds introduced in Lemma 1, Section 6. Firstly we compare ∆N (µ)
to ∆impr
(µ) := (∆N,U (µ)2 + ∆N,P (µ)2 )1/2 in Figure 10. Although the estimate ∆impr
(µ) is slightly
N
N
sharper it is still highly overestimating the error.
We expect to get a further improvement, if we use ∆N,U (µ) as error indicator, which was defined
in (20) as an estimate for keN,U (µ)kXU . This is confirmed by the results shown in Figure 11. We were
able to improve the sharpness of the estimate by a factor of at least 1 · 102 .
8 Conclusions and perspectives
We presented a framework for RB approximation of the parametrized Stokes-Darcy system, picking up
recent developments in the field of RB methods for homogeneous domain decomposition problems. Our
algorithm allows for a heterogeneous problem structure and different ways of choosing interface modes
on the inner boundary. The a-posteriori error estimation for inf-sup stable saddle point problems has
been extended to the setting and a problem using the Babuˇska framework was revealed.
The strength of the method is that, even in a high-dimensional parameter domain, we need only
few global snapshots. In most realistic cases the subproblems will depend on fewer parameters than
the global problem. So high dimensions in the parameter space can be better treated. Especially, large
RB approximation for the coupled Stokes-Darcy system
3
19
8
10
10
k
k
N
N
X
impr k
k
(µ ) / ||e (µ )||
N
N
X
∆ (µk) / ||e (µk)||
N,U
N,U
X
∆ (µ ) / ||e (µ )||
∆
2
10
7
10
U
1
10
||e
10
6
10
k
(µ )||
N,U
X
U
0
k
∆N,U(µ )
5
10
−1
10
4
10
−2
10
−3
10
0
3
20
40
60
k: test parameter index
80
100
10
0
20
40
60
k: test parameter index
80
100
Fig. 11 Left: Evaluation of the error keN,U (µk )kXU and the estimate ∆N,U (µk ) at 100 randomly chosen test
parameters µk , k = 1, . . . , 100. Right: Corresponding effectivities of ∆N (µk ), ∆impr
(µk ) and ∆N,U (µk ).
N
parts of the basis generation process can be performed in parallel, leading to important advantages
and improvements.
Furthermore, there is a wide range of possible applications, including biological systems, environmental problems and food technologies. An extension of the method to FSI problems is foreseen.
Regarding future improvements, the conjunction of snapshot based basis generation and basis enrichment with Fourier functions has to be investigated more extensively.
Acknowledgements The authors would like to thank the German Research Foundation (DFG) for financial
support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University
of Stuttgart and the Baden-W¨
urttemberg Stiftung gGmbH. G. Rozza acknowledges the NOFYSAS program:
New Opportunities for Young Scientists at SISSA, Trieste, Italy.
References
1. Almroth, B.O., Stern, P., Brogan, F.A.: Automatic choice of global shape functions in structural
analysis. AIAA Journal 16, 525–528 (1978)
2. Barrault, M., Maday, Y., Nguyen, N., Patera, A.: An “empirical interpolation” method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus de
l’Acad´emie des Sciences, Series I 339, 667–672 (2004)
3. Brezzi, F.: On the existence, uniqueness and approximation of saddle point problems arising from
Lagrangian multipliers. R.A.I.R.O. 8, 129–151 (1974)
4. Brezzi, F., Fortin, M.: Mixed and Hybrid Finite Element Methods. Springer Verlag Berlin Heidelberg New York (1991)
5. Dihlmann, M.A., Haasdonk, B.: Certified PDE-constrained parameter optimization using reduced
basis surrogate models for evolution problems. Tech. rep., University of Stuttgart (2013). Submitted to Computational Optimization and Applications
6. Discacciati, M.: Domain decomposition methods for the coupling of surface and groundwater flows.
Ph.D. thesis, EPFL, Lausanne (2004)
7. Discacciati, M., Miglio, E., Quarteroni, A.: Mathematical and numerical models for coupling surface
and groundwater flows. Applied Numerical Mathematics 43, 57–74 (2002)
8. Eftang, J.L., Patera, A.T.: Port reduction in parametrized component static condensation: approximation and a posteriori error estimation. International Journal for Numerical Methods in
Engineering. (2013). DOI 10.1002/nme.4543
9. Gerner, A.L., Veroy, K.: Certified reduced basis methods for parametrized saddle point problems.
SIAM Journal on Scientific Computing 34(5), A2812–A2836 (2012)
20
Immanuel Maiera et al.
10. Grepl, M.A., Patera, A.T.: A posteriori error bounds for reduced-basis approximations of
parametrized parabolic partial differential equations. Mathematical Modelling and Numerical
Analysis 39, 157–181 (2005)
11. Haasdonk, B., Ohlberger, M.: Reduced basis method for finite volume approximations of
parametrized linear evolution equations. Mathematical Modelling and Numerical Analysis 42,
277–302 (2008)
12. Huynh, D.B.P., Knezevic, D.J., Patera, A.T.: A static condensation reduced basis element method:
Approximation and a posteriori error estimation. Mathematical Modelling and Numerical Analysis
47(1), 213–251 (2013)
13. Huynh, D.B.P., Rozza, G., Sen, S., Patera, A.T.: A successive constraint linear optimization method
for lower bounds of parametric coercivity and inf-sup stability constants. Comptes Rendus de
l’Acad´emie des Sciences, Series I 345, 473–478 (2007)
14. Iapichino, L.: Reduced basis method for the solution of parametrized PDEs in repetitive and
complex networks with application to CFD. Ph.D. thesis, EPFL, Lausanne (2012)
15. Iapichino, L., Quarteroni, A., Rozza, G.: A reduced basis hybrid method for the coupling of
parametrized domains represented by fluidic networks. Computer Methods in Applied Mechanics
and Engineering 221–222, 62–82 (2012)
16. Layton, W.J., Schieweck, F., Yotov, I.: Coupling fluid flow with porous media flow. SIAM Journal
on Numerical Analysis 40(6), 2195–2218 (2003)
17. Løvgren, A.E., Maday, Y., Rønquist, E.M.: A reduced basis element method for the steady Stokes
problem. Mathematical Modelling and Numerical Analysis 40, 529–552 (2006)
18. Maday, Y., Rønquist, E.M.: A reduced-basis element method. Journal of Scientific Computing 17,
447–459 (2002)
19. Maday, Y., Rønquist, E.M.: The reduced basis element method: Application to a thermal fin
problem. Journal of Scientific Computing 26, 240–258 (2004)
20. Maier, I., Haasdonk, B.: A Dirichlet-Neumann reduced basis method for homogeneous domain
decomposition problems. Applied Numerical Mathematics (2013). Accepted
21. Manzoni, A.: Reduced models for optimal control, shape optimization and inverse problems in
haemodynamics. Ph.D. thesis, EPFL, Lausanne (2012)
22. Oliveira, I.B., Patera, A.T.: Reduced-basis techniques for rapid reliable optimization of systems
described by affinely parametrized coercive elliptic partial differential equations. Optimization and
Engineering 8, 43–65 (2007)
23. Patera, A.T., Rozza, G.: Reduced Basis Approximation and A Posteriori Error Estimation for
Parametrized Partial Differential Equations. MIT Pappalardo Graduate Monographs in Mechanical
Engineering (2007). http://augustine.mit.edu
24. Prud’homme, C., Rovas, D.V., Veroy, K., Machiels, L., Maday, Y., Patera, A.T., Turinici, G.:
Reliable real-time solution of parametrized partial differential equations: Reduced-basis output
bound methods. Journal of Fluids Engineering 124, 70–80 (2002)
25. Quarteroni, A., Rozza, G., Manzoni, A.: Certified reduced basis approximation for parametrized
partial differential equations and applications. Journal of Mathematics in Industry 1:3 (2011).
DOI 10.1186/2190-5983-1-3
26. Quarteroni, A., Valli, A.: Domain Decomposition Methods for Partial Differential Equations. Oxford University Press (1999)
27. Rozza, G.: Shape design by optimal flow control and reduced basis techniques: applications to
bypass configurations in haemodynamics. Ph.D. thesis, EPFL, Lausanne (2005)
28. Rozza, G., Huynh, D.B.P., Manzoni, A.: Reduced basis approximation and a posteriori error estimation for Stokes flow in parametrized geometries: roles of the inf-sup stability constants. Numerische Mathematik 125, 115–152 (2013)
29. Rozza, G., Huynh, D.B.P., Patera, A.T.: Reduced basis approximation and a posteriori error
estimation for affinely parametrized elliptic coercive partial differential equations. Archives of
Computational Methods in Engineering 15, 229–275 (2008)
30. Rozza, G., Veroy, K.: On the stability of the reduced basis method for Stokes equations in
parametrized domains. Computer methods in applied mechanics and engineering 196, 1244–1260
(2007)
31. Toselli, A., Widlund, O.: Domain Decomposition Methods - Algorithms and Theory. SpringerVerlag Berlin Heidelberg New York (2005)
RB approximation for the coupled Stokes-Darcy system
21
32. Vallagh´e, S., Patera, A.T.: The static condensation reduced basis element method for a mixed-mean
conjugate heat exchanger model. SIAM Journal on Scientific Computing (2013). Revised
33. Veroy, K., Prud’homme, C., Rovas, D.V., Patera, A.T.: A posteriori error bounds for reduced-basis
approximation of parametrized noncoercive and nonlinear elliptic partial differential equations. In:
Proceedings of 16th AIAA computational fluid dynamics conference (2003). Paper 2003-3847
34. Xu, J., Zikatanov, L.: Some observations on Babuˇska and Brezzi theories. Numerische Mathematik
94, 195–202 (2003)

Download Report