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

preconditioners

parent b72cfbf8
......@@ -3,21 +3,27 @@ from xfem.lsetcurv import *
import matplotlib.pyplot as plt
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import math
from functools import reduce
def refineMesh(mesh, phi, order, hs=10 ** 6):
if order>2:
hs=10**8
with TaskManager():
mesh.SetDeformation(None)
lsetp1 = GridFunction(H1(mesh, order=order))
print("lsetp1 interpolation..")
InterpolateToP1(phi, lsetp1)
print("lsetp1 refinement..")
RefineAtLevelSet(lsetp1)
mesh.Refine()
print("lsetmeshadap...")
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
print("calc deform...")
deformation = lsetmeshadap.CalcDeformation(phi)
mesh.SetDeformation(deformation)
print("set deform..")
if order>1:
mesh.SetDeformation(deformation)
lset_approx = lsetmeshadap.lset_p1
return lset_approx, deformation
......@@ -116,12 +122,117 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
else:
plt.show()
def computeEV(phi, mesh, lset_approx, deformation, order=2):
vlams = []
def computeEV(phi, mesh, lset_approx, deformation, order=2):
with TaskManager():
print("starting ev comp")
vlams = []
Vh = VectorH1(mesh, order=order, dirichlet=[], low_order_space=True)
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
blocks = []
freedofs = VhG.FreeDofs()
class SymmetricGS(BaseMatrix):
def __init__ (self, smoother):
super(SymmetricGS, self).__init__()
self.smoother = smoother
def Mult (self, x, y):
y[:] = 0.0
self.smoother.Smooth(y, x)
self.smoother.SmoothBack(y,x)
def Height (self):
return self.smoother.height
def Width (self):
return self.smoother.height
vertexdofs = BitArray(Vh.ndof)
vertexdofs[:] = False
print("doing stuff for precond")
"""
for v in mesh.vertices:
# vdofs = set(d for d in VhG.GetDofNrs(el) if freedofs[d])
for d in Vh.GetDofNrs(v):
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]
def blockcr(FEspace):
for v in FEspace.mesh.vectices:
vdofs = set()
for el in FEspace.mesh[v].elements:
vdofs |= set(d for d in FEspace.GetDofNrs(el)
if freedofs[d])
# 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)
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
pre1 = Preconditioner(a2, "direct", inverse="pardiso")
a2.Assemble()
coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
# pre = Projector(VhG.FreeDofs(), True)
# pre1 = 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")
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)
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)
......@@ -158,13 +269,18 @@ 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
pre = Preconditioner(a2, "direct")
a2.Assemble()
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
def computeEV2(Vh, VhG, a2, b, pre):
Vh.Update()
VhG.Update()
print("updated..")
a2.Assemble()
b.Assemble()
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True,maxit = 20)
\ No newline at end of file
......@@ -263,18 +263,22 @@ if __name__ == "__main__":
for k in range(3):
vlams.append([])
for i in range(n_cut_ref + 1):
print("refinement is: ", i)
if i==0:
# a,b,pre, Vh, VhG = getevproblem(mesh, phi, lset_approx, deformation, order)
pass
if i != 0:
lset_approx, deformation = refineMesh(mesh, phi, order)
with TaskManager():
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
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):
if math.log(abs(1 - vlams[k][1]), 2) <-(order+1)*(i):
print("corresponding log: ")
print(math.log(abs(1 - vlams[k][1]), 2))
print("condition: ", -(order+1)*(i+1))
......@@ -287,12 +291,13 @@ if __name__ == "__main__":
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))
print("corresponding log: ")
print(math.log(abs(1 - temp), 2))
print("condition: ", -(2*order+1)*(i))
if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i):
symmetries += 1
if i!=0:
if True and i!=0:
print("starting solve..")
print("symmetries: ", symmetries)
with TaskManager():
......
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