AN EFFICIENT SOLVER FOR MIXED PROBLEMS

YUNRONG ZHU

1. Problem Description

We present here the augmented Lagrangian method for solving systems arising from mixed finite element discretization of elliptic boundary value problem:

Δp = f in Ω, p|∂Ω = 0.
(1.1)

The aim is to show that implementing an efficient iterative method for the resulting indefinite linear system reduces to designing an efficient method for the solution of an auxiliary nearly singular H(div) problem.

The model problem (1.1) can be casted into the following mixed problem by introducing the variable u := ∇p: Find (u,p) ∈ H(div) × L2(Ω) such that

{ (u,v) + (p,divv) = 0, ∀v ∈ H(div) (divu,q) = (f,q), ∀q ∈ L2 (Ω ).
(1.2)

Given a conforming triangulation Th, the mixed finite element method is to solve the model problem (1.2) in the finite element spaces: V h(div) ⊂ H(div) and V h(0) ⊂ L2(Ω). That is,to find (uh,ph) ∈ V h(div) × V h(0) such that

{ (uh,vh) + (ph,divvh) = 0, ∀vh ∈ Vh(div) (divuh,qh) = (f,qh), ∀qh ∈ Vh(0).
(1.3)

Here we restrict ourselves to the lowest order Raviart-Thomas spaces for V h(div).

The mixed finite element method (1.3) results in the following linear system:

[A B*][u] [0] B 0 p = f .
(1.4)

It is not difficult to see that A is the mass matrix of the Raviart-Thomas element and B is the matrix representations of div*.

2. Augmented Lagrangian Method

The augmented Lagrangian method solves the following equivalent problem to (1.4) by the Uzawa method:

[ ][ ] [ ] A + ϵ-1B*B B * u ϵ-1B*f B 0 p = f .
(2.1)

The iteration reads: Given (u(k),p(k)), the new iterate (u(k+1),p(k+1)) is obtained by solving the following systems:

{ (A + ϵ-1B*B )u(k+1) = ϵ- 1B*f - B *p(k), p(k+1) = p(k) - ϵ-1(f - Bu (k+1)).
(2.2)

This algorithm has the following convergence.

Theorem 2.1. Let (u(0),p(0)) be a given initial guess and for k ≥ 1, let (u(k),p(k)) be the iterates obtained via the augmented Lagrangian method. Then the following estimates hold:

( ϵ )k ∥p- p(k)∥0,Ω ≤ ------ ∥p - p(0)∥0,Ω, ϵ+ λ0 ( ) (k) √- (k) √ - --ϵ---k (0) ∥u- u ∥A ≤ ϵ∥p- p ∥0,Ω ≤ ϵ ϵ+ λ0 ∥p - p ∥0,Ω,
where λ0 is the minimum eigenvalue of S = BA-1B*.

According to this theorem, the iteration procedure (2.2) converges very fast to the solution of (1.3) for small ϵ. However, one needs to solve a nearly singular H(div) system

(ϵA + B *B)u(k+1) = B*f - ϵB *p(k).
(2.3)

Thus, an efficient and robust H(div) solver will result in an optimal iterative method for the saddle point problem (1.3). For the nearly singular H(div) system (2.3), we used auxiliary AMG solver.

3. A Model Problem and Sample Mesh Data

Here, we solve equation (1.1) in a unit cube [0,1]3 with structured tetrahedra mesh. The right hand side f = 1. Each cube is partitioned into 6 tetrahedrons.

PIC PIC

Three sample data set of the mesh are given in mesh*.mat files. The data files *.mat are the corresponding stiffness and mass matrices corresponding to the mesh*.mat files.
  • mesh8.mat 8.mat: The data files for the structure mesh with 8 intervals in each direction.
  • mesh17.mat 17.mat: The data files for the structure mesh with 17 intervals in each direction.
  • mesh26.mat 26.mat: The data files for the structure mesh with 26 intervals in each direction.
Click here for a full description of the source codes, to download of the codes and data files, as well as the performance of the algorithms.