Commit 82d779d5 authored by Tilman Aleman's avatar Tilman Aleman
Browse files

check precond

parent c9696c5e
......@@ -7,6 +7,9 @@ import scipy
import time
import logging as lg
from functools import reduce
import scipy.sparse as sp
import matplotlib.pyplot as plt
from ngsolve.la import EigenValues_Preconditioner
def refineMesh(mesh, phi, order, hs=10 ** 6):
if order>2:
hs=10**8
......@@ -128,11 +131,11 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1, sca=2):
def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
with TaskManager():
condense=True
#condense=True
print("starting ev comp")
vlams = []
if order<3:
if order<6:
precond = "bddc"
else:
precond = "direct"
......@@ -184,11 +187,11 @@ 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]
#blocks1 = [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]
......@@ -251,21 +254,14 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)).Compile() * dX
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))).Compile() * ds
pre1 = Preconditioner(a2, precond, inverse="pardiso")
pre1 = Preconditioner(a2, precond, inverse="pardiso", coarsetype="h1amg")
#pre1 = Preconditioner(a2, "bddc", inverse="pardiso")
# prejac = Preconditioner(a2, "local", inverse="pardiso")
# pre = Projector(VhG.FreeDofs(), True)
#
if not condense:
pre = pre1#+prejac
else:
pre=pre1
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat * u2, Pmat * v2).Compile() * ds
b += (rhob * (grad(u2) * n) * (grad(v2) * n)).Compile() * dX
......@@ -273,20 +269,46 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2, condense=True):
b.Assemble()
start = time.time()
a2.Assemble()
#vblocks = [VhG.GetDofNrs(vertex) for el in VhG.Elements() for vertex in el.vertices ]
#eblocks = [VhG.GetDofNrs(edge) for edge in VhG.mesh.edges]
#fblocks = [VhG.GetDofNrs(face) for face in VhG.mesh.faces]
#vinv = a2.mat.CreateSmoother()
#einv = a2.mat.CreateBlockSmoother(eblocks)
#finv = a2.mat.CreateBlockSmoother(fblocks)
print("smoother")
#prejac = vinv#finv@einv@vinv
if not condense:
pre = pre1#+prejac
#pre=prejac
else:
pre=pre1#+prejac
#pre=prejac
rows,cols,vals = a2.mat.COO()
A1 = sp.csr_matrix((vals,(rows,cols)))
plt.spy(A1)
#plt.show()
#coarsepre = Preconditioner(a2, "h1amg", inverse="pardiso")
#coarsepre = a2.mat.Inverse(freedofs=vertexdofs, inverse="pardiso")
lams = EigenValues_Preconditioner(a2.mat, pre1.mat)#+coarsepre)
print("Preconditioned EVs: ", lams[-1]/lams[0])
#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 = 10)
ev, evecs = PINVIT2(a2, b.mat, pre, VhG, 4, GramSchmidt=True,maxit = 1)
else:
ev, evecs =solvers.PINVIT(a2.mat, b.mat, pre, 4, GramSchmidt=True, maxit=10)
ev, evecs =solvers.LOBPCG(a2.mat, b.mat, pre, 4, maxit=40)
return ev, evecs, VhG, a2, pre, ci, start
def PINVIT2(a, matm, pre,fes, num=1, maxit=20, printrates=True, GramSchmidt=False):
......
......@@ -6,6 +6,7 @@ from netgen.csg import CSGeometry, OrthoBrick, Pnt
from helperfunctions import *
import csv
from time import perf_counter
# -------------------------------- PARAMETERS ---------------------------------
......@@ -101,10 +102,17 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, VhG, a,pre, ci, ste
f.Assemble()
# pre=Preconditioner(a,"direct", inverse="pardiso")
# a.Assemble()
if condense:
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
if True:
start = perf_counter()
solvers.GMRes(A=a.mat, b=f.vec,x=gfu.vec,pre=pre,maxsteps=10000, tol=1e-7)
#solvers.CG(mat=a.mat, rhs=f.vec, pre=pre,sol=gfu.vec, maxsteps=10000, tol=1e-7)
print("Iterative: " + str(perf_counter()-start))
#else:
start=perf_counter()
solvers.BVP(a, f, gfu, inverse="pardiso", needsassembling=False)
#a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
print("Direct: " + str(perf_counter()-start))
kf = GridFunction(VhG)
print(vlam)
......
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