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{split} \left( \begin{array}{ccc} A & B^T \\ B & 0 \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} f \\ g \end{array} \right). \end{split}\]

Let \(\widehat A\) be a preconditioner to \(A\) scaled such that

\[ \widehat A < A \leq \alpha \widehat A \]

Thanks to the scaling, the matrix \(A - \widehat A\) is positive definite. The spectrum of the preconditioned system satisfies \(\sigma(\widehat A^{-1} A) \subset (1, \alpha]\).

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

\[\begin{split} \left( \begin{array}{ccc} \widehat A^{-1} A & \widehat A^{-1} B^T \\ B & 0 \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} \widehat A^{-1} f \\ g \end{array} \right). \end{split}\]

If \(\widehat 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{split} \left( \begin{array}{ccc} \widehat A^{-1} A & \widehat A^{-1} B^T \\ B \widehat A^{-1} A - B & B \widehat A^{-1} B^T \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} \widehat A^{-1} f \\ B \widehat A^{-1} f - g \end{array} \right). \end{split}\]

finally, we multiply the first row by the positive definite matrix \(A-\widehat A\) to obtain

\[\begin{split} \left( \begin{array}{ccc} (A-\widehat A) \widehat A^{-1} A & (A-\widehat A) \hat A^{-1} B^T \\ B \widehat A^{-1} (A - \hat A) & B \widehat A^{-1} B^T \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} (A-\widehat A) \widehat A^{-1} f \\ B \widehat A^{-1} f -g \end{array} \right) \end{split}\]

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{split} \left( \begin{array}{ccc} A-\widehat A & 0 \\ 0 & I \end{array} \right) \left( \begin{array}{ccc} \widehat I & 0 \\ B & -I \end{array} \right) \left( \begin{array}{ccc} \widehat A^{-1} & 0 \\ 0 & I \end{array} \right) \end{split}\]

Theorem

Let

\[\begin{split} {\mathcal K} = \left( \begin{array}{cc} (A-\widehat A) \widehat A^{-1} A & (A-\widehat A) \hat A^{-1} B^T \\ B \widehat A^{-1} (A - \hat A) & B \widehat A^{-1} B^T \end{array} \right) \end{split}\]

and

\[\begin{split} {\mathcal C} = \left( \begin{array}{cc} A - \widehat A & \\ 0 & B \widehat A^{-1} B^T \end{array} \right) \end{split}\]

Then there hold the spectral estimates

\[ \gamma_1 {\mathcal C} \leq {\mathcal K} \leq \gamma_2 {\mathcal C} \]

with some \(\gamma_2 < 2\) and \(\gamma_1 = O(\alpha^{-1})\).

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

\[\begin{align*} \left( \begin{array}{c} u \\ p \end{array} \right) {\mathcal K} \left( \begin{array}{c} u \\ p \end{array} \right) & = u^T (A-\widehat A) \widehat A^{-1} A u + 2 p^T B \widehat A^{-1} (A - \widehat A) u + p^T B \widehat A^{-1} B^T p \\ &\leq u^T (A-\widehat A) \widehat A^{-1} A u + \gamma u^T (A-\widehat A) \widehat A^{-1} (A - \widehat A) u \\ & + (1 + \gamma^{-1}) p^T B \widehat A^{-1} B^T p \end{align*}\]

For the \(u\) term term we estimate

\[ u^T (A - \widehat A) \widehat A^{-1} (A + \gamma (A-\widehat A)) u \leq (\alpha + \gamma \alpha - \gamma) u^T (A - \widehat A) u \]

and for the \(p\)-term we already have

\[ (1 + \gamma^{-1}) p^T B \widehat A^{-1} B^T p \leq (1 + \gamma^{-1}) p^T B \widehat A^{-1} B^T p \]

Now we determine the optimal \(\gamma\) from the equation

\[ \alpha + \gamma \alpha - \gamma = 1 + \gamma^{-1} \]

which leads to \(\gamma = \frac{1}{2} + \sqrt{ \frac{1}{4} + \frac{1}{\alpha-1} }\) and

\[ {\mathcal K} \leq 1 + \left( \frac{1}{2} + \sqrt{ \frac{1}{4} + \frac{1}{\alpha-1} } \right)^{-1} {\mathcal C} \leq 2 {\mathcal C} \]

The lower estimate is proven with a similarly (exercise). \(\Box\)

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