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 =  42885 + 14295 = 57180

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.5351907179280806 / 0.9991090979878732
Iteration 0 err= 0.6436410377222933     
Iteration 1 err= 0.7798068324142924     
Iteration 2 err= 0.5555064772141202     
Iteration 3 err= 0.49641340527791505     
Iteration 4 err= 0.5750625624154669     
Iteration 5 err= 0.6429009382013575     
Iteration 6 err= 0.7697892376140483     
Iteration 7 err= 0.7468571547958558     
Iteration 8 err= 0.7122867512112588     
Iteration 9 err= 0.6623820069786397     
Iteration 10 err= 0.6407358970496791     
Iteration 11 err= 0.6931166001430553     
Iteration 12 err= 0.6158906452421831     

Iteration 13 err= 0.5091804322895684     
Iteration 14 err= 0.46805838639251396     
Iteration 15 err= 0.4694473517289632     
Iteration 16 err= 0.44070796770231957     
Iteration 17 err= 0.447726682447315     
Iteration 18 err= 0.4998365934711988     
Iteration 19 err= 0.537582560329931     
Iteration 20 err= 0.4852288855364436     
Iteration 21 err= 0.48139180779231205     
Iteration 22 err= 0.5721489190674225     
Iteration 23 err= 0.605660910933646     
Iteration 24 err= 0.5854217137271133     
Iteration 25 err= 0.5494180805538459     
Iteration 26 err= 0.4717884743361036     
Iteration 27 err= 0.47876614737461604     
Iteration 28 err= 0.5115677106735579     
Iteration 29 err= 0.5018653937270958     
Iteration 30 err= 0.4892571677313947     
Iteration 31 err= 0.49328861656413536     
Iteration 32 err= 0.541402347784054     
Iteration 33 err= 0.5698205322280904     
Iteration 34 err= 0.5583008183997115     
Iteration 35 err= 0.5276053690062187     
Iteration 36 err= 0.5093415369119706     
Iteration 37 err= 0.511111850290832     
Iteration 38 err= 0.47890980891721874     
Iteration 39 err= 0.47385558730919713     
Iteration 40 err= 0.4982068560393076     
Iteration 41 err= 0.47184759741850163     
Iteration 42 err= 0.4734552791329646     
Iteration 43 err= 0.46244517270533586     
Iteration 44 err= 0.4401447202156156     
Iteration 45 err= 0.39661523008733546     
Iteration 46 err= 0.40627936640106477     
Iteration 47 err= 0.42270823305084637     
Iteration 48 err= 0.3800606247083981     
Iteration 49 err= 0.3230579778971618     
Iteration 50 err= 0.2685983052034653     
Iteration 51 err= 0.24608271354431752     
Iteration 52 err= 0.2275008646683409     
Iteration 53 err= 0.21305203828655556     
Iteration 54 err= 0.18826909910976689     
Iteration 55 err= 0.1574403189034756     
Iteration 56 err= 0.131550247632397     

Iteration 57 err= 0.11363569956955362     
Iteration 58 err= 0.09501356347090993     
Iteration 59 err= 0.07701505821027804     
Iteration 60 err= 0.06171535126289979     
Iteration 61 err= 0.04895455644348039     
Iteration 62 err= 0.03855420586578525     
Iteration 63 err= 0.03433761925971192     
Iteration 64 err= 0.031287141945907315     
Iteration 65 err= 0.025136313109144295     
Iteration 66 err= 0.019138996028011966     
Iteration 67 err= 0.01706969320560678     
Iteration 68 err= 0.01335176192309377     
Iteration 69 err= 0.009951851185014305     
Iteration 70 err= 0.008269181261284158     

Iteration 71 err= 0.007123744957512292     
Iteration 72 err= 0.005512547236650239     
Iteration 73 err= 0.00403363254669762     
Iteration 74 err= 0.003052604737458166     
Iteration 75 err= 0.0023546215756527567     
Iteration 76 err= 0.001910962952258433     
Iteration 77 err= 0.00160593968330534     
Iteration 78 err= 0.0013219197823532323     
Iteration 79 err= 0.0009853199809872698     
Iteration 80 err= 0.0007519050616378379     
Iteration 81 err= 0.0005717987825239775     
Iteration 82 err= 0.00047030129086439546     
Iteration 83 err= 0.0004190629862700898     
Iteration 84 err= 0.0003466379679389567     
Iteration 85 err= 0.00026374703981760683     
Iteration 86 err= 0.00020448418540804043     
Iteration 87 err= 0.00015152910757983864     
Iteration 88 err= 0.00012324283630135637     
Iteration 89 err= 0.00010763971287415228     
Iteration 90 err= 9.685603105124634e-05     
Iteration 91 err= 8.649182470197881e-05     
Iteration 92 err= 7.274486316168028e-05     
Iteration 93 err= 5.746222029499832e-05     
Iteration 94 err= 5.192246254730332e-05     
Iteration 95 err= 5.1117461698994634e-05     
Iteration 96 err= 4.5019959544900016e-05     
Iteration 97 err= 3.543717305689242e-05     
Iteration 98 err= 2.881591938873914e-05     
Iteration 99 err= 2.268292518770427e-05     
Iteration 100 err= 1.7842879089085803e-05     
Iteration 101 err= 1.4303409937383303e-05     
Iteration 102 err= 1.1333891123426938e-05     
Iteration 103 err= 9.061800736747913e-06     
Iteration 104 err= 6.7984238434970315e-06     
Iteration 105 err= 5.127415011170147e-06     
Iteration 106 err= 3.806073602314382e-06     
Iteration 107 err= 2.965069997955484e-06     
Iteration 108 err= 2.336622347149322e-06     
Iteration 109 err= 1.8212418866229577e-06     
Iteration 110 err= 1.4148879700089398e-06     
Iteration 111 err= 1.0872066763062005e-06     
Iteration 112 err= 8.725857078140307e-07     
Iteration 113 err= 6.88034921532866e-07     
Iteration 114 err= 4.879617533871071e-07     

Iteration 115 err= 3.7134171894078226e-07     
Iteration 116 err= 2.9535367255708207e-07     
Iteration 117 err= 2.2627233486001648e-07     
Iteration 118 err= 1.7052810543646857e-07     
Iteration 119 err= 1.2560669246782775e-07     
Iteration 120 err= 9.672661383016375e-08     
Iteration 121 err= 7.893866872024923e-08     
Iteration 122 err= 6.595710562085434e-08     
Iteration 123 err= 5.342697926024504e-08     
Iteration 124 err= 4.4027163610840136e-08     
Iteration 125 err= 3.754347320962985e-08     
Iteration 126 err= 3.233914952443855e-08     
Iteration 127 err= 2.881936860927959e-08     
Iteration 128 err= 2.4131522693192584e-08     
Iteration 129 err= 1.958579643410052e-08     

Iteration 130 err= 1.7088207598677408e-08     
Iteration 131 err= 1.520178769399592e-08     
Iteration 132 err= 1.203616984682676e-08     
Iteration 133 err= 8.718233009633121e-09     
Iteration 134 err= 6.464133001404661e-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
)
Contents