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

small changes

parent 085c7cec
for ORDER in 1 2 3
for ORDER in 1 2 3 4 5
do
for C in 1 1.1 1.5
for C in 1 1.0005 1.001
do
./killing.py --refinements 5 -c $C --order $ORDER --plot --saveplot
......@@ -8,7 +8,7 @@ do
done
done
for ORDER in 1 2 3
for ORDER in 1 2 3 4 5
do
for C in "1" "0.5" "0.1"
......
......@@ -150,8 +150,9 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
vertexdofs = BitArray(Vh.ndof)
vertexdofs[:] = False
print("doing stuff for precond")
"""
print("doing stuff for precond")
for v in mesh.vertices:
......@@ -166,9 +167,10 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
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(VhG.GetDofNrs(el)) for el in VhG.Elements(VOL)]
# 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]
"""
def blockcr(FEspace):
for v in FEspace.mesh.vectices:
vdofs = set()
......@@ -179,6 +181,7 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
# blocks = pool.map(partial(calcstuff, mesh=mesh), range(len(mesh.vertices)))
vertexdofs &= VhG.FreeDofs()
"""
"""
print("done with stuff for precond")
n = Normalize(grad(lset_approx))
ntilde = GridFunction(VhG)
......@@ -213,24 +216,26 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
a2 += (Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
pre1 = Preconditioner(a2, "direct", inverse="pardiso")
# pre1 = Preconditioner(a2, "bddc", inverse="pardiso")
# print(reduce(lambda x,y:x and y, vertexdofs and VhG.FreeDofs()))
a2.Assemble()
coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
# coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
# pre = Projector(VhG.FreeDofs(), True)
# pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
# pre2 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
# pre1 = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
pre = pre1#+coarsepre
# pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
# pre = pre1+coarsepre#+pre2+coarsepre
#pre=SymmetricGS(pre1)
pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
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,maxit = 15)
return solvers.PINVIT(a2.mat, b.mat, pre, 5, GramSchmidt=True,maxit = 10)
def getevproblem(mesh, phi, lset_approx, deformation, order=2):
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
......
......@@ -275,9 +275,11 @@ if __name__ == "__main__":
lset_approx, deformation = refineMesh(mesh, phi, order)
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
print("eigenvalues:", vlam)
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
for k in range(3):
vlams[k].append(vlam[k])
print("vlamk: ", vlam[k])
if i < 2 and i != 0:
symmetries = 0
for k in range(3):
......@@ -300,7 +302,7 @@ if __name__ == "__main__":
print(math.log(abs(1 - temp), 2))
print("condition: ", -(2*order+1)*(i+1))
print("log10: ", math.log(abs(1 - temp), 10))
if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i+1)+1:
if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i+1)+3:
symmetries += 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