Commit 8143b898 authored by Tilman Alemán's avatar Tilman Alemán
Browse files

working convergence modulo killing

parent 22fa5570
import sys
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
from ngsolve.internal import *
from xfem import *
from xfem.lsetcurv import *
from math import pi
import scipy.sparse as sp
from scipy.sparse.linalg import eigs
import numpy as np
import time
# -------------------------------- PARAMETERS ---------------------------------
# Interpolation for the rate-of-strain tensor when computing the r.h.s.
# Introduces an interpolation error, but significantly improves performance
interpol = True
# Mesh diameter
maxh = 0.6
# Geometry
geom = "ellipsoid"
# Polynomial order of FE space
order = 2
hs=10**8
# Problem parameters
reac_cf = 1
diff_cf = 1
# Ellipsoid parameter
c = 1
if c==1:
maxk = 3
else:
maxk=1
# Geometry
#SetHeapSize(10**9)
#SetNumThreads(48)
with TaskManager():
cube = CSGeometry()
if geom == "ellipsoid":
phi = Norm(CoefficientFunction((x , y, z/c))) - 1
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
elif geom == "decocube":
phi = (x ** 2 + y ** 2 - 4) ** 2 + (y ** 2 - 1) ** 2 + (y ** 2 + z ** 2 - 4) ** 2 + (x ** 2 - 1) ** 2 + (
x ** 2 + z ** 2 - 4) ** 2 + (z ** 2 - 1) ** 2 - 13
cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
# Preliminary refinements. As assembling the r.h.s. is very expensive without interpolation, keep it small in that case.
# When interpolating, getting good results requires more refinements with complex geometries, e.g. the decocube
n_cut_ref = int(sys.argv[1])
# Class to compute the mesh transformation needed for higher order accuracy
# * order: order of the mesh deformation function
# * threshold: barrier for maximum deformation (to ensure shape regularity)
if geom == "ellipsoid":
kvfs = [CoefficientFunction((y, -x, 0)),CF((z,0,-x)), CF((0,z,-y))]
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
deformation = lsetmeshadap.CalcDeformation(phi)
lset_approx = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
# Background FESpaces
# Helper Functions for tangential projection and gradients of the exact solution
def coef_grad(u):
dirs = {1: [x], 2: [x, y], 3: [x, y, z]}
if u.dim == 1:
return CF(tuple([u.Diff(r) for r in dirs[mesh.dim]]))
else:
return CF(tuple([u[i].Diff(r) for i in range(u.dim) for r in dirs[mesh.dim]]), dims=(u.dim, mesh.dim))
def grad(u):
if type(u) in [ProxyFunction, GridFunction]:
return ngsolve.grad(u)
else:
return coef_grad(u)
def Pmats(n):
return Id(3) - OuterProduct(n, n)
# Coefficients / parameters:
# Exact tangential projection
#.Compile(realcompile=True)
def doKilling(mesh, lset_approx, deformation, step):
Ps = Pmats(Normalize(grad(phi)))
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
n = Normalize(grad(lset_approx))
ntilde = Interpolate(Normalize(grad(phi)), VhG)
h = specialcf.mesh_size
elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
cuthmax = max([(2*vol)**(1/2) for vol in elvol])
print("hmax:")
print(cuthmax)
alpha=3
# epsilon =(cuthmax)**(alpha)
epsilon=h**alpha
# print("Epsilon:")
# print(epsilon)
# epsilon=1
eta = 1000.0 / (h)**2
rhou = 1.0 / h
rhop = h
# Measure on surface
ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
# Measure on the bulk around the surface
dX = dx(definedonelements=ba_IF, deformation=deformation)
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = ngsolve.grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
if step>4:
comp=True
else:
comp=False
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat)
u2, v2 = VhG.TnT()
a2 = BilinearForm(VhG, symmetric=True)
a2 += InnerProduct(E(u2),E(v2)) * ds
# a2 +=(Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
a2 += (eta * ((u2 * n) * (v2 * n))) * ds
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(u2,v2)*dX
b.Assemble()
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 0)
egfu = GridFunction(VhG)
egfu.vec.data = evecs.vecs[0].data
print("is it eigen?")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat)**2,mesh=mesh))
fun = GridFunction(VhG)
kf = GridFunction(VhG)
l2errs = []
print(len(vlam))
for k in range(3):
lam = vlam[k]
print(k)
print("lam: "+ str(abs(lam)))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
evecnorm = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5))
kvfnorm = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kvfs[k],kvfs[k]), mesh=mesh, order=5))
scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k)/evecnorm,kvfs[k]/kvfnorm), mesh=mesh, order=5)
sgn = np.sign(scal)
l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(sgn*kvfs[k]/kvfnorm-evecs.MDComponent(k)/evecnorm)**2, mesh=mesh, order=5)/evecnorm))
return l2errs, evecs.MDComponent(0), evecs.MDComponent(1), evecs.MDComponent(2)
l2errors = []
for i in range(n_cut_ref+1):
if i!=0:
mesh.SetDeformation(None)
lsetp1 = GridFunction(H1(mesh, order=2))
InterpolateToP1(phi, lsetp1)
RefineAtLevelSet(lsetp1)
mesh.Refine()
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
deformation = lsetmeshadap.CalcDeformation(phi)
lset_approx = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
print("starting solve..")
l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
l2errors.append(l2err)
print("L2 errors:")
print(l2errors)
if False:# and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
# input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
# with TaskManager():
if maxk==1:
vtk = VTKOutput(ma=mesh,
coefs=[lset_approx, ev1, kvfs[0]],
names=["P1-levelset", "ev", "kf"],
filename="killing", subdivision=0, legacy = False)
vtk.Do()
else:
vtk = VTKOutput(ma=mesh,
coefs=[lset_approx, ev1, ev2, ev3, kvfs[0], kvfs[1], kvfs[2]],
names=["P1-levelset", "ev1", "ev2", "ev3", "kf1", "kf2", "kf3"],
filename="killing", subdivision=0, legacy = False)
vtk.Do()
......@@ -41,15 +41,23 @@ maxh = 0.6
geom = "ellipsoid"
# Polynomial order of FE space
order = 2
if order==2:
hs = 10**6
else:
hs= 10**8
# Problem parameters
reac_cf = 1
diff_cf = 1
# Ellipsoid parameter
c = 1
if c==1:
maxk = 3
else:
maxk=1
# Geometry
#SetHeapSize(10**9)
SetNumThreads(16)
#SetNumThreads(48)
with TaskManager():
cube = CSGeometry()
if geom == "ellipsoid":
......@@ -67,12 +75,12 @@ with TaskManager():
# When interpolating, getting good results requires more refinements with complex geometries, e.g. the decocube
n_cut_ref = int(sys.argv[1])
# Class to compute the mesh transformation needed for higher order accuracy
# * order: order of the mesh deformation function
# * threshold: barrier for maximum deformation (to ensure shape regularity)
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
deformation = lsetmeshadap.CalcDeformation(phi)
lset_approx = lsetmeshadap.lset_p1
......@@ -104,7 +112,7 @@ with TaskManager():
# Coefficients / parameters:
# Exact tangential projection
Ps = Pmats(Normalize(grad(phi))).Compile(realcompile=True)
Ps = Pmats(Normalize(grad(phi)))#.Compile(realcompile=True)
# define solution and right-hand side
......@@ -120,8 +128,9 @@ with TaskManager():
}
pSol = (x * y ** 3 + z * (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)) / ((x ** 2 + y ** 2 + z ** 2) ** 2)
uSol = Ps* CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))"""
uSol = (Ps* CF((-z**2,x,y))).Compile(realcompile=True)
uSol = (Ps* CF((-z**2,x,y)))#.Compile(realcompile=True)
pSol = CF(x*y**3+z)
kvfs = [CoefficientFunction((y, -x, 0)),CF((z,0,-x)), CF((0,z,-y))]
elif geom == "decocube":
uSol = Ps*CoefficientFunction((-z ** 2, y, x))
......@@ -131,7 +140,7 @@ with TaskManager():
cf=x * y ** 3 + z, mesh=mesh))
weingex = grad(Normalize(grad(phi))).Compile(realcompile=True)
# Helper Functions for rate-of-strain-tensor
def eps(u):
return Ps * Sym(grad(u)) * Ps
......@@ -149,7 +158,8 @@ with TaskManager():
N = 3
r = Grad(u) * Ps
return CF(tuple([Trace(r[N*i:N*(i+1),0:N]) for i in range(N)]))
def solveStokes(mesh, lset_approx, deformation, step):
def solveStokes(mesh, lset_approx, deformation, step, uSol):
uSol = (Ps* CF((-z**2,x,y)))
Vh = VectorH1(mesh, order=order, dirichlet=[])
Qh = H1(mesh, order=order - 1, dirichlet=[])
......@@ -163,17 +173,18 @@ with TaskManager():
gfu = GridFunction(X)
n = Normalize(grad(lset_approx))
ntilde = Interpolate(grad(phi), VhG)
ntilde = Interpolate(Normalize(grad(phi)), VhG)
h = specialcf.mesh_size
elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
cuthmax = max([(2*vol)**(1/2) for vol in elvol])
alpha=2
epsilon =(cuthmax)**(alpha)
# elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
# cuthmax = max([(2*vol)**(1/2) for vol in elvol])
alpha=3
# epsilon =(cuthmax)**(alpha)
epsilon=h**alpha
print("Epsilon:")
print(epsilon)
# epsilon=1
eta = 100.0 / (h * h)
eta = 1000.0 / (h)**2
rhou = 1.0 / h
rhop = h
......@@ -188,13 +199,14 @@ with TaskManager():
if interpol:
tmpeps = GridFunction(Lpmat)
tmpeps.Set(eps(uSol) ,definedonelements=ba_IF)
tmpdiv = GridFunction(Lpvec)
tmpdiv.Set(Ps*-divG(tmpeps), definedonelements=ba_IF)
tmpeps.Set(eps(uSol)-(uSol*Normalize(grad(phi)))*weingex ,definedonelements=ba_IF)
# tmpdiv = GridFunction(Lpvec)
# tmpdiv.Set(Ps*-divG(tmpeps), definedonelements=ba_IF)
# tmpuSol = GridFunction(L2(mesh, order=order+1, dim=3))
# tmpuSol.Set(uSol, definedonelements=ba_IF)
tmptotal = GridFunction(Lpvec)
tmptotal.Set(tmpdiv+Ps*grad(pSol), definedonelements=ba_IF)
# tmptotal.Set(tmpdiv+Ps*grad(pSol), definedonelements=ba_IF)
tmptotal.Set(Ps*(-divG(tmpeps)+grad(pSol)), definedonelements=ba_IF)
return tmptotal
else:
......@@ -204,8 +216,10 @@ with TaskManager():
# Weingarten mappings
weing = grad(n)
weingex = grad(Normalize(grad(phi)))
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = ngsolve.grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
L21 = L2(mesh, order=order+1, dim=9)
L22 = L2(mesh, order=order+1, dim=3)
......@@ -222,7 +236,7 @@ with TaskManager():
def E(u):
return (Pmat * Sym(grad(u)) * Pmat)#.Compile(realcompile=comp)# - (u * n) * weing
return (Pmat * Sym(grad(u)) * Pmat)- (u * n) * weing
a = BilinearForm(X, symmetric=True)
......@@ -233,14 +247,14 @@ with TaskManager():
a += InnerProduct(u, Pmat * grad(q)) * ds
a += InnerProduct(v, Pmat * grad(p)) * ds
a += -rhop * InnerProduct(n * grad(p), n * grad(q)) * dX
a.Assemble()
u2, v2 = VhG.TnT()
a2 = BilinearForm(VhG, symmetric=True)
a2 += InnerProduct(E(u2),E(v2)).Compile(realcompile=comp) * ds
# a2 +=epsilon* (Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
# a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
a2.Assemble()
......@@ -250,7 +264,7 @@ with TaskManager():
b.Assemble()
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 0)
# vlam, evecs = solvers.PINVIT(a2.mat, b.mat,pre=IdentityMatrix(), num=3)
"""
......@@ -263,52 +277,83 @@ with TaskManager():
"""
# R.h.s. linear form
f = LinearForm(X)
f += rhs1 * v * ds
f += -rhs2 * q * ds
f.Assemble()
a.Assemble()
"""f2 = LinearForm(X)
f2 += rhs1ex * v * ds
f2 += -rhs2 * q * ds
gfu2 = GridFunction(X)
f2.Assemble()"""
egfu = GridFunction(VhG)
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
print("removing killing VFs")
print(vlam)
fun = GridFunction(VhG)
kf = GridFunction(VhG)
for k, lam in enumerate(vlam):
print("lam: "+ str(abs(lam)))
print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if False and k==0:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
kf.vec.data =evecs.vecs[0].data
evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
print("removing number: "+str(k))
integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh, order=5)/evecnorm
print("scalar product before removing killing field: " + str(integral))
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh)
gfu.components[0].vec.data+= integral*evecs.vecs[k]
# print(gfu.components[0].vec.data-kf.vec.data)
print("scalar product with killing field: ")
"""
for k in range(3):
print(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF},
cf = InnerProduct(gfu.components[0], evecs.MDComponent(k)), mesh=mesh, order=10))
"""
#gfu2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso")*f2.vec
l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
uerr = gfu.components[0] - uSol
print("l2error : ", l2err)
return l2err, uerr, gfu, kf
if True:
print("removing killing VFs")
fun = GridFunction(VhG)
kf = GridFunction(VhG)
print(vlam)
for k, lam in enumerate(vlam):
print("lam: "+ str(abs(lam)))
egfu.vec.data = evecs.vecs[k].data
print("is it eigen?")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat)**2,mesh=mesh))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if False and k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
print("removing number: "+str(k))
evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh, order=5)/evecnorm
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh)
gfu.components[0].vec.data+= integral*evecs.vecs[k]
# print(gfu.components[0].vec.data-kf.vec.data)
kf2 = GridFunction(VhG)
for k in range(maxk):
# kf2.Set(kvfs[k],definedonelements=ba_IF)
# kfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kf2,kf2), mesh=mesh, order=5)
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu2.components[0],kf2), mesh=mesh, order=5)/kfnorm
kfexnorm = Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
uSol = uSol - Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(uSol, kvfs[k]), mesh=mesh, order=5) / kfexnorm * kvfs[k]
# gfu2.components[0].vec.data += integral*kf2.vec
kf.Set(uSol,definedonelements = ba_IF)
kf2.vec.data = gfu.components[0].vec.data-kf.vec.data
for k in range(maxk):
kf.Set(kvfs[k],definedonelements=ba_IF)
integral = -Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kf2, kvfs[k]), mesh=mesh, order=5)/Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
kf2.vec.data += integral*kf.vec
# gfu.components[0].vec.data +=-Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
# cf=InnerProduct(gfu.components[0], kvfs[k]), mesh=mesh, order=5) / kfexnorm * kf.vec
#gfu2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso")*f2.vec
l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
l2errnokf = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kf2,kf2) , mesh=mesh))
uerr = gfu.components[0] - uSol
print("l2error : ", l2err)
print("error without killing error: ", l2errnokf)
return l2err, l2errnokf, kf2, gfu, uSol
else:
l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
return l2err, None, None, None, None, None
l2errors = []
l2errkf = []
for i in range(n_cut_ref+1):
if i!=0:
mesh.SetDeformation(None)
......@@ -318,7 +363,7 @@ with TaskManager():
mesh.Refine()
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
deformation = lsetmeshadap.CalcDeformation(phi)
......@@ -328,17 +373,22 @@ with TaskManager():
print("starting solve..")
l2err, uerr, uh, kf = solveStokes(mesh, lset_approx, deformation, step=i)
l2err, l2errnokf, uerr, uh, uSolnokf = solveStokes(mesh, lset_approx, deformation, step=i, uSol=uSol)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
print(l2errors)
print("error with no killing field: ")
print(l2errkf)
if False:# and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
# input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
# with TaskManager():
vtk = VTKOutput(ma=mesh,
coefs=[lset_approx, uh.components[0],kf, uerr],
names=["P1-levelset", "uh","KF", "uerr"],
filename="stokestracefem3d", subdivision=2, legacy = False)
coefs=[lset_approx, uh.components[0], uerr, uSolnokf],
names=["P1-levelset", "uh", "uerr", "uex"],
filename="stokestracefem3d", subdivision=0, legacy = False)
vtk.Do()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment