A reduced basis element method for the steady Stokes problem: Application

In this thesis we consider the reduced basis element method for approximating the solution of parameter dependent problems described by partial differential equations. In particular we focus on fluid flow in pipes, bifurcations and hierarchical systems, where the flow is described by the steady Stokes equations, or the steady Navier-Stokes equations. The thesis consists of four papers and this introduction.The reduced basis element method is different from traditional reduced basis methods in that it combines these methods with domain decomposition. A given geometry is decomposed into building blocks with the same topology as a few reference domains, e.g. a rectangle and a reference bifurcation. Relative to each reference domain we precompute and store basis functions found on a preselected set of deformations of the respective reference domains. A reduced basis solution is then found by mapping the basis functions from their respective reference domains to corresponding domains in the domain decomposition. A local approximation of the “true” solution on one domain is found using the basis functions belonging to that specific domain, and the global approximation is found by “gluing” the local approximations together with constraints across domain interfaces. Geometries where building blocks of the same topology are reused many times are attractive candidates for the reduced basis element method. When there is only one domain in the geometry, the reduced basis element method is seen as a traditional reduced basis method where the geometry of the domain is one of the independent parameters.In the first part of the introduction we present the reduced basis method and the reduced basis element method. In the second part of the introduction we motivate the work done in the thesis, and give a summary of the papers.


Introduction
The reduced basis element method is a new approach for approximating the solution of problems described by partial differential equations.The method takes its roots in domain decomposition methods and in reduced basis discretizations.For a given parameter dependent problem: Find u ∈ X such that the computational effort needed to find an approximate discrete solution often makes the problem unsuitable for repetitive solves or real time control.Similarily, if (1) represents a very complex system, the resolution requirements may be so severe that even a single approximative solution may be hard to obtain.The idea behind reduced basis methods is to precompute several solutions of (1), {u i } N i=1 , corresponding to a preselected set of parameter values, S N = {µ i } N i=1 .
If the resolution of each u i is represented by N , then N N .These precomputed solutions (or "snapshots") are then used as a basis for the solution space of (1) to find the reduced basis solution for a generic µ where the coefficients α i (µ) are determined through a Galerkin method.The error between the reduced basis solution u N and the high resolution solution u N depends on the quality of the reduced basis space and on the underlying regularity of u N with respect to µ.As long as N is small, the work needed to find u N is negligible.Examples of reduced basis methods following this computational approach include the Proper Orthogonal Decomposition [5], Centroidal Voronoi Tessellations [7], and Output Bound methods [11].
In the reduced basis element method we consider the geometry of the computational domain as the generic parameter.The domain is decomposed into smaller blocks, all of which can be considered to be deformations of a few reference shapes.Associated with each reference shape are precomputed solutions for different deformations of the shapes.The precomputed solutions are mapped from the reference shapes to the different blocks of the decomposed domain, and the solution on each block is found as a linear combination of the mapped precomputed solutions.The solutions on the different blocks are glued together using Lagrange multipliers.We will in this work focus on hierarchical flow systems, which can be decomposed into pipes and bifurcations.We limit ourselves to the steady Stokes equations in the modeling of the flow through such systems, and in the next section we describe how the geometry enters the equations as a parameter.We use spectral elements in the modeling, but the method applies to other discretization techniques as well.

The steady Stokes problem
We consider here the two-dimensional steady Stokes equations where u = (u 1 , u 2 ) is the velocity field, p is the pressure, f = ( f 1 , f 2 ) is a prescribed volumetric body force, and ν is the fluid viscosity; see [1].For all the problems studied in this paper, this model will suffice.
The domain Ω has an inflow boundary Γ in , an outflow boundary Γ out , and wall boundaries Γ w .On this domain we introduce the velocity space where v t is the tangential velocity component.In addition, we have the Neumann type boundary conditions given by specifying σ n = ν ∂u n ∂n − p to be σ in n = −1 along Γ in and σ out n = 0 along Γ out ; here, u n is the normal velocity component and ∂/∂n denotes the derivative in the outward normal direction.For all the problems solved in this study, the exact solution of (4) satisfies ∂u n ∂n = 0 along Γ in and Γ out , which implies that the Neumann conditions correspond to specifying the pressure along the inflow and outflow boundaries (in a weak sense).With the given boundary conditions, we define the pressure space to be In order to solve the steady Stokes equations we define the bilinear forms and consider the weak form: Find u ∈ X (Ω) and p ∈ M(Ω) such that where For all the problems considered in this work, the body force f will be zero.
To ensure a unique solution of the steady Stokes problem (9), the coercivity condition (11) and the inf-sup condition inf must be satisfied; see [2] and [6].These conditions are fulfilled for our particular Stokes problem.

Discretization
We consider here a computational domain Ω which represents a two-dimensional flow system as depicted in Figure 1.We assume that the domain can be decomposed as a union of E non-overlapping subdomains Ω e , e = 1, . . ., E, with each subdomain representing a deformed square.Each deformed square is again a regular one-to-one deformation, φ e , of the reference square Ω = (−1, 1) 2 , i.e., In our case, each subdomain will be considered to be a single spectral element; see [12].Let P n ( Ω) be the space of all functions which are polynomials of degree less than or equal to n in each spatial direction on Ω.
For v e = v | Ωe , the discrete space for the velocity is then taken to be , v e • φ e ∈ (P N ( Ω)) 2 , e =1, ..., E}, (14) while the discrete space for the pressure is The bases for X N (Ω) and M N (Ω) are conveniently expressed in terms of the reference variables ξ and η.As a basis for X N (Ω) we use a nodal basis through the tensor-product Gauss-Lobatto Legendre (GLL) points, while the basis for M N (Ω) is a nodal basis through the tensor-product Gauss-Legendre (GL) points; see [12] and [13].Specifically, we write where i (ξ) refers to a one-dimensional N -th order Lagrangian interpolant through the GLL points ξ m , m = 0, ..., N ; here, i (ξ m ) j (ξ n ) = δ im δ jn for a given point (ξ m , ξ n ) in the underlying tensor-product GLL grid.
In a similar fashion, we write where ˜ i (ξ) refers to a one-dimensional (N −2)th order Lagrangian interpolant through the (interior) GL points ζ m , m = 0, ..., N −2; here, for a given point in the tensor-product GL grid.
Based on (13), where we have defined the global mapping Ω = Φ( Ω), we express the bilinear forms ( 7) and ( 8) in terms of the reference variables ξ and η as The elemental contributions to these sums are where J e is the Jacobian of φ e , and J e its determinant.The operator ∇ = [ ∂ ∂ξ , ∂ ∂η ] T .This gives us the following discrete system: Find u N ∈ X N (Ω) and where a N , b N and l N refer to integration of the bilinear and linear forms using Gauss-type quadrature.We thus see that the geometry enters the equations as a parameter via the mapping Φ, or more specifically, through the elemental mappings φ e , e = 1, . . ., E.

The reduced basis
We now define the reduced basis solution spaces X N (Ω) ⊂ X N (Ω) and M N (Ω) ⊂ M N (Ω).Our objective is to find a unique reduced basis solution As before the coercivity of a(•, •; Φ) holds for all v ∈ X N (Ω), since it is a subset of X N (Ω).The inf-sup condition (12), however, depends strongly on X N (Ω) and M N (Ω).
In a hierarchical flow system we differ between building blocks with pipe structure, and building blocks with bifurcation structure.The precomputation of basis functions is done separately for the two types of building blocks.By grouping the spectral elements of E k spectral elements, we may write where k = 1, ..., K 1 indicates the pipe blocks, and Each building block can again be expressed as for k = 1, . . ., K 1 + K 2 , and with φ k e ( Ω) denoting the mapping of the reference square to each individual spectral element in the building block.Since each block B k may consist of several spectral elements, we have, in general, K 1 + K 2 < E. We will also denote the restriction of a field v to a block B k as v k .
In this work, a pipe building block consists of a single spectral element, i.e., E k = 1, k = 1, . . ., K 1 .In order to generate the basis functions to be used on a pipe, we solve the steady Stokes problem (22) with the same possible boundary conditions as on {B k } K 1 k=1 for a preselected set of deformations of the reference domain, i=1 .This is achieved by solving the steady Stokes problem on a deformed pipe comprising three spectral elements.The restriction of the solution to each of these three elements will thus take care of the three possible types of pipe segments (inflow, interior, outflow) in the hierarchical flow system.The resulting velocity fields, {u i } 3N 1 i=1 , are then mapped to the reference domain Ω by the Piola transformation In this way, the direction of the velocity relative to the geometry is preserved, and by construction b(u i , q; Φ i ) = b( ûi , q • Φ i ; I), where I is the identity mapping.Hence, each precomputed velocity field ûi is also incompressible on Ω.A similar result holds when the inverse Piola transformation is applied in order to map the velocity from the reference domain to The pressure is a scalar field, and is mapped from Only the N 1 fields with the proper boundary conditions are used as basis functions on each B k .
A building block which represents a bifurcation must necessarily comprise several spectral elements.In this work, we use six elements to build each bifurcation; see Figure 2. Hence, E k = 6 for k = K 1 +1, . . ., K 1 +K 2 in (25).The basis functions associated with bifurcations are constructed by solving the steady Stokes equations on a preselected set of deformed bifurca- , where B is a reference bifurcation and F i is a regular mapping.To take care of the different boundary conditions needed in the hierarchical flow system, basis functions are constructed with pipes added to either the inflow boundary, or the outflow boundaries of the bifurcations, or both.Only the restriction of the solutions to the bifurcation blocks are used as basis functions.The resulting velocity solutions, {u i } 3N 2  i=1 , are again mapped to the reference domain Ω by the Piola transformation.In contrast to the solutions found on pipes, all velocity (and pressure) solutions comprise six segments, u i = {u ie } 6 e=1 , and they are mapped to the reference domain one element at a time.On each we use the N 2 basis functions found on deformed bifurcations with the right boundary conditions, mapped to the generic domain by the inverse Piola transformation, Now, on each of the six elements of a bifurcation the Jacobian is smooth and continuous, but across internal interfaces in a bifurcation it is not.To ensure that ũk i is a continuous function we must define both the preselected bifurcations {B i } N 2 i=1 , and the actual bifurcations in the hierarchical system {B k } K 2 k=1 , as C 1 deformations of the same reference bifurcation B. A method for achieving this is presented in [10].The pressure solutions on the deformed bifurcations also comprise six elemental contributions, and on the generic bifurcation B k , these are evaluated as pk We define the spaces for which we know that the inf-sup condition is not fulfilled since the velocity fields in Y 0 N (Ω) are all divergence free.The index N denotes the dimension of Y 0 N (Ω) and M N (Ω), and may be expressed as Recall that N 1 is the number of precomputed basis functions for the K 1 pipe blocks, and N 2 is the number of precomputed basis functions for the K 2 bifurcation blocks.We have that N 1 N and N 2 N , and both are independent of the number of spatial dimensions.We also need to enforce a continuity condition across the block interfaces in Ω, Γ kl = B k T B l .We try to minimize the jump across these block interfaces by introducing the constraints and where n is the unit normal vector of Γ kl , t is the unit tangential vector, and W n k,l and W t k,l are spaces of low order polynomials defined on Γ kl .In [9] it is shown that the order of these polynomial spaces can be used to control the jump across the interfaces for multielement pipes.We thus define the reduced basis velocity space and remark that X 0 N (Ω) X N (Ω) due to the jump across the block interfaces.If we are not interested in the pressure on the generic domain, we may solve the problem: Find In this case the inf-sup condition is insignificant.If we also want to find the pressure however, the reduced basis velocity space X N (Ω) must be enriched.This is due to the fact that the space X 0 N is spanned by divergence free basis functions.We define the enriched space as ), where X e N consists of velocity fields constructed in order to guarantee the infsup condition, together with the constraints in (31) and (32); see [9] and [10] for details on how to construct these velocity fields.The inf-sup condition ( 12) is then fulfilled for the spaces M N (Ω) and X N (Ω), and we may solve (23) to find both the velocity and the pressure simultaneously involving a system of size 3N.Alternatively, we could solve the two separate N-sized problems (34) and for the velocity and pressure, respectively.Note that neither X 0 N (Ω) nor X e N (Ω) is a subset of X N (Ω), and that both methods are non-conforming.

A posteriori error estimation
Since the approximation abilities of the reduced basis method strongly depends on the quality of the precomputed basis functions, we have no a priori knowledge of how well the reduced basis solution for a generic parameter will approximate the actual solution.To get an estimate of how good our solution is we need a posteriori error estimation.Based on the theory developed in [17], and following the strategy of [18], the lower and upper output bounds, s − (u N ) and s + (u N ), for the compliant output was constructed in [9] in the single block case, where For a diffusion operator on the reference domain Ω, where v and w are functions on B, and g(Φ) is a positive function depending on the mapping Φ : Ω → Ω, g(Φ) is chosen such that for some positive real constant α 0 .The reconstructed error e is then defined as the field that satisfies where XN (Ω) = {v • Φ ∈ (P N ( Ω)) 2 , v | Γw = 0}.For the bounds defined by s − (u N ) = l(u N ), and (40) we then get In the non-conforming setting when Ω consists of several blocks, the work is on-going.

Output driven reduction
When we generate the reduced basis, the number of basis functions quickly increases when more parameters are introduced.Potentially we could end up with more than thousand basis functions, which would make the method rather costly and impractical.Fortunately, the basis functions typically contain much redundant information.Different post-processing techniques [5,7,11] may be applied to reduce the number of basis functions needed, while preserving the approximation capabilities of the generated basis.
We follow the method presented in [19], where the output bound gap developed in Section 4 is used to reorder the basis functions, such that the error in the output of interest, s(u), is minimized.We also recall that S N is the set of preselected parameter values.Adapted to geometric parameters, we proceed as follows, separately for the different block structures, i.e., pipe and bifurcation.
Offline we choose an arbitrary parameter µ 1 = µ i ∈ S N , with corresponding geometry B i and basis functions u i , u e i , and p i .These basis functions are saved as v 1 , v e 1 , and q 1 , and they span the spaces and calculate where u N 1 is the resulting reduced basis velocity.The basis functions corresponding to µ 2 are saved as v 2 , v e 2 , and q 2 , and together with v 1 , v e 1 , and q 1 they span the spaces X N 2 and M N 2 .We denote S N 2 = {µ 1 , µ 2 } and repeat the process above for all µ j ∈ S N \S N 2 . In a recursive manner we thus choose µ i with corresponding velocity and pressure basis functions until the maximum bound gap reaches a predefined level.
In the online computation of generic solutions, we then start with µ 1 and its corresponding basis functions.We solve for the reduced basis solution, and calculate the bound gap.If the bound gap is larger than a specified limit, we include the basis functions corresponding to the next parameter in S N .The bound gap limit in the online case has to be larger than the bound gap limit used to sort the basis functions.If the number of basis functions in S N is not too large we may also include all of them to find the reduced basis solution in a nonadaptive fashion.Since the a posteriori analysis for the multi-block case is missing, this is what is done when solving the hierarchical flow system in Figure 1.

Parameterizing the geometries
In ( 20) and ( 21) we see how the geometric mapping Φ enters the steady Stokes equations in a natural way.
To control different instantiations of this mapping we define Φ( Ω; µ), where µ ∈ D ⊂ R P .The P elements in µ are parameters which describe, for example, the length, thickness and opening angle of a bifurcation block.Given µ, a pipe or bifurcation block is constructed by defining its outer edges according to µ, and then the internal nodes are found using a Gordon-Hall algorithm.Finally, for multi-element blocks, all the internal nodes are adjusted according to the smoothing process described above.After all the points are found, the blocks may be rotated to any desired orientation.We note that all corners in the blocks are right angles.
Once the final values of all the nodal points have been computed, the Jacobian, J , of the mapping Φ( Ω; µ) and its determinant, J, are calculated and stored for each node.If we again study ( 20) and ( 21), we see that these quantities appear nonlinearly in the equations.Thus we have nonlinear parameter dependence, which is fundamentally different from most reduced basis applications.The equations are not affine in their parameter dependence either; see [17] for affine, linear parameter dependence, and [3] for non-affine parameter dependence.
Since we are only interested in giving a proof of concept, we choose P = 2 for the bifurcations, and let µ = (µ 1 , µ 2 ).We let the first parameter, µ 1 , define the difference in length between the upper leg of the bifurcation and the lower leg, i.e., µ 1 = L l − L u ; see Figure 3.The second parameter, µ 2 , is taken to be the difference in the opening angle of two legs of the bi- The outline of the bifurcation is defined through its corner points and the length of the body relative to the length of the legs.The opening angle is adjusted by a rigid body rotation of the two cornerpoints of the upper leg around the centerpoint of the inflow boundary.Before any rotation, the difference in length between the two legs is defined by setting the x-coordinates of the corner points of the upper leg, (both are the same before the rotation).The non-linear edges of the bifurcation are constructed such that they pass through the corner points and the common edge point of the two elements sharing an edge, and such that they are perpendicular to the inflow and outflow boundaries.The upper and lower edges are fourth order polynomials, while the edge connecting the two legs is, for optimal flexibility, constructed by the use of cubic splines.
We let the first parameter denote the rotation of the outflow boundary relative to the inflow boundary, see Figure 4 for µ 1 = 0 and µ 1 = − π 2 .The second and third parameters are used to define the length of the inflow and outflow boundaries, and the fourth parameter defines the fluctuation of the wall boundaries.

Numerical examples
We now present some examples of the method applied to different geometric structures.The first is a pipe consisting of three elements, where the solution is found as a linear combination on each element, glued together with Lagrange multipliers.The second is a six-element bifurcation, where the solution is found as a linear combination of global basis functions.The third structure is a hierarchical system consisting of one pipe and three bifurcations.The solution is now found as a linear combination on each block structure, i.e. pipe or six-element bifurcation, and glued together Figure 4: The deformed pipes used to construct the basis functions for the pipe blocks.
Table 1: The reduced basis error on a generic multiblock pipe with three blocks.N = 3N 1 is the total number of degrees-of-freedom in the reduced basis spaces X 0 N , X e N , and M N .N 1 is the number of basis geometries used to generate the basis functions.
with Lagrange multipliers across the block interfaces.The final structure is a "bypass" system with three pipe blocks and two bifurcation blocks.

Pipes
We consider the eight geometries in Figure 4 to be pipes with inflow boundary along the left vertical edge, and outflow boundary along the opposite edge.Each pipe is decomposed into three sub-domains, all of which are regular one-to-one deformations of the reference square Ω = (−1, 1) 2 , and the restriction to each sub-domain of the steady Stokes solutions found on these geometries are stored on Ω.In addition we store the reflection of the solutions across the ξ-axis in Ω.This accounts to solving the steady Stokes equations on the reflection of the geometries across the xaxis.Since the first geometry is symmetric, we thus end up with 15 precomputed solutions for each subdomain.We compute the associated enriched velocity solutions, and solve (23) when Ω is taken to be a generic deformed pipe, decomposed into three subdomains.Since we only have N 1 = 15 basis functions in this case, we do not apply the selection algorithm described in Section 4.1.When we use cubic Lagrange multipliers in both the normal and tangential direction to glue the solution together across the block interfaces, the error of the re- 1.4   duced basis solution for an increasing number of basis functions is as presented in Table 1.

Bifurcations
We consider bifurcations characterized by the length and angle of the upper leg relative to the length and angle of the lower leg.In the tensor product parameter space generated by eight relative lengths and eight relative angles, we generate 64 bifurcations.We precompute the steady Stokes solutions on these bifurcations, and store them on Ω. Again we compute the associated enriched velocity solutions, but before we find the reduced basis solution we apply the selection algorithm described earlier.
The resulting errors in velocity and pressure are presented in Table 2, and we see that the convergence is very good.It is better than the convergence seen for a multi-element pipe in Table 1, both because we don't have any consistency error from the element interfaces, and because the basis bifurcations span the generic bifurcation better than the deformed pipes represent the generic pipe in the previous example.In this single-block case we may apply the a posteriori error analysis described in Section 4, and we compute both the upper and the lower bound gaps.We do this without using the selection algorithm, and the bound  gaps converge as shown in Figure 5.Even without the selection algorithm the convergence is exponential.

Hierarchical flow system
An example of a multi-block domain comprising both pipe blocks and bifurcation blocks, is the complex flow system shown in Figure 1.To precompute the basis solutions, we use the same geometries for both pipes and bifurcations as described above.For the pipes we only use the restrictions to the inflow element, while we for the bifurcations precompute the solutions by adding pipe elements to the inflow and outflow boundaries in order to get the right boundary conditions.Only the restrictions of the solutions to the bifurcation block are stored and used as basis solutions.For the pipe block we use all 15 precomputed solutions, while we for the bifurcation blocks again use the selection process to limit the number of precomputed solutions to 30.To glue the blocks together across block interfaces, we again use Lagrange multipliers.Since each bifurcation block consists of two elements on the interface to an adjacent block, we now use linear Lagrange multipliers defined on one half of the interface.In Table 3 we see how the errors in velocity and pressure behave as the number of basis functions increases.

A "bypass"
As the final example we combine both block structures in the bypass system shown in Figure 6.Here the upper branch illustrates the effect of a clogged vein, while the lower branch is the bypass-vein.To model this domain with the reduced basis element method, we use snapshot solutions computed on three-  is the total number of degrees-of-freedom in each of the reduced basis spaces X 0 N , X e N , and M N .N 1 is the number of basis geometries used to generate the basis functions on the pipe block, N 2 is the number of basis functions used on the bifurcation blocks.domain pipes to generate the basis functions for the pipe blocks.The restriction of the snapshot solutions to each of the three sub-domains are now used as basis functions on their respective pipe block in the bypass system.As basis functions for the bifurcation blocks we use the same basis functions that were used on the hierarchical flow system in the previous example.In this case we have two more block-interfaces compared to the hierarchical flow system, each contributing eight constraints on the reduced basis velocity solution u N (2 constraints in each spatial direction for each half of one interface).We see in Table 4 that the error convergence is good, but if too few basis functions are used we get spurious pressure modes due to the severe constraints on the reduced basis velocity space X N (Ω).In Figure 6 we present a contour plot of the error in the reduced basis pressure solution p N when N 1 = 15 and N 2 = 30.Most of the error is located around the pipe block modeling the clogged vein.Compared to the deformed pipes in Figure 4, used to generate the basis functions for the pipe blocks, this pipe block differs

Future work
We have seen how the reduced basis element method works on the steady Stokes problem when the geometry is considered to be a parameter.In a forthcoming paper we will consider the steady Navier-Stokes equations, and theory for the a posteriori error estimation in the multi-block case.In addition we will incorporate the non-affine theory of [3] in order to do more of the necessary computations in the precomputation stage.Other issues to investigate include the extension to time-dependent problems, possibly with moving boundaries, and extension to three dimensional domains.

Figure 1 :
Figure 1: A hierarchical flow system with one pipeblock and three bifurcation-blocks.

Figure 2 :
Figure 2: The reference bifurcation B constructed from six elements.

Figure 3 :
Figure 3: The parameters used to define the bifurcations.

The number of basis functions log 10 Figure 5 :
Figure 5: The bound gaps when varying two geometric parameters on a single bifurcation.

Figure 6 :
Figure 6: The bypass with three pipe blocks and two bifurcation blocks.

Figure 7 :
Figure 7: The contour of the error in the reduced basis pressure solution p N when N 1 = 15 and N 2 = 30.

Table 2 :
The reduced basis error on a single bifurcation.N is the total number of degrees-of-freedom in the reduced basis spaces X 0 N , X e N , and M N .

Table 3 :
The error in the reduced basis steady Stokes solution on a multi-block system corresponding to Figure1.N = N 1 + 3N 2 is the total number of degrees-offreedom in the reduced basis spaces X 0 N , X e N , and M N .N 1 is the number of basis geometries used to generate the basis functions on the pipe block, N 2 is the number of basis functions used on the bifurcation blocks.

Table 4 :
The error in the reduced basis steady Stokes solution on a multi-block bypass with three pipe blocks and two bifurcation blocks.N = 3N 1 + 2N 2