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.7798068324142925     
Iteration 2 err= 0.5555064772141202     
Iteration 3 err= 0.49641340527791494     
Iteration 4 err= 0.5750625624154668     
Iteration 5 err= 0.6429009382013572     
Iteration 6 err= 0.7697892376140484     
Iteration 7 err= 0.746857154795856     
Iteration 8 err= 0.7122867512112594     
Iteration 9 err= 0.6623820069786396     
Iteration 10 err= 0.640735897049679     
Iteration 11 err= 0.6931166001430554     
Iteration 12 err= 0.6158906452421831     
Iteration 13 err= 0.5091804322895683     
Iteration 14 err= 0.46805838639251396     
Iteration 15 err= 0.4694473517289633     
Iteration 16 err= 0.4407079677023196     
Iteration 17 err= 0.44772668244731517     
Iteration 18 err= 0.49983659347119885     

Iteration 19 err= 0.5375825603299311     
Iteration 20 err= 0.4852288855364437     
Iteration 21 err= 0.4813918077923122     
Iteration 22 err= 0.5721489190674226     
Iteration 23 err= 0.6056609109336458     
Iteration 24 err= 0.5854217137271128     
Iteration 25 err= 0.5494180805538462     
Iteration 26 err= 0.47178847433610366     
Iteration 27 err= 0.4787661473746158     
Iteration 28 err= 0.5115677106735578     
Iteration 29 err= 0.5018653937270958     
Iteration 30 err= 0.4892571677313946     
Iteration 31 err= 0.49328861656413503     
Iteration 32 err= 0.5414023477840532     
Iteration 33 err= 0.5698205322280896     
Iteration 34 err= 0.5583008183997124     
Iteration 35 err= 0.5276053690062192     
Iteration 36 err= 0.5093415369119707     
Iteration 37 err= 0.5111118502908318     
Iteration 38 err= 0.47890980891721857     
Iteration 39 err= 0.4738555873091972     
Iteration 40 err= 0.4982068560393078     
Iteration 41 err= 0.4718475974185015     
Iteration 42 err= 0.4734552791329644     
Iteration 43 err= 0.4624451727053352     
Iteration 44 err= 0.44014472021561263     
Iteration 45 err= 0.39661523008733074     
Iteration 46 err= 0.4062793664010531     
Iteration 47 err= 0.4227082330508443     
Iteration 48 err= 0.3800606247084043     
Iteration 49 err= 0.32305797789716384     
Iteration 50 err= 0.2685983052034662     
Iteration 51 err= 0.24608271354431885     
Iteration 52 err= 0.22750086466833944     
Iteration 53 err= 0.213052038286529     
Iteration 54 err= 0.18826909910956494     
Iteration 55 err= 0.15744031890263702     
Iteration 56 err= 0.13155024762954337     
Iteration 57 err= 0.11363569955725918     
Iteration 58 err= 0.09501356338789099     
Iteration 59 err= 0.07701505755418783     
Iteration 60 err= 0.0617153469649411     
Iteration 61 err= 0.04895453265526377     
Iteration 62 err= 0.038554084863854594     
Iteration 63 err= 0.034336759864460414     
Iteration 64 err= 0.031282901910129314     
Iteration 65 err= 0.02511716708304179     

Iteration 66 err= 0.01902990907238441     
Iteration 67 err= 0.016905886922391825     
Iteration 68 err= 0.013410131700390916     
Iteration 69 err= 0.00997017444201796     
Iteration 70 err= 0.008272435142292846     
Iteration 71 err= 0.007124204184060756     
Iteration 72 err= 0.005512605964568883     
Iteration 73 err= 0.004033641508869306     
Iteration 74 err= 0.003052606137507243     
Iteration 75 err= 0.002354621255304202     
Iteration 76 err= 0.0019109602136377112     
Iteration 77 err= 0.001605927526242757     
Iteration 78 err= 0.0013218634821901898     
Iteration 79 err= 0.0009850239936529722     
Iteration 80 err= 0.0007505272614215356     
Iteration 81 err= 0.0005664988672406867     
Iteration 82 err= 0.0004481563778307417     
Iteration 83 err= 0.0003683591541516672     
Iteration 84 err= 0.00031613120772595193     
Iteration 85 err= 0.00025627491173273496     
Iteration 86 err= 0.00020877197015145796     

Iteration 87 err= 0.00015638674875803304     
Iteration 88 err= 0.00012527114671885965     
Iteration 89 err= 0.00011049005595070225     
Iteration 90 err= 0.00010037698805559284     
Iteration 91 err= 8.915009078973712e-05     
Iteration 92 err= 7.78921848459878e-05     
Iteration 93 err= 6.30060861463638e-05     
Iteration 94 err= 5.388132858834472e-05     
Iteration 95 err= 4.949917023037078e-05     
Iteration 96 err= 4.318092372257852e-05     
Iteration 97 err= 3.5029740396925735e-05     
Iteration 98 err= 2.8641141845073934e-05     
Iteration 99 err= 2.2619051654786746e-05     
Iteration 100 err= 1.782643701462638e-05     
Iteration 101 err= 1.4300413219922837e-05     
Iteration 102 err= 1.1333341008150001e-05     
Iteration 103 err= 9.06170537599474e-06     
Iteration 104 err= 6.798376489857956e-06     
Iteration 105 err= 5.127405290308138e-06     
Iteration 106 err= 3.806078332005516e-06     
Iteration 107 err= 2.965111878682875e-06     
Iteration 108 err= 2.336868613444752e-06     
Iteration 109 err= 1.8225810824308396e-06     
Iteration 110 err= 1.4198834918437845e-06     
Iteration 111 err= 1.1086102924423612e-06     
Iteration 112 err= 9.143867460876161e-07     
Iteration 113 err= 6.804460593213838e-07     
Iteration 114 err= 4.833967255992788e-07     
Iteration 115 err= 3.7079784087804903e-07     
Iteration 116 err= 2.953032958185205e-07     
Iteration 117 err= 2.262672390985073e-07     
Iteration 118 err= 1.70527380521108e-07     
Iteration 119 err= 1.2560656520965702e-07     
Iteration 120 err= 9.672660877897361e-08     
Iteration 121 err= 7.893869789625041e-08     
Iteration 122 err= 6.595715338358314e-08     
Iteration 123 err= 5.3426980327368506e-08     
Iteration 124 err= 4.402701727941045e-08     
Iteration 125 err= 3.7542608875622183e-08     
Iteration 126 err= 3.233309048020016e-08     
Iteration 127 err= 2.87934633779183e-08     
Iteration 128 err= 2.4072864224623353e-08     
Iteration 129 err= 1.9348208134853554e-08     
Iteration 130 err= 1.6086746829439893e-08     
Iteration 131 err= 1.3047828643779448e-08     
Iteration 132 err= 1.0470394666941598e-08     
Iteration 133 err= 8.917496024870079e-09     
Iteration 134 err= 6.751007573778622e-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