L2 Error of pressure blows up outside the unit sphere

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

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

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 / (2*mu[0])*(x**2+y**2) apos = 1/ (2*mu[1])*(x**2+y**2) gammaf=0 vel_exact = [a*CoefficientFunction((-y, x ,0)) for a in [aneg, apos]] pres_exact =[x,x] # define RHS coef_g = [CoefficientFunction((1+4*y, -4*x,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))

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