Commit 5dbd127b authored by Tilman Alemán's avatar Tilman Alemán
Browse files

more or less working

parent 8fc38709
for ORDER in 1 2 3 4 5
for ORDER in 2 3 4
do
for C in 1 1.0005 1.001
......@@ -7,13 +7,13 @@ do
done
done
#
#for ORDER in 2 3 4
#do
#for C in "1" "0.5" "0.1"
for ORDER in 1 2 3 4 5
do
for C in "1" "0.5" "0.1"
do
./killing.py --geometry arrow --refinements 5 -c $C --order $ORDER --plot --saveplot
#do
#./killing.py --geometry arrow --refinements 5 -c $C --order $ORDER --plot --saveplot
done
done
\ No newline at end of file
#done
#done
\ No newline at end of file
......@@ -128,8 +128,9 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
with TaskManager():
#condense=False
condense=True
print("starting ev comp")
vlams = []
if order<3:
precond = "bddc"
......@@ -141,17 +142,9 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
"""
n = VhG.ndof
markdof = [0]*n
for el in VhG.Elements():
for d in VhG.GetDofNrs(el):
markdof[d]+=1
#dofs = VhG.GetDofs()
for k in range(VhG.ndof):
if markdof[k]<2:
VhG.SetCouplingType(k, COUPLING_TYPE.LOCAL_DOF)
"""
# VhG.SetCouplingType(IntRange(localdofs),COUPLING_TYPE.LOCAL_DOF)
blocks = []
freedofs = VhG.FreeDofs()
......@@ -173,9 +166,10 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
vertexdofs = BitArray(VhG.ndof)
vertexdofs[:] = False
"""
print(len(vertexdofs))
print(VhG.ndof)
print("doing stuff for precond")
"""
for v in mesh.vertices:
......@@ -187,14 +181,28 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
# vdofs = set(d for d in VhG.GetDofNrs(el) if freedofs[d])
for d in VhG.GetDofNrs(v):
vertexdofs[d] = 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 = [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
markdof = [0]*n
"""
for el in VhG.Elements():
for d in VhG.GetDofNrs(el):
markdof[d]+=1
#dofs = VhG.GetDofs()
for k in range(VhG.ndof):
if markdof[k]<2:
VhG.SetCouplingType(k, COUPLING_TYPE.LOCAL_DOF)
"""
def blockcr(FEspace):
for v in FEspace.mesh.vectices:
vdofs = set()
......@@ -236,89 +244,50 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
return ((Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing)#.Compile()
a2 = BilinearForm(VhG, symmetric=True, condense=condense, store_inner=condense)
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
a2 += InnerProduct(E(u2), E(v2)).Compile() * ds
# a2 += (Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)).Compile() * dX
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))).Compile() * ds
pre1 = Preconditioner(a2, precond, inverse="sparsecholesky")
pre1 = Preconditioner(a2, precond, inverse="pardiso")
# prejac = Preconditioner(a2, "local", inverse="pardiso")
# coarsepre = a2.mat.Inverse(vertexdofs, 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
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 += InnerProduct(Pmat * u2, Pmat * v2).Compile() * ds
b += (rhob * (grad(u2) * n) * (grad(v2) * n)).Compile() * dX
b += ((u2 * n) * (v2 * n)).Compile() * ds
b.Assemble()
start = time.time()
a2.Assemble()
# pre=a2.mat.Inverse(freedofs=VhG.FreeDofs(True), inverse="pardiso") + coarsepre
print("assemble time: ", time.time()-start)
if condense:
return *PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 10), start
ev, evecs = PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 5)
else:
return *solvers.PINVIT(a2.mat, b.mat, pre, 4, GramSchmidt=True, maxit=10), start
def getevproblem(mesh, phi, lset_approx, deformation, order=2):
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, "local")
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
return a2, b, pre, Vh, VhG
ev, evecs =solvers.PINVIT(a2.mat, b.mat, pre, 4, GramSchmidt=True, maxit=5)
return ev, evecs, VhG, a2, pre, ci, start
def PINVIT2(a, matm, pre,fes, num=1, maxit=20, printrates=True, GramSchmidt=False):
"""preconditioned inverse iteration"""
mata = (IdentityMatrix()-a.harmonic_extension_trans)@(a.mat + a.inner_matrix)@(IdentityMatrix()-a.harmonic_extension_trans)
matainv = (IdentityMatrix()+a.harmonic_extension)@(pre + a.inner_solve)@(IdentityMatrix()-a.harmonic_extension)
mata = (IdentityMatrix()-a.harmonic_extension_trans)@(a.mat + a.inner_matrix)@(IdentityMatrix()-a.harmonic_extension)
matainv = (IdentityMatrix()+a.harmonic_extension)@(pre + a.inner_solve)@(IdentityMatrix()-a.harmonic_extension_trans)
r = mata.CreateRowVector()
......@@ -353,3 +322,5 @@ def PINVIT2(a, matm, pre,fes, num=1, maxit=20, printrates=True, GramSchmidt=Fals
uvecs[:] = vecs * Matrix(evec[:,0:num])
return lams, uvecs
......@@ -9,7 +9,7 @@ from helperfunctions import *
# -------------------------------- PARAMETERS ---------------------------------
def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsym = 0, maxsymex = 0, condense=False):
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)))
......@@ -23,9 +23,9 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
X = VhG
......@@ -80,12 +80,13 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
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
a += epsilon * (Pmat * u) * (Pmat * v) * ds
a += (eta * ((u * ntilde) * (v * ntilde))) * ds
a += rhou * ((grad(u) * n) * (grad(v) * n)) * dX
"""
L21 = L2(mesh, order=order + 1, dim=9)
L22 = L2(mesh, order=order + 1, dim=3)
lpmat = Restrict(L21, ba_IF)
......@@ -110,10 +111,10 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
f.Assemble()
pre=Preconditioner(a,"direct", inverse="pardiso")
a.Assemble()
# pre=Preconditioner(a,"direct", inverse="pardiso")
# a.Assemble()
if condense:
solvers.BVP(bf=a, lf=f,gf=gfu,pre=pre)
solvers.BVP(bf=a, lf=f,gf=gfu,pre=pre, needsassembling=False, inverse="pardiso", )
else:
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
kf = GridFunction(VhG)
......@@ -124,15 +125,14 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
# gfu2.components[0].vec.data += integral*kf2.vec
kf.Set(uSol, definedonelements=ba_IF)
kf2.Set(kvfs[k], definedonelements=ba_IF)
for k in range(0):
kf.Set(kvfs[k], definedonelements=ba_IF)
kf.Set(uSol-gfu, definedonelements=ba_IF)
kf2 = evecinter
for k in range(maxsym):
#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
cf=InnerProduct(kf2[k], kf), mesh=mesh)
kf.vec.data += integral * kf2[k].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
......@@ -147,7 +147,7 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
l2err = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(gfu - uSol) ** 2, mesh=mesh))
l2errnokf = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=InnerProduct(kf2, kf2), mesh=mesh))
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=InnerProduct(kf, kf), mesh=mesh))
uerr = gfu - uSol
return l2err, l2errnokf, uerr, gfu, uSol
......@@ -170,10 +170,10 @@ if __name__ == "__main__":
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=1)
parser.add_argument("--log", dest="loglevel", type=str, required=False, default="DEBUG")
parser.add_argument("--log", dest="loglevel", type=str, required=False, default="NOTSET")
args = parser.parse_args()
numeric_level = getattr(lg, args.loglevel.upper())
lg.basicConfig(level=numeric_level)
#lg.basicConfig(level=numeric_level)
interpol = True
# Mesh diameter
maxh = 0.6*args.radius
......@@ -280,7 +280,7 @@ if __name__ == "__main__":
if i != 0:
lset_approx, deformation = refineMesh(mesh, phi, order)
vlam, ev, start = computeEV(phi,mesh,lset_approx, deformation, order, condense=condense)
vlam, ev, VhG, a,pre, ci,start = computeEV(phi,mesh,lset_approx, deformation, order, condense=condense)
print("time: ", time.time()-start)
print("condense: ", condense)
print("eigenvalues:", vlam)
......@@ -292,9 +292,9 @@ if __name__ == "__main__":
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):
if math.log(abs( vlams[k][1]), 2) <-(order+1)*(i):
print("corresponding log: ")
print(math.log(abs(1 - vlams[k][1]), 2))
print(math.log(abs( vlams[k][1]), 2))
print("condition: ", -(order+1)*(i+1))
symmetries += 1
......@@ -307,10 +307,10 @@ if __name__ == "__main__":
(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(1 - temp), 2))
print(math.log(abs( temp), 2))
print("condition: "+ str(-(2*order+1)*(i+1)))
print("log10: "+ str(math.log(abs(1 - temp), 10)))
if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i+1)+3:
print("log10: "+ str(math.log(abs( temp), 10)))
if False or math.log(abs( temp), 2) < -(2*order+1) * (i+1)+3:
symmetries += 1
......@@ -318,7 +318,7 @@ if __name__ == "__main__":
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, condense=condense)
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)
......@@ -329,7 +329,7 @@ if __name__ == "__main__":
if args.plot:
if geom=="ellipsoid" and c==1:
plot_errors(l2errors, l2errkf, "circle", 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)
else:
plot_errors(l2errors, l2errkf, "ellipsoid", args.saveplot, order=order, refs=n_cut_ref, c=c)
......
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