merging

Author: rckirby

demos/patch/hcurl_riesz_star.py.rst

@@ -1,5 +1,5 @@
-Using patch relaxation for H(div)
-=================================
+Using patch relaxation for H(curl)
+==================================
 
 Contributed by `Robert Kirby <https://sites.baylor.edu/robert_kirby/>`_
 and `Pablo Brubeck <https://www.maths.ox.ac.uk/people/pablo.brubeckmartinez/>`_.
@@ -13,7 +13,7 @@ Here, we demonstrate how to do this in the latter case.::
   mh = MeshHierarchy(mesh, 3)
   mesh = mh[-1]
 
-We consider the Riesz map on H(div), discretized with lowest order
+We consider the Riesz map on H(curl), discretized with lowest order
 Nedelec elements.  We force the system with a random right-hand side and
 impose homogeneous Dirichlet boundary conditions::
 
@@ -45,13 +45,13 @@ patches yield a robust method.::
           "ksp_type": "cg",
           "pc_type": "mg",
           "mg_levels": {
-               "ksp_type": "chebyshev",
-               "ksp_max_it": 1,
-                **relax
+              "ksp_type": "chebyshev",
+              "ksp_max_it": 1,
+              **relax
           },
           "mg_coarse": {
               "ksp_type": "preonly",
-               "pc_type": "cholesky"
+              "pc_type": "cholesky"
           }
       }
 
@@ -65,27 +65,30 @@ patches yield a robust method.::
       }
 
 Hiptmair proposed a finer space decomposition for Nedelec elements using edge
-patches and vertex patches on the gradient of a Lagrange space. The python type
-preconditioner :class:`~.HiptmairPC` automatically sets up a two-level method
-using the auxiliary Lagrange space in a multigrid hierarchy. ::
+patches on the original Nedelec space and vertex patches on the gradient of a Lagrange space. The python type
+preconditioner :class:`~.HiptmairPC` automatically sets up an additive two-level method
+using the auxiliary Lagrange space in a multigrid hierarchy. Therefore, the overall multigrid relaxation composes the edge patches with the auxiliary space relaxation. For the latter, the residual on each level is restricted from the dual of H(curl) into the dual of H1 via the adjoint of the gradient, where a vertex patch relaxation is applied to obtain a correction that is prolonged from H1 into H(curl) via the gradient. ::
 
 
   def hiptmair_params():
       return {
-         "pc_type": "python",
-         "pc_python_type": "firedrake.HiptmairPC",
-         "hiptmair_mg_levels": asm_params(1),
-         "hiptmair_mg_coarse": asm_params(0),
+          "pc_type": "python",
+          "pc_python_type": "firedrake.HiptmairPC",
+          "hiptmair_mg_coarse": asm_params(0),
+          "hiptmair_mg_levels": asm_params(1),
+          "hiptmair_mg_levels_ksp_type": "richardson",
+          "hiptmair_mg_levels_ksp_max_it": 1,
+          "hiptmair_mg_coarse_ksp_type": "preonly",
       }
 
 
 Now, for each parameter choice, we report the iteration count for the Riesz map
-over a range of meshes.  We see that the auxiliary space approach gives lower
-iteration counts than vertex patches, while being cheaper to invert.::
+over a range of meshes.  We see that vertex patches approach give lower
+iteration counts than the Hiptmair approach, but they are more expensive.::
 
   names = {
-     "Vertex Star": mg_params(asm_params(0)),
-     "Hiptmair": mg_params(hiptmair_params()),
+      "Vertex Star": mg_params(asm_params(0)),
+      "Hiptmair": mg_params(hiptmair_params()),
   }
 
   for name, parameters in names.items():
@@ -105,16 +108,16 @@ For vertex patches, we expect output like,
   3        16
 ======== ============
 
-and with Hiptmair (edge patches + vertex patches on gradients of H1)
+and with Hiptmair (edge patches + vertex patches on gradients of Lagrange)
 
 ======== ============
  Level    Iterations
 ======== ============
-  1        10
-  2        12
-  3        13
+  1        18
+  2        20
+  3        21
 ======== ============
 
 and additional mesh refinement will lead to these numbers leveling off.
 
-A runnable python version of this demo can be found :demo:`here<hdiv_riesz_star.py>`.
+A runnable python version of this demo can be found :demo:`here<hcurl_riesz_star.py>`.

demos/patch/hdiv_riesz_star.py.rst

@@ -39,22 +39,28 @@ Having done both :class:`~.ASMStarPC` and :class:`~.PatchPC` in other demos, her
 Arnold, Falk, and Winther show that either vertex (`construct_dim=0`) or edge patches (`construct_dim=1`)  will be acceptable in three dimensions.::
 
 
-  def asm_params(construct_dim):
+  def mg_params(relax):
       return {
           "ksp_type": "cg",
-	  "pc_type": "mg",
-	  "mg_levels": {
-	      "ksp_type": "chebyshev",
-	      "ksp_max_it": 1,
-	      "pc_type": "python",
-              "pc_python_type": "firedrake.ASMStarPC",
-              "pc_star_construct_dim": construct_dim,
-              "pc_star_backend_type": "tinyasm"
-	  },
-	  "mg_coarse": {
-	      "ksp_type": "preonly",
-	      "pc_type": "cholesky",
-	  }
+          "pc_type": "mg",
+          "mg_levels": {
+              "ksp_type": "chebyshev",
+              "ksp_max_it": 1,
+              **relax
+          },
+          "mg_coarse": {
+              "ksp_type": "preonly",
+              "pc_type": "cholesky"
+          }
+      }
+
+
+  def asm_params(construct_dim):
+      return {
+          "pc_type": "python",
+          "pc_python_type": "firedrake.ASMStarPC",
+          "pc_star_construct_dim": construct_dim,
+          "pc_star_backend_type": "tinyasm"
       }
 
 Now, for each parameter choice, we report the iteration count for the Riesz map
@@ -63,10 +69,11 @@ edge patches, but they are more expensive.::
 
 
   for cdim in (0, 1):
+      params = mg_params(asm_params(cdim))
       print(f"Relaxation with patches of dimension {cdim}")
       print("Level | Iterations")
       for lvl, msh in enumerate(mh[1:], start=1):
-          its = run_solve(msh, asm_params(cdim))
+          its = run_solve(msh, params)
           print(f"{lvl}     | {its}")
 
 For vertex patches, we expect output like,

demos/patch/poisson_mg_patches.py.rst

@@ -11,7 +11,7 @@ based on local patch-based decompositions through two different paths.
 This demonstration illustrates basic usage of these methods for the Poisson
 problem.  Here, multigrid with point Jacobi relaxation works, but the iteration
 count degrades with polynomial degree, while vertex star patches give
-degree-indendent iteration counts.  Here, all the degrees of freedom in the cells
+degree-independent iteration counts.  Here, all the degrees of freedom in the cells
 (but not on the boundary) around each vertex are included in a patch
 
 .. figure:: star.svg
@@ -63,7 +63,7 @@ on the coarsest level of the multigrid hierarchy.::
   }
 
 When we use a matrix-free method, there will not be an assembled matrix to factor
-This forces the matrix to be assembled.::
+on the coarse level. This forces the matrix to be assembled.::
 
   assembled_lu = {
       "ksp_type": "preonly",
@@ -103,8 +103,8 @@ degree increases.::
 
   jacobi_relax = mg_params({"pc_type": "jacobi"}, mat_type="matfree")
 
-These options specify an additive Schwarz relaxation through PatchPC.
-PatchPC builds the patch operators by assembling the bilineary form over
+These options specify an additive Schwarz relaxation through :class:`~.PatchPC`.
+:class:`~.PatchPC` builds the patch operators by assembling the bilineary form over
 each subdomain.  Hence, it does not require the global stiffness
 matrix to be assembled.
 These options tell the patch mechanism to use vertex star patches, storing
@@ -128,7 +128,7 @@ in dense format.::
           "sub_pc_type": "lu"}},
       mat_type="matfree")
 
-ASMStarPC, on the other hand, does no re-discretization, but extracts the
+:class:`~.ASMStarPC`, on the other hand, does no re-discretization, but extracts the
 patch operators for each patch from the already-assembled global stiffness matrix.
 
 
@@ -144,7 +144,7 @@ direct or Krylov method on each one.::
 
 Now, for each parameter choice, we report the iteration count for the Poisson problem
 over a range of polynomial degrees.  We see that the Jacobi relaxation leads to growth
-in iteration count, while both PatchPC and ASMStarPC do not.  Mathematically, the two
+in iteration count, while both :class:`~.PatchPC` and :class:`~.ASMStarPC` do not.  Mathematically, the two
 latter options do the same operations, just via different code paths.::
 
   names = {"Jacobi": jacobi_relax,
@@ -173,7 +173,7 @@ For Jacobi, we expect output such as
    7         19
 ======== ================
 
-While for either PatchPC or ASMStarPC, we expect
+While for either :class:`~.PatchPC` or :class:`~.ASMStarPC`, we expect
 
 ======== ================
  Degree    Iterations

demos/patch/stokes_vanka_patches.py.rst

@@ -14,7 +14,7 @@ of the patch but not pressures.
    :align: center
 
 In practice, we arrive at mesh-independent multigrid convergence using these relaxation.
-We can construct Vanka patches either through firedrake.PatchPC, in which the bilinear form
+We can construct Vanka patches either through :class:`~.PatchPC`, in which the bilinear form
 is assembled on each vertex patch, or through firedrake.ASMVankaPC, in which the patch
 operators are extracted from the globally assembled stiffness matrix.::
 
@@ -62,7 +62,7 @@ on the coarsest level of the multigrid hierarchy.::
   }
 
 When we use a matrix-free method, there will not be an assembled matrix to factor
-This forces the matrix to be assembled.::
+on the coarse level. This forces the matrix to be assembled.::
 
   assembled_ldlt = {
       "ksp_type": "preonly",
@@ -94,8 +94,8 @@ relaxation options and matrix assembled type.::
       }
 
 
-These options specify an additive Schwarz relaxation through PatchPC.
-PatchPC builds the patch operators by assembling the bilineary form over
+These options specify an additive Schwarz relaxation through :class:`~.PatchPC`.
+:class:`~.PatchPC` builds the patch operators by assembling the bilineary form over
 each subdomain.  Hence, it does not require the global stiffness
 matrix to be assembled.  These are quite similar to the options used in
 <poisson_mg_patches.py>::
@@ -115,7 +115,7 @@ matrix to be assembled.  These are quite similar to the options used in
            "pc_patch_precompute_element_tensors": None}},
       mat_type="matfree")
 
-ASMStarPC, on the other hand, does no re-discretization, but extracts the
+:class:`~.ASMStarPC`, on the other hand, does no re-discretization, but extracts the
 patch operators for each patch from the already-assembled global stiffness matrix.::
 
   asm_relax = mg_params(
@@ -133,7 +133,7 @@ direct or Krylov method on each one.
 
 Now, for each parameter choice, we report the iteration count for the Poisson problem
 over a range of polynomial degrees.  We see that the Jacobi relaxation leads to growth
-in iteration count, while both PatchPC and ASMStarPC do not.  Mathematically, the two
+in iteration count, while both :class:`~.PatchPC` and :class:`~.ASMStarPC` do not.  Mathematically, the two
 latter options do the same operations, just via different code paths.::
 
   names = {"ASM Vanka": asm_relax,