started killing field conv studies

statstokes.py

+"""
+In this example we solve a surface stokes problem with a
+consistent penalty method following [1]
+
+Used features:
+--------------
+* Higher order geometry approximation, cf. jupyter tutorial `lsetint`
+
+* Restricted finite element space to condense the system to active dofs,
+  cf. jupyter tutorial `basics`
+
+* Visualization: The visualization of the solution is most convenient
+  with paraview and the generated vtk file.
+
+Literature:
+-----------
+[1] T. Jankuhn, A. Reusken, Higher order Trace Finite Element Methods for the Surface Stokes Equation.
+    arXiv preprint arXiv:1909.08327, 2019.
+
+"""
+
+# ------------------------------ LOAD LIBRARIES -------------------------------
 import sys
 
 from netgen.csg import CSGeometry, OrthoBrick, Pnt
@@ -6,84 +28,59 @@ from ngsolve.internal import *
 from xfem import *
 from xfem.lsetcurv import *
 from math import pi
-import datetime
-begin_time=datetime.datetime.now()
+import scipy.sparse as sp
+from scipy.sparse.linalg import eigs
+import time
+# -------------------------------- PARAMETERS ---------------------------------
+# Interpolation for the rate-of-strain tensor when computing the r.h.s.
+# Introduces an interpolation error, but significantly improves performance
+interpol = True
+# Mesh diameter
+maxh = 0.6
+# Geometry
+geom = "ellipsoid"
+# Polynomial order of FE space
+order = 2
+
+# Problem parameters
+reac_cf = 1
+diff_cf = 1
+# Ellipsoid parameter
+c = 1.5
+# Geometry
+#SetHeapSize(10**10)
 with TaskManager():
-#if True:
-    # -------------------------------- PARAMETERS ---------------------------------
-    # Mesh diameter
-    maxh = 0.6
-    geom = "circle"
-    # Polynomial order of FE space
-    order = 2
-
-    # Problem parameters
-    reac_cf = 1
-    diff_cf = 1
-    c=1.3
-    # Geometry
     cube = CSGeometry()
-    if geom == "circle":
-        phi = Norm(CoefficientFunction((x/c, y, z))) - 1
+    if geom == "ellipsoid":
+        phi = Norm(CoefficientFunction((x , y, z/c))) - 1
 
         cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
         mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
-    elif geom=="decocube":
-        phi=(x**2+y**2-4)**2+(y**2-1)**2+(y**2+z**2-4)**2+(x**2-1)**2+(x**2+z**2-4)**2+(z**2-1)**2-13
-        cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
-        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
-    elif geom=="torus":
-        phi = sqrt(z**2+(sqrt(x**2+y**2)-1)**2)-0.5
+    elif geom == "decocube":
+        phi = (x ** 2 + y ** 2 - 4) ** 2 + (y ** 2 - 1) ** 2 + (y ** 2 + z ** 2 - 4) ** 2 + (x ** 2 - 1) ** 2 + (
+                    x ** 2 + z ** 2 - 4) ** 2 + (z ** 2 - 1) ** 2 - 13
         cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
         mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
-    n_cut_ref = 3
-    for i in range(n_cut_ref):
-        lsetp1 = GridFunction(H1(mesh, order=2))
-        InterpolateToP1(phi, lsetp1)
-        RefineAtLevelSet(lsetp1)
-        mesh.Refine()
-    lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order)
-    deformation = lsetmeshadap.CalcDeformation(phi)
-    lset_approx = lsetmeshadap.lset_p1
 
-    mesh.SetDeformation(deformation)
+    # Preliminary refinements. As assembling the r.h.s. is very expensive without interpolation, keep it small in that case.
+    # When interpolating, getting good results requires more refinements with complex geometries, e.g. the decocube
+    n_cut_ref = int(sys.argv[1])
+
+    
     # Class to compute the mesh transformation needed for higher order accuracy
     #  * order: order of the mesh deformation function
     #  * threshold: barrier for maximum deformation (to ensure shape regularity)
 
-    # Background FESpace
-    Vh = VectorH1(mesh, order=order, dirichlet=[])
-    Qh = H1(mesh, order=order-1, dirichlet = [])
-    ci = CutInfo(mesh, lset_approx)
-    ba_IF = ci.GetElementsOfType(IF)
-    VhG = Restrict(Vh, ba_IF)
-    L2G = Restrict(Qh, ba_IF)
-    print("pre product space")
-    X=VhG*L2G
-    gfu = GridFunction(X)
-
-    # Coefficients / parameters:
-
-    def coeffscalgrad(n):
-        temp = []
-        var = [x,y,z]
-        for v in var:
-            temp.append(n.Diff(v))
-        return CoefficientFunction(tuple(temp))
-    ntilde = Interpolate(coeffscalgrad(phi),VhG)
+    lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
+    deformation = lsetmeshadap.CalcDeformation(phi)
+    lset_approx = lsetmeshadap.lset_p1
 
-    # Tangential projection
-    n = Normalize(grad(lset_approx))
-#    n = Normalize(coeffscalgrad(lset_approx2))#/Norm(coeffscalgrad(lset_approx2))
-    n.Compile()
-    print(n.dim)
-    h = specialcf.mesh_size
+    mesh.SetDeformation(deformation)
 
-    eta = 100.0 / (h * h)
-    rhou = 1.0 / h
-    rhop = h
+# Background FESpaces
 
 
+    # Helper Functions for tangential projection and gradients of the exact solution
     def coef_grad(u):
         dirs = {1: [x], 2: [x, y], 3: [x, y, z]}
         if u.dim == 1:
@@ -103,11 +100,15 @@ with TaskManager():
         return Id(3) - OuterProduct(n, n)
 
 
+    # Coefficients / parameters:
+
+    # Exact tangential projection
     Ps = Pmats(Normalize(grad(phi)))
-    # define solution and right-hand side
 
+    # define solution and right-hand side
 
-    if geom=="circle":
+    if geom == "ellipsoid":
+        """
         functions = {
             "extvsol1": ((-y - z) * x + y * y + z * z),
             "extvsol2": ((-x - z) * y + x * x + z * z),
@@ -116,29 +117,21 @@ with TaskManager():
             "rhs2": ((-x - z) * y + x * x + z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
             "rhs3": ((-x - y) * z + x * x + y * y) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
         }
-        extpsol = (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"]))
-        pSol = CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
-        extpsol = CoefficientFunction(0)
-#        extpsol = CoefficientFunction((x*y**3+z*(x**2+y**2+z**2)**(3/2))/(x**2+y**2+z**2)**2)
-#        rhs1 = CoefficientFunction((functions["rhs1"], functions["rhs2"], functions["rhs3"]))
-    elif geom=="decocube":
-        uSol=Ps*CoefficientFunction((-z**2,y,x))
-        pSol = CoefficientFunction(x*y**3+z-1/Integrate({"levelset": lset_approx, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh)*Integrate({"levelset": lset_approx, "domain_type":IF},cf=x*y**3+z,mesh=mesh))
-        extpsol = pSol
-        extpsol.Compile()
-#        extpsol = CoefficientFunction(0)
-    elif geom=="torus":
-        rhs1 = CoefficientFunction((0,0,0))
-        rhs2 = CoefficientFunction(x+y)
-    u,p = X.TrialFunction()
-    v, q = X.TestFunction()
-    # Measure on surface
-    ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
-    # Measure on the bulk around the surface
-    dx = dx(definedonelements=ba_IF, deformation=deformation)
+        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))
+        pSol = CF(x*y**3+z)
+
+    elif geom == "decocube":
+        uSol = Ps*CoefficientFunction((-z ** 2, y, x))
+        pSol = CoefficientFunction(
+            x * y ** 3 + z - 1 / Integrate({"levelset": lset_approx, "domain_type": IF}, cf=CoefficientFunction(1),
+                                        mesh=mesh) * Integrate({"levelset": lset_approx, "domain_type": IF},
+                                                                cf=x * y ** 3 + z, mesh=mesh))
 
 
+
+    # Helper Functions for rate-of-strain-tensor
     def eps(u):
         return Ps * Sym(grad(u)) * Ps
 
@@ -147,70 +140,191 @@ with TaskManager():
         if u.dim == 3:
             return Trace(grad(u) * Ps)
         if u.dim == 9:
-            N = 3
-            divGi = [divG(u[i, :]) for i in range(N)]
-            return CF(tuple([divGi[i] for i in range(N)]))
-
-    weing = grad(n)
-    weingex = grad(Normalize(grad(phi)))
-#    weing.Compile()
-    print("start rhs")
-    if geom!="torus":
-    #    rhs1 = -Pmat*divtens(tangsymgrad(uSol)-(uSol*n)*weing)+Pmat*uSol+Pmat*CoefficientFunction((extpsol.Diff(x),extpsol.Diff(y),extpsol.Diff(z)))
-        rhs1 = Ps*-divG(eps(uSol)-(uSol*Normalize(grad(phi)))*weingex)+uSol+grad(extpsol)
-
- #       rhs1.Compile()
-        rhs2 = Trace(Ps*grad(uSol))
- #       rhs2.Compile()
-        print("end rhs")
-#    rhs2 = 0
-    # bilinear forms:
-    Pmat=Pmats(n)
-    def E(u):
-        return Pmat*Sym(grad(u))*Pmat-(u*n)*weing
-
-    a = BilinearForm(X, symmetric=True)
-#    a += (InnerProduct(Pmat * Sym(grad(u)) * Pmat-(u*n)*weing, Pmat * Sym(grad(v)) * Pmat-(v*n)*weing)) * ds
-#    a += (InnerProduct(Pmat * Sym(grad(u)) * Pmat, Pmat * Sym(grad(v)) * Pmat)) * ds
-    a += InnerProduct(E(u), E(v))*ds
-    a += (Pmat*u) * (Pmat*v) * ds
-    a += (eta * ((u * ntilde) * (v * ntilde))) * ds
-    a += rhou * ((grad(u) * n) * (grad(v) * n)) * dx
-    a += InnerProduct(u,Pmat*grad(q))*ds
-    a += InnerProduct(v, Pmat*grad(p)) * ds
-    a += -rhop*InnerProduct(n*grad(p),n*grad(q))*dx
-    print("preassembling")
-    a.Assemble()
-    print("assembled a")
-    f = LinearForm(X)
-    f += rhs1 * v * ds
-    f += -rhs2*q*ds
-    f.Assemble()
-    print("assembled f, inverting")
-    gfu.vec.data = a.mat.Inverse(freedofs = X.FreeDofs(), inverse="pardiso") * f.vec
-    if geom != "torus":
-        print("computing error")
-        print("l2error : ", sqrt(Integrate((gfu.components[0] - uSol)**2 * ds, mesh=mesh)))
-        uerr = (gfu.components[0]-uSol)
-
-    Draw(deformation, mesh, "deformation")
-    Draw(gfu, mesh, "u")
-    print(sys.argv[1])
-    if (sys.argv[1]=='1'):
-#        input("Continue (press enter) to create a VTK-Output to tracefem3d.vtk")
-
-        print("writing output")
-        if geom != "torus":
-            vtk = VTKOutput(ma=mesh,
-                        coefs=[lset_approx, gfu.components[0], uSol, uerr],
-                        names=["P1-levelset", "uh", "uex", "uerr"],
-                        filename="stokespenalty3d", subdivision=2)
+            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)]))
+    def solveStokes(mesh, lset_approx, deformation, step):
+        Vh = VectorH1(mesh, order=order, dirichlet=[])
+        Qh = H1(mesh, order=order - 1, dirichlet=[])
+
+
+        ci = CutInfo(mesh, lset_approx)
+        ba_IF = ci.GetElementsOfType(IF)
+        VhG = Restrict(Vh, ba_IF)
+        L2G = Restrict(Qh, ba_IF)
+
+        X = VhG * L2G
+        gfu = GridFunction(X)
+
+        n = Normalize(grad(lset_approx))
+        ntilde = Interpolate(grad(phi), VhG)
+
+        h = specialcf.mesh_size
+        elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
+        cuthmax = max([(6*vol)**(1/3) for vol in elvol])
+        alpha=3/2
+        epsilon =(cuthmax)**(alpha)
+        print("Epsilon:")
+        print(epsilon)
+#        epsilon=1
+        eta = 100.0 / (h * h)
+        rhou = 1.0 / h
+        rhop = h
+
+        # Measure on surface
+        ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
+        # Measure on the bulk around the surface
+        dX = dx(definedonelements=ba_IF, deformation=deformation)
+
+        u, p = X.TrialFunction()
+        v, q = X.TestFunction()
+        def rhsfromsol(uSol, pSol, interpol=True, Lpmat = None, Lpvec = None):
+            if interpol:
+
+                tmpeps = GridFunction(Lpmat)
+                tmpeps.Set(eps(uSol) ,definedonelements=ba_IF)
+                tmpdiv = GridFunction(Lpvec)
+                tmpdiv.Set(Ps*-divG(tmpeps), definedonelements=ba_IF)
+#                tmpuSol = GridFunction(L2(mesh, order=order+1, dim=3))
+#                tmpuSol.Set(uSol, definedonelements=ba_IF)
+                tmptotal = GridFunction(Lpvec)
+                tmptotal.Set(tmpdiv+Ps*grad(pSol), definedonelements=ba_IF)
+
+                return tmptotal
+            else:
+
+                return Ps * (-divG(eps(uSol).Compile())) +Ps* grad(pSol).Compile()
+
+        # Weingarten mappings
+
+
+        weing = grad(n)
+        weingex = grad(Normalize(grad(phi)))
+#        rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
+        L21 = L2(mesh, order=order+1, dim=9)
+        L22 = L2(mesh, order=order+1, dim=3)
+        lpmat = Restrict(L21, ba_IF)
+        lpvec = Restrict(L22, ba_IF)
+        if step>4:
+            comp=True
         else:
-            vtk = VTKOutput(ma=mesh,
-                            coefs=[lset_approx, gfu.components[0]],
-                            names=["P1-levelset", "uh"],
-                            filename="stokespenalty3d", subdivision=2)
+            comp=False
+        rhs1 = rhsfromsol(uSol, pSol, interpol, lpmat, lpvec).Compile(realcompile=comp)
+        rhs2 = Trace(Ps * grad(uSol)).Compile(realcompile=comp)
+
+        Pmat = Pmats(n)
+
+
+        def E(u):
+            return Pmat * Sym(grad(u)) * Pmat# - (u * n) * weing
+
+
+        a = BilinearForm(X, symmetric=True)
+        a += InnerProduct(E(u), E(v)).Compile(realcompile=comp) * ds
+        a +=epsilon* (Pmat * u) * (Pmat * v) * ds
+        a += (eta * ((u * ntilde) * (v * ntilde))) * ds
+        a += rhou * ((grad(u) * n) * (grad(v) * n)) * dX
+        a += InnerProduct(u, Pmat * grad(q)) * ds
+        a += InnerProduct(v, Pmat * grad(p)) * ds
+        a += -rhop * InnerProduct(n * grad(p), n * grad(q)) * dX
+        
+        a.Assemble()
+        u2, v2 = VhG.TnT()
+        a2 = BilinearForm(VhG, symmetric=True)
+
+        a2 += InnerProduct(E(u2),E(v2)).Compile(realcompile=comp) * ds
+#        a2 +=epsilon* (Pmat * u2) * (Pmat * v2) * ds
+        a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
+        a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
+
+        a2.Assemble()
+
+        b = BilinearForm(VhG, symmetric=True)
+        b += InnerProduct(u2,v2)*dX
+        b.Assemble()
+
+        evecs = GridFunction(VhG, multidim=3)
+        
+        vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 0)
+#        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
+
+        
+
+        f = LinearForm(X)
+        f += rhs1 * v * ds
+        f += -rhs2 * q * ds
+
+        f.Assemble()
+
+        """f2 = LinearForm(X)
+        f2 += rhs1ex * v * ds
+        f2 += -rhs2 * q * ds
+        gfu2 = GridFunction(X)
+        f2.Assemble()"""
+        gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
+        print("removing killing VFs")
+        print(vlam)
+        fun = GridFunction(VhG)
+        for k, lam in enumerate(vlam):
+            print("lam: "+ str(abs(lam)))
+            print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
+            if abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
+                
+                print("removing number: "+str(k))
+                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]
+        #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))
+        print("l2error : ", l2err)
+        return l2err
+
+
+
+    l2errors = []
+    for i in range(n_cut_ref+1):
+        if i!=0:
+            mesh.SetDeformation(None)
+            lsetp1 = GridFunction(H1(mesh, order=2))
+            InterpolateToP1(phi, lsetp1)
+            RefineAtLevelSet(lsetp1)
+
+            mesh.Refine()
+
+            lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
+
+            deformation = lsetmeshadap.CalcDeformation(phi)
+
+            lset_approx = lsetmeshadap.lset_p1
+
+            mesh.SetDeformation(deformation)
+
+
+        print("starting solve..")
+        l2errors.append(solveStokes(mesh, lset_approx, deformation, step=i))
+    print(l2errors)
+if len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
+#    input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
 
-        vtk.Do()
-print(datetime.datetime.now()-begin_time)
+#    with TaskManager():
+    pass
+    """vtk = VTKOutput(ma=mesh,
+                    coefs=[lset_approx, gfu.components[0], uSol, uerr],
+                    names=["P1-levelset", "uh", "uex", "uerr"],
+                    filename="stokestracefem3d", subdivision=2)
 
+    vtk.Do()"""