Hdiv for Stokes, error at boundary

Hello,

I am trying to implement the Stokes problem for Discontinuous Galerkin (and eventually HDG) H(div) elements. The code is posted below. I do not however get the correct order of convergence. When I plot the norm of the error it is clear that something is wrong with the boundary integrals, but I cannot figure out what. I have attached the plot of the norm of the error.

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), (0.99, 0.99), bcs = ("top", "left", "bottom", "right"))




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))

l22_list = []


def l2error(refinement = 0):
h0 = 0.1
order = 2
alpha = 10

mesh = Mesh(geo.GenerateMesh(maxh = h0*2**(-refinement), quad_dominated = False))
#Draw(mesh)

Wh = HDiv(mesh, order = order, dirichlet="top|left|bottom|right", dgjumps = True BDM = True) #Velocity space
Qh = L2(mesh, order = 0) #Pressure space

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

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

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

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

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)

#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,v), BND, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(grad(v)*n,u), BND, skeleton = True) #symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(u,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

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)

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

err = vel-gfuex

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

l22 = Integrate(eu, mesh)
l22_list.append(l22)
print("h = ", h0*2**(-refinement),", L2^2 error in velocity: ", l22)

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

for i in range(2,3):
l2error(i)


Best regards,

Alvar

Hi Alvar,

You use Dirichlet-Flags for the Hdiv-FESpace, but don't set the Dirichlet-Values. Hence, your discrete solution will have zero values at the dirichlet dofs of you Hdiv space, i.e. u·n is zero on all boundaries. Your Nitsche-type terms are not acting here as u·n and v·n vanishes with the Dirichlet-Flags. Either remove the dirichlet-flags from the space or (preferably) use gfu.components[0].Set(uex) to set the correct Dirichlet values.

Further:
· You may want to use "Grad" instead of "grad" (the latter is not the Jacobi matrix for hdiv elements).
· As you are normal-continuous, you may want to add a tangential projection for all jump and boundary terms of u and v.


Best,
Christoph