diff --git a/docs/src/PETScPC/poisson.py.rst b/docs/src/PETScPC/poisson.py.rst index 32f79fe..def6580 100755 --- a/docs/src/PETScPC/poisson.py.rst +++ b/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. :: diff --git a/docs/src/PETScPC/stokes.py.rst b/docs/src/PETScPC/stokes.py.rst index 3b7fb2c..de32860 100755 --- a/docs/src/PETScPC/stokes.py.rst +++ b/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@(b.mat.T@mpre.mat@b.mat)@apre, None], [None, mGpre.mat] ] ) gfu.vec.data[:] = 0; gfp.vec.data[:] = 0; diff --git a/ngsPETSc/pc.py b/ngsPETSc/pc.py index 4dafc03..108fa62 100644 --- a/ngsPETSc/pc.py +++ b/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):