93. 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.9999999999999996 / 0.9999999999999996
Iteration 0 err= 2.275231485672021
Iteration 1 err= 2.4595346348299874
Iteration 2 err= 2.5609251911988244
Iteration 3 err= 1.8286400837497736
Iteration 4 err= 1.6430809405007067
Iteration 5 err= 1.1384217332188646
Iteration 6 err= 1.5534199211755717
Iteration 7 err= 1.5442590512095804
Iteration 8 err= 1.3360563648468742
Iteration 9 err= 1.641873140604289
Iteration 10 err= 1.1878304221931175
Iteration 11 err= 1.3135666330705387
Iteration 12 err= 0.799942696280966
Iteration 13 err= 0.7081697695522882
Iteration 14 err= 0.35466296023698746
Iteration 15 err= 0.2543844210629822
Iteration 16 err= 0.1103447064010411
Iteration 17 err= 0.07206149199653222
Iteration 18 err= 0.03817866085104358
Iteration 19 err= 0.015478243833255394
Iteration 20 err= 0.012593235496691812
Iteration 21 err= 0.006569583458348788
Iteration 22 err= 0.006232827417991425
Iteration 23 err= 0.0027951959113932845
Iteration 24 err= 0.002276858549194963
Iteration 25 err= 0.0012372842216012837
Iteration 26 err= 0.0006433262689863725
Iteration 27 err= 0.00039986779835627306
Iteration 28 err= 0.00017058117698420488
Iteration 29 err= 0.00011106809555478819
Iteration 30 err= 4.03388714456461e-05
Iteration 31 err= 3.037560275714898e-05
Iteration 32 err= 1.3750164072862934e-05
Iteration 33 err= 9.285383502906159e-06
Iteration 34 err= 4.17350599374367e-06
Iteration 35 err= 2.345158547203137e-06
Iteration 36 err= 1.4667974311645842e-06
Iteration 37 err= 5.157853297005392e-07
Iteration 38 err= 3.9420255804794306e-07
Iteration 39 err= 1.673587031423844e-07
Iteration 40 err= 8.788545999950837e-08
Iteration 41 err= 4.191758153920234e-08
Iteration 42 err= 2.2520155196774335e-08
Iteration 43 err= 1.4224051803555327e-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