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
Let \(\widehat A\) be a preconditioner to \(A\) scaled such that
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
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):
finally, we multiply the first row by the positive definite matrix \(A-\widehat A\) to obtain
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
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
For the \(u\) term term we estimate
and for the \(p\)-term we already have
Now we determine the optimal \(\gamma\) from the equation
which leads to \(\gamma = \frac{1}{2} + \sqrt{ \frac{1}{4} + \frac{1}{\alpha-1} }\) and
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