Commit a26642fc authored by tilman's avatar tilman
Browse files

working static condensation

parent 1d749739
......@@ -3,6 +3,9 @@ from xfem.lsetcurv import *
import matplotlib.pyplot as plt
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import math
import scipy
import time
import logging as lg
from functools import reduce
def refineMesh(mesh, phi, order, hs=10 ** 6):
if order>2:
......@@ -123,16 +126,33 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
else:
plt.show()
def computeEV(phi, mesh, lset_approx, deformation, order=2):
def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
with TaskManager():
#condense=False
print("starting ev comp")
vlams = []
if order<3:
precond = "bddc"
else:
precond = "direct"
start=time.time()
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)
"""
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()
class SymmetricGS(BaseMatrix):
......@@ -171,12 +191,8 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
# blocks.append(vdofs)
<<<<<<< HEAD
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)]
>>>>>>> ea4519effe74102893ea15789e6c4633f4f61bbf
# 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):
......@@ -218,15 +234,17 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
def E(u):
return ((Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing)#.Compile()
a2 = BilinearForm(VhG, symmetric=True)
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
pre1 = Preconditioner(a2, "direct", inverse="sparsecholesky")
pre1 = Preconditioner(a2, precond, inverse="sparsecholesky")
# prejac = Preconditioner(a2, "local", inverse="pardiso")
a2.Assemble()
# coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
......@@ -234,15 +252,22 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
# 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)
pre = pre1#+prejac
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.Assemble()
a2.Assemble()
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True,maxit = 20)
if condense:
return *PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 10), start
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)
......@@ -289,10 +314,42 @@ def getevproblem(mesh, phi, lset_approx, deformation, order=2):
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,maxit = 20)
\ No newline at end of file
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)
r = mata.CreateRowVector()
uvecs = MultiVector(r, num)
vecs = MultiVector(r, 2*num)
# hv = MultiVector(r, 2*num)
for v in vecs[0:num]:
v.SetRandom()
uvecs[:] = matainv * vecs[0:num]
print("condensed")
lams = Vector(num * [1])
for i in range(maxit):
vecs[0:num] = mata * uvecs - (matm * uvecs).Scale (lams)
vecs[num:2*num] = matainv * vecs[0:num]
vecs[0:num] = uvecs
vecs.Orthogonalize(matm)
# hv[:] = mata * vecs
# asmall = InnerProduct (vecs, hv)
# hv[:] = matm * vecs
# msmall = InnerProduct (vecs, hv)
asmall = InnerProduct (vecs, mata * vecs)
msmall = InnerProduct (vecs, matm * vecs)
ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)
lams = Vector(ev[0:num])
if printrates:
print (i, ":", list(lams))
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):
def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsym = 0, maxsymex = 0, condense=False):
uSol = (Ps * CF((-z ** 2, x, y)))
uSol2 = (Ps * CF((-z ** 2, x, y)))
......@@ -81,7 +81,7 @@ 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=False)
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
......@@ -110,9 +110,12 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
f.Assemble()
pre=Preconditioner(a,"direct", inverse="pardiso")
a.Assemble()
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
if condense:
solvers.BVP(bf=a, lf=f,gf=gfu,pre=pre)
else:
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
kf = GridFunction(VhG)
print(vlam)
......@@ -167,8 +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")
args = parser.parse_args()
numeric_level = getattr(lg, args.loglevel.upper())
lg.basicConfig(level=numeric_level)
interpol = True
# Mesh diameter
maxh = 0.6*args.radius
......@@ -264,6 +269,7 @@ if __name__ == "__main__":
l2errkf = []
vlams = []
symmetries = 0
condense=True
for k in range(3):
vlams.append([])
for i in range(n_cut_ref + 1):
......@@ -274,7 +280,9 @@ if __name__ == "__main__":
if i != 0:
lset_approx, deformation = refineMesh(mesh, phi, order)
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
vlam, ev, start = computeEV(phi,mesh,lset_approx, deformation, order, condense=condense)
print("time: ", time.time()-start)
print("condense: ", condense)
print("eigenvalues:", vlam)
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
for k in range(3):
......@@ -300,8 +308,8 @@ if __name__ == "__main__":
# 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("condition: ", -(2*order+1)*(i+1))
print("log10: ", math.log(abs(1 - temp), 10))
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:
symmetries += 1
......@@ -310,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)
l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, vlam, ev, step=i, alpha=args.alpha, maxsym = symmetries, maxsymex = maxk, condense=condense)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
......
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