MIXED-HYBRID FINITE ELEMENT APPROXIMATION OF THE POTENTIAL FLUID FLOW PROBLEM

J. MARYŠKA, M. ROZLOŽNÍK and M. TŮMA^{∗}

Institute of Computer Science, Academy of Sciences of the Czech Republic, Pod vodárenskou věží 2, 182 07 Prague 8, Czech Republic

Abstract: In the paper a mixed-hybrid approximation of the potential fluid flow problem based on prismatic discretization of the domain is presented. Trilateral prismatic elements with vertical faces and nonparallel bases suitable for the modelling of real geological circumstances are considered. The set of linearly independent vector basis functions is defined and existence and uniqueness of the approximate solution from the resulting symmetric indefinite system are examined. Possible approaches to the solution of the discretized system are discussed.

Keywords: potential flow problem in porous media, mixed-hybrid formulation, general prismatic elements, symmetric indefinite linear systems.

AMS classification: 65N30, 65K10, 73C99

1 Introduction

Solution of the underground water flow problem in real conditions must reflect complex geological structure of sedimented minerals. Layers of stratified rocks with substantially different physical properties must be modelled using an appropriate discretization of the geological region. These geological characteristics can be correspondingly described by the mixed-hybrid finite element method using trilateral prismatic elements with vertical faces and generally nonparallel bases.

An outline of the paper is as follows. First, we introduce the mathematical formulation of the problem. In Section 2, we consider the mixed-hybrid formulation. Finite-dimensional approximation, existence and uniqueness of the approximate solution are derived in Section 3. Finally, different approaches to the solution of the discretized linear system with a symmetric indefinite matrix are discussed and promising ways are proposed.

Let Ω be a bounded domain in R^{3} with a Lipschitz continuous boundary ∂Ω. Potential fluid flow in saturated porous media can be modelled by the velocity u given by Darcy’s law

u = −A^{−1}∇p,  (1.1)

∗This work was supported by the GA CR under grant 201/93/0067.

where p is the piezometric potential (fluid pressure) and A^{−1} is the symmetric and uniformly positive definite second rank tensor of hydraulic permeability of porous medium, i.e., there exists a positive constant α_{0} such that

α_{0}∥ξ∥_{2}^{2} ≤ (A^{−1}(x)ξ, ξ)

holds for all ξ ∈ R^{3} and almost all x ∈ Ω. Further we assume [A^{−1}(x)]_{ij} ∈ L^{∞}(Ω) for all i, j ∈ {1, 2, 3}. Consider also the continuity equation for the incompressible flow

∇ · u = q,  (1.2)

where q represents the density of potential sources in the medium. The boundary conditions are given by

p = p_{D} on ∂Ω_{D},  (1.3)
u · n = −A^{−1}∇p · n = u_{N} on ∂Ω_{N},  (1.4)

where ∂Ω = \overline{∂Ω_{D}} ∪ \overline{∂Ω_{N}} are such that ∂Ω_{D} ∩ ∂Ω_{N} = ∅ and n is the outward normal vector defined (almost everywhere) on the boundary ∂Ω.

A remark on the notation. We denote by L^{2}(Ω) the Lebesgue space defined as

L^{2}(Ω) = {ϕ : Ω → R; ∫_{Ω} |ϕ|^{2}dx < ∞}

with the scalar product (ϕ_{1}, ϕ_{2})_{0,Ω} = ∫_{Ω} ϕ_{1}ϕ_{2}dx and the standard norm ∥ϕ∥_{0,Ω} = (ϕ, ϕ)0,Ω^{\frac{1}{2}} . Further we denote by L^{2}(Ω) the Lebesgue space of vector valued functions v, where the components v_{i}, i = 1, …, 3, belong to the space L^{2}(Ω), and consider the Sobolev space

H^{1}(Ω) = {φ ∈ L^{2}(Ω); ∇φ ∈ L^{2}(Ω)}

with the scalar product (φ_{1}, φ_{2})_{1,Ω} = ∫_{Ω} [φ_{1}φ_{2}+∇φ_{1}·∇φ_{2}]dx and the norm ∥φ∥_{1,Ω} = (φ, φ)1,Ω^{\frac{1}{2}} . We introduce the space of vector valued functions

H(div, Ω) = {v ∈ L^{2}(Ω); ∇ · v ∈ L^{2}(Ω)}

with the norm defined as ∥ v ∥_{div,Ω} = (∥ v ∥_{0,Ω}^{2} + ∥ ∇ · v ∥_{0,Ω}^{2} ) ^{\frac{1}{2}} . We will also introduce the bilinear form < ϕ, μ >_{∂Ω} = ∫_{∂Ω} ϕμdS, where ϕ and μ are functions from L^{2}(∂Ω).

2 Mixed-hybrid formulation of the problem

Denote the collection of subdomains of the domain Ω by E_{h} and the collection of faces of subdomains e ∈ E_{h} which are not adjacent to the boundary ∂Ω_{D} by Γ_{h} = ∪_{e∈E_{h}}∂e − ∂Ω_{D}, where h is the discretization parameter (see [4])

h = max_{e∈E_{h}}{diam e}.

Denote the restriction of any function on the subdomain e ∈ E_{h} by the superscript e, i.e. ϕ^{e} = ϕ|_{e}. Let us introduce the functional spaces defined on E_{h} and Γ_{h}. Let H(div, E_{h}) be the space of square integrable vector functions v ∈ L^{2}(Ω), whose divergences are square integrable on every subdomain e ∈ E_{h}, i.e.

H(div, E_{h}) = {v ∈ L^{2}(Ω); ∇ · v^{e} ∈ L^{2}(e), ∀e ∈ E_{h}}  (2.1)

with the norm given as

∥ v ∥_{div,E_{h}}= [∥ v ∥_{0,Ω}^{2} + ∑_{e∈E_{h}} ∥ ∇ · v^{e} ∥_{0,e}^{2}]^{\frac{1}{2}} .  (2.2)

We also consider the space of traces

HD^{\frac{1}{2}}(Γ_{h}) = \{μ : Γ_{h} → R; ∃φ ∈ H_{D}^{1}(Ω), μ = γ_{h}φ\},  (2.3)

where the space H_{D}^{1}(Ω) is defined as H_{D}^{1}(Ω) = {φ ∈ H^{1}(Ω); γφ = 0 on ∂Ω_{D}} and γφ = φ|_{∂Ω} is the trace of the function φ ∈ H^{1}(Ω) on the boundary ∂Ω; γ_{h}φ = φ|_{Γ_{h}} is the trace of the

function φ ∈ H^{1}(Ω) on the structure of faces Γ_{h}. The space HD^{\frac{1}{2}}(Γ_{h}) is equipped with the norm

∥ μ ∥_{\frac{1}{2} ,Γ_{h}} = inf_{φ∈H_{D}^{1}(Ω)} \{|φ|_{1,Ω} ; γ_{h}φ = μ on Γ_{h}\},  (2.4)

where |φ|_{1,Ω} denotes the seminorm |φ|_{1,Ω} = (∇φ, ∇φ)0,Ω^{\frac{1}{2}} .

Thus, the mixed-hybrid formulation of the problem (1.1), (1.2) with the boundary conditions (1.3), (1.4) and the discretization E_{h} of the domain Ω can be stated as follows (see also [12], [9]):

F ind (u, p, λ) ∈ H(div, E_{h}) × L^{2}(Ω) × HD^{\frac{1}{2}}(Γ_{h}) such that :

∑_{e∈E_{h}}{(A^{−1}u^{e}, v^{e})_{0,e} − (p^{e}, ∇ · v^{e})_{0,e}+ < λ^{e}, n^{e} · v^{e} >_{∂e∩Γ_{h}} } =  (2.5)

= ∑_{e∈E_{h}} < p_{D}^{e}, n^{e} · v^{e} >_{∂e∩∂Ω_{D}} , ∀v ∈ H(div, E_{h});

− ∑_{e∈E_{h}}(∇ · u^{e}, ϕ^{e})_{0,e} = − ∑_{e∈E_{h}}(q^{e}, ϕ^{e})_{0,e}, ∀ϕ^{e} ∈ L_{2}(Ω);  (2.6)
∑_{e∈E_{h}} < n^{e} · u^{e}, μ^{e} >_{∂e}= ∑_{e∈E_{h}} < u_{N}^{e}, μ^{e} >_{∂e∩∂Ω_{N}} , ∀μ ∈ HD^{\frac{1}{2}}(Γ_{h}).  (2.7)

3 Discretization of the domain and finite-dimensional approximation

In this section we introduce the discretization of the domain Ω and the lowest order finitedimensional approximation of (2.5)-(2.7).

Assume from now that the domain Ω is a polyhedron and is subdivided into a collection of subdomains such that every subdomain is a trilateral prism with suitably chosen vertices x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6},

\begin{matrix}x_{1} = (x_{1}, y_{1}, z_{1}) , & x_{2} = (x_{2}, y_{2}, z_{2}) , & x_{3} = (x_{3}, y_{3}, z_{3}) , \\ x_{4} = (x_{1}, y_{1}, z_{4}) , & x_{5} = (x_{2}, y_{2}, z_{5}) , & x_{6} = (x_{3}, y_{3}, z_{6}) . \end{matrix}

We allow also elements, that approximate boundary of the domain, such that for some i ∈ {1, 2, 3} is x_{i} = x_{i+3} or there exist i, j ∈ {1, 2, 3}, i ≠ j, such that x_{i} = x_{i+3} and x_{j} = x_{j+3} .

[figure: x_{1}; x_{3}; x_{4}; x_{6}; x_{2}; x_{5}; n_{1}^{e}; n_{2}^{e}; n_{4}^{e}; n_{5}^{e}; n_{3}^{e}]

FIG 3.1

We will assume that the mesh obtained is strongly regular, i.e., there exists a positive constant ζ independent of the mesh size h such that d′/d ≥ ζ holds for every element edges d, d′ from the decomposition E_{h}. Denote by η the minimum angle of the triangulation obtained as a horizontal projection of the prismatic mesh.

The velocity function u will be approximated by vector functions linear on every element e ∈ E_{h}. We define the Raviart-Thomas space RT^{0}(e),

RT^{0}(e) = {v^{e}; v^{e}(x) = ∑_{i=1}^{5}ν_{j}v_{j}(x), x ∈ e},  (3.1)

with linearly independent basis functions v_{j}^{e}, j = 1, ..., 5, of the form

v_{1}^{e} = k_{1}^{e} [ \begin{matrix}0 \\ 0 \\ x_{3} − α_{13}^{e} \end{matrix} ] , v_{2}^{e} = k_{2}^{e} [ \begin{matrix}0 \\ 0 \\ x_{3} − α_{23}^{e} \end{matrix} ] ,  (3.2)

v_{3}^{e} = k_{3}^{e} [ \begin{matrix}x_{1} − α_{31}^{e} \\ x_{2} − α_{32}^{e} \\ β_{3}^{e}x_{3} − α_{33}^{e} \end{matrix} ] , v_{4}^{e} = k_{4}^{e} [ \begin{matrix}x_{1} − α_{41}^{e} \\ x_{2} − α_{42}^{e} \\ β_{4}^{e}x_{3} − α_{43}^{e} \end{matrix} ] , v_{5}^{e} = k_{5}^{e} [ \begin{matrix}x_{1} − α_{51}^{e} \\ x_{2} − α_{52}^{e} \\ β_{5}^{e}x_{3} − α_{53}^{e} \end{matrix} ]

such that

F_{j}(v_{i}^{e}) = ∫_{f_{j}^{e}} n_{j}^{e} · v_{i}^{e}dS = δ_{ij}, i, j = 1, ..., 5.  (3.3)

Here f_{j}^{e} denotes the jth face of the element e and n_{j}^{e} its outward normal vector (with respect to the element e).

LEMMA 3.1: The system of functional equations

F_{j}(v_{i}^{e}) := ∫_{f_{j}^{e}} n_{j}^{e} · v_{i}^{e}dS = δ_{ij}, i, j = 1, …, 5 ,

generates the unique set of basis vector functions of the form (3.2).

PROOF: Substituting (3.2) into (3.3) we obtain

∫_{f_{j}^{e}} k_{i}^{e} [n_{j1}^{e}x_{1} + n_{j2}^{e}x_{2} + n_{j3}^{e}β_{i}^{e}x_{3} − n_{j1}^{e}α_{i1}^{e} − n_{j2}^{e}α_{i2}^{e} − n_{j3}^{e}α_{i3}^{e}]dS = δ_{ij}.  (3.4)

For i ≠ j we get the condition

n_{j1}^{e}α_{i1}^{e} + n_{j2}^{e}α_{i2}^{e} + n_{j3}^{e}α_{i3}^{e} − n_{j3}^{e}x_{T 3}(f_{j}^{e})β_{i}^{e} = n_{j1}^{e}x_{T 1}(f_{j}^{e}) + n_{j2}^{e}x_{T 2}(f_{j}^{e}),

where x_{T} (f_{j}^{e}) denotes the centre of gravity of the jth face. Thus for i = 1, 2 we have

α_{i3}^{e} = x_{T 3}(f_{j}^{e}), i = 1, j = 2 or i = 2, j = 1.  (3.5)

From (3.4) we obtain

k_{i}^{e} = {n_{i3}^{e}[x_{T 3}(f_{i}^{e}) − x_{T 3}(f_{j}^{e})]|f_{i}^{e}|}^{−1},  (3.6)

where |f_{i}^{e}| denotes the area of the ith face, so |f_{i}^{e}| > 0. For i = 3, 4, 5 we obtain the system of equations for the unknowns α_{i1}^{e}, α_{i2}^{e}, α_{i3}^{e} and β_{i}^{e};

n_{j1}^{e}α_{i1}^{e} + n_{j2}^{e}α_{i2}^{e} + n_{j3}^{e}α_{i3}^{e} − n_{j3}^{e}x_{T 3}(f_{j}^{e})β_{i}^{e} = n_{j1}^{e}x_{T 1}(f_{j}^{e}) + n_{j2}^{e}x_{T 2}(f_{j}^{e}),
j = 1, 2 ,  (3.7)
n_{j1}^{e}α_{i1}^{e} + n_{j2}^{e}α_{i2}^{e} = n_{j1}^{e}x_{T 1}(f_{j}^{e}) + n_{j2}^{e}x_{T 2}(f_{j}^{e}),
j = 3, 4, 5; j ≠ i.  (3.8)

Arbitrary two vectors from the set {n_{3}^{e}, n_{4}^{e}, n_{5}^{e}} are linearly independent. Hence, the system (3.8) has a unique solution a_{i(1,2)}^{e} = (α_{i1}^{e}, α_{i2}^{e}). Consider the subdeterminant

n_{13}^{e}n_{23}^{e} | \begin{matrix}1 & −x_{T 3}(f_{1}^{e}) \\ 1 & −x_{T 3}(f_{2}^{e}) \end{matrix} | = n_{13}^{e}n_{23}^{e}[x_{T 3}(f_{1}^{e}) − x_{T 3}(f_{2}^{e})] ≠ 0.

Then substituting for α_{i1}^{e} and α_{i2}^{e}, the system (3.7) has the unique solution (α_{i3}^{e}, β_{i}^{e}). For k_{i}^{e} we have

k_{i}^{e}[n_{i}^{e} · x_{i}(f_{i}^{e}) − n_{i}^{e} · a_{i(1,2)}^{e}]|f_{i}^{e}| = 1.  (3.9)

The coefficient k_{i}^{e} is invariant of translation. We use a translation of the element e to element ê such that x_{3} = fê ; j ∩ fê ; i for some j ≠ i. Then the right-hand side of the system (3.8) is equal to zero, i.e. nê ; j1x_{T 1}(fê ; j) + nê ; j2x_{T 2}(fê ; j) = 0, and so aê ; i(1,2) = (0, 0). Estimate then

nê ; i · x_{i}(fê ; i) ≥ |nê ; i||x_{i}(fê ; i)| sinη > 0.

Therefore

k_{i}^{e} = kê ; i = [nê ; i · x_{i}(fê ; i)]^{−1}|fê ; i|^{−1} > 0.

□

We also define

RT _{−1}^{0}(E_{h}) = {v_{h} ∈ L^{2}(Ω); v_{h}^{e} ∈ RT^{0}(e), ∀e ∈ E_{h}},  (3.10)

the space that consists of vector functions linear on every element. We note that in the case of nonparallel bases these functions are not continuous across the interelement boundaries Γ_{h}. Denote the space of scalar functions constant on the element e by M^{0}(e) with the basis function of the form:

ϕ_{h}^{e}(x) = 1, x ∈ e; ϕ_{h}^{e}(x) = 0, x ∉ e.

Then we introduce the space

M_{−1}^{0}(E_{h}) = {ϕ_{h} ∈ L^{2}(Ω); ϕ_{h}^{e} ∈ M^{0}(e), ∀e ∈ E_{h}}  (3.11)

which consists of elementwise constant functions that will approximate the piezometric potential p. Let M^{0}(f) be the space of constant functions on the interelement face f ∈ Γ_{h} and define the space M_{−1}^{0}(Γ_{h}),

M_{−1}^{0}(Γ_{h}) = {μ_{h} : Γ_{h} → R; μ_{h}^{f} ∈ M^{0}(f), ∀f ∈ Γ_{h}},  (3.12)

which consists of functions constant on every face from Γ_{h}. Further let p_{D,h}, u_{N,h} be the functions from M_{−1}^{0}(∂Ω) which approximate the functions p_{D} and u_{N} given in the boundary conditions and which satisfy

∫_{f}(p_{D,h} − p_{D})dS = 0; ∀f ∈ ∂Ω_{D},  (3.13)

∫_{f}(u_{N,h} − u_{N})dS = 0; ∀f ∈ ∂Ω_{N}.  (3.14)

Then the Raviart-Thomas approximation of the mixed-hybrid formulation for the problem reads as follows (we refer also to [9] or [10]):

F ind (u_{h}, p_{h}, λ_{h}) ∈ RT _{−1}^{0}(E_{h}) × M_{−1}^{0}(E_{h}) × M_{−1}^{0}(Γ_{h}) such that

∑_{e∈E_{h}}{(A^{−1}u_{h}, v_{h})_{0,e} − (p_{h}, ∇ · v_{h})_{0,e}+ < λ_{h}, n^{e} · v_{h} >_{∂e∩Γ_{h}} } =

=< p_{D,h}, n^{e} · v_{h} >_{∂e}; ∀v_{h} ∈ RT _{−1}^{0}(E_{h}).  (3.15)
− ∑_{e∈E_{h}}(∇ · u_{h}, ϕ_{h})_{0,e} = −(q_{h}, ϕ_{h})_{0,Ω} ; ∀ϕ_{h} ∈ M_{−1}^{0}(E_{h}).  (3.16)
∑_{e∈E_{h}} < n^{e} · u_{h}, μ_{h} >_{∂e}=< u_{N,h}, μ_{h} >_{∂Ω_{N}} ; ∀μ_{h} ∈ M_{−1,D}^{0}(Γ_{h}).  (3.17)

Next, the system of linear equations equivalent to the problem (3.15) - (3.17) will be derived. Let e_{j} ∈ E_{h}, j = 1, …, J, be the numbered system of prismatic elements. On every element we have defined the five-dimensional space RT^{0}(e) of linear vector functions v_{i}^{e}, i = 1, …, 5. The finite-dimensional space RT _{−1}^{0}(E_{h}) is then spanned by 5 × J linearly independent basis functions v_{i}, i = 1, …, I = 5 × J. Let f_{k}, k = 1, …, K, be the numbered system of interelement faces from Γ_{h}. By the introduced approximation, the functions u_{h}, p_{h} and λ_{h} belong to spaces RT _{−1}^{0}(E_{h}), M_{−1}^{0}(E_{h}) and M_{−1}^{0}(Γ_{h}), respectively, and so they can be expanded in the form

u_{h}(x) = ∑_{i∈I} ũ_{i}v_{i}(x), p_{h}(x) = ∑_{j∈J} p̃_{j}ϕ_{j}(x), x ∈ Ω,

λ_{h}(x) = ∑_{k∈K} λ̃_{k}μ_{k}(x), x ∈ Γ_{h}.

We put u = (ũ_{1}, ..., ũ_{I})^{T} , p = (p̃_{1}, ..., p̃_{J})^{T} , λ = (λ̃_{1}, ..., λ̃_{K})^{T} and

A_{ij} = (A^{−1}v_{i}, v_{j})_{0,Ω} ; i = 1, …, I, j = 1, …, I,  (3.18)
B_{ij} = −(∇ · v_{i}, 1)_{0,e_{j}}; i = 1, …, I, j = 1, …, J,  (3.19)
C_{ik} =< n_{k} · v_{i}, 1 >_{f_{k}}; i = 1, …, I, k = 1, …, K.  (3.20)

Here n_{k} is the outward normal vector to the face f_{k} with respect to the element corresponding to the support of the function v_{i}. Then we compute

[q_{1}]_{i} = − < p_{D,h}, n_{i} · v_{i} >_{∂Ω_{D}} ; i = 1, …, I,  (3.21)
[q_{2}]_{j} = −(q, 1)_{0,e_{j}}; j = 1, …, J,  (3.22)
[q_{3}]_{k} = < u_{N,h}, 1 >_{f_{k}}; k = 1, …, K.  (3.23)

Substituting u_{h}, p_{h} and λ_{h} into (3.15) - (3.17), we can now write the system of linear equations

( \begin{matrix}A & B & C \\ B^{T} \\ C^{T} \end{matrix} ) ( \begin{matrix}u \\ p \\ λ \end{matrix} ) = ( \begin{matrix}q_{1} \\ q_{2} \\ q_{3} \end{matrix} ) .  (3.24)

LEMMA 3.2: Assuming ∂Ω_{D} ≠ ∅, the matrix (B C) ∈ R^{I,J+K} defined in (3.19) and (3.20) has the full column rank, i.e. rank (B C) = J + K.

PROOF: Suppose Bp + Cλ = 0 holds for some vectors p and λ. Then v^{T} (Bp + Cλ) = 0, ∀v ∈ R^{I}. Equivalently, we have

− ∑_{e∈E_{h}}(∇ · v, p)_{0,e} + ∑_{e∈E_{h}} < n^{e} · v, λ >_{∂e∩Γ_{h}} = 0, ∀v ∈ RT _{−1}^{0}(E_{h}) .  (3.25)

Applying the Green formula to the first term in (3.25), we have for all v ∈ RT _{−1}^{0}(E_{h})

− ∑_{e∈E_{h}} < n^{e} · v, p >_{∂e} + ∑_{e∈E_{h}} < n^{e} · v, λ >_{∂e∩Γ_{h}} =

= ∑_{e∈E_{h}} < n^{e} · v, λ − p >_{∂e∩Γ_{h}} − ∑_{e∈E_{h}} < n^{e} · v, p >_{∂e∩∂Ω_{D}} = 0.  (3.26)

Because the Dirichlet boundary condition is defined at least on one face from ∂Ω_{D}, there exists ē ∈ E_{h} such that ∂ē∩∂Ω_{D} ≠ ∅. Then for some v̄, ē ⊂ supp(v̄) we have < nē ·v, p >_{∂ē∩∂Ω_{D}} = 0 and so p = 0 on ē. Consequently p = 0 also on some faces belonging to Γ_{h}. Since the first sum in (3.26) implies λ = p on Γ_{h} we get p = 0 on all e ∈ E_{h} and λ = 0 on all f ∈ Γ_{h}. □

4 Solution of the discretized linear system

A considerable interest has been devoted to the solution of the set of linear equations (3.24) in recent years. These systems arise frequently, e.g. from mixed finite element or finite difference discretizations of the Stokes equations in computational fluid dynamics or other second order elliptic problems.

We will briefly recall some possible ways for the solution of (3.24) and particularly we will concentrate on some approaches which we consider to be promising.

Consider first the Uzawa-like approach (see [6], [5]), often advocated as an efficient solution technique. These algorithms are in fact variants of some classical iterative schemes applied to the system of linear equations. As an example, the inexact Uzawa scheme for (3.24) uses the splitting matrix of the form

M = ( \begin{matrix}Q_{A} \\ B^{T} & \frac{1}{α}Q_{MB} \\ C^{T} & & \frac{1}{β} Q_{MC} \end{matrix} ) ,  (4.1)

where Q_{A} is an approximation of the matrix A; Q_{MB} and Q_{MC} are some preconditioning matrices and α, β are fixed parameters. This leads to the iterative scheme:

u_{i+1} = u_{i} + Q_{A}^{−1}(q_{1} − (Au_{i} + Bp_{i} + Cλ_{i})) ,
p_{i+1} = p_{i} + αQ_{MB}^{−1}(q_{2} − B^{T} u_{i+1} ) ,
λ_{i+1} = λ_{i} + βQ_{MC}^{−1}(q_{3} − C^{T} u_{i+1} ) .

Note that in the “exact” case (Q_{A} = A) this is a first order Richardson iteration with the two fixed parameters applied to the Schur complement system for the unknowns (p^{T} , λ^{T})^{T} preconditioned by the matrix diag(Q_{MB}, Q_{MC}).

Since the Uzawa-like algorithms are stationary iterative methods it is only natural to apply the standard and powerful nonstationary conjugate gradient method for the solution of systems with a symmetric positive definite matrix (see [8]). We have two possible ways to reduce the problem (3.24) to subproblems with a symmetric and positive definite matrix. We can use a partial substitution for the unknown u to obtain the Schur complement system for the unknowns (p^{T} , λ^{T})^{T}. This approach was discussed in [9]. Similarly, two successive substitutions were used in [2], where Schur complement system for the unknown vector λ was solved.

Although a block diagonal block matrix can be easily invertible in scalar computing environment, another efficient strategy may be based on the solution of the global system (3.24) by the preconditioned conjugate gradient method. The system (3.24) is, however, symmetric indefinite. Motivated by [3] and [5], suppose that Q_{A} satisfies

1 < α_{1} ≤ \frac{(v, Av)}{(v, Q_{A}, v)} ≤ α_{2}

and let Q_{T} be an approximation of the matrix

T = ( \begin{matrix}(A − Q_{A})Q_{A}^{−1}A & (A − Q_{A})Q_{A}^{−1}B \\ B^{T} Q_{A}^{−1}(A − Q_{A}) & B^{T} Q_{A}^{−1}B \end{matrix} )

such that

1 < β_{1} ≤ \frac{(v, T v)}{(v, Q_{T} , v)} ≤ β_{2}.

Then the product of matrices

( \begin{matrix}A − Q_{A} \\ & I \\ & & I \end{matrix} ) ( \begin{matrix}Q_{A}^{−1} \\ B^{T} Q_{A}^{−1} & −I \\ C^{T} Q_{A}^{−1} & & −I \end{matrix} ) ( \begin{matrix}A & B & C \\ B^{T} \\ C^{T} \end{matrix} )  (4.2)

is a symmetric positive definite matrix (see [3]). Premultiplying (4.2) by the matrix

\binom{Q_{T}^{−1}}{C^{T} Q_{A}^{−1} (A − Q_{A} B)Q_{T}^{−1} −I} ,

we obtain a symmetric positive matrix with respect to the inner product

[( \begin{matrix}x \\ p \\ λ \end{matrix} ) , ( x̃ ; p̃ ; λ̃ )] = (T ( \begin{matrix}x \\ p \end{matrix} ) , ( x̃ ; p̃ )) − (Q_{T} ( \begin{matrix}x \\ p \end{matrix} ) , ( x̃ ; p̃ )) + (λ, λ̃).

Consequently, the matrix of the system (3.24) is symmetrizable and, following Hageman-Young [7], we can apply the conjugate gradient method.

Another conjugate gradient-type method, which can be applied also to a symmetric indefinite system is the MINRES method presented in [13]. This strategy, based on the preconditioned MINRES scheme, has been tested in [14] and [15], where different types of preconditioners were investigated. Our preconditioner is based on the incomplete Bunch-Parlett decomposition of (3.24), which is obtained from the left-looking algorithm based on the directed graph model (see [1]) to get a structure of the rows of the Bunch-Parlett factor L.

Detailed results of numerical experiments and a comparison of different approaches used in our underground water flow applications will be published in the forthcoming paper [11].

References

[1] A.V. Aho, J.E. Hopcroft, J.D. Ullman: Data Structures and Algorithms, Addison-Wesley, Reading, MA, 1983.

[2] L. Bergamaschi, S. Mantica, F. Saleri: Mixed finite element approximation of Darcy’s law in porous media, manuscript, 1994.

[3] J.H. Bramble, J.E. Pasciak: A preconditioning technique for indefinite systems resulting from mixed approximations of elliptic problems, Math. Comp. 50 (1988), 1–17.

[4] P.G. Ciarlet: The Finite Element Method for Elliptic Problems, North-Holland Publishing Company, Amsterdam, 1978.

[5] H.C. Elman: Multigrid and Krylov subspace methods for the discrete Stokes equations, Tech. Rep. UMIACS-TR-94-76, IACS, University of Maryland, 1994.

[6] H.C. Elman, G.H. Golub: Inexact and preconditioned Uzawa algorithms for saddle point problems, Tech. Rep. UMIACS-TR-93-41, IACS, University of Maryland, 1993.

[7] L.A.Hageman, D.M.Young: Applied Iterative Methods, Academic Press, New York, 1981.

[8] M.R. Hestenes, E.Stiefel: Methods of conjugate gradients for solving linear systems, J. Res. Nat. Bureau of Standards 49 (1952), 409–436.

[9] E.F. Kaasschieter, A.J.M. Huijben: Mixed-hybrid finite elements and streamline computation for the potential flow problem, Report PN-90-02-A, TNO Institute of Applied Geoscience, Delft, 1990.

[10] J.Maryška, D.Frydrych: Mixed-hybrid method finite elements for the potential flow problem, Proceedings of Mathematics, Liberec, 1993.

[11] J. Maryška, M. Rozložník, M. Tůma: Numerical experiments with the mixed-hybrid approximation in the underground water flow modelling, in preparation.

[12] J.T. Oden, J.K. Lee: Dual-mixed hybrid finite element method for second-order elliptic problems, In: Mathematical Aspects of Finite Element Methods, Lecture Notes in Mathematics 606, Ed. I. Galliani, E. Magenes, Springer-Verlag, Berlin, 1977, 275-291.

[13] C.C. Paige, M.A.Saunders: Solution of sparse indefinite systems of linear equations, SIAM J. Num. Anal. 12 (1975), 617–629.

[14] A. Ramage, A.J. Wathen: Iterative solution techniques for the Stokes and Navier-Stokes equations, Int. J. Numer. Methods Fluids 19 (1994), 67–83.

[15] T. Rusten, R. Winther: A preconditioner for saddle point problems, in: T. Rusten: Iterative Methods for Mixed Finite Element Systems, PhD. Thesis, University of Oslo,