3. Schäfer-Turek Flow benchmark

https://wwwold.mathematik.tu-dortmund.de/lsiii/cms/papers/SchaeferTurek1996.pdf

Using MCS mixed methods

Thesis P. Lederer, Lederer-Gopalakrishnan-Schöberl, Lederer-Gopalakrishnan-Kogler-Schöberl,

Velocity \(u \in BDM_k \subset H(\operatorname{div})\), pressure \(p \in P^{k-1}\), stress \(\sigma\) such that

\[ \operatorname{div} \sigma \in H(\operatorname{div})^\ast \]

Thus \(\left< \operatorname{div} \sigma , u \right>\) is well defined.

Then \(\sigma \in H(\operatorname{curl} \operatorname{div})\), what means \(\sigma_{nt}\) is continuous.

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from time import time
import ipywidgets as widgets

order = 2

box = Box((0,0,0), (2.5,0.41,0.41))
box.faces.name="wall"
box.faces.Min(X).name="inlet"
box.faces.Max(X).name="outlet"
cyl = Cylinder((0.5,0.2,0), Z, h=0.41,r=0.05)
cyl.faces.name="cyl"
cyl.faces.maxh=0.03
shape = box-cyl

mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.1)).Curve(order)
print ("ne =", mesh.GetNE(VOL))
Draw (mesh);

ne = 5401

from NavierStokesIPCS_MCS import NavierStokes

# original benchmark:
um=2.25
timestep = 2e-3/um

# Reynolds 1000
# um = 20   
# timestep = 1e-3/um

uin = (um*16*y*(0.41-y)*z*(0.41-z)/0.41**4,0, 0)

with TaskManager():
    solver = NavierStokes(mesh, nu=1e-3, inflow="inlet", outflow="outlet", wall="wall|cyl",
                          uin=uin, timestep=timestep, order=order, verbose=1)                    

ndof X = 201672  =  132036 + 69636

H1AMG: level 0, num_edges = 7552, nv = 1347

Xproj.ndof= 253270

with TaskManager(): 
    solver.SolveInitial()

clipping = { "function" : True,  "pnt" : (2.5,0.2,0.2), "vec" : (0,0,-1) }

scene = Draw (solver.velocity, clipping=clipping, order=order, min=0, max=um);
scenep = Draw (solver.pressure, mesh, clipping=clipping, order=order, draw_surf=False);

tw = widgets.Text(value='t = 0')
perstep = widgets.Text(value="per step: 0")
display(tw, perstep)

tend = 0.1
t = 0
ts = time()
steps = 0
with TaskManager(): # pajetrace=10**9):
    while t < tend:
        t += timestep
        steps+=1
        solver.DoTimeStep()
        tw.value = "t = " + str(t)
        perstep.value = "per step: " + str( (time()-ts)/ steps)
        scene.Redraw()
        scenep.Redraw()