The Bramble-Pasciak Transformation

91. The Bramble-Pasciak Transformation#

In this lecture we derive the BP-Transfromation, see Bramble-Pasciak: A preconditioning technique for indefinite systems resulting from mixed approximations of elliptic problems, 1988

For an application of the BP-transformation see application.

We consider the saddle point system

\begin{array}{r} (\begin{array}{ccc} A & B^{T} \\ B & 0 \end{array}) (\begin{array}{c} u \\ p \end{array}) = (\begin{array}{c} f \\ g \end{array}) . \end{array}

Let $\hat{A}$ be a preconditioner to $A$ scaled such that

\hat{A} < A \leq α \hat{A}

Thanks to the scaling, the matrix $A - \hat{A}$ is positive definite. The spectrum of the preconditioned system satisfies $σ ({\hat{A}}^{- 1} A) \subset (1, α]$ .

We follow the idea of block Gaussian elimination. Instead of multiplying the first row by $A^{- 1}$ , we multiply by the preconditioning action ${\hat{A}}^{- 1}$ and obtain

\begin{array}{r} (\begin{array}{ccc} {\hat{A}}^{- 1} A & {\hat{A}}^{- 1} B^{T} \\ B & 0 \end{array}) (\begin{array}{c} u \\ p \end{array}) = (\begin{array}{c} {\hat{A}}^{- 1} f \\ g \end{array}) . \end{array}

If $\hat{A}$ is a good preconditioner, the upper-left block is close to the identity matrix. The next step in Gaussian elimination is to zero the lower-left block. We subtrace $B$ times the first from from the second row (and change sign):

\begin{array}{r} (\begin{array}{ccc} {\hat{A}}^{- 1} A & {\hat{A}}^{- 1} B^{T} \\ B {\hat{A}}^{- 1} A - B & B {\hat{A}}^{- 1} B^{T} \end{array}) (\begin{array}{c} u \\ p \end{array}) = (\begin{array}{c} {\hat{A}}^{- 1} f \\ B {\hat{A}}^{- 1} f - g \end{array}) . \end{array}

finally, we multiply the first row by the positive definite matrix $A - \hat{A}$ to obtain

\begin{array}{r} (\begin{array}{ccc} (A - \hat{A}) {\hat{A}}^{- 1} A & (A - \hat{A}) {\hat{A}}^{- 1} B^{T} \\ B {\hat{A}}^{- 1} (A - \hat{A}) & B {\hat{A}}^{- 1} B^{T} \end{array}) (\begin{array}{c} u \\ p \end{array}) = (\begin{array}{c} (A - \hat{A}) {\hat{A}}^{- 1} f \\ B {\hat{A}}^{- 1} f - g \end{array}) \end{array}

This matrix is symmetric. The next theorem shows that it is positive definite as well.

All transformation together can be combined to the transformation matrix

\begin{array}{r} (\begin{array}{ccc} A - \hat{A} & 0 \\ 0 & I \end{array}) (\begin{array}{ccc} \hat{I} & 0 \\ B & - I \end{array}) (\begin{array}{ccc} {\hat{A}}^{- 1} & 0 \\ 0 & I \end{array}) \end{array}

Theorem

Let

$\begin{array}{r} K = (\begin{array}{cc} (A - \hat{A}) {\hat{A}}^{- 1} A & (A - \hat{A}) {\hat{A}}^{- 1} B^{T} \\ B {\hat{A}}^{- 1} (A - \hat{A}) & B {\hat{A}}^{- 1} B^{T} \end{array}) \end{array}$

and

$\begin{array}{r} C = (\begin{array}{cc} A - \hat{A} \\ 0 & B {\hat{A}}^{- 1} B^{T} \end{array}) \end{array}$

Then there hold the spectral estimates

$γ_{1} C \leq K \leq γ_{2} C$

with some $γ_{2} < 2$ and $γ_{1} = O (α^{- 1})$ .

Proof: We expand the quadratic form of $K$ , and use Young’s inequality with some $γ \in R^{+}$ for the second term

\begin{aligned} (\begin{array}{c} u \\ p \end{array}) K (\begin{array}{c} u \\ p \end{array}) & = u^{T} (A - \hat{A}) {\hat{A}}^{- 1} A u + 2 p^{T} B {\hat{A}}^{- 1} (A - \hat{A}) u + p^{T} B {\hat{A}}^{- 1} B^{T} p \\ \leq u^{T} (A - \hat{A}) {\hat{A}}^{- 1} A u + γ u^{T} (A - \hat{A}) {\hat{A}}^{- 1} (A - \hat{A}) u \\ + (1 + γ^{- 1}) p^{T} B {\hat{A}}^{- 1} B^{T} p \end{aligned}

For the $u$ term term we estimate

u^{T} (A - \hat{A}) {\hat{A}}^{- 1} (A + γ (A - \hat{A})) u \leq (α + γ α - γ) u^{T} (A - \hat{A}) u

and for the $p$ -term we already have

(1 + γ^{- 1}) p^{T} B {\hat{A}}^{- 1} B^{T} p \leq (1 + γ^{- 1}) p^{T} B {\hat{A}}^{- 1} B^{T} p

Now we determine the optimal $γ$ from the equation

α + γ α - γ = 1 + γ^{- 1}

which leads to $γ = \frac{1}{2} + \sqrt{\frac{1}{4} + \frac{1}{α - 1}}$ and

K \leq 1 + {(\frac{1}{2} + \sqrt{\frac{1}{4} + \frac{1}{α - 1}})}^{- 1} C \leq 2 C

The lower estimate is proven with a similarly (exercise). $◻$

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
mesh = Mesh( OCCGeometry(shape,dim=2).GenerateMesh(maxh=0.05)).Curve(3)

V = VectorH1(mesh, order=2, dirichlet="wall|inlet|cyl")
V.SetOrder(TRIG,3)
V.Update()
Q = L2(mesh, order=1)

u,v = V.TnT()
p,q = Q.TnT()

bfa = BilinearForm(InnerProduct(grad(u),grad(v))*dx).Assemble()
bfb = BilinearForm(div(u)*q*dx).Assemble()

# prea = Preconditioner(bfa, "local")
prea = Preconditioner(bfa, "direct")
bfa.Assemble()

bfschur = BilinearForm(p*q*dx, diagonal=True).Assemble()
preschur = bfschur.mat.Inverse()

from ngsolve.la import EigenValues_Preconditioner
def BramblePasciakCG(A, B, f, g, preA, preS, maxit=1000, tol=1e-8):
    lam = EigenValues_Preconditioner(A,preA)
    print ("lammin/lammax = ", lam[0], '/', lam[-1])
    preA = 1.2/lam[0]*preA   # scaling
    
    x = BlockVector([f.CreateVector(), g.CreateVector()])
    w,r,p,ap = [x.CreateVector() for i in range(4)]
    
    # r.data = b
    # p.data = pre*r
    pru = (preA * f).Evaluate()
    r[0].data = A*pru - f
    r[1].data = B*pru - g
    p[0].data = pru    
    p[1].data = preS*r[1]
    
    wrn = InnerProduct(r,p)
    err0 = sqrt(wrn)
        
    x[:] = 0
    for it in range(maxit):
        # ap.data = A * p
        hv = (A * p[0] + B.T * p[1]).Evaluate()
        papu = (preA * hv).Evaluate()
        ap[0].data = A * papu - hv
        ap[1].data = B * (papu - p[0])
        
        pap = InnerProduct(p, ap)
        wr = wrn
        alpha = wr / pap
        
        x += alpha * p
        r -= alpha * ap
        pru -= alpha * papu

        # w.data = pre*r
        w[0].data = pru
        w[1].data = preS * r[1]
        
        wrn = InnerProduct(w, r)
        err = sqrt(wrn)
        print ("Iteration",it,"err=",err)
        if err < tol * err0: break
            
        beta = wrn / wr
        
        p *= beta
        p.data += w 
    
    return x[0], x[1]

gfu = GridFunction(V)
gfp = GridFunction(Q)

uin = (1.5*4*y*(0.41-y)/(0.41*0.41), 0)
gfu.Set(uin, definedon=mesh.Boundaries("inlet"))

resf = (-bfa.mat * gfu.vec).Evaluate()
resg = (-bfb.mat * gfu.vec).Evaluate()

sol = BramblePasciakCG (A=bfa.mat, B=bfb.mat, f=resf, g=resg, preA=prea, preS=preschur)

gfu.vec.data += sol[0]
gfp.vec.data += sol[1]

Draw (gfu)
Draw (gfp);

lammin/lammax =  0.9999999999999993 / 0.9999999999999993
Iteration 0 err= 2.2693827572253804
Iteration 1 err= 2.4521077551150774
Iteration 2 err= 2.556456083826408
Iteration 3 err= 1.8399291261764208
Iteration 4 err= 1.6397979567617955
Iteration 5 err= 1.1227484703156794
Iteration 6 err= 1.571878152242305
Iteration 7 err= 1.578665958559203
Iteration 8 err= 1.3114977115398425
Iteration 9 err= 1.6181368587478409
Iteration 10 err= 1.2159255177579773
Iteration 11 err= 1.3063314760812088
Iteration 12 err= 0.7753694641455832
Iteration 13 err= 0.7037791513458759
Iteration 14 err= 0.3542094017002498
Iteration 15 err= 0.2582316265611039
Iteration 16 err= 0.10956896174344156
Iteration 17 err= 0.06822155257283025
Iteration 18 err= 0.03867743230879999
Iteration 19 err= 0.01574233336548143
Iteration 20 err= 0.012556324138378401
Iteration 21 err= 0.006268014060966591
Iteration 22 err= 0.00579803836835753
Iteration 23 err= 0.0027214089686471436
Iteration 24 err= 0.0022470117450723555
Iteration 25 err= 0.0012135653881944293
Iteration 26 err= 0.0006815596893261466
Iteration 27 err= 0.00042485165355397067
Iteration 28 err= 0.00021280037638335268
Iteration 29 err= 0.0001487969376679355
Iteration 30 err= 5.871898147569218e-05
Iteration 31 err= 4.6508768563981184e-05
Iteration 32 err= 1.9521726909671235e-05
Iteration 33 err= 1.2977358264383947e-05
Iteration 34 err= 4.66811773135793e-06
Iteration 35 err= 2.562726802531528e-06
Iteration 36 err= 1.518172758996038e-06
Iteration 37 err= 5.437213226838117e-07
Iteration 38 err= 4.0101979369744123e-07
Iteration 39 err= 1.7427297251387739e-07
Iteration 40 err= 1.0091029727105706e-07
Iteration 41 err= 4.781165681543078e-08
Iteration 42 err= 2.8420679358690194e-08
Iteration 43 err= 1.5742754062254667e-08

for convenience of the reader, we started from this version of the usual preconditioned cg-iteration:

def CG(A,b,pre, maxit=200, tol=1e-8):
    x = b.CreateVector()
    w = b.CreateVector()
    r = b.CreateVector()
    p = b.CreateVector()
    ap = b.CreateVector()

    r.data = b
    p.data = pre*r
    wrn = InnerProduct(r,p)
    err0 = sqrt(wrn)
    x[:] = 0
    
    for it in range(maxit):
        ap.data = A * p

        pap = InnerProduct(p, ap)
        wr = wrn
        alpha = wr / pap
        
        x.data += alpha * p
        r.data -= alpha * ap
        w.data = pre*r
        
        wrn = InnerProduct(w, r)
        err = sqrt(wrn)
        print ("Iteration",it,"err=",err)
        if err < tol * err0: break
            
        beta = wrn / wr
        
        p *= beta
        p.data += w 
        
    return x