Refining mesh - Stokes problem

Hi,

I am working on the stationary Stokes problem. My advisor told me that the stability of Taylor-hood elements for the Stokes problem when using triangular elements requires that each element only has maximum one side on the boundary. Hence, we thought to use a large mesh that satisfies this, and then refine it. However, the error of the solution diverges with more refinements. Any idea's? The code is found below.

Best regards,

Alvar

from ngsolve import *
from ngsolve.internal import *
# basic xfem functionality
# basic geometry features (for the background mesh)
from netgen.geom2d import SplineGeometry
import numpy as np
import netgen.gui

from numpy import pi

import scipy
import scipy.sparse as sp

geo = SplineGeometry()
geo.AddRectangle( (0, 0), (1, 1), bcs = ("top", "left", "bottom", "right"))


#uex = CoefficientFunction((2*pi*sin(pi*x)**2*cos(pi*y)*sin(pi*y), -2*pi*sin(pi*x)*cos(pi*x)*sin(pi*y)**2))
#force = CoefficientFunction((12*pi**3*sin(pi*x)**2*sin(pi*y)*cos(pi*y) - pi*sin(pi*x)*cos(pi*y) - 4*pi**3*sin(pi*y)*cos(pi*x)**2*cos(pi*y),-12*pi**3*sin(pi*x)*sin(pi*y)**2*cos(pi*x) + 4*pi**3*sin(pi*x)*cos(pi*x)*cos(pi*y)**2 - pi*sin(pi*y)*cos(pi*x)))
#pex = CoefficientFunction(cos(pi*x)*cos(pi*y))


nu = 1
#uex = CoefficientFunction((sin(pi*y), -cos(pi*x)))
#force = CoefficientFunction((-pi*sin(pi*x)*cos(pi*y) + nu * pi**2*sin(pi*y), -pi*sin(pi*y)*cos(pi*x) - nu * pi**2*cos(pi*x)))
#pex = CoefficientFunction(cos(pi*x)*cos(pi*y))

#uex = CoefficientFunction((sin(y) + cos(y), 0))
#force = CoefficientFunction((sin(y) + cos(y) + cos(x),0))
#pex = sin(x)

uex = CoefficientFunction((y*(1-y), 0))
force = CoefficientFunction((3,0))
pex = CoefficientFunction(x-1/2)

l22_list = []
l2p_list = []

def my_abs(mesh, v, h0):
a = 0

x = np.arange(0,1+h0, h0)
for i in range(len(x)):
for j in range(len(x)):
p1 = np.array([x])
p2 = np.array([x[j]])
mip = mesh(p1, p2)
val = abs(v(mip))
a = max(a, val)
return a



def l2error(refinement = 0, structured = False, quad = False):

order = 2
alpha = 10
h0 = 0.3

if structured:
mesh = Mesh(geo.GenerateMesh(maxh = h0, quad_dominated = quad))
for i in range(refinement):
mesh.Refine()
#Draw(mesh)
else:
mesh = Mesh(geo.GenerateMesh(maxh = h0*2**(-refinement), quad_dominated = quad))


#Wh = HDiv(mesh, order = order+1, dgjumps = True, BDM = True) #Velocity space
Wh = VectorH1(mesh, order = order, dgjumps = True) #Velocity space

Qh = L2(mesh, order = order-1) #Pressure space
Nh = NumberSpace(mesh) #Lagrange multiplier

Vh = FESpace([Wh,Qh,Nh], dgjumps = True)

u,p,r = Vh.TrialFunction()
v,q,s = Vh.TestFunction()

n = specialcf.normal(2)
h = specialcf.mesh_size

tang = lambda v:v #If no tangential projections are wanted
#tang = lambda v: v-(v*n)*n #If tangential projections are wanted

jump_u = tang(u-u.Other())
jump_v = tang(v-v.Other())

flux_u = 0.5*(grad(u) + grad(u.Other()))*n
flux_v = 0.5*(grad(v) + grad(v.Other()))*n


a = BilinearForm(Vh)

#Volume integrals
a += SymbolicBFI(nu*InnerProduct(grad(u), grad(v)))
a += SymbolicBFI(-div(u)*q) #consistent symmetrizing term
a += SymbolicBFI(-div(v)*p)

a += SymbolicBFI(p*s)
a += SymbolicBFI(q*r)

#Interior facet integrals
a += SymbolicBFI(-nu*InnerProduct(flux_u,jump_v), VOL, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(flux_v,jump_u), VOL, skeleton = True) #consistent symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(jump_u,jump_v), VOL, skeleton=True) #consistent stabilizing term

#Boundary facet integrals

a += SymbolicBFI(-nu*InnerProduct(Grad(u)*n,tang(v)), BND, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(Grad(v)*n,tang(u)), BND, skeleton = True) #symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(tang(u),tang(v)), BND, skeleton=True) #stabilizing term
a += SymbolicBFI(v*n*p, BND, skeleton=True)
a += SymbolicBFI(u*n*q, BND, skeleton=True) #symmetrizing term

#Right-hand side

f = LinearForm(Vh)
f += SymbolicLFI(force*v)
f += SymbolicLFI(-nu*InnerProduct(Grad(v)*n,uex),BND, skeleton = True) #symmetrizing term
f += SymbolicLFI(h**(-1) * alpha * order**2 * InnerProduct(uex,v), BND, skeleton=True) #stabilizing term
f += SymbolicLFI(uex*n*q, BND, skeleton = True) #symmetrizing term
f += SymbolicLFI(pex*s)

a.Assemble()
f.Assemble()


gfu = GridFunction(Vh)

gfuex = GridFunction(Wh)
gfuex.Set(uex)
gfpex = GridFunction(Qh)
gfpex.Set(pex)

a_inv = a.mat.Inverse(freedofs = Vh.FreeDofs(), inverse="umfpack")

gfu.vec.data = a_inv*f.vec
vel = gfu.components[0]
pres = gfu.components[1]
#Draw(pres, mesh, "p")
#Draw(vel, mesh, "u")
#Draw(div(vel), mesh, "div(u)")
#Draw(Norm(vel), mesh, "|u|")


dv = grad(vel)
dg = grad(gfuex)

divu = div(vel)

#eu = sum((dv-dg)**2)
eu = (vel - gfuex)**2
ep = (pres - gfpex)**2

err = vel-gfuex

Draw(Norm(err), mesh, "|e|")

l22 = Integrate(eu, mesh)
l2p = Integrate(ep, mesh)

l22_list.append(l22)
l2p_list.append(l2p)

print("\nRefinement level ", refinement)
print("L2 error in velocity: ", sqrt(l22))
print("L2 error in pressure: ", sqrt(l2p))
print("div(u)-norm", Integrate(divu**2, mesh))

if len(l22_list)>1:
print("EOCu: ", (np.log(sqrt(l22_list[-2]))-(np.log(sqrt(l22_list[-1]))))/np.log(2))
print("EOCp: ", (np.log(sqrt(l2p_list[-2]))-(np.log(sqrt(l2p_list[-1]))))/np.log(2))

for i in range(3,4):
l2error(i, structured = True, quad = False)

I have solved it, I had to switch from L2 to H1(dgjumps=true) in the pressure space. But why should this be needed? Isnt L2-space same as H1(dgjumps=true)?

It's not - just look at fes.ndof

The dgjumps=True does not change your FESpace. It only changes the non-zero entries to be reserved in the matrix. In NGSolve the sparsity pattern is set up before the actual assembling of the element entries. Typically, the sparsity is precomputed based on overlapping support of basis functions. If two basis functions have support on the same element, a non-zero entry is reserved in the sparse matrix. For DG-type couplings this is however not enough. If you are on a facet and want to involve all functions of the two adjacent elements, some of these may not overlap, but may still interact due to some jump/derivative-interaction. To make sure to increase the sparsity pattern we introduced the dgjumps-Flag. If it is set to true all basis functions that are linked (through the neighboring volume element) to the same facet, get a non-zero entry now.

I hope that helps (and that I didn't forget to discuss some corner cases...).

So, L2 is not the same as H1(dgjumps=True). However, L2 is the same as Discontinuous(H1), but this is a different story...

Best,
Christoph

Thank you for your help. I will stay away from L2 in the future