Gradient recovery by L2-projection (#13)

* Start recovery module
* Add unit tests for recover_gradient_l2

---------

Co-authored-by: Stephan Kramer <s.kramer@imperial.ac.uk>

Author: joewallwork

adapt_common/__init__.py

@@ -1,3 +1,4 @@
 from adapt_common.norms import *  # noqa
+from adapt_common.recovery import *  # noqa
 from adapt_common.reduction import *  # noqa
 from adapt_common.utility import *  # noqa

adapt_common/recovery.py

@@ -0,0 +1,53 @@
+"""Module containing functions for recovering derivatives of Functions."""
+
+import firedrake as fd
+import ufl
+from firedrake.petsc import PETSc
+
+__all__ = ["recover_gradient_l2"]
+
+
+@PETSc.Log.EventDecorator()
+def recover_gradient_l2(f, target_space=None):
+    r"""Recover the gradient of a scalar or vector field using :math:`L^2` projection.
+
+    :arg f: the scalar field whose derivatives we seek to recover
+    :type f: :class:`firedrake.function.Function`
+    :kwarg mesh: the underlying mesh
+    :type mesh: :class:`firedrake.mesh.MeshGeometry`
+    :kwarg target_space: the vector-valued function space to recover the gradient in
+    :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`
+    :returns: recovered gradient
+    :rtype: :class:`firedrake.function.Function`
+    """
+    if target_space is None:
+        if not isinstance(f, fd.Function):
+            val_err = (
+                "If a target space is not provided then the input must be a Function."
+            )
+            raise ValueError(val_err)
+        source_degree = f.ufl_element().degree()
+        if source_degree <= 0:
+            val_err = "Input Function must be at least degree 1."
+            raise ValueError(val_err)
+        target_degree = source_degree - 1
+        mesh = f.function_space().mesh()
+        rank = len(f.function_space().value_shape)
+        family = "DG" if target_degree == 0 else "CG"
+        if rank == 0:
+            target_space = fd.VectorFunctionSpace(mesh, family, target_degree)
+        elif rank == 1:
+            target_space = fd.TensorFunctionSpace(
+                mesh,
+                family,
+                target_degree,
+                shape=(f.function_space().value_size, mesh.geometric_dimension),
+            )
+        else:
+            val_err = (
+                "L2 projection can only be used to compute gradients of scalar or"
+                f" vector Functions, not Functions of rank {rank}."
+            )
+            raise ValueError(val_err)
+
+    return fd.project(ufl.grad(f), target_space)

test/test_recovery_l2.py

@@ -0,0 +1,119 @@
+"""Unit tests for derivative recovery based on L2-projection."""
+
+import firedrake as fd
+import pytest
+
+from adapt_common.recovery import recover_gradient_l2
+
+
+@pytest.fixture(params=[1, 2, 3])
+def dim(request):
+    """Set the spatial dimension."""
+    return request.param
+
+
+@pytest.fixture
+def mesh(dim):
+    """Create a uniform simplex mesh."""
+    n = 4
+    return {
+        1: fd.UnitIntervalMesh(n),
+        2: fd.UnitSquareMesh(n, n),
+        3: fd.UnitCubeMesh(n, n, n),
+    }[dim]
+
+
+def test_recover_gradient_p2_scalar(mesh):
+    """Test gradient recovery for a P2 scalar field."""
+    # Define a scalar function in P2 space
+    f = fd.Function(fd.FunctionSpace(mesh, "CG", 2))
+    x = fd.SpatialCoordinate(mesh)
+    f.interpolate(sum([0.5 * xi**2 for xi in x]))
+
+    # Recovery its gradient using L2 projection
+    grad_f = recover_gradient_l2(f)
+
+    # Check the function space of the recovered gradient
+    element = grad_f.function_space().ufl_element()
+    assert element.family() == "Lagrange"
+    assert element.degree() == 1
+    assert element.num_sub_elements == mesh.geometric_dimension
+
+    # Verify the accuracy of the recovered gradient
+    expected = x
+    assert fd.errornorm(expected, grad_f, norm_type="L2") == pytest.approx(0, abs=1e-8)
+
+
+def test_recover_gradient_p2_vector(mesh):
+    """Test gradient recovery for a P2 vector field."""
+    # Define a vector function in P2 space
+    f = fd.Function(fd.VectorFunctionSpace(mesh, "CG", 2))
+    x = fd.SpatialCoordinate(mesh)
+    f.interpolate(fd.as_vector([0.5 * xi**2 for xi in fd.SpatialCoordinate(mesh)]))
+
+    # Recovery its gradient using L2 projection
+    grad_f = recover_gradient_l2(f)
+
+    # Check the function space of the recovered gradient
+    element = grad_f.function_space().ufl_element()
+    assert element.family() == "Lagrange"
+    assert element.degree() == 1
+    assert element.num_sub_elements == mesh.geometric_dimension**2
+
+    # Verify the accuracy of the recovered gradient
+    expected = fd.Function(fd.TensorFunctionSpace(mesh, "CG", 1))
+    expected.interpolate(
+        fd.as_tensor(
+            [
+                [xi if i == j else 0 for j in range(mesh.geometric_dimension)]
+                for i, xi in enumerate(x)
+            ]
+        )
+    )
+    assert fd.errornorm(expected, grad_f, norm_type="L2") == pytest.approx(0, abs=1e-8)
+
+
+def test_recover_gradient_p1_scalar(mesh):
+    """Test gradient recovery for a P1 scalar field."""
+    # Define a scalar function in P1 space
+    f = fd.Function(fd.FunctionSpace(mesh, "CG", 1))
+    x = fd.SpatialCoordinate(mesh)
+    f.interpolate(sum(x))
+
+    # Recovery its gradient using L2 projection
+    grad_f = recover_gradient_l2(f)
+
+    # Check the function space of the recovered gradient
+    element = grad_f.function_space().ufl_element()
+    assert element.family() == "Discontinuous Lagrange"
+    assert element.degree() == 0
+    assert element.num_sub_elements == mesh.geometric_dimension
+
+    # Verify the accuracy of the recovered gradient
+    expected = fd.Function(grad_f.function_space()).assign(1.0)
+    assert fd.errornorm(expected, grad_f, norm_type="L2") == pytest.approx(0, abs=1e-8)
+
+
+def test_recover_gradient_invalid_input():
+    """Test that an error is raised for invalid input."""
+    val_err = "If a target space is not provided then the input must be a Function."
+    with pytest.raises(ValueError, match=val_err):
+        recover_gradient_l2("not_a_function")
+
+
+def test_recover_gradient_degree_error(mesh):
+    """Test that an error is raised for degree below 1."""
+    val_err = "Input Function must be at least degree 1."
+    with pytest.raises(ValueError, match=val_err):
+        recover_gradient_l2(fd.Function(fd.FunctionSpace(mesh, "DG", 0)))
+
+
+def test_recover_gradient_rank_error(mesh):
+    """Test that an error is raised for unsupported function ranks."""
+    f = fd.Function(fd.TensorFunctionSpace(mesh, "CG", 2))
+    val_err = (
+        "L2 projection can only be used to compute gradients of scalar or vector"
+        " Functions, not Functions of rank 2."
+    )
+    with pytest.raises(ValueError, match=val_err):
+        recover_gradient_l2(f)