Commit 54411e7a authored by tilman's avatar tilman
Browse files

Merge remote-tracking branch 'origin/main'

# Conflicts:
#	killing.py
parents a04dcd64 6f6ec326
......@@ -25,7 +25,7 @@ hs=10**8
reac_cf = 1
diff_cf = 1
# Ellipsoid parameter
c = 1
c = 1.5
if c==1:
maxk = 3
else:
......@@ -64,6 +64,7 @@ with TaskManager():
mesh.SetDeformation(deformation)
# Background FESpaces
......@@ -105,17 +106,14 @@ with TaskManager():
ntilde.Set(Normalize(grad(phi)), definedonelements=ba_IF)
h = specialcf.mesh_size
elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
cuthmax = max([(2*vol)**(1/2) for vol in elvol])
print("hmax:")
print(cuthmax)
alpha=3
# epsilon =(cuthmax)**(alpha)
epsilon=h**alpha
# print("Epsilon:")
# print(epsilon)
# epsilon=1
eta = 1.0 / (h)**2
eta = 1000.0 / (h)**2
rhou = 1.0 / h
rhop = h
......@@ -128,7 +126,7 @@ with TaskManager():
ninter = GridFunction(VhG)
ninter.Set(Normalize(grad(phi)), definedonelements=ba_IF)
ninter.Set(n, definedonelements=ba_IF)
weing = ngsolve.grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
......@@ -140,7 +138,6 @@ with TaskManager():
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat)- (u * n) * weing
......@@ -155,36 +152,104 @@ with TaskManager():
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
a2.Assemble()
# a2inv = a2.Inverse(dofs=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()
def bh(u,v):
ret = InnerProduct(Pmat*u, Pmat*v)
ret += rhob*(Grad(u)*n)*(Grad(v)*n)
ret += ((u * n) * (v * n))
return ret
def l2sca(u,v):
return Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(u,v), mesh=mesh)
def l2norm(u):
return l2sca(u,u)
def orthproj(f, kfs):
res = GridFunction(VhG)
res.vec.data = f.vec.data
for k in range(maxk):
sca = l2sca(f,kfs[k])#/l2norm(kfs[k])
res.vec.data += -sca*kfs[k].vec
return res
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 1)
rows, cols, vals = a2.mat.COO()
A = sp.csr_matrix((vals, (rows, cols)))
rows2, cols2, vals2 = b.mat.COO()
B = sp.csr_matrix((vals2, (rows2, cols2)))
spevs ,spevecs = eigs(A,3,M=B,sigma=1,which='LM')
print("sp evals: ", spevs)
funs = []
print(len(spevecs[:,0]))
evecs2 = GridFunction(VhG, multidim=3)
evecs3 = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs2.vecs, 1)
evecinter = []
for k in range(3):
fun = GridFunction(VhG)
evc = spevecs[:,k]
fun.vec.FV().NumPy()[:] = evc
funs.append(fun)
print(len(fun.vec.FV().NumPy()[:]))
print("what does scipy say:")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(evecs.MDComponent(0)-fun)**2, mesh=mesh))
kf = GridFunction(VhG)
# kf.vec.data = spevecs
evecs3.vecs[k].data =b.mat* evecs2.vecs[k]
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].Set(evecs.MDComponent(k), definedonelements=ba_IF)
print("evecsgrad: ")
# print(l2norm(Grad(evecinter[k])))
print("IP: ", InnerProduct(evecs.vecs[1],evecs2.vecs[0]))
if False:
sizeA = a2.mat.height
print("sizeA: ", sizeA)
atemp = a2.mat.CreateMatrix()
atemp.AsVector().data = a2.mat.AsVector() -b.mat.AsVector()
# atemp.AsVector().data = m.mat.AsVector()
ainv = atemp.Inverse(VhG.FreeDofs(), inverse="pardiso")
x1 = ninter.vec.CreateVector()
y1 = ninter.vec.CreateVector()
def A_inv(v):
for j in range(sizeA):
y1[j] = v[j]
x1.data = ainv*y1
for j in range(sizeA):
v[j] = x1[j]
return v
Ainv = sp.linalg.LinearOperator( (sizeA,sizeA), A_inv)
rowsA,colsA,valsA = a2.mat.COO()
Amattemp = sp.csr_matrix((valsA,(rowsA,colsA)))
Amat = Amattemp + Amattemp.transpose() - sp.diags(Amattemp.diagonal(),0)
rowsM,colsM,valsM = b.mat.COO()
Mmattemp = sp.csr_matrix((valsM,(rowsM,colsM)))
Mmat = Mmattemp + Mmattemp.transpose() - sp.diags(Mmattemp.diagonal(),0)
vals, vecs = sp.linalg.eigsh(Amat, k=3, M=Mmat, sigma=1, which='LM', OPinv=Ainv)
ui = []
for k in range(maxk):
ui.append(GridFunction(VhG))
ui[k].vec.FV().NumPy()[:] = vecs[:,k]
evecs.vecs[k].data = ui[0].vec.data
print("evecsnorms: ")
# print(l2norm(evecs.MDComponent(k)))
egfu = GridFunction(VhG)
egfu.vec.data = evecs.vecs[0].data
def divG(u):
if u.dim == 3:
return Trace(grad(u) * Ps)
if u.dim == 9:
if u.dims[0] == 3:
N = 3
divGi = [divG(u[i, :]) for i in range(N)]
return CF(tuple([divGi[i] for i in range(N)]))
else:
N = 3
r = Grad(u) * Ps
return CF(tuple([Trace(r[N*i:N*(i+1),0:N]) for i in range(N)]))
# print("is killing? ", l2norm(divG(evecs2.MDComponent(0))-evecs2.MDComponent(0)))
fun = GridFunction(VhG)
kf = GridFunction(VhG)
l2errs = []
print(len(vlam))
temp = GridFunction(VhG)
......@@ -198,26 +263,28 @@ with TaskManager():
print("lam: "+ str(abs(lam)))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
# temp.Set( Pmat*evecs.MDComponent(k), definedonelements=ba_IF)
# temp.vec.data = evecs.vecs[k].data
temp.vec.data = funs[k].vec.data
temp.vec.data = evecs.vecs[k].data
for j in range(maxk):
if True:
kvfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kvfs[j],kvfs[j]), mesh=mesh, order=5)
# kvfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(kvfs[j],kvfs[j]), mesh=mesh, order=5)
scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),kvfs[j]), mesh=mesh, order=5)
scal = Integrate(levelset_domain = {"levelset": phi, "domain_type": IF},
cf = InnerProduct(funs[k], kvfs[j]), mesh = mesh, order = 5)
temp.vec.data += -scal*temparr[k].vec
scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs2.MDComponent(k),kvfs[j]), mesh=mesh, order=5)/kvfnorm
# scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(evecs.MDComponent(k),kvfs[j]), mesh=mesh, order=5)/kvfnorm
temp.vec.data += -scal*temparr[j].vec
# l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(temp)**2, mesh=mesh, order=5)))
l2errs.append(
sqrt(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF}, cf = (temp) ** 2, mesh = mesh, order = 5)))
# print("scipy ev problem solution")
# print(sqrt(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF}, cf = (temp) ** 2, mesh = mesh)))
l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(temp)**2, mesh=mesh, order=5)))
for k in range(3):
lam = vlam[k]
print(k)
print("lam: "+ str(abs(lam)))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
print("reverse filtering: ",sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(orthproj(temparr[k],evecinter))**2, mesh=mesh, order=5)))
return l2errs, evecs.MDComponent(0), evecs.MDComponent(1), evecs.MDComponent(2)
l2errors = []
......@@ -241,9 +308,7 @@ with TaskManager():
print("starting solve..")
# l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step = i)
l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
l2errors.append(l2err)
print("L2 errors:")
......
......@@ -173,8 +173,9 @@ with TaskManager():
gfu = GridFunction(X)
n = Normalize(grad(lset_approx))
ntilde = Interpolate(Normalize(grad(phi)), VhG)
ntilde = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi)), definedonelements=ba_IF)
# n=ntilde
h = specialcf.mesh_size
# elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
# cuthmax = max([(2*vol)**(1/2) for vol in elvol])
......@@ -266,9 +267,12 @@ with TaskManager():
b.Assemble()
evecs2 = GridFunction(VhG, multidim=3)
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 1)
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()
......@@ -305,13 +309,13 @@ with TaskManager():
egfu.vec.data = evecs.vecs[k].data
print("is it eigen?")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat)**2,mesh=mesh))
# 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)
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
......
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