Do-nothing boundary for Stokes behaving weirdly

Hi,

I am implementing Stokes equation on surfaces, from Lehrenfeld-Lederer-Schöberl-2020. My code is posted below. I am now trying to have an inflow boundary and a outflow/do-nothing boundary. The solution is a constant parabolic velocity-profile in the x-direction, and a linear pressure in the x-direction [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_x.png[/attachment]. However, when I solve the problem something strange happens. At the do-nothing boundary I get a sinusoidal profile in the y-direction, and the x-direction is also slightly perturbed. I have attached plots. I suspect this has to do with how I set the boundary-conditions on the do-nothing boundary, but I believe I have followed LLS-2020 virtually to the point.

I am very thankful for any possible help/ideas.

Best regards,

Alvar


Code:

### First the definition of the solver:
import re
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 netgen.csg import *
from netgen.meshing import MeshingStep

from numpy import pi

def SurfaceStokes_HdivHDG(mesh, force, g, order_u, order_p, order_g, bcdict, exactsoldict = None, mass=False, refinements = 0, pLagrange = False):



"""
Mesh is a predefined surface mesh with boundary conditions set
force is the right-hand side f in the PDE
g is the divergence of the velocity field

order_u, order_p and order_g are the orders for the velocity u, the pressure p and the geometric mapping.
bcdict is a dicitonary of the form

bcdict = {'dir': [], 'donothing': []}

where the lists should have the entries of the correct DOFS and the corresponding u-value.
Any of the entries can be replaced by 'None'.
exactsoldict is also a dictionary of the form

exactsoldict = {'velocity': uex, 'pressure': pex}

where pex can be set to 'None' if not needed/known. Any of the entries can be replaced by 'None'.

mass is a boolean which determines whether the Stokes problem showuld include a mass matrix,
necessary for uniqueness if the surface is closed.
refinements is the number of times the mesh should be refined before calculating, necessary for convergence tests.
"""
# Extract eventual solutions
if exactsoldict != None:
uex = exactsoldict
pex = exactsoldict
else:
uex = None
pex = None

# Refine the mesh and give it the correct geometric order
for i in range(refinements):
mesh.Refine()
mesh.Curve(order_g)



# Defines the mixed Finite element space. Here Wh is the velocity space, What is the facet space defined for HDG,
# Qh is the pressure space and Nh is the real line, needed for Lagrange multiplier to p, for uniqueness. They are
# put together in the mixed space Vh.
Wh = HDivSurface(mesh, order = order_u)
What = HCurl(mesh, order = order_u, flags={"orderface": 0})
Qh = SurfaceL2(mesh, order=order_p)
if pLagrange:
Nh = NumberSpace(mesh)
Vh = FESpace([Wh, What, Qh, Nh])
else:
Vh = FESpace([Wh, What, Qh])




if bcdict != None:
print("Here")

dirichlet = ""

for key in bcdict:
if key != "outflow":
dirichlet += "|"+key

dirichlet = dirichlet[1:]
print(dirichlet)

mark_boundary_dofs([Vh.FreeDofs(False), Vh.FreeDofs(True)], "inflow|wall", Vh, clear=True)
else:
SetFreeDofs(Vh)




gamma = 4*order_u*(order_u+1)


print("Step 1 done")

# Defines the trial- and testfunctions
if pLagrange:
u, uhat, p, r = Vh.TrialFunction()
v, vhat, q, s = Vh.TestFunction()
else:
u, uhat, p = Vh.TrialFunction()
v, vhat, q = Vh.TestFunction()

uhat = uhat.Trace()
vhat = vhat.Trace()



# Sets up the normal n, the face-tangential vector tau and the tangential co-normal mu
n = specialcf.normal(3)
tau = specialcf.tangential(3)
mu = Cross(n,tau)
# Defines the tangential projection, needed for HDG as the Hdiv space is normal continuous.
tang = lambda v: v-(v*mu)*mu

# Defines the projection matrix proj, the mesh size h and the stability parameter gamma.
nmat = CoefficientFunction(n, dims = (3,1))
proj = Id(3) - nmat*nmat.trans
h = specialcf.mesh_size



# Defines the symmetric gradients epsu and epsv, AKA the displacement matrices
gradu = proj * CoefficientFunction((grad(u),), dims = (3,3)).trans * proj # full gradient
gradv = proj * CoefficientFunction((grad(v),), dims = (3,3)).trans * proj # full gradient

epsu = 0.5 * (gradu + gradu.trans)
epsv = 0.5 * (gradv + gradv.trans)

print("Step 1.1 done")

# Defines the Stokes bilinear form

a = BilinearForm(Vh)
a += SymbolicBFI(2*InnerProduct(epsu,epsv), BND)
a += SymbolicBFI(-2*InnerProduct(epsu*mu, tang(v.Trace()-vhat)), BND, element_boundary = True)
a += SymbolicBFI(-2*InnerProduct(epsv*mu, tang(u.Trace()-uhat)), BND, element_boundary = True) #Symmetrizing
a += SymbolicBFI(gamma/h * InnerProduct(tang(v.Trace()-vhat), tang(u.Trace()-uhat)), BND, element_boundary = True) #Stabilizing


print("Step 1.2 done")

a += SymbolicBFI(- div(u.Trace())*q.Trace(), BND) #Enforces divergence condition
a += SymbolicBFI(- div(v.Trace())*p.Trace(), BND)



print("Step 1.3 done")
if pLagrange:
a += SymbolicBFI(s.Trace()*p.Trace(), BND) #Uniqueness of pressure through Lagrange multiplier
a += SymbolicBFI(r.Trace()*q.Trace(), BND) #Symmetrizing

print("Step 1.4 done")

# Adds mass form if necessary
if mass:
a += SymbolicBFI(u.Trace()*v.Trace(), BND)

print("Step 1.5 done")

a.Assemble()

print("Step 2 done")

#Defines linear form for right-hand side

f = LinearForm(Vh)
f += SymbolicLFI( proj * force * v.Trace(), BND, bonus_intorder = order_u)
f += SymbolicLFI(-g*q.Trace(), BND) #Enforces divergence condition
if pex != None and pLagrange:
f += SymbolicLFI (s.Trace()*pex, BND) #Ensures the pressure has the correct mean value

f.Assemble()

print("Step 3 done")

# Homogenize the problem

gfu = GridFunction(Vh)

res = gfu.vec.CreateVector()
res.data = f.vec



if bcdict != None:

dof_in = BitArray(Vh.components[0].FreeDofs())
dof_in.Clear()
mark_boundary_dofs(dof_in, "inflow", Vh.components[0], clear = False)

u_in = bcdict


mass = BilinearForm(Vh.components[0])
mass += SymbolicBFI((Vh.components[0].TestFunction().Trace()*mu) * (Vh.components[0].TrialFunction().Trace()*mu), BND, element_boundary = True)
mass.Assemble()


mass_rhs = LinearForm(Vh.components[0])
mass_rhs += SymbolicLFI((Vh.components[0].TestFunction().Trace()*mu) * (u_in * mu), BND, element_boundary = True)
mass_rhs.Assemble()

gfu.components[0].vec.data = mass.mat.Inverse(dof_in)*mass_rhs.vec


#Draw(gfu.components[0], mesh, 'u_in')


res -= a.mat*gfu.vec
else:
SetFreeDofs(Vh)



# Inverts the matrix and solves the problem
ainv = a.mat.Inverse(Vh.FreeDofs(), inverse = "umfpack")
gfu.vec.data += ainv * res

uh = gfu.components[0]
ph = gfu.components[2]

# Calculates eventual L2-errors of the solution

if uex != None:
l2u = sqrt(Integrate((uex-uh)**2, mesh, BND))
else:
l2u = None
if pex != None:
l2p = sqrt(Integrate((pex-ph)**2, mesh, BND))
else:
l2p = None

return {'velocity': uh, 'pressure': ph, 'l2u': l2u, 'l2p': l2p}

def mark_boundary_dofs(bitarrays, boundary, Vh, component=None, clear = False):
# bitarray should be a bitarray for the DOFS of the space Vh
if type(bitarrays) != list:
bitarrays = [bitarrays]

# Define what part of the space you want to mark the dofs for (if applicable), and
# start counting dofs at the correct spot

offset = 0
if component!= None:
X = Vh.components[component]
if components>0:
for i in range(component):
offset += Vh.components.ndof
else:
X = Vh

# Searches through all the relevant dofs, check if they match the correct type of boundary, and either
# clears or sets the dof.
for el in X.Elements(BBND):
bnd_name = mesh.GetBBoundaries()[el.index]
match = re.search(boundary, bnd_name)
if match != None:
for dofs in el.dofs:
for ba in bitarrays:
if clear:
ba.Clear(dofs+offset)
else:
ba.Set(dofs+offset)

def SetFreeDofs(V):
for el in V.components[0].Elements(BBND):
for dofs in el.dofs:
V.FreeDofs().Clear(dofs)
V.FreeDofs(True).Clear(dofs)
for el in V.components[1].Elements(BBND):
for dofs in el.dofs:
V.FreeDofs().Clear(dofs+V.components[0].ndof)
V.FreeDofs(True).Clear(dofs+V.components[0].ndof)



### And now the specific example

import netgen.gui

geo = CSGeometry()
bottom = Plane (Pnt(0,0,0), Vec(0,0, 1) )

surface = SplineSurface(bottom)
pts = [(0,0,0),(0,1,0),(1,1,0),(1,0,0)]
geopts = [surface.AddPoint(*p) for p in pts]
for p1,p2,bc in [(0,1,"inflow"), (1, 2,"wall"),(2,3,"outflow"),(3,0,"wall")]:
surface.AddSegment(geopts[p1],geopts[p2],bc)
geo.AddSplineSurface(surface)

ngmesh = geo.GenerateMesh(maxh=1)
mesh = Mesh(ngmesh)
for i in range(3):
mesh.Refine()




uex = CoefficientFunction(( y*(1 - y), 0, 0 )).Compile(True)
g = CoefficientFunction( 0 ).Compile(True)
pex = CoefficientFunction( 1 - x ).Compile(True)
force = CoefficientFunction(( 1.00000000000000, 0, 0 )).Compile(True)

uin = uex
bcdict = {'inflow': uin, 'wall': None, 'outflow': None}

exactsoldict = {'velocity': uex, 'pressure': pex}

solutiondict = SurfaceStokes_HdivHDG(mesh, force, g, order_u=2, order_p=1, order_g=2, bcdict=bcdict, exactsoldict = exactsoldict, refinements = 0, mass=False, pLagrange = False)


import netgen.gui
uh = solutiondict
ph = solutiondict


Draw(uex, mesh, 'uex')
Draw(uh, mesh, 'u')
Draw(Norm(uh), mesh, '|u|')
Draw(Norm(uh - uex), mesh, '|uh - uex|')
Draw(Norm(uex), mesh, '|uex|')
Draw(pex, mesh, 'pe [attachment=undefined]p.png[/attachment]x')
Draw(ph, mesh, 'p')
Draw(g - div(uh), mesh, 'g')
print(solutiondict)
print(solutiondict)