Commit b72cfbf8 authored by Tilman Aleman's avatar Tilman Aleman
Browse files

Merge branch 'main' of git.rwth-aachen.de:tilmandiego.aleman/reproduction

parents 9ce434cc faf8e6b9
for ORDER in 1 2 3
do
for C in 1 1.1 1.5
do
./killing.py --refinements 5 -c $C --order $ORDER --plot --saveplot
done
done
for ORDER in 1 2 3
do
for C in "1" "0.5" "0.1"
do
./killing.py --geometry arrow --refinements 5 -c $C --order $ORDER --plot --saveplot
done
done
\ No newline at end of file
......@@ -4,9 +4,11 @@ import matplotlib.pyplot as plt
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import math
def refineMesh(mesh, phi, order, hs=10 ** 6):
if order>2:
hs=10**8
with TaskManager():
mesh.SetDeformation(None)
lsetp1 = GridFunction(H1(mesh, order=2))
lsetp1 = GridFunction(H1(mesh, order=order))
InterpolateToP1(phi, lsetp1)
RefineAtLevelSet(lsetp1)
......@@ -91,99 +93,78 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
hvals = []
ticks = range(1,len(l2errs)+1)
for k in range(len(l2errs)):
hvals.append((1 / 2 ** 3) ** k)
plt.semilogy(l2errs, label=geom+" with KF removal")
plt.semilogy(l2errskf, label=geom + " without KF removal")
plt.semilogy(hvals, label="O(h^"+ str(order+1)+")")
plt.ylabel("L2 error")
hvals.append((1 / 2 ** (order+1)) ** k)
hvals2 = []
for k in range(len(l2errs)):
hvals2.append((1 / 2 ** order) ** k)
plt.semilogy(ticks,l2errs, label=geom+" with KF removal")
# plt.semilogy(l2errskf, label=geom + " without KF removal")
plt.semilogy(ticks,hvals, label=r'$O(h^'+ str(order+1)+r')$')
plt.semilogy(ticks,hvals2, label=r'$O(h^'+ str(order)+r')$')
plt.ylabel(r'$L2$ error')
plt.xlabel("Refinement")
plt.title("Convergence removing KF")
plt.legend()
plt.xticks(range(len(l2errs)))
plt.xticks(ticks)
if saveplot:
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
else:
plt.show()
def computelimit(phi, order=2, iters=5, radius=1):
maxh = 0.6
cube = CSGeometry()
cubesize = 2*radius
with TaskManager():
cube.Add(OrthoBrick(Pnt(-cubesize, -cubesize, -cubesize), Pnt(cubesize, cubesize, cubesize)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(phi)
lset_approx = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
vlams = []
for k in range(3):
vlams.append([])
for k in range(max(iters,3)):
if k!=0:
lset_approx, deformation = refineMesh(mesh, phi, order)
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 = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
h = specialcf.mesh_size
# epsilon=1
eta = 1 / (h) ** 2
rhou = 1.0 / h
rhob = 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)
u2 = VhG.TrialFunction()
v2 = VhG.TestFunction()
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = grad(ninter)
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
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 * ntilde) * (v2 * ntilde))) * ds
pre = Preconditioner(a2, "direct")
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat * u2, Pmat * v2) * ds
b += rhob * (grad(u2) * n) * (grad(v2) * n) * dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
vlam, evecs = solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
for j in range(3):
vlams[j].append(vlam[j])
result = []
for j in range(3):
print(len(vlams[j]))
temp=0
temp += vlams[j][-1]
temp -= (vlams[j][-1]-vlams[j][-2])**2/((vlams[j][-1]-vlams[j][-2])-(vlams[j][-2]-vlams[j][-3]))
result.append(math.log(abs(temp-1),10)<math.log(abs(temp-vlams[j][-1]),10))
print("convergence for eigval ", j)
print("logdiff: ", math.log(abs(temp-vlams[j][-1]),10))
print("logex, no accel: ", math.log(abs(1-vlams[j][-1]),10))
print("logex: ", math.log(abs(temp-1),10))
print("logdiff, no accel: ",math.log(abs(vlams[j][-2]-vlams[j][-1]),10))
return result
\ No newline at end of file
def computeEV(phi, mesh, lset_approx, deformation, order=2):
vlams = []
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 = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
h = specialcf.mesh_size
# epsilon=1
eta = 1 / (h) ** 2
rhou = 1.0 / h
rhob = 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)
u2 = VhG.TrialFunction()
v2 = VhG.TestFunction()
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = grad(ninter)
Pmat = Pmats(n)#.Compile()
def E(u):
return ((Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing)#.Compile()
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 * ntilde) * (v2 * ntilde))) * ds
pre = Preconditioner(a2, "direct")
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat * u2, Pmat * v2) * ds
b += rhob * (grad(u2) * n) * (grad(v2) * n) * dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
......@@ -9,11 +9,11 @@ from helperfunctions import *
# -------------------------------- PARAMETERS ---------------------------------
def getKilling(mesh, lset_approx, deformation, step, alpha=3):
def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsym = 0, maxsymex = 0):
uSol = (Ps * CF((-z ** 2, x, y)))
uSol2 = (Ps * CF((-z ** 2, x, y)))
for k in range(1):
for k in range(maxsymex):
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},
......@@ -35,10 +35,13 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
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])
# 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])
# epsilon =(cuthmax)**(alpha)
alpha=order+1
epsilon = h ** alpha
if False:
epsilon=0
print("Epsilon:")
print(epsilon)
# epsilon=1
......@@ -87,23 +90,7 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
tmptotal = GridFunction(lpvec)
tmptotal.Set(Ps * (-divG(tmpeps, Ps, mesh)), definedonelements=ba_IF)
rhs1 = tmptotal
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 * ntilde) * (v2 * ntilde))) * ds
pre = Preconditioner(a2, "direct")
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat * u2, Pmat * v2) * ds
b += rhob * (grad(u2) * n) * (grad(v2) * n) * dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
vlam, evecs = solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
evecinter = []
for k in range(3):
......@@ -143,13 +130,13 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
# 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
for k in range(3):
for k in range(maxsym):
if True and (vlam[k]-1)+2*(vlam[k]-1)*(cuthmax**3-2*cuthmax**(5/2))<cuthmax**(3):
print("removed ", k)
print("condition: ", (vlam[k]-1)+2*(vlam[k]-1)*(cuthmax**3-2*cuthmax**(5/2)))
sca = l2sca(gfu, evecinter[k], phi, mesh)
gfu.vec.data += -sca * evecinter[k].vec
print("removed ", k)
sca = l2sca(gfu, evecinter[k], phi, mesh)
gfu.vec.data += -sca * evecinter[k].vec
l2err = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(gfu - uSol) ** 2, mesh=mesh))
l2errnokf = sqrt(
......@@ -175,7 +162,7 @@ if __name__ == "__main__":
parser.add_argument("--plot", action='store_true')
parser.add_argument("--saveplot", action='store_true')
parser.add_argument("--alpha", dest="alpha", type=float, required=False, default=3)
parser.add_argument("--radius", dest="radius", type=float, required=False, default=0.5)
parser.add_argument("--radius", dest="radius", type=float, required=False, default=1)
args = parser.parse_args()
interpol = True
......@@ -199,14 +186,14 @@ if __name__ == "__main__":
elif c>1 and geom=="ellipsoid":
maxk = 1
elif geom=="arrow":
maxk=1
maxk=0
# Geometry
# SetHeapSize(10**9)
# SetNumThreads(48)
cube = CSGeometry()
if geom == "ellipsoid":
phi = Norm(CoefficientFunction((x, y, z / c))) - args.radius
phi = (Norm(CoefficientFunction((x, y, z / c))) - args.radius).Compile(realcompile=True)
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
......@@ -271,20 +258,52 @@ if __name__ == "__main__":
weingex = grad(Normalize(grad(phi, mesh)), mesh).Compile(realcompile=True)
l2errors = []
l2errkf = []
vlams = []
symmetries = 0
for k in range(3):
vlams.append([])
for i in range(n_cut_ref + 1):
if i != 0:
lset_approx, deformation = refineMesh(mesh, phi, order)
print("starting solve..")
with TaskManager():
l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, step=i, alpha=args.alpha)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
for k in range(3):
vlams[k].append(vlam[k])
if i < 2 and i != 0:
symmetries = 0
for k in range(3):
# if math.log(abs(1 - vlams[k][1]), 10) < math.log(abs(vlams[k][0] - vlams[k][1]), 10):
if math.log(abs(1 - vlams[k][1]), 2) <-(order+1)*(i+1):
print("corresponding log: ")
print(math.log(abs(1 - vlams[k][1]), 2))
print("condition: ", -(order+1)*(i+1))
symmetries += 1
elif i >= 2:
symmetries = 0
for j in range(3):
temp = 0
temp += vlams[j][-1]
temp -= (vlams[j][-1] - vlams[j][-2]) ** 2 / (
(vlams[j][-1] - vlams[j][-2]) - (vlams[j][-2] - vlams[j][-3]))
# if math.log(abs(1 - vlams[j][-1]), 10) < math.log(abs(temp - vlams[j][-1]), 10):
if math.log(abs(1 - temp), 2) < -(order + 1) * (i+2):
print("corresponding log: ")
print(math.log(abs(1 - temp), 2))
print("condition: ", -(order+1)*(i+1))
symmetries += 1
if i!=0:
print("starting solve..")
print("symmetries: ", symmetries)
with TaskManager():
l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, vlam, ev, step=i, alpha=args.alpha, maxsym = symmetries, maxsymex = maxk)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
print(l2errors)
print("error with no killing field: ")
print(l2errkf)
print("Limits: ", computelimit(phi, radius=args.radius))
if args.plot:
if geom=="ellipsoid" and c==1:
......
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