test not removing kvf

statstokes.py

@@ -46,9 +46,10 @@ order = 2
 reac_cf = 1
 diff_cf = 1
 # Ellipsoid parameter
-c = 1.1
+c = 1
 # Geometry
-#SetHeapSize(10**10)
+#SetHeapSize(10**9)
+SetNumThreads(16)
 with TaskManager():
     cube = CSGeometry()
     if geom == "ellipsoid":
@@ -103,7 +104,7 @@ with TaskManager():
     # Coefficients / parameters:
 
     # Exact tangential projection
-    Ps = Pmats(Normalize(grad(phi)))
+    Ps = Pmats(Normalize(grad(phi))).Compile(realcompile=True)
 
     # define solution and right-hand side
 
@@ -119,7 +120,7 @@ with TaskManager():
         }
         pSol = (x * y ** 3 + z * (x ** 2 + y ** 2 + z ** 2) ** (3 / 2)) / ((x ** 2 + y ** 2 + z ** 2) ** 2)
         uSol = Ps* CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))"""
-        uSol = Ps* CF((-z**2,x,y))
+        uSol = (Ps* CF((-z**2,x,y))).Compile(realcompile=True)
         pSol = CF(x*y**3+z)
 
     elif geom == "decocube":
@@ -221,7 +222,7 @@ with TaskManager():
 
 
         def E(u):
-            return Pmat * Sym(grad(u)) * Pmat# - (u * n) * weing
+            return (Pmat * Sym(grad(u)) * Pmat)#.Compile(realcompile=comp)# - (u * n) * weing
 
 
         a = BilinearForm(X, symmetric=True)
@@ -283,7 +284,7 @@ with TaskManager():
         for k, lam in enumerate(vlam):
             print("lam: "+ str(abs(lam)))
             print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
-            if k==0:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
+            if False and k==0:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
                 kf.vec.data =evecs.vecs[0].data
                 evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
                 print("removing number: "+str(k))
@@ -294,10 +295,11 @@ with TaskManager():
                 gfu.components[0].vec.data+=  integral*evecs.vecs[k]
 #                print(gfu.components[0].vec.data-kf.vec.data)
                 print("scalar product with killing field: ")
-
+        """
         for k in range(3):
             print(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF},
                             cf = InnerProduct(gfu.components[0], evecs.MDComponent(k)), mesh=mesh, order=10))
+        """
         #gfu2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso")*f2.vec
         l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
         uerr = gfu.components[0] - uSol