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.4595346348299865
Iteration 2 err= 2.560925191198827
Iteration 3 err= 1.8286400837497732
Iteration 4 err= 1.6430809405007099
Iteration 5 err= 1.1384217332188638
Iteration 6 err= 1.5534199211755728
Iteration 7 err= 1.5442590512095782
Iteration 8 err= 1.3360563648468766
Iteration 9 err= 1.641873140604289
Iteration 10 err= 1.1878304221931186
Iteration 11 err= 1.313566633070539
Iteration 12 err= 0.7999426962809673
Iteration 13 err= 0.7081697695522884
Iteration 14 err= 0.35466296023698846
Iteration 15 err= 0.2543844210629829
Iteration 16 err= 0.1103447064010416
Iteration 17 err= 0.0720614919965319
Iteration 18 err= 0.038178660851043135
Iteration 19 err= 0.015478243833255205
Iteration 20 err= 0.012593235496691791
Iteration 21 err= 0.00656958345834857
Iteration 22 err= 0.006232827417991075
Iteration 23 err= 0.002795195911393293
Iteration 24 err= 0.002276858549195023
Iteration 25 err= 0.001237284221601205
Iteration 26 err= 0.0006433262689863229
Iteration 27 err= 0.00039986779835634
Iteration 28 err= 0.00017058117698422878
Iteration 29 err= 0.00011106809555473569
Iteration 30 err= 4.033887144580721e-05
Iteration 31 err= 3.0375602757399626e-05
Iteration 32 err= 1.3750164072671311e-05
Iteration 33 err= 9.285383502670899e-06
Iteration 34 err= 4.173505994090854e-06
Iteration 35 err= 2.3451585474356334e-06
Iteration 36 err= 1.4667974309083101e-06
Iteration 37 err= 5.157853295269396e-07
Iteration 38 err= 3.9420255786777766e-07
Iteration 39 err= 1.6735870303018525e-07
Iteration 40 err= 8.788546025741023e-08
Iteration 41 err= 4.1917581918502326e-08
Iteration 42 err= 2.252015537650467e-08
Iteration 43 err= 1.4224051767640286e-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