PCPatch: validate mesh overlap (#4684)

* PCPatch: validate mesh overlap

Author: pbrubeck

demos/patch/hcurl_riesz_star.py.rst

@@ -5,11 +5,17 @@ Contributed by `Robert Kirby <https://sites.baylor.edu/robert_kirby/>`_
 and `Pablo Brubeck <https://www.maths.ox.ac.uk/people/pablo.brubeckmartinez/>`_.
 
 Multigrid in H(div) and H(curl) also requires relaxation based on topological patches.
-Here, we demonstrate how to do this in the latter case. ::
+Here, we demonstrate how to do this in the latter case.
+
+We start by importing firedrake and setting up a :func:`.MeshHierarchy` and the
+exact solution and forcing data. Crucially, the meshes must have an overlapping
+parallel domain decomposition that supports the vertex star patches. This is set
+via the `distribution_parameters` kwarg of the :func:`.Mesh` constructor. ::
 
   from firedrake import *
 
-  base = UnitCubeMesh(2, 2, 2)
+  dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
+  base = UnitCubeMesh(2, 2, 2, distribution_parameters=dparams)
   mh = MeshHierarchy(base, 3)
   mesh = mh[-1]
 

demos/patch/hdiv_riesz_star.py.rst

@@ -5,11 +5,17 @@ Contributed by `Robert Kirby <https://sites.baylor.edu/robert_kirby/>`_
 and `Pablo Brubeck <https://www.maths.ox.ac.uk/people/pablo.brubeckmartinez/>`_.
 
 Multigrid in H(div) and H(curl) also requires relaxation based on topological patches.
-Here, we demonstrate how to do this in the former case. ::
+Here, we demonstrate how to do this in the former case.
+
+We start by importing firedrake and setting up a :func:`.MeshHierarchy` and the
+exact solution and forcing data. Crucially, the meshes must have an overlapping
+parallel domain decomposition that supports the vertex star patches. This is set
+via the `distribution_parameters` kwarg of the :func:`.Mesh` constructor. ::
 
   from firedrake import *
 
-  base = UnitCubeMesh(2, 2, 2)
+  dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
+  base = UnitCubeMesh(2, 2, 2, distribution_parameters=dparams)
   mh = MeshHierarchy(base, 3)
   mesh = mh[-1]
 

demos/patch/poisson_mg_patches.py.rst

@@ -4,7 +4,7 @@ Using patch relaxation for multigrid
 Contributed by `Robert Kirby <https://sites.baylor.edu/robert_kirby/>`_
 and `Pablo Brubeck <https://www.maths.ox.ac.uk/people/pablo.brubeckmartinez/>`_.
 
-Simple relaxation like point Jacobi are not optimal or even suitable
+Simple multigrid relaxation methods like point Jacobi are not optimal or even suitable
 smoothers for all applications.  Firedrake supports additive Schwarz methods
 based on local patch-based decompositions through two different paths.
 
@@ -20,12 +20,15 @@ degree-independent iteration counts.  Here, all the degrees of freedom in the ce
 For many problems, point Jacobi is even worse, and patches are required even to
 get a convergent method.  We refer the reader to other demos.
 
-We start by importing firedrake and setting up a mesh hierarchy and the
-exact solution and forcing data. ::
+We start by importing firedrake and setting up a :func:`.MeshHierarchy` and the
+exact solution and forcing data. Crucially, the meshes must have an overlapping
+parallel domain decomposition that supports the vertex star patches. This is set
+via the `distribution_parameters` kwarg of the :func:`.Mesh` constructor. ::
 
   from firedrake import *
 
-  base = UnitSquareMesh(4, 4)
+  dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
+  base = UnitSquareMesh(4, 4, distribution_parameters=dparams)
   mh = MeshHierarchy(base, 1)
   mesh = mh[-1]
 

demos/patch/stokes_vanka_patches.py.rst

@@ -16,11 +16,17 @@ of the patch but not pressures.
 In practice, we arrive at mesh-independent multigrid convergence using these relaxation.
 We can construct Vanka patches either through :class:`~.PatchPC`, in which the bilinear form
 is assembled on each vertex patch, or through :class:`~.ASMVankaPC`, in which the patch
-operators are extracted from the globally assembled stiffness matrix. ::
+operators are extracted from the globally assembled stiffness matrix.
+
+We start by importing firedrake and setting up a :func:`.MeshHierarchy` and the
+exact solution and forcing data. Crucially, the meshes must have an overlapping
+parallel domain decomposition that supports the Vanka patches. This is set
+via the `distribution_parameters` kwarg of the :func:`.Mesh` constructor. ::
 
   from firedrake import *
 
-  base = UnitSquareMesh(4, 4)
+  dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 2)}
+  base = UnitSquareMesh(4, 4, distribution_parameters=dparams)
   mh = MeshHierarchy(base, 3)
   mesh = mh[-1]
 

docs/source/preconditioning.rst

@@ -69,6 +69,19 @@ multiplicatively within an MPI rank and additively between ranks.
    entities into lines or planes (useful for advection-dominated
    problems).
 
+.. note::
+   The additive Schwarz preconditioners listed here construct patches around
+   mesh entities.  Crucially, the mesh must have an overlapping parallel domain
+   decomposition that supports the patches. This is set via the
+   `distribution_parameters` kwarg of the :func:`.Mesh` constructor.  For
+   instance, vertex-star patches require ::
+
+      distribution_parameters["overlap_type"] = (DistributedMeshOverlapType.VERTEX, 1)
+
+   while Vanka patches require ::
+
+      distribution_parameters["overlap_type"] = (DistributedMeshOverlapType.VERTEX, 2)
+
 Multigrid methods
 =================
 

firedrake/preconditioners/asm.py

@@ -4,6 +4,7 @@
 from firedrake.preconditioners.base import PCBase
 from firedrake.petsc import PETSc
 from firedrake.dmhooks import get_function_space
+from firedrake.mesh import DistributedMeshOverlapType
 from firedrake.logging import warning
 from tinyasm import _tinyasm as tinyasm
 from mpi4py import MPI
@@ -160,6 +161,8 @@ def get_patches(self, V):
         opts = PETSc.Options(self.prefix)
         depth = opts.getInt("construct_dim", default=0)
         ordering = opts.getString("mat_ordering_type", default="natural")
+        validate_overlap(mesh, depth, "star")
+
         # Accessing .indices causes the allocation of a global array,
         # so we need to cache these for efficiency
         V_local_ises_indices = tuple(iset.indices for iset in V.dof_dset.local_ises)
@@ -233,8 +236,11 @@ def get_patches(self, V):
         ises = []
         if depth != -1:
             (start, end) = mesh_dm.getDepthStratum(depth)
+            patch_dim = depth
         else:
             (start, end) = mesh_dm.getHeightStratum(height)
+            patch_dim = mesh_dm.getDimension() - height
+        validate_overlap(mesh, patch_dim, "vanka")
 
         for seed in range(start, end):
             # Only build patches over owned DoFs
@@ -453,6 +459,7 @@ def get_patches(self, V):
             else:
                 continue
 
+            validate_overlap(mesh, base_depth, "star")
             start, end = mesh_dm.getDepthStratum(base_depth)
             for seed in range(start, end):
                 # Only build patches over owned DoFs
@@ -508,3 +515,26 @@ def get_patches(self, V):
                     iset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF)
                     ises.append(iset)
         return ises
+
+
+def validate_overlap(mesh, patch_dim, patch_type):
+    if patch_type == "python":
+        return
+    patch_depth = {"pardecomp": 0, "star": 1, "vanka": 2}[patch_type]
+
+    tdim = mesh.topology_dm.getDimension()
+    overlap_entity, overlap_depth = mesh._distribution_parameters["overlap_type"]
+    overlap_dim = {
+        DistributedMeshOverlapType.VERTEX: 0,
+        DistributedMeshOverlapType.FACET: tdim-1,
+        DistributedMeshOverlapType.NONE: tdim,
+    }[overlap_entity]
+
+    if mesh.comm.size > 1:
+        if overlap_dim > patch_dim:
+            patch_entity = {0: "vertex", 1: "edge", 2: "face", tdim: "cell"}[patch_dim]
+            warning(f"{overlap_entity} does not support {patch_entity}-patches. "
+                    "Did you forget to set overlap_type in your mesh's distribution_parameters?")
+        if overlap_depth < patch_depth:
+            warning(f"Mesh overlap depth of {overlap_depth} does not support {patch_type}-patches. "
+                    "Did you forget to set overlap_type in your mesh's distribution_parameters?")

firedrake/preconditioners/patch.py

@@ -1,4 +1,5 @@
 from firedrake.preconditioners.base import PCBase, SNESBase, PCSNESBase
+from firedrake.preconditioners.asm import validate_overlap
 from firedrake.petsc import PETSc
 from firedrake.cython.patchimpl import set_patch_residual, set_patch_jacobian
 from firedrake.solving_utils import _SNESContext
@@ -772,15 +773,24 @@ def initialize(self, obj):
         if mesh.cell_set._extruded:
             raise NotImplementedError("Not implemented on extruded meshes")
 
-        if "overlap_type" not in mesh._distribution_parameters:
-            if mesh.comm.size > 1:
-                # Want to do
-                # warnings.warn("You almost surely want to set an overlap_type in your mesh's distribution_parameters.")
-                # but doesn't warn!
-                PETSc.Sys.Print("Warning: you almost surely want to set an overlap_type in your mesh's distribution_parameters.")
+        # Validate the mesh overlap
+        prefix = (obj.getOptionsPrefix() or "") + "patch_"
+        opts = PETSc.Options(prefix)
+        petsc_prefix = self._petsc_prefix
+        patch_type = opts.getString(f"{petsc_prefix}construct_type")
+        patch_dim = opts.getInt(f"{petsc_prefix}construct_dim", -1)
+        patch_codim = opts.getInt(f"{petsc_prefix}construct_codim", -1)
+        if patch_dim != -1:
+            assert patch_codim == -1, "Cannot set both dim and codim"
+        elif patch_codim != -1:
+            assert patch_dim == -1, "Cannot set both dim and codim"
+            patch_dim = self.plex.getDimension() - patch_codim
+        else:
+            patch_dim = 0
+        validate_overlap(mesh, patch_dim, patch_type)
 
         patch = obj.__class__().create(comm=mesh.comm)
-        patch.setOptionsPrefix((obj.getOptionsPrefix() or "") + "patch_")
+        patch.setOptionsPrefix(prefix)
         self.configure_patch(patch, obj)
         patch.setType("patch")
 
@@ -937,6 +947,8 @@ def view(self, pc, viewer=None):
 
 class PatchPC(PCBase, PatchBase):
 
+    _petsc_prefix = "pc_patch_"
+
     def configure_patch(self, patch, pc):
         (A, P) = pc.getOperators()
         patch.setOperators(A, P)
@@ -949,6 +961,9 @@ def applyTranspose(self, pc, x, y):
 
 
 class PatchSNES(SNESBase, PatchBase):
+
+    _petsc_prefix = "snes_patch_"
+
     def configure_patch(self, patch, snes):
         patch.setTolerances(max_it=1)
         patch.setConvergenceTest("skip")