7. Navier-Stokes with MCS

Benchmark from Schäfer and Turek

MCS … Mass conserving Mixed Stress formulation, Phd thesis Philip Lederer

Idea: TDNNS upside down

vector variable \(u\): normal component continuous -> Raviart Thomas elements

matrix variable \(\sigma\): tangential component of normal vector continuous

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

order = 2
nu = 1e-3

rec = Rectangle(2,0.41).Face()
circ = Circle((0.2,0.2),0.05).Face()
circ.edges.maxh=0.01
shape = rec-circ
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
mesh = shape.GenerateMesh(maxh=0.05, dim=2).Curve(order)

print ("ne = ", mesh.GetNE(VOL))
Draw (mesh);

ne =  919

V = HDiv(mesh, order=order, dirichlet="inlet|wall", highest_order_dc=True)
Vhat = TangentialFacetFESpace(mesh, order=order-1, dirichlet="inlet|wall")
Sigma = Discontinuous(HCurlDiv(mesh, order = order-1, orderinner=order))
S = MatrixValued(L2(mesh, order=order-1), skewsymmetric=True)
Sigma = Compress(PrivateSpace(Sigma))
S = Compress(PrivateSpace(S))
Q = L2(mesh, order=order-1)
X = V*Vhat*Sigma*S*Q

u, uhat, sigma, W, p = X.TrialFunction()
v, vhat, tau, R, q = X.TestFunction()

dS = dx(element_boundary=True, bonus_intorder=2)
n = specialcf.normal(mesh.dim)
def tang(u): return u-(u*n)*n
        
stokesA = -0.5/nu * InnerProduct(sigma,tau) * dx + \
           (div(sigma)*v+div(tau)*u) * dx + \
           (InnerProduct(W,tau) + InnerProduct(R,sigma)) * dx + \
           -(((sigma*n)*n) * (v*n) + ((tau*n)*n )* (u*n)) * dS + \
           (-(sigma*n)*tang(vhat) - (tau*n)*tang(uhat)) * dS

stokes = stokesA + div(u)*q*dx + div(v)*p*dx - 1e-8*p*q*dx
astokes = BilinearForm (stokes, eliminate_hidden = True).Assemble()
f = LinearForm(X)

gfu = GridFunction(X)
umean = 1.5

uin = CF( (umean*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
gfu.components[1].Set(uin, definedon=mesh.Boundaries("inlet"))

inv_stokes = astokes.mat.Inverse(X.FreeDofs())

res = f.vec - astokes.mat*gfu.vec
gfu.vec.data += inv_stokes * res

Draw (gfu.components[0], mesh);

tau = 0.001 # timestep

mstar = BilinearForm(u*v*dx+tau*stokes, eliminate_hidden = True).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

# conv = BilinearForm(X, nonassemble=True)
# conv += (Grad(u) * u) * v * dx


conv = BilinearForm(X, nonlinear_matrix_free_bdb=True)
conv += -InnerProduct(grad(v)*u, u) * dx(bonus_intorder=order)

un = u.Operator("normalcomponent")
uhat_t = uhat.Operator("tangentialcomponent")
vn = v.Operator("normalcomponent")
conv += IfPos(un, un * u * v, 2 * un * un * vn + un * ((2*uhat_t-u) * v) ) * dx(element_boundary = True, bonus_intorder=order)
conv.Assemble()
conv_operator = conv.mat

u1, v1 = V.TnT()
uh1, vh1 = Vhat.TnT()

massfacet = BilinearForm(tang(uh1)*tang(vh1)*dx(element_boundary=True)).Assemble()
mixedmassfacet = BilinearForm(u1*tang(vh1)*dx(element_boundary=True)).Assemble()
invfacet = massfacet.mat.CreateBlockSmoother(Vhat.CreateSmoothingBlocks(blocktype="facet"))

eu, euhat, es, ew, ep = X.embeddings
ru, ruhat, rw, rw, rp = X.restrictions
avg_op = eu@ru + euhat@invfacet@mixedmassfacet.mat@ru

t = 0; i = 0
tend = 10
vel = gfu.components[0]
scene = Draw (gfu.components[0], mesh, min=0, max=2, autoscale=False)

with TaskManager():
    while t < tend:
        res = astokes.mat*gfu.vec + conv_operator@avg_op * gfu.vec
        gfu.vec.data -= tau * inv * res    
        V.Average(gfu.components[0].vec)
        t = t + tau; i = i + 1
        if i%10 == 0: 
            scene.Redraw()