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

removing KF works if dim known

parent 6f6ec326
......@@ -51,7 +51,7 @@ reac_cf = 1
diff_cf = 1
# Ellipsoid parameter
c = 1
if c==1:
if c==1 or c<1.01:
maxk = 3
else:
maxk=1
......@@ -257,7 +257,7 @@ with TaskManager():
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, "bddc")
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
......@@ -267,24 +267,16 @@ with TaskManager():
b.Assemble()
evecs2 = GridFunction(VhG, multidim=3)
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs2.vecs, 1)
for k in range(3):
evecs.vecs[k].data = b.mat*evecs2.vecs[k]
# vlam, evecs = solvers.PINVIT(a2.mat, b.mat,pre=IdentityMatrix(), num=3)
"""
rows,cols,vals = a.mat.COO()
A = sp.csr_matrix((vals,(rows,cols)))
print("computing eigenvectors..")
evecs, evs = eigs(A, k=3, sigma=0, which='SM',tol=10**-6, maxiter=1000)
print("done with eigenvectors")
critevals = [evs<cuthmax**Epsilon-2*cuthmax**2]
"""
# R.h.s. linear form
vlam, evecs = solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
evecinter = []
for k in range(3):
# evn = sqrt(Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(evecs3.MDComponent(k)) ** 2,
# mesh=mesh))
# evecs.vecs[k].data = evecs3.vecs[k].data
evecinter.append(GridFunction(VhG))
evecinter[k].vec.data = evecs[k].data
f = LinearForm(X)
f += rhs1 * v * ds
......@@ -292,11 +284,7 @@ with TaskManager():
f.Assemble()
a.Assemble()
"""f2 = LinearForm(X)
f2 += rhs1ex * v * ds
f2 += -rhs2 * q * ds
gfu2 = GridFunction(X)
f2.Assemble()"""
egfu = GridFunction(VhG)
gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
if True:
......@@ -307,20 +295,20 @@ with TaskManager():
for k, lam in enumerate(vlam):
print("lam: "+ str(abs(lam)))
egfu.vec.data = evecs.vecs[k].data
egfu.vec.data = evecinter[k].vec.data
print("is it eigen?")
# print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat-egfu)**2,mesh=mesh))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if True and k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
print("removing number: "+str(k))
evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecinter[k],evecinter[k]), mesh=mesh, order=5)
print("Evecnorm: ", evecnorm)
integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh, order=5)/evecnorm
integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecinter[k]), mesh=mesh, order=5)/evecnorm
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh)
gfu.components[0].vec.data+= integral*evecs.vecs[k]
gfu.components[0].vec.data+= integral*evecinter[k].vec
# print(gfu.components[0].vec.data-kf.vec.data)
kf2 = GridFunction(VhG)
......
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