From 65504742dd9bd7d0639f99d611390ccecec7959a Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Fri, 31 May 2024 12:19:02 +0100 Subject: [PATCH] wip on the docs Signed-off-by: Umberto Zerbinati --- docs/src/PETScPC/poisson.py | 131 ++++++++++++++++++ docs/src/{PETSc PC => PETScPC}/poisson.py.rst | 9 +- docs/src/index.rst | 11 +- 3 files changed, 147 insertions(+), 4 deletions(-) create mode 100644 docs/src/PETScPC/poisson.py rename docs/src/{PETSc PC => PETScPC}/poisson.py.rst (95%) diff --git a/docs/src/PETScPC/poisson.py b/docs/src/PETScPC/poisson.py new file mode 100644 index 0000000..ccbbb94 --- /dev/null +++ b/docs/src/PETScPC/poisson.py @@ -0,0 +1,131 @@ +# Preconditioning the Poisson problem +# ===================================== +# +# In this tutorial, we explore using `PETSc PC` as a preconditioner inside NGSolve preconditioning infrastructure. +# We will focus our attention on the Poisson problem, and we will consider different preconditioning strategies. +# Not all the preconditioning strategies are equally effective for all the discretizations of the Poisson problem, and this demo is intended to provide a starting point for the exploration of the preconditioning strategies rather than providing a definitive answer. +# We begin by creating a discretisation of the Poisson problem using H1 elements, in particular, we consider the usual variational formulation +# +# .. math:: +# +# \text{find } u\in H^1_0(\Omega) \text{ s.t. } a(u,v) := \int_{\Omega} \nabla u\cdot \nabla v \; d\vec{x} = L(v) := \int_{\Omega} fv\; d\vec{x}\qquad v\in H^1_0(\Omega). +# +# Such a discretisation can easily be constructed using NGSolve as follows: :: + +from ngsolve import * +import netgen.gui +import netgen.meshing as ngm +from mpi4py.MPI import COMM_WORLD + +if COMM_WORLD.rank == 0: + mesh = Mesh(unit_square.GenerateMesh(maxh=0.1).Distribute(COMM_WORLD)) +else: + mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD)) +for _ in range(5): + mesh.Refine() +order = 4 +fes = H1(mesh, order=order, dirichlet="left|right|top|bottom") +print("Number of degrees of freedom: ", fes.ndof) +u,v = fes.TnT() +a = BilinearForm(grad(u)*grad(v)*dx) +f = LinearForm(fes) +f += 32 * (y*(1-y)+x*(1-x)) * v * dx +a.Assemble() +f.Assemble() + +# We now construct an NGSolve preconditioner wrapping a `PETSc PC`, we begin considering an Algebraic MultiGrid preconditioner constructed using `HYPRE`. +# In this tutorial, we will use the Krylov solver implemented inside NGSolve to solve the linear system. :: + +from ngsPETSc.pc import * +from ngsolve.krylovspace import CG +pre = Preconditioner(a, "PETScPC", pc_type="hypre") +gfu = GridFunction(fes) +print("-------------------|HYPRE p={}|-------------------".format(order)) +gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pre.mat, printrates=True) +Draw(gfu) + +# We see that the HYPRE preconditioner is quite effective for the Poisson problem discretised using linear elements, but it is not as effective for higher-order elements. +# .. list-table:: Title +# :widths: 25 25 50 +# :header-rows: 1 +# +# * - Preconditioner +# - p=1 +# - p=2 +# - p=3 +# * - HYPRE +# - 15 (8.49e-13) +# - 100 (4.81e-8) +# - 100 (3.60e-9) +# * - Two Level Additive Schwarz +# - 59 (1.74e-12) +# - 58 (2.01e-12) +# - 59 (1.72e-12) +# +# We can use PETSc preconditioner as one of the building blocks of a more complex preconditioner. For example, we can construct an additive Schwarz preconditioner. +# In this case, we will use as fine space correction, the inverse of the local matrices associated with the patch of a vertex. :: + +def VertexPatchBlocks(mesh, fes): + blocks = [] + freedofs = fes.FreeDofs() + for v in mesh.vertices: + vdofs = set() + for el in mesh[v].elements: + vdofs |= set(d for d in fes.GetDofNrs(el) + if freedofs[d]) + blocks.append(vdofs) + return blocks + +blocks = VertexPatchBlocks(mesh, fes) +blockjac = a.mat.CreateBlockSmoother(blocks) + +# 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. :: + +def VertexDofs(mesh, fes): + vertexdofs = BitArray(fes.ndof) + vertexdofs[:] = False + for v in mesh.vertices: + for d in fes.GetDofNrs(v): + vertexdofs[d] = True + vertexdofs &= fes.FreeDofs() + return vertexdofs + +vertexdofs = VertexDofs(mesh, fes) +preCoarse = PETScPreconditioner(a.mat, vertexdofs, solverParameters={"pc_type": "hypre"}) +pretwo = preCoarse + blockjac +print("-------------------|Additive Schwarz p={}|-------------------".format(order)) +gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pretwo, printrates=True) + + +# We can also use the PETSc preconditioner as an auxiliary space preconditioner. +# Let us consdier the disctinuous Galerkin discretisation of the Poisson problem. +# +# fesDG = L2(mesh, order=3, dgjumps=True) +# u,v = fesDG.TnT() +# aDG = BilinearForm(fesDG) +# jump_u = u-u.Other(); jump_v = v-v.Other() +# n = specialcf.normal(2) +# mean_dudn = 0.5*n * (grad(u)+grad(u.Other())) +# mean_dvdn = 0.5*n * (grad(v)+grad(v.Other())) +# alpha = 4 +# h = specialcf.mesh_size +# aDG = BilinearForm(fesDG) +# aDG += grad(u)*grad(v) * dx +# aDG += alpha*3**2/h*jump_u*jump_v * dx(skeleton=True) +# aDG += alpha*3**2/h*u*v * ds(skeleton=True) +# aDG += (-mean_dudn*jump_v -mean_dvdn*jump_u)*dx(skeleton=True) +# aDG += (-n*grad(u)*v-n*grad(v)*u)*ds(skeleton=True) +# fDG = LinearForm(fesDG) +# fDG += 1*v * dx +# aDG.Assemble() +# fDG.Assemble() +# +# We can now use the PETSc PC assembled for the confroming Poisson problem as an auxiliary space preconditioner for the DG discretisation. +# +# from ngsPETSc import pc +# smoother = Preconditioner(aDG, "PETScPC", pc_type="sor") +# transform = fes.ConvertL2Operator(fesDG) +# preDG = transform @ pre.mat @ transform.T + smoother.mat +# gfuDG = GridFunction(fesDG) +# gfuDG.vec.data = CG(aDG.mat, rhs=fDG.vec, pre=preDG, printrates=True) +# Draw(gfuDG) \ No newline at end of file diff --git a/docs/src/PETSc PC/poisson.py.rst b/docs/src/PETScPC/poisson.py.rst similarity index 95% rename from docs/src/PETSc PC/poisson.py.rst rename to docs/src/PETScPC/poisson.py.rst index 033c65b..a996889 100755 --- a/docs/src/PETSc PC/poisson.py.rst +++ b/docs/src/PETScPC/poisson.py.rst @@ -45,6 +45,7 @@ In this tutorial, we will use the Krylov solver implemented inside NGSolve to so Draw(gfu) We see that the HYPRE preconditioner is quite effective for the Poisson problem discretised using linear elements, but it is not as effective for higher-order elements. + .. list-table:: Preconditioner performance :widths: 25 25 50 :header-rows: 1 @@ -53,6 +54,7 @@ We see that the HYPRE preconditioner is quite effective for the Poisson problem - p=1 - p=2 - p=3 + - p=4 * - HYPRE - 15 (8.49e-13) - 100 (4.81e-8) @@ -93,6 +95,8 @@ We now isolate the degrees of freedom associated with the vertices and construct print("-------------------|Additive Schwarz p={}|-------------------".format(order)) gfu.vec.data = CG(a.mat, rhs=f.vec, pre=pretwo, printrates=True) +We can see that the two-level additive Schwarz preconditioner where the coarse space correction is performed using HYPRE is more effective than using just HYPRE preconditioner for higher-order elements. + .. list-table:: Preconditioner performance :widths: 25 25 50 :header-rows: 1 @@ -101,6 +105,7 @@ We now isolate the degrees of freedom associated with the vertices and construct - p=1 - p=2 - p=3 + - p=4 * - HYPRE - 15 (8.49e-13) - 100 (4.81e-8) @@ -114,7 +119,7 @@ We now isolate the degrees of freedom associated with the vertices and construct We can also use the PETSc preconditioner as an auxiliary space preconditioner. -Let us consdier the disctinuous Galerkin discretisation of the Poisson problem. +Let us consdier the disctinuous Galerkin discretisation of the Poisson problem. :: fesDG = L2(mesh, order=3, dgjumps=True) u,v = fesDG.TnT() @@ -136,7 +141,7 @@ Let us consdier the disctinuous Galerkin discretisation of the Poisson problem. aDG.Assemble() fDG.Assemble() -We can now use the PETSc PC assembled for the confroming Poisson problem as an auxiliary space preconditioner for the DG discretisation. +We can now use the PETSc PC assembled for the confroming Poisson problem as an auxiliary space preconditioner for the DG discretisation. :: from ngsPETSc import pc smoother = Preconditioner(aDG, "PETScPC", pc_type="sor") diff --git a/docs/src/index.rst b/docs/src/index.rst index 5633e1c..894c3f7 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -7,10 +7,17 @@ Welcome to ngsPETSc's documentation! ==================================== .. toctree:: - :maxdepth: 2 - :caption: Contents: + :maxdepth: 1 + :caption: Getting Started install + +.. toctree:: + :maxdepth: 1 + :caption: PETSc PC + + PETScPC/poisson.py +