Enable mass lumping for projection operator (#115)

Closes #95.

This PR allows the projection operator (and its adjoint) to be mass
lumped. For P1 spaces, this implies that extrema are preserved following
the transfer.

While developing these changes, I noticed #113 and #114, so the
associated functionality and tests are turned off here.

animate/interpolation.py

@@ -73,20 +73,77 @@ def transfer(source, target_space, transfer_method="project", **kwargs):
 
 @PETSc.Log.EventDecorator()
 def interpolate(source, target_space, **kwargs):
-    """
-    A wrapper for :func:`transfer` with ``transfer_method="interpolate"``.
+    r"""
+    Overload function :func:`firedrake.__future__.interpolate` to account for the case
+    of two mixed function spaces defined on different meshes and for the adjoint
+    interpolation operator when applied to :class:`firedrake.cofunction.Cofunction`\s.
+
+    :arg source: the function to be transferred
+    :type source: :class:`firedrake.function.Function` or
+        :class:`firedrake.cofunction.Cofunction`
+    :arg target_space: the function space which we seek to transfer onto, or the
+        function or cofunction to use as the target
+    :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`,
+        :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction`
+    :returns: the transferred function
+    :rtype: :class:`firedrake.function.Function` or
+        :class:`firedrake.cofunction.Cofunction`
+
+    Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate`
     """
     return transfer(source, target_space, transfer_method="interpolate", **kwargs)
 
 
 @PETSc.Log.EventDecorator()
 def project(source, target_space, **kwargs):
-    """
-    A wrapper for :func:`transfer` with ``transfer_method="interpolate"``.
+    r"""
+    Overload function :func:`firedrake.projection.project` to account for the case of
+    two mixed function spaces defined on different meshes and for the adjoint
+    projection operator when applied to :class:`firedrake.cofunction.Cofunction`\s.
+
+    For details on the approach for achieving boundedness through mass lumping and
+    post-processing, see :cite:`FPP+:2009`.
+
+    :arg source: the function to be transferred
+    :type source: :class:`firedrake.function.Function` or
+        :class:`firedrake.cofunction.Cofunction`
+    :arg target_space: the function space which we seek to transfer onto, or the
+        function or cofunction to use as the target
+    :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`,
+        :class:`firedrake.function.Function` or :class:`firedrake.cofunction.Cofunction`
+    :returns: the transferred function
+    :rtype: :class:`firedrake.function.Function` or
+        :class:`firedrake.cofunction.Cofunction`
+    :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness
+    :type bounded: :class:`bool`
+
+    Extra keyword arguments are passed to :func:`firedrake.projection.project`.
     """
     return transfer(source, target_space, transfer_method="project", **kwargs)
 
 
+# TODO: Reimplement by introducing a LumpedSupermeshProjector subclass of
+#       firedrake.projection.SupermeshProjector (#123)
+# TODO: Implement minimal diffusion correction (#124)
+def _supermesh_project(source, target, bounded=False):
+    Vs = source.function_space()
+    Vt = target.function_space()
+    element_t = Vt.ufl_element()
+    if bounded and (element_t.family(), element_t.degree()) != ("Lagrange", 1):
+        raise ValueError("Mass lumping is not recommended for spaces other than P1.")
+
+    # Create a linear system using the lumped mass matrix for the target space
+    mixed_mass = assemble_mixed_mass_matrix(Vs, Vt)
+    ksp = petsc4py.KSP().create()
+    ksp.setOperators(assemble_mass_matrix(Vt, lumped=bounded))
+
+    # Solve the linear system
+    with source.dat.vec_ro as s, target.dat.vec_wo as t:
+        rhs = t.copy()
+        mixed_mass.mult(s, rhs)
+        ksp.solve(rhs, t)
+
+
 @PETSc.Log.EventDecorator()
 def _transfer_forward(source, target, transfer_method, **kwargs):
     """
@@ -102,12 +159,17 @@ def _transfer_forward(source, target, transfer_method, **kwargs):
     :kwarg transfer_method: the method to use for the transfer. Options are
         "interpolate" (default) and "project"
     :type transfer_method: str
+    :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness
+        (project only)
+    :type bounded: :class:`bool`
     :returns: the transferred Function
     :rtype: :class:`firedrake.function.Function`
 
     Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` or
         :func:`firedrake.projection.project`.
     """
+    is_project = transfer_method == "project"
+    bounded = is_project and kwargs.pop("bounded", False)
     Vs = source.function_space()
     Vt = target.function_space()
     _validate_matching_spaces(Vs, Vt)
@@ -117,7 +179,10 @@ def _transfer_forward(source, target, transfer_method, **kwargs):
             if transfer_method == "interpolate":
                 t.interpolate(s, **kwargs)
             elif transfer_method == "project":
-                t.project(s, **kwargs)
+                if bounded:
+                    _supermesh_project(s, t, bounded=True)
+                else:
+                    t.project(s, **kwargs)
             else:
                 raise ValueError(
                     f"Invalid transfer method: {transfer_method}."
@@ -127,7 +192,10 @@ def _transfer_forward(source, target, transfer_method, **kwargs):
         if transfer_method == "interpolate":
             target.interpolate(source, **kwargs)
         elif transfer_method == "project":
-            target.project(source, **kwargs)
+            if bounded:
+                _supermesh_project(source, target, bounded=True)
+            else:
+                target.project(source, **kwargs)
         else:
             raise ValueError(
                 f"Invalid transfer method: {transfer_method}."
@@ -148,12 +216,17 @@ def _transfer_adjoint(target_b, source_b, transfer_method, **kwargs):
     :kwarg transfer_method: the method to use for the transfer. Options are
         "interpolate" (default) and "project"
     :type transfer_method: str
+    :kwarg bounded: apply mass lumping to the mass matrix to ensure boundedness
+        (project only)
+    :type bounded: :class:`bool`
     :returns: the transferred Cofunction
     :rtype: :class:`firedrake.cofunction.Cofunction`
 
     Extra keyword arguments are passed to :func:`firedrake.__future__.interpolate` or
         :func:`firedrake.projection.project`.
     """
+    is_project = transfer_method == "project"
+    bounded = is_project and kwargs.pop("bounded", False)
 
     # Map to Functions to apply the adjoint transfer
     if not isinstance(target_b, firedrake.Function):
@@ -178,10 +251,12 @@ def _transfer_adjoint(target_b, source_b, transfer_method, **kwargs):
     # Apply adjoint transfer operator to each component
     for i, (t_b, s_b) in enumerate(zip(target_b_split, source_b_split)):
         if transfer_method == "interpolate":
-            s_b.interpolate(t_b, **kwargs)
+            raise NotImplementedError(
+                "Adjoint of interpolation operator not implemented."
+            )  # TODO (#113)
         elif transfer_method == "project":
             ksp = petsc4py.KSP().create()
-            ksp.setOperators(assemble_mass_matrix(t_b.function_space()))
+            ksp.setOperators(assemble_mass_matrix(t_b.function_space(), lumped=bounded))
             mixed_mass = assemble_mixed_mass_matrix(Vt[i], Vs[i])
             with t_b.dat.vec_ro as tb, s_b.dat.vec_wo as sb:
                 residual = tb.copy()

animate/utility.py

@@ -5,6 +5,7 @@
 from collections import OrderedDict
 
 import firedrake
+import firedrake.function as ffunc
 import firedrake.mesh as fmesh
 import ufl
 from firedrake.__future__ import interpolate
@@ -43,7 +44,7 @@ def Mesh(arg, **kwargs):
 
     # Facet area
     boundary_markers = sorted(mesh.exterior_facets.unique_markers)
-    one = firedrake.Function(P1).assign(1.0)
+    one = ffunc.Function(P1).assign(1.0)
     bnd_len = OrderedDict(
         {i: firedrake.assemble(one * ufl.ds(int(i))) for i in boundary_markers}
     )
@@ -202,16 +203,16 @@ def errornorm(u, uh, norm_type="L2", boundary=False, **kwargs):
         u = cofunction2function(u)
     if isinstance(uh, firedrake.Cofunction):
         uh = cofunction2function(uh)
-    if not isinstance(uh, firedrake.Function):
+    if not isinstance(uh, ffunc.Function):
         raise TypeError(f"uh should be a Function, is a '{type(uh)}'.")
     if norm_type[0] == "l":
-        if not isinstance(u, firedrake.Function):
+        if not isinstance(u, ffunc.Function):
             raise TypeError(f"u should be a Function, is a '{type(u)}'.")
 
     if len(u.ufl_shape) != len(uh.ufl_shape):
         raise RuntimeError("Mismatching rank between u and uh.")
 
-    if isinstance(u, firedrake.Function):
+    if isinstance(u, ffunc.Function):
         degree_u = u.function_space().ufl_element().degree()
         degree_uh = uh.function_space().ufl_element().degree()
         if degree_uh > degree_u:
@@ -243,14 +244,16 @@ def errornorm(u, uh, norm_type="L2", boundary=False, **kwargs):
 
 
 @PETSc.Log.EventDecorator()
-def assemble_mass_matrix(space, norm_type="L2"):
+def assemble_mass_matrix(space, norm_type="L2", lumped=False):
     """
     Assemble a mass matrix associated with some finite element space and norm.
 
     :arg space: function space to build the mass matrix with
     :type space: :class:`firedrake.functionspaceimpl.FunctionSpace`
     :kwarg norm_type: the type norm to build the mass matrix with
     :type norm_type: :class:`str`
+    :kwarg lumped: if `True`, mass lumping is applied
+    :type lumped: :class:`bool`
     :returns: the corresponding mass matrix
     :rtype: petsc4py.PETSc.Mat
     """
@@ -265,7 +268,14 @@ def assemble_mass_matrix(space, norm_type="L2"):
         )
     else:
         raise ValueError(f"Norm type '{norm_type}' not recognised.")
-    return firedrake.assemble(lhs).petscmat
+    mass_matrix = firedrake.assemble(lhs).petscmat
+    if not lumped:
+        return mass_matrix
+    rhs = ffunc.Function(space).assign(1.0)
+    product = ffunc.Function(space)
+    with rhs.dat.vec_ro as b, product.dat.vec as x:
+        mass_matrix.mult(b, x)
+        return mass_matrix.createDiagonal(x)
 
 
 def cofunction2function(cofunc):
@@ -275,7 +285,7 @@ def cofunction2function(cofunc):
     :returns: a function with the same underyling data
     :rtype: :class:`firedrake.function.Function`
     """
-    func = firedrake.Function(cofunc.function_space().dual())
+    func = ffunc.Function(cofunc.function_space().dual())
     if isinstance(func.dat.data_with_halos, tuple):
         for i, arr in enumerate(func.dat.data_with_halos):
             arr[:] = cofunc.dat.data_with_halos[i]

demos/combining_metrics.py

@@ -67,10 +67,10 @@
 axes[1].set_aspect("equal")
 axes[1].set_title("Mesh based on metric2")
 fig.show()
-fig.savefig("combining_two_metrics.jpg")
+fig.savefig("combining_metrics-inputs.jpg")
 
 #
-# .. figure:: combining_two_metrics.jpg
+# .. figure:: combining_metrics-inputs.jpg
 #    :figwidth: 90%
 #    :align: center
 #
@@ -88,10 +88,10 @@
 axes.set_aspect("equal")
 axes.set_title("Mesh based on intersected metric")
 fig.show()
-fig.savefig("combining_intersection.jpg")
+fig.savefig("combining_metrics-intersection.jpg")
 
 #
-# .. figure:: combining_intersection.jpg
+# .. figure:: combining_metrics-intersection.jpg
 #    :figwidth: 90%
 #    :align: center
 #
@@ -112,12 +112,14 @@
 axes.set_aspect("equal")
 axes.set_title("Mesh based on averaged metric")
 fig.show()
-fig.savefig("combining_averaging.jpg")
+fig.savefig("combining_metrics-averaging.jpg")
 
 #
-# .. figure:: combining_averaging.jpg
+# .. figure:: combining_metrics-averaging.jpg
 #    :figwidth: 90%
 #    :align: center
 #
 # the resolution in, e.g., the region :math:`x<0.5, y<0.3` is now based on an
 # average of :math:`hm=0.02` and :math:`hc=0.1`, i.e. an edge length of 0.06.
+#
+# This demo can also be accessed as a `Python script <combining_metrics.py>`__.