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.5351907179280804 / 0.9991090979878732
Iteration 0 err= 0.6436410377222936     
Iteration 1 err= 0.7798068324142925     
Iteration 2 err= 0.5555064772141202     
Iteration 3 err= 0.49641340527791505     
Iteration 4 err= 0.5750625624154668     
Iteration 5 err= 0.6429009382013573     
Iteration 6 err= 0.769789237614048     
Iteration 7 err= 0.7468571547958556     
Iteration 8 err= 0.7122867512112587     
Iteration 9 err= 0.6623820069786395     
Iteration 10 err= 0.640735897049679     
Iteration 11 err= 0.693116600143055     
Iteration 12 err= 0.6158906452421823     
Iteration 13 err= 0.5091804322895681     
Iteration 14 err= 0.46805838639251396     
Iteration 15 err= 0.46944735172896307     
Iteration 16 err= 0.44070796770231935     
Iteration 17 err= 0.44772668244731506     

Iteration 18 err= 0.4998365934711986     
Iteration 19 err= 0.5375825603299305     
Iteration 20 err= 0.48522888553644306     
Iteration 21 err= 0.4813918077923119     
Iteration 22 err= 0.5721489190674222     
Iteration 23 err= 0.6056609109336452     
Iteration 24 err= 0.5854217137271122     
Iteration 25 err= 0.5494180805538457     
Iteration 26 err= 0.47178847433610355     
Iteration 27 err= 0.47876614737461587     
Iteration 28 err= 0.5115677106735574     
Iteration 29 err= 0.5018653937270953     
Iteration 30 err= 0.48925716773139405     
Iteration 31 err= 0.49328861656413475     
Iteration 32 err= 0.5414023477840532     
Iteration 33 err= 0.5698205322280894     
Iteration 34 err= 0.5583008183997115     
Iteration 35 err= 0.5276053690062186     
Iteration 36 err= 0.5093415369119704     
Iteration 37 err= 0.5111118502908312     
Iteration 38 err= 0.47890980891721796     
Iteration 39 err= 0.47385558730919675     
Iteration 40 err= 0.49820685603930714     
Iteration 41 err= 0.47184759741850063     
Iteration 42 err= 0.4734552791329634     
Iteration 43 err= 0.46244517270533386     
Iteration 44 err= 0.4401447202156105     
Iteration 45 err= 0.3966152300873285     
Iteration 46 err= 0.4062793664010494     
Iteration 47 err= 0.4227082330508433     
Iteration 48 err= 0.38006062470840507     
Iteration 49 err= 0.3230579778971632     
Iteration 50 err= 0.26859830520346545     
Iteration 51 err= 0.24608271354431807     
Iteration 52 err= 0.2275008646683361     
Iteration 53 err= 0.21305203828649594     
Iteration 54 err= 0.18826909910932568     
Iteration 55 err= 0.15744031890165874     
Iteration 56 err= 0.13155024762622888     
Iteration 57 err= 0.11363569954298637     
Iteration 58 err= 0.0950135632915174     
Iteration 59 err= 0.07701505679258615     
Iteration 60 err= 0.061715341975890646     
Iteration 61 err= 0.04895450504222048     
Iteration 62 err= 0.03855394440440356     
Iteration 63 err= 0.034335762121178684     
Iteration 64 err= 0.03127797364435805     

Iteration 65 err= 0.025094792684768726     
Iteration 66 err= 0.018898362872886387     
Iteration 67 err= 0.01665843502011037     
Iteration 68 err= 0.013482613849355406     
Iteration 69 err= 0.010000452880549284     
Iteration 70 err= 0.008278026513663442     
Iteration 71 err= 0.007124998307657447     
Iteration 72 err= 0.005512707703550439     
Iteration 73 err= 0.004033657573714523     
Iteration 74 err= 0.003052610711025041     
Iteration 75 err= 0.002354629561659234     
Iteration 76 err= 0.0019109939798426945     
Iteration 77 err= 0.0016060602990358973     
Iteration 78 err= 0.0013224488931256992     
Iteration 79 err= 0.0009880142862928824     
Iteration 80 err= 0.0007636492609903181     
Iteration 81 err= 0.0006061362876643011     
Iteration 82 err= 0.0005209707988575197     
Iteration 83 err= 0.00041345697521742914     
Iteration 84 err= 0.00033258071386250255     

Iteration 85 err= 0.00026166187356345443     
Iteration 86 err= 0.00020429279565955152     
Iteration 87 err= 0.00015171687950284852     
Iteration 88 err= 0.00012383596129459672     
Iteration 89 err= 0.00010918952836765955     
Iteration 90 err= 9.930285789021855e-05     
Iteration 91 err= 8.894241686032606e-05     
Iteration 92 err= 7.765262912485472e-05     
Iteration 93 err= 6.311079386252812e-05     
Iteration 94 err= 5.413026605870057e-05     
Iteration 95 err= 4.958622225418996e-05     
Iteration 96 err= 4.319467850313892e-05     
Iteration 97 err= 3.503132785517901e-05     
Iteration 98 err= 2.8641592784745954e-05     
Iteration 99 err= 2.2619190243663354e-05     
Iteration 100 err= 1.7826471067589977e-05     
Iteration 101 err= 1.4300419360771886e-05     
Iteration 102 err= 1.1333342114528846e-05     
Iteration 103 err= 9.061705463222328e-06     
Iteration 104 err= 6.798376134010088e-06     
Iteration 105 err= 5.127402845896933e-06     
Iteration 106 err= 3.806064207533803e-06     
Iteration 107 err= 2.96502618785763e-06     
Iteration 108 err= 2.3363697530492775e-06     
Iteration 109 err= 1.819863724581574e-06     
Iteration 110 err= 1.4096183395184093e-06     
Iteration 111 err= 1.0615475984891518e-06     
Iteration 112 err= 7.707596515618936e-07     
Iteration 113 err= 5.678939623099133e-07     
Iteration 114 err= 4.772703684176031e-07     
Iteration 115 err= 3.998963965934087e-07     
Iteration 116 err= 3.0088163131911137e-07     
Iteration 117 err= 2.2692166549564326e-07     
Iteration 118 err= 1.7062856284833003e-07     
Iteration 119 err= 1.2562579190728693e-07     
Iteration 120 err= 9.672949757657273e-08     
Iteration 121 err= 7.894075687284693e-08     
Iteration 122 err= 6.596598669513807e-08     
Iteration 123 err= 5.345575685599375e-08     
Iteration 124 err= 4.413072846655022e-08     
Iteration 125 err= 3.802210370225776e-08     
Iteration 126 err= 3.478383785614803e-08     
Iteration 127 err= 3.189332410423553e-08     
Iteration 128 err= 2.5200471200452652e-08     
Iteration 129 err= 2.150388496342634e-08     
Iteration 130 err= 1.7847365240445305e-08     
Iteration 131 err= 1.4473137259214665e-08     

Iteration 132 err= 1.158533469490621e-08     
Iteration 133 err= 8.655748526369781e-09     
Iteration 134 err= 6.4593823524986765e-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