Skip to content
GitLab
Commit f6b10e15 authored by Tilman Alemán's avatar Tilman Alemán
Browse files

csv output

parent 5dbd127b
Hide whitespace changes
Inline Side-by-side
helper.sh
View file @ f6b10e15
for ORDER in 2 3 4
for SCA in 2
do
for C in 1 1.0005 1.001
for ORDER in 2 3
do
for C in 1 1.1 1.001
do
./killing.py --refinements 5 -c $C --order $ORDER --plot --saveplot
./killing.py --refinements 5 -c $C --order $ORDER --plot --saveplot -l $SCA
done
done
for ORDER in 4
do
for C in 1 1.1 1.001
do
./killing.py --refinements 4 -c $C --order $ORDER --plot --saveplot -l $SCA
done
done
done
#
......
helperfunctions.py
View file @ f6b10e15
......@@ -98,7 +98,7 @@ def removekf(VhG, f, kvs, vlams, phi, mesh, maxk=3, tol=0, autofilter=True, inpl
sca = l2sca(f.components[0], kvs[k], phi, mesh)
f.components[0].vec.data += -sca * kvs[k].vec
return
def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1, sca=2):
hvals = []
......@@ -111,7 +111,7 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
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, l2errskf, label=geom + " modulo discrete KF")
plt.semilogy(ticks,hvals, label=r'$O(h^'+ str(order+1)+r')$')
plt.semilogy(ticks,hvals2, label=r'$O(h^'+ str(order)+r')$')
......@@ -122,7 +122,7 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
plt.legend()
plt.xticks(ticks)
if saveplot:
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+"_sca"+str(sca)+".png")
else:
plt.show()
......@@ -137,7 +137,7 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
else:
precond = "direct"
start=time.time()
Vh = VectorH1(mesh, order=order, dirichlet=[], low_order_space=True)
Vh = VectorH1(mesh, order=order, dirichlet=[], low_order_space=True, wb_fulledges=True)
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
......@@ -149,7 +149,7 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
blocks = []
freedofs = VhG.FreeDofs()
class SymmetricGS(BaseMatrix):
def __init__ (self, smoother, inv):
def __init__ (self, smoother):
super(SymmetricGS, self).__init__()
self.k=0
self.smoother = smoother
......@@ -184,12 +184,13 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
if d>=0:
vertexdofs[d] = True
"""
# blocks.append(vdofs)
"""
blocks = [set(d for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
"""
#blocks = [set(d for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
"""
# blocks = [reduce(lambda x,y: x or y,[set(d for d in VhG.GetDofNrs(el) if freedofs[d])for el in mesh[v].elements]) for v in mesh.vertices]
"""
n = VhG.ndof
......@@ -211,7 +212,7 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
if freedofs[d])
# blocks = pool.map(partial(calcstuff, mesh=mesh), range(len(mesh.vertices)))
vertexdofs &= VhG.FreeDofs()
#vertexdofs &= VhG.FreeDofs()
"""
"""
print("done with stuff for precond")
......@@ -250,20 +251,19 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))).Compile() * ds
pre1 = Preconditioner(a2, precond, inverse="pardiso")
#pre1 = Preconditioner(a2, "bddc", inverse="pardiso")
# prejac = Preconditioner(a2, "local", inverse="pardiso")
# pre = Projector(VhG.FreeDofs(), True)
# pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
# inv = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="sparsecholesky")
# prejac = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True), inv)
#
if not condense:
pre = pre1#+prejac
else:
pre=pre1
pre=pre1
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat * u2, Pmat * v2).Compile() * ds
......@@ -272,15 +272,20 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
b.Assemble()
start = time.time()
a2.Assemble()
#coarsepre = a2.mat.Inverse(freedofs=vertexdofs, inverse="pardiso")
#prejac = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
#prejac = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
#coarsepre = Preconditioner(a2, "multigrid", inverse="pardiso")
#pre = pre + coarsepre
# pre=a2.mat.Inverse(freedofs=VhG.FreeDofs(True), inverse="pardiso") + coarsepre
print("assemble time: ", time.time()-start)
if condense:
ev, evecs = PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 5)
ev, evecs = PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 10)
else:
ev, evecs =solvers.PINVIT(a2.mat, b.mat, pre, 4, GramSchmidt=True, maxit=5)
ev, evecs =solvers.PINVIT(a2.mat, b.mat, pre, 4, GramSchmidt=True, maxit=10)
return ev, evecs, VhG, a2, pre, ci, start
def PINVIT2(a, matm, pre,fes, num=1, maxit=20, printrates=True, GramSchmidt=False):
......
killing.py
View file @ f6b10e15
......@@ -5,22 +5,22 @@ import argparse
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from helperfunctions import *
import csv
# -------------------------------- PARAMETERS ---------------------------------
def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, step,alpha=3, maxsym = 0, maxsymex = 0, condense=False):
uSol = (Ps * CF((-z ** 2, x, y)))
uSol2 = (Ps * CF((-z ** 2, x, y)))
for k in range(maxsymex):
print("removing something?")
print("Removing ex. KF")
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]
Vh = VectorH1(mesh, order=order, dirichlet=[])
......@@ -32,10 +32,10 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
gfu = GridFunction(X)
n = Normalize(grad(lset_approx))
ntilde = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# ntilde = GridFunction(VhG)
# 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])
# epsilon =(cuthmax)**(alpha)
......@@ -43,17 +43,7 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
# print(len(elvol))
# print("hmax:")
# print(cuthmax)
alpha=order+1
epsilon = h ** alpha
if True:
epsilon=0
print("Epsilon:")
print(epsilon)
# epsilon=1
eta = 1 / (h) ** 2
rhou = 1.0 / h
rhop = h
rhob = h
# Measure on surface
ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
# Measure on the bulk around the surface
......@@ -64,9 +54,7 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
......@@ -76,10 +64,7 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
comp = False
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
"""
a = BilinearForm(X, symmetric=True, condense=condense)
a += InnerProduct(E(u), E(v)).Compile(realcompile=comp) * ds
......@@ -114,7 +99,7 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
# pre=Preconditioner(a,"direct", inverse="pardiso")
# a.Assemble()
if condense:
solvers.BVP(bf=a, lf=f,gf=gfu,pre=pre, needsassembling=False, inverse="pardiso", )
solvers.BVP(bf=a, lf=f,gf=gfu,pre=pre, needsassembling=False, inverse="pardiso",maxsteps=10000)
else:
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
kf = GridFunction(VhG)
......@@ -171,7 +156,15 @@ if __name__ == "__main__":
parser.add_argument("--alpha", dest="alpha", type=float, required=False, default=3)
parser.add_argument("--radius", dest="radius", type=float, required=False, default=1)
parser.add_argument("--log", dest="loglevel", type=str, required=False, default="NOTSET")
parser.add_argument("-l", dest="sca", type=float, required=False, default=2)
args = parser.parse_args()
fname = "data_order"+str(args.order)+"_refs"+str(args.n_cut_ref)+"_c"+str(args.c)+"_symsca"+str(args.sca)
with open(fname, 'w') as f:
f.write("Order: "+str(args.order)+", ")
with open(fname, 'a') as f:
f.write("Refinemets: " + str(args.n_cut_ref)+", ")
f.write("Geometry param: "+ str(args.c) + "\n")
f.write("Scaling for symmetry: "+ str(args.sca))
numeric_level = getattr(lg, args.loglevel.upper())
#lg.basicConfig(level=numeric_level)
interpol = True
......@@ -215,9 +208,6 @@ if __name__ == "__main__":
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
# Class to compute the mesh transformation needed for higher order accuracy
......@@ -270,6 +260,14 @@ if __name__ == "__main__":
vlams = []
symmetries = 0
condense=True
fnameeig = fname + "eig.csv"
header = ["Refinement"]
for k in range(1,5):
header += ["\\lambda_"+str(k)+"(h)"]
with open(fnameeig, "w", newline='') as f:
writer = csv.DictWriter(f, delimiter="&", fieldnames = header)
writer.writeheader()
for k in range(3):
vlams.append([])
for i in range(n_cut_ref + 1):
......@@ -284,6 +282,16 @@ if __name__ == "__main__":
print("time: ", time.time()-start)
print("condense: ", condense)
print("eigenvalues:", vlam)
with open(fnameeig, "a", newline='') as f:
values = [i]+[lam for lam in vlam]
dictionary = dict(zip(header, values))
writer = csv.DictWriter(f, delimiter=",", fieldnames = header)
writer.writerow(dictionary)
with open(fname, "a") as f:
f.write("eigenvalues for refinement "+ str(i)+":\n")
f.write(str(vlam) + "\n")
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
for k in range(3):
vlams[k].append(vlam[k])
......@@ -300,20 +308,24 @@ if __name__ == "__main__":
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):
print("corresponding log: ")
print(math.log(abs( temp), 2))
print("condition: "+ str(-(2*order+1)*(i+1)))
print("log10: "+ str(math.log(abs( temp), 10)))
if False or math.log(abs( temp), 2) < -(2*order+1) * (i+1)+3:
with open(fname, "a") as f:
f.write("Condition for discrete KF: "+ str(-(args.sca*order+1) * (i+1)+3)+":\n")
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):
print("corresponding log: ")
print(math.log(abs( temp), 2))
print("condition: "+ str(-(args.sca*order+1)*(i+1)))
print("log10: "+ str(math.log(abs( temp), 10)))
f.write("Actual value: "+str(math.log(abs( temp), 2)) + "\n")
if False or math.log(abs( temp), 2) < -(args.sca*order+1) * (i+1)+3:
symmetries += 1
symmetries += 1
f.write("Symmetries: " + str(symmetries)+"\n")
if True and i!=0:
print("starting solve..")
print("symmetries: ", symmetries)
......@@ -321,6 +333,17 @@ if __name__ == "__main__":
l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, vlam, ev, VhG, a, pre, ci, step=i, alpha=args.alpha, maxsym = symmetries, maxsymex = maxk, condense=condense)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
with open(fname, "a") as f:
f.write("l2 errors: \n")
f.write(str(l2errors)+ "\n")
f.write("l2 errors for LaTeX: \n")
f.write(" & ".join([str(l2) for l2 in l2errors]) + "\n")
f.write("l2 errors without discrete KF: \n")
f.write(str(l2errkf)+ "\n")
f.write("l2 errors without discrete KF for LaTeX: \n")
f.write(" & ".join([str(l2) for l2 in l2errkf])+ "\n")
print(l2errors)
print("error with no killing field: ")
......@@ -329,10 +352,10 @@ if __name__ == "__main__":
if args.plot:
if geom=="ellipsoid" and c==1:
plot_errors(l2errors, l2errkf, "sphere", args.saveplot, order=order, refs=n_cut_ref, c=c)
plot_errors(l2errors, l2errkf, "sphere", args.saveplot, order=order, refs=n_cut_ref, c=c, sca=args.sca)
else:
plot_errors(l2errors, l2errkf, "ellipsoid", args.saveplot, order=order, refs=n_cut_ref, c=c)
plot_errors(l2errors, l2errkf, "ellipsoid", args.saveplot, order=order, refs=n_cut_ref, c=c, sca=args.sca)
if vtkout: # and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
# input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
......
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