Wraped PETSc PCASM

Signed-off-by: Umberto Zerbinati <[email protected]>

docs/src/PETScPC/poisson.py.rst

@@ -18,12 +18,13 @@ Such a discretisation can easily be constructed using NGSolve as follows: ::
    from mpi4py.MPI import COMM_WORLD
 
    if COMM_WORLD.rank == 0:
-      mesh = Mesh(unit_square.GenerateMesh(maxh=0.1).Distribute(COMM_WORLD))
+      ngmesh = unit_square.GenerateMesh(maxh=0.1)
+      for _ in range(4):
+         ngmesh.Refine()
+      mesh = Mesh(ngmesh.Distribute(COMM_WORLD))
    else:
       mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
-   for _ in range(5):
-      mesh.Refine()
-   order = 4
+   order = 3
    fes = H1(mesh, order=order, dirichlet="left|right|top|bottom")
    print("Number of degrees of freedom: ", fes.ndof)
    u,v = fes.TnT()
@@ -70,7 +71,11 @@ In this case, we will use as fine space correction, the inverse of the local mat
     return blocks
 
    blocks = VertexPatchBlocks(mesh, fes)
-   blockjac = a.mat.CreateBlockSmoother(blocks)
+   dofs = BitArray(fes.ndof); dofs[:] = True
+   blockjac = PETScPreconditioner(a.mat, dofs, blocks=blocks,
+                                  solverParameters={"pc_type": "asm",
+                                                    "sub_ksp_type": "preonly",
+                                                    "sub_pc_type": "lu"})  
 
 We now isolate the degrees of freedom associated with the vertices and construct a two-level additive Schwarz preconditioner, where the coarse space correction is the inverse of the local matrices associated with the vertices. ::
 

docs/src/PETScPC/stokes.py.rst

@@ -33,8 +33,8 @@ Such a discretization can easily be constructed using NGSolve as follows: ::
    else:
       mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
    nu = Parameter(1.0)
-   V = VectorH1(mesh, order=2, dirichlet="wall|inlet|cyl", autoupdate=True)
-   Q = L2(mesh, order=0, autoupdate=True)
+   V = VectorH1(mesh, order=4, dirichlet="wall|inlet|cyl", autoupdate=True)
+   Q = L2(mesh, order=2, autoupdate=True)
    u,v = V.TnT(); p,q = Q.TnT()
    a = BilinearForm(nu*InnerProduct(Grad(u),Grad(v))*dx)
    a.Assemble()
@@ -172,7 +172,7 @@ This is not ideal for large problems, and we can use a `Hypre` preconditioner fo
 Our first attempt at using a `HYPRE` preconditioner for the Laplacian block did not converge. This is because the top left block of the saddle point problem now contains the augmentation term, which has a very large kernel.
 It is well known that algebraic multi-grid methods do not work well with indefinite problems, and this is what we are observing here. ::
 
-Let uss us consider an alternative approach to the augmented Lagrangian formulation. We begin by constructing the augmented Lagrangian formulation in more numerical linear algebra terms, i.e. ::
+Let us consider an alternative approach to the augmented Lagrangian formulation. We begin by constructing the augmented Lagrangian formulation in more numerical linear algebra terms, i.e. ::
 
    d = BilinearForm((1/gamma)*p*q*dx)
    d.Assemble()
@@ -221,7 +221,6 @@ In fact, since we know that the augmentation block has a lower rank than the Lap
                    printrates=True, initialize=False)
    Draw(gfu)
 
-   
    C = BlockMatrix( [ [apre + apre@([email protected]@b.mat)@apre, None], [None, mGpre.mat] ] )
 
    gfu.vec.data[:] = 0; gfp.vec.data[:] = 0;

ngsPETSc/pc.py

@@ -25,8 +25,8 @@ class PETScPreconditioner(BaseMatrix):
     MKL sparse: mklaij or CUDA: aijcusparse
 
     '''
-    def __init__(self, mat, freeDofs, solverParameters=None, optionsPrefix=None, nullspace=None,
-                 matType="aij"):
+    def __init__(self, mat, freeDofs, solverParameters=None, optionsPrefix="", nullspace=None,
+                 matType="aij", blocks=None):
         BaseMatrix.__init__(self)
         if hasattr(solverParameters, "ToDict"):
             solverParameters = solverParameters.ToDict()
@@ -39,18 +39,38 @@ def __init__(self, mat, freeDofs, solverParameters=None, optionsPrefix=None, nul
         petscMat = Matrix(self.ngsMat, (dofs, freeDofs, None), matType).mat
         if nullspace is not None:
             if nullspace.near:
-                petscMat.setNearNullSpace(nullspace.nullspace)
+                petscMat.mat.setNearNullSpace(nullspace.nullspace)
         self.petscPreconditioner = PETSc.PC().create(comm=petscMat.getComm())
         self.petscPreconditioner.setOperators(petscMat)
         options_object = PETSc.Options()
         if solverParameters is not None:
             for optName, optValue in solverParameters.items():
-                if optName not in ["matType"]:
-                    options_object[optName] = optValue
+                if "sub_pc_type"==optName and blocks is not None:
+                    options_object["sub_pc_type"] = "lu"
+                if "sub_pc_factor_mat_ordering_type"==optName and blocks is not None:
+                    options_object["sub_pc_factor_mat_ordering_type"] = "natural"
+                options_object[optName] = optValue
         self.petscPreconditioner.setOptionsPrefix(optionsPrefix)
         self.petscPreconditioner.setFromOptions()
-        self.petscPreconditioner.setUp()
-        #self.petscPreconditioner.view()
+        self.asmpc = None
+        if blocks is not None:
+            if len (blocks) == 0:
+                self.ises = [PETSc.IS().createGeneral([0], comm=PETSc.COMM_SELF)]
+            else:
+                lgmap = petscMat.getLGMap()[0] #TODO: Fiox this only works for rc
+                self.ises = [lgmap.applyIS(PETSc.IS().createGeneral(list(block),
+                             comm=PETSc.COMM_SELF)) for block in blocks]
+            self.asmpc = PETSc.PC().create(comm=self.petscPreconditioner.getComm())
+            self.asmpc.incrementTabLevel(1, parent=self.petscPreconditioner)
+            self.asmpc.setOptionsPrefix(optionsPrefix)
+            self.asmpc.setFromOptions()
+            self.asmpc.setOperators(*self.petscPreconditioner.getOperators())
+            self.asmpc.setType(self.asmpc.Type.ASM)
+            self.asmpc.setASMType(PETSc.PC.ASMType.BASIC)
+            self.asmpc.setASMLocalSubdomains(len(self.ises), self.ises)
+            self.asmpc.setUp()
+        else:
+            self.petscPreconditioner.setUp()
         self.petscVecX, self.petscVecY = petscMat.createVecs()
     def Shape(self):
         '''
@@ -76,7 +96,10 @@ def Mult(self,x,y):
 
         '''
         self.vecMap.petscVec(x,self.petscVecX)
-        self.petscPreconditioner.apply(self.petscVecX, self.petscVecY)
+        if self.asmpc is not None:
+            self.asmpc.apply(self.petscVecX, self.petscVecY)
+        else:
+            self.petscPreconditioner.apply(self.petscVecX, self.petscVecY)
         self.vecMap.ngsVec(self.petscVecY, y)
 
     def MultTrans(self,x,y):
@@ -87,7 +110,10 @@ def MultTrans(self,x,y):
 
         '''
         self.vecMap.petscVec(x,self.petscVecX)
-        self.petscPreconditioner.applyTranspose(self.petscVecX, self.petscVecY)
+        if self.asmpc is not None:
+            self.asmpc.applyTranspose(self.petscVecX, self.petscVecY)
+        else:
+            self.petscPreconditioner.applyTranspose(self.petscVecX, self.petscVecY)
         self.vecMap.ngsVec(self.petscVecY, y)
 
 def createPETScPreconditioner(mat, freeDofs, solverParameters):