Contents

Example solved by the Bramble-Pasciak Transformation

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

dim = 3

if dim==2:
    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)
else:
    box = Box((0,0,0), (2,0.41,0.41))
    box.faces.name="wall"
    box.faces.Min(X).name="inlet"
    box.faces.Max(X).name="outlet"
    cyl = Cylinder((0.2,0,0.2), Y, h=0.41,r=0.05)
    cyl.faces.name="cyl"
    cyl.faces.maxh=0.02
    shape = box-cyl
    Draw (shape)
    
    mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.03))
    Draw(mesh)
V = VectorH1(mesh, order=1, dirichlet="wall|inlet|cyl")
V1 = H1(mesh, order=1, dirichlet="wall|inlet|cyl")
Q = H1(mesh, order=1)
print ("ndof = ", V.ndof,'+',Q.ndof,'=',V.ndof+Q.ndof)

u,v = V.TnT()
u1,v1 = V1.TnT()
p,q = Q.TnT()

h = specialcf.mesh_size

with TaskManager():
    bfa1 = BilinearForm(InnerProduct(grad(u1),grad(v1))*dx)
    bfb = BilinearForm(div(u)*q*dx).Assemble()
    bfc = BilinearForm(h*h*grad(p)*grad(q)*dx).Assemble()

    # prea1 = Preconditioner(bfa1, "direct", inverse="sparsecholesky")
    prea1 = Preconditioner(bfa1, "h1amg")
    bfa1.Assemble()

mata = sum( [Ri.T@bfa1.mat@Ri for Ri in V.restrictions] )
prea = sum( [Ei@prea1@Ei.T for Ei in V.embeddings])    
    
with TaskManager():
    bfschur = BilinearForm(p*q*dx, diagonal=True).Assemble()
    preschur = bfschur.mat.Inverse()
ndof =  43101 + 14367 = 57468
from ngsolve.krylovspace import BramblePasciakCG

gfu = GridFunction(V)
gfp = GridFunction(Q)

if mesh.dim==2:
    uin = (1.5*4*y*(0.41-y)/(0.41*0.41), 0)
else:
    uin = (1.5*4*y*(0.41-y)/(0.41*0.41)*z*(0.41-z)/0.41**2,0, 0)

gfu.Set(uin, definedon=mesh.Boundaries("inlet"))

resf = (-mata * gfu.vec).Evaluate()
resg = (-bfb.mat * gfu.vec).Evaluate()

with TaskManager():
    sol = BramblePasciakCG (A=mata, B=bfb.mat, C=bfc.mat, f=resf, g=resg, \
                            preA=prea, preS=preschur, printrates='\r')

gfu.vec.data += sol[0]
gfp.vec.data += sol[1]
lammin/lammax =  0.5358648222961435 / 0.999105875654954
Iteration 0 err= 0.6422574868811712     
Iteration 1 err= 0.7626795389173298     
Iteration 2 err= 0.5401318134887052     
Iteration 3 err= 0.4892925374272816     
Iteration 4 err= 0.5491577368906568     
Iteration 5 err= 0.7166280513160017     
Iteration 6 err= 0.7917614043413376     
Iteration 7 err= 0.7672069371787942     
Iteration 8 err= 0.7287649730899978     
Iteration 9 err= 0.6792916466594386     
Iteration 10 err= 0.6612431188635448     
Iteration 11 err= 0.7175647832785991     
Iteration 12 err= 0.6508451302837635     
Iteration 13 err= 0.5203140321780023     
Iteration 14 err= 0.506146812166123     
Iteration 15 err= 0.4814853770420756     
Iteration 16 err= 0.4628668141158896     
Iteration 17 err= 0.44579395591044385     
Iteration 18 err= 0.4757002448871478     
Iteration 19 err= 0.5355291443173381     
Iteration 20 err= 0.5025194394231776     
Iteration 21 err= 0.45942323241555777     
Iteration 22 err= 0.525126310870192     
Iteration 23 err= 0.5937316187541024     
Iteration 24 err= 0.562975421337546     
Iteration 25 err= 0.500161551012948     
Iteration 26 err= 0.46733689671044293     
Iteration 27 err= 0.4883456689663639     
Iteration 28 err= 0.5017360038189036     
Iteration 29 err= 0.5107434241281112     
Iteration 30 err= 0.49754848815741537     
Iteration 31 err= 0.5061729904018872     
Iteration 32 err= 0.5036485771930268     
Iteration 33 err= 0.5205385000271844     
Iteration 34 err= 0.5178148357977969     
Iteration 35 err= 0.5253742720462282     
Iteration 36 err= 0.5175925277867022     
Iteration 37 err= 0.4947798856049582     
Iteration 38 err= 0.46524770740059573     
Iteration 39 err= 0.4954387686690926     
Iteration 40 err= 0.5073268053984671     
Iteration 41 err= 0.48361914822908236     
Iteration 42 err= 0.4695364680634301     
Iteration 43 err= 0.4734783166384601     
Iteration 44 err= 0.44194258384137874     
Iteration 45 err= 0.3943624005154073     
Iteration 46 err= 0.4014077862310191     
Iteration 47 err= 0.40869482932348294     
Iteration 48 err= 0.36691584808979344     
Iteration 49 err= 0.306458230175725     
Iteration 50 err= 0.26349648762311184     
Iteration 51 err= 0.24345692145641926     
Iteration 52 err= 0.23383221647466643     
Iteration 53 err= 0.22096721750816906     
Iteration 54 err= 0.18164767317754468     
Iteration 55 err= 0.14633377404439676     
Iteration 56 err= 0.12721915379325038     
Iteration 57 err= 0.11432269168627315     
Iteration 58 err= 0.09896024590430351     
Iteration 59 err= 0.08856644126791419     
Iteration 60 err= 0.06794013863864694     
Iteration 61 err= 0.05151361257467961     
Iteration 62 err= 0.04169918663782892     
Iteration 63 err= 0.03581022491124209     
Iteration 64 err= 0.033092458683837936     
Iteration 65 err= 0.027330382914969233     
Iteration 66 err= 0.02190029435793318     
Iteration 67 err= 0.01810148105552501     
Iteration 68 err= 0.013687196333319095     
Iteration 69 err= 0.010935017223374172     
Iteration 70 err= 0.00896326148211464     
Iteration 71 err= 0.007386609067897706     
Iteration 72 err= 0.005865492231097758     
Iteration 73 err= 0.004362376462786854     
Iteration 74 err= 0.0031658432800965625     
Iteration 75 err= 0.0025620873641519076     
Iteration 76 err= 0.0020542294451275636     
Iteration 77 err= 0.0016652265452639387     
Iteration 78 err= 0.0013787070789876878     
Iteration 79 err= 0.0010534824949476145     
Iteration 80 err= 0.000778573187902241     
Iteration 81 err= 0.0006150757220387196     
Iteration 82 err= 0.0004994633776232117     
Iteration 83 err= 0.000392910803733831     
Iteration 84 err= 0.0003221989357829235     
Iteration 85 err= 0.0002561028796699069     
Iteration 86 err= 0.00019586587106754536     
Iteration 87 err= 0.00014753949588548073     
Iteration 88 err= 0.00012008137501909425     
Iteration 89 err= 0.00010508880471855321     
Iteration 90 err= 9.574788922448143e-05     
Iteration 91 err= 8.234688003079029e-05     
Iteration 92 err= 7.131205540565637e-05     
Iteration 93 err= 5.936126119878902e-05     
Iteration 94 err= 5.154452942584167e-05     
Iteration 95 err= 4.598740783393944e-05     
Iteration 96 err= 4.1805232351868455e-05     
Iteration 97 err= 3.6219441622724636e-05     
Iteration 98 err= 2.7675200674893923e-05     
Iteration 99 err= 2.1114130220237352e-05     
Iteration 100 err= 1.7152489425908266e-05     
Iteration 101 err= 1.403677861313346e-05     
Iteration 102 err= 1.1372572476463817e-05     
Iteration 103 err= 8.522085007612375e-06     
Iteration 104 err= 6.414533876461114e-06     
Iteration 105 err= 4.881612294618145e-06     
Iteration 106 err= 3.835846559330044e-06     
Iteration 107 err= 3.1778203806989535e-06     
Iteration 108 err= 2.624427968030782e-06     
Iteration 109 err= 2.009650710200033e-06     
Iteration 110 err= 1.5158049796646495e-06     
Iteration 111 err= 1.1228818020100263e-06     
Iteration 112 err= 8.204149766520061e-07     
Iteration 113 err= 5.991154259597335e-07     
Iteration 114 err= 4.6507274819498893e-07     
Iteration 115 err= 3.7257899444650243e-07     
Iteration 116 err= 2.9708934938237364e-07     
Iteration 117 err= 2.326153482571042e-07     
Iteration 118 err= 1.8926141872257927e-07     
Iteration 119 err= 1.4544420525127089e-07     
Iteration 120 err= 1.1009905789328344e-07     
Iteration 121 err= 9.129552785510817e-08     
Iteration 122 err= 7.968513132372621e-08     
Iteration 123 err= 6.77199207211961e-08     
Iteration 124 err= 5.590067442127978e-08     
Iteration 125 err= 4.763754794604962e-08     
Iteration 126 err= 4.1846458610721273e-08     
Iteration 127 err= 3.523365070590343e-08     
Iteration 128 err= 2.9378462200803565e-08     
Iteration 129 err= 2.412803916563848e-08     
Iteration 130 err= 1.8652483209788057e-08     
Iteration 131 err= 1.4712663597729102e-08     
Iteration 132 err= 1.179741413526433e-08     
Iteration 133 err= 8.935021255972776e-09     
Iteration 134 err= 6.775044661780795e-09     
Draw (Norm(gfu), mesh, order=1, clipping={"y":1,"z":0, "function":True})
Draw (gfp, order=1);
help (BramblePasciakCG)
Help on function BramblePasciakCG in module ngsolve.krylovspace:

BramblePasciakCG(A, B, C, f, g, preA, preS, maxit=1000, tol=1e-08, printrates=False)