L2 Error of pressure blows up outside the unit sphere

Hi,

I am working on a two-phase stationary Stokes problem. I am trying to modify the geometry in stokesxfem.py in the py_tutorial folder to a unit sphere in the 333 cube scenario. I modified the setting for geometry

cube = CSGeometry() cube.Add (OrthoBrick(Pnt(-1.5,-1.5,-1.5), Pnt(1.5,1.5,1.5))) mesh = Mesh(cube.GenerateMesh(maxh=0.5, quad_dominated=False)) names = [ "bottom", "top", "left", "right", "front", "back" ] for i, name in enumerate(names): mesh.ngmesh.SetBCName(i, name) phi = Norm(CoefficientFunction((x,y,z)))**2-1 lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=10.5, discontinuous_qn=True) deformation = lsetmeshadap.CalcDeformation(phi) lsetp1 = lsetmeshadap.lset_p1 ci = CutInfo(mesh, lsetp1) d = mesh.dim

and added a part to perform [tex]p^+{new}=p^+{old}-\frac{\int_{\Omega^+}p_{old}\ d\mathscr{H}^{3}}{|\Omega^+|} [/tex], where [tex]\Omega^+[/tex] is the cube minus (unit ball).

const_eliminate=GridFunction(WhG) p_eliminator=const_eliminate.components[1] ones_pressure = [1, 1] area=[Integrate(lset_doms[i], gfp.components[i], mesh=mesh)/Integrate(lset_doms[i],1, mesh=mesh) for i in [0, 1]] p_eliminator.components[0].Set(area[0]) p_eliminator.components[1].Set(area[1]) L0_gfp=[gfp.components[i]-p_eliminator.components[i] for i in [0, 1]]
When I print out the result, L2 error of velocity and L2 error of pressure within the unit sphere seems fine. But L2 error of pressure outside unit sphere seems blows up. Here is the result I am getting,
For maxh=0.5
L2 Error of velocity: 0.013408092886879695
L2 Error of pressure within unit sphere: 0.018069442530407473
L2 Error of pressure outside unit sphere: 9289563.362909812
For maxh=0.25
L2 Error of velocity: 0.0014547447652420093
L2 Error of pressure within unit sphere: 0.003051588344220756
L2 Error of pressure outside unit sphere: 148769.50294472638
I plot the picture of pressure I got, it seems p^+ blows up near the corner of the cube.
[attachment=undefined]Screen Shot 2020-11-30 at 8.53.29 PM.png[/attachment]
Am I missing some settings for geometry or FEM space? Thank you very much.
Here I attached my code:

[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from ngsolve import *
from xfem import *

visualization stuff

from ngsolve.internal import *

basic xfem functionality

basic geometry features (for the background mesh)

from netgen.csg import CSGeometry, OrthoBrick, Pnt

pi

from math import pi

For LevelSetAdaptationMachinery

from xfem.lsetcurv import *

3D: unit sphere configuration

cube = CSGeometry()
cube.Add (OrthoBrick(Pnt(-1.5,-1.5,-1.5), Pnt(1.5,1.5,1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=0.25, quad_dominated=False))
names = [ “bottom”, “top”, “left”, “right”, “front”, “back” ]
for i, name in enumerate(names):
mesh.ngmesh.SetBCName(i, name)
phi = Norm(CoefficientFunction((x,y,z)))**2-1
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=10.5, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
ci = CutInfo(mesh, lsetp1)
d = mesh.dim
mu=[1.0,1.0]

set analytic problem

aneg = 1 / (2mu[0])(x2+y2)
apos = 1/ (2mu[1])(x2+y2)
gammaf=0
vel_exact = [a*CoefficientFunction((-y, x ,0)) for a in [aneg, apos]]
pres_exact =[x,x]

define RHS

coef_g = [CoefficientFunction((1+4y, -4x,0))for i in [0, 1]]

stabilization parameter for ghost-penalty

gamma_stab_vel = 0.05 # if set to zero: no GP-stabilization for velocity
gamma_stab_pres = 0.05

order=2

stabilization parameter for Nitsche

lambda_nitsche = 0.5 * (mu[0] + mu[1]) * 20 * order * order

Discretization

Vhbase = VectorH1(mesh, order=order, dirichlet=“bottom|top|left|right|front|back”)
Qhbase = H1(mesh, order=order - 1)

Vhneg = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))
Vhpos = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASPOS)))
Qhneg = Compress(Qhbase, GetDofsOfElements(Qhbase, ci.GetElementsOfType(HASNEG)))
Qhpos = Compress(Qhbase, GetDofsOfElements(Qhbase, ci.GetElementsOfType(HASPOS)))

WhG = FESpace([Vhneg * Vhpos, Qhneg * Qhpos, NumberSpace(mesh)], flags={“dgjumps”: True})

gfup = GridFunction(WhG)
gfu, gfp, gfn = gfup.components
gfu_neg, gfu_pos = gfu.components
gfp_neg, gfp_pos = gfp.components

h = specialcf.mesh_size

lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
lset_if = {“levelset”: lsetp1, “domain_type”: IF}

element, facet and dof marking w.r.t. boundary approximation with lsetp1:

ba_facets = [GetFacetsWithNeighborTypes(mesh, a=ci.GetElementsOfType(HASNEG), b=ci.GetElementsOfType(IF)),
GetFacetsWithNeighborTypes(mesh, a=ci.GetElementsOfType(HASPOS), b=ci.GetElementsOfType(IF))]

n_lset = 1.0 / Norm(grad(lsetp1)) * grad(lsetp1)
kappa = [CutRatioGF(ci), 1.0 - CutRatioGF(ci)]

a = BilinearForm(WhG, symmetric=False)
f = LinearForm(WhG)

u, p, n = WhG.TrialFunction()
v, q, m = WhG.TestFunction()

some helper expressions:

def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)

def sigma(i, u, p):
return - 2 * mu[i] * eps(u[i]) + p[i] * Id(mesh.dim)

def average_flux(u, p):
return sum([kappa[i] * sigma(i, u, p) * n_lset for i in range(2)])

def average_inv(u):
return sum([kappa[1 - i] * u[i] for i in range(2)])

def jump(u):
return u[0] - u[1]

Stokes variational formulation:

for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2 * mu[i] * InnerProduct(eps(u[i]), eps(v[i])))
a += SymbolicBFI(lset_doms[i], form=- div(u[i]) * q[i] - div(v[i]) * p[i])
f += SymbolicLFI(lset_doms[i], form=coef_g[i] * v[i])

constraining the pressure integral

a += SymbolicBFI(lset_neg, form=n * q[0] + m * p[0])

Nitsche parts:

a += SymbolicBFI(lset_if, form=average_flux(u, p) * jump(v))
a += SymbolicBFI(lset_if, form=average_flux(v, p) * jump(u))
a += SymbolicBFI(lset_if, form=lambda_nitsche / h * jump(u) * jump(v))

surface tension:

f += SymbolicLFI(lset_if, form=-gammaf * average_inv(v) * n_lset)

ghost penalty terms:

for i in range(2):
if gamma_stab_vel > 0:
a += SymbolicFacetPatchBFI(form=gamma_stab_vel / h ** 2 * (u[i] - u[i].Other()) * (v[i] - v[i].Other()),
skeleton=False, definedonelements=ba_facets[i])
if gamma_stab_pres > 0:
a += SymbolicFacetPatchBFI(form=- gamma_stab_pres * (p[i] - p[i].Other()) * (q[i] - q[i].Other()),
skeleton=False, definedonelements=ba_facets[i])

apply mesh adaptation

mesh.SetDeformation(deformation)

a.Assemble()
f.Assemble()

gfu.components[0].Set(vel_exact[0])

gfu.components[1].Set(vel_exact[1],BND)

Solve linear system

f.vec.data -= a.mat * gfup.vec
gfup.vec.data += a.mat.Inverse(WhG.FreeDofs()) * f.vec

enforce p\in L_0^2, by setting L0_gfp^+=gfp^±\integrate(gfp^+)/(\integrate(\Omega^+))

const_eliminate=GridFunction(WhG)
p_eliminator=const_eliminate.components[1]
ones_pressure = [1, 1]
area=[Integrate(lset_doms[i], gfp.components[i], mesh=mesh)/Integrate(lset_doms[i],1, mesh=mesh) for i in [0, 1]]
p_eliminator.components[0].Set(area[0])
p_eliminator.components[1].Set(area[1])
L0_gfp=[gfp.components[i]-p_eliminator.components[i] for i in [0, 1]]

#measure the error

vl2error = sqrt(sum([Integrate(lset_doms[i], Norm(gfu.components[i] - vel_exact[i]) ** 2, mesh=mesh) for i in [0, 1]]))
pl2error_inner =sqrt(sum([Integrate(lset_doms[i], (L0_gfp[i] - pres_exact[i]) ** 2, mesh=mesh) for i in [0]]))
pl2error_outer =sqrt(sum([Integrate(lset_doms[i], (L0_gfp[i] - pres_exact[i]) ** 2, mesh=mesh) for i in [1]]))

print(“L2 Error of velocity: {0}”.format(vl2error))
print(“L2 Error of pressure within unit sphere: {0}”.format(pl2error_inner))
print(“L2 Error of pressure outside unit sphere: {0}”.format(pl2error_outer))

[/code]

Hi Ryanleaf,

Great that you are using ngsxfem!
Away from the interface you are using a plain standard Taylor-Hood element. However, in the case of Dirichlet bnd. conds. for the velocity, on elements where there is no interior vertex, the Taylor-Hood element is known to be not inf-sup stable. This is the case for some of your corners.

You can fix this by refining the problematic elements. Simply add this after your mesh generation:

vs = [0 for i in range(len(mesh.vertices))] for el in mesh.Elements(): mesh.SetRefinementFlag(el,False) for v in el.vertices: vs[v.nr] += 1 for el in mesh.Elements(): for v in el.vertices: if vs[v.nr] == 1: mesh.SetRefinementFlag(el,True) mesh.Refine()

Best,
Christoph

Hi Christoph,

Thanks a lot for your help.
This part fixed the issue. I was very confused about this problem.
Thank you very much.

Have a nice day,
Ryan