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.9999999999999993 / 0.9999999999999993
Iteration 0 err= 2.27523148567202
Iteration 1 err= 2.459534634829986
Iteration 2 err= 2.5609251911988267
Iteration 3 err= 1.828640083749773
Iteration 4 err= 1.6430809405007087
Iteration 5 err= 1.1384217332188629
Iteration 6 err= 1.5534199211755715
Iteration 7 err= 1.5442590512095775
Iteration 8 err= 1.3360563648468748
Iteration 9 err= 1.6418731406042861
Iteration 10 err= 1.187830422193117
Iteration 11 err= 1.3135666330705376
Iteration 12 err= 0.7999426962809655
Iteration 13 err= 0.7081697695522868
Iteration 14 err= 0.35466296023698773
Iteration 15 err= 0.2543844210629828
Iteration 16 err= 0.11034470640104128
Iteration 17 err= 0.07206149199653164
Iteration 18 err= 0.03817866085104304
Iteration 19 err= 0.015478243833255145
Iteration 20 err= 0.012593235496691717
Iteration 21 err= 0.006569583458348432
Iteration 22 err= 0.006232827417990785
Iteration 23 err= 0.0027951959113931696
Iteration 24 err= 0.0022768585491949834
Iteration 25 err= 0.0012372842216013314
Iteration 26 err= 0.0006433262689864749
Iteration 27 err= 0.00039986779835608397
Iteration 28 err= 0.00017058117698389257
Iteration 29 err= 0.00011106809555467636
Iteration 30 err= 4.0338871445907486e-05
Iteration 31 err= 3.0375602757739977e-05
Iteration 32 err= 1.3750164072902666e-05
Iteration 33 err= 9.285383502534869e-06
Iteration 34 err= 4.173505993943293e-06
Iteration 35 err= 2.3451585476408004e-06
Iteration 36 err= 1.4667974311573793e-06
Iteration 37 err= 5.157853296439992e-07
Iteration 38 err= 3.9420255790013686e-07
Iteration 39 err= 1.6735870298436908e-07
Iteration 40 err= 8.788546016467882e-08
Iteration 41 err= 4.191758161776519e-08
Iteration 42 err= 2.2520155287756295e-08
Iteration 43 err= 1.422405219316397e-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