Merge branch 'master' into indiamai/new_def_integration

demos/poisson/poisson_mixed.py.rst

@@ -171,5 +171,5 @@ Don't forget to show the image::
     warning("Cannot show figure. Error msg '%s'" % e)
 
 This demo is based on the corresponding `DOLFIN mixed Poisson demo
-<http://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/mixed-poisson/python/documentation.html>`__
+<https://olddocs.fenicsproject.org/dolfin/1.3.0/python/demo/documented/mixed-poisson/python/documentation.html>`__
 and can be found as a script in :demo:`poisson_mixed.py <poisson_mixed.py>`.

docs/source/team.ini

@@ -49,6 +49,7 @@ Patrick E. Farrell: https://pefarrell.org, https://pefarrell.org/images/profile.
 Pablo D. Brubeck: https://www.maths.ox.ac.uk/people/pablo.brubeckmartinez, https://avatars.githubusercontent.com/u/6486130
 India Marsden: https://www.maths.ox.ac.uk/people/india.marsden, https://pefarrell.org/images/marsden.jpg
 Daiane I. Dolci: https://www.imperial.ac.uk/people/d.dolci, https://avatars.githubusercontent.com/u/63597005?v=4
+Joshua Hope-Collins: https://www.imperial.ac.uk/people/joshua.hope-collins13
 
 [inactive-team]
 Sophia Vorderwuelbecke: https://www.imperial.ac.uk/people/s.vorderwuelbecke18

firedrake/adjoint/__init__.py

@@ -21,10 +21,11 @@
                             pause_annotation, continue_annotation, \
                             stop_annotating, annotate_tape  # noqa F401
 from pyadjoint.reduced_functional import ReducedFunctional  # noqa F401
+from pyadjoint.checkpointing import disk_checkpointing_callback  # noqa F401
 from firedrake.adjoint_utils.checkpointing import \
     enable_disk_checkpointing, pause_disk_checkpointing, \
     continue_disk_checkpointing, stop_disk_checkpointing, \
-    checkpointable_mesh   # noqa F401
+    checkpointable_mesh  # noqa F401
 from firedrake.adjoint_utils import get_solve_blocks  # noqa F401
 
 from pyadjoint.verification import taylor_test, taylor_to_dict  # noqa F401

firedrake/adjoint_utils/checkpointing.py

@@ -1,5 +1,5 @@
 """A module providing support for disk checkpointing of the adjoint tape."""
-from pyadjoint import get_working_tape, OverloadedType
+from pyadjoint import get_working_tape, OverloadedType, disk_checkpointing_callback
 from pyadjoint.tape import TapePackageData
 from pyop2.mpi import COMM_WORLD
 import tempfile
@@ -10,6 +10,8 @@
 from numbers import Number
 _enable_disk_checkpoint = False
 _checkpoint_init_data = False
+disk_checkpointing_callback["firedrake"] = "Please call enable_disk_checkpointing() "\
+    "before checkpointing on the disk."
 
 __all__ = ["enable_disk_checkpointing", "disk_checkpointing",
            "pause_disk_checkpointing", "continue_disk_checkpointing",
@@ -204,6 +206,12 @@ def restore_from_checkpoint(self, state):
         self.init_checkpoint_file = state["init"]
         self.current_checkpoint_file = state["current"]
 
+    def continue_checkpointing(self):
+        continue_disk_checkpointing()
+
+    def pause_checkpointing(self):
+        pause_disk_checkpointing()
+
 
 def checkpointable_mesh(mesh):
     """Write a mesh to disk and read it back.
@@ -251,7 +259,7 @@ def restore(self):
         pass
 
 
-class CheckpointFunction(CheckpointBase):
+class CheckpointFunction(CheckpointBase, OverloadedType):
     """Metadata for a Function checkpointed to disk.
 
     An object of this class replaces the :class:`~firedrake.Function` stored as
@@ -304,6 +312,9 @@ def restore(self):
         return type(function)(function.function_space(),
                               function.dat, name=self.name(), count=self.count)
 
+    def _ad_restore_at_checkpoint(self, checkpoint):
+        return checkpoint.restore()
+
 
 def maybe_disk_checkpoint(function):
     """Checkpoint a Function to disk if disk checkpointing is active."""

scripts/firedrake-install

@@ -2073,6 +2073,7 @@ with environment(**compiler_env):
 
     if options["netgen"]:
         packages += ["ngsPETSc"]
+        run_pip(["install", "-U", "ngsPETSc"])
 
 with pipargs("--no-deps"):
     if options["opencascade"]:

tests/regression/test_assemble.py

@@ -334,3 +334,13 @@ def test_assemble_sparsity_diagonal_entries_for_bc():
     A = assemble(inner(u[1], v[0]) * dx, bcs=[bc], mat_type="nest")
     # Make sure that diagonals are allocated.
     assert np.all(A.M.sparsity[1][1].nnz == np.ones(4, dtype=IntType))
+
+
+@pytest.mark.skipcomplex
+def test_assemble_power_zero_minmax():
+    mesh = UnitSquareMesh(1, 1)
+    V = FunctionSpace(mesh, "CG", 1)
+    f = Function(V).assign(1.)
+    g = Function(V).assign(2.)
+    assert assemble(zero()**min_value(f, g) * dx) == 0.0
+    assert assemble(zero()**max_value(f, g) * dx) == 0.0

tests/regression/test_interpolate_cross_mesh.py

@@ -215,7 +215,7 @@ def parameters(request):
         expected = np.asarray(
             [np.outer(coords[i], coords[i]) for i in range(len(coords))]
         )
-        V_src = FunctionSpace(m_src, "Regge", 2)  # Not point evaluation nodes
+        V_src = FunctionSpace(m_src, "Regge", 2, variant="point")  # Not point evaluation nodes
         V_dest = TensorFunctionSpace(m_dest, "CG", 4)
         V_dest_2 = TensorFunctionSpace(m_dest, "DQ", 2)
     elif request.param == "spheresphere":

tests/regression/test_stress_els.py → tests/regression/test_stress_elements.py

@@ -7,19 +7,21 @@
 convergence_orders = lambda x: -np.log2(relative_magnitudes(x))
 
 
-@pytest.fixture(scope='module', params=["conforming", "nonconforming", "macro"])
+@pytest.fixture(params=["macro", "nonconforming", "conforming", "high-order"])
 def stress_element(request):
-    if request.param == "conforming":
-        return FiniteElement("AWc", triangle, 3)
+    if request.param == "macro":
+        return FiniteElement("JM", triangle, 1)
     elif request.param == "nonconforming":
         return FiniteElement("AWnc", triangle, 2)
-    elif request.param == "macro":
-        return FiniteElement("JM", triangle, 1)
+    elif request.param == "conforming":
+        return FiniteElement("AWc", triangle, 3)
+    elif request.param == "high-order":
+        return FiniteElement("HZ", triangle, 3)
     else:
         raise ValueError("Unknown family")
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture
 def mesh_hierarchy(request):
     N_base = 2
     mesh = UnitSquareMesh(N_base, N_base)
@@ -64,7 +66,10 @@ def epsilon(u):
     l2_u = []
     l2_sigma = []
     l2_div_sigma = []
-    displacement_element = VectorElement("DG", cell=mesh.ufl_cell(), degree=1, variant="integral")
+    sdegree = stress_element.degree()
+    sfamily = stress_element.family()
+    vdegree = sdegree-1 if sfamily == "Hu-Zhang" else 1
+    displacement_element = VectorElement("DG", cell=mesh.ufl_cell(), degree=vdegree, variant="integral")
     element = MixedElement([stress_element, displacement_element])
     for msh in mesh_hierarchy[1:]:
         x, y = SpatialCoordinate(msh)
@@ -126,20 +131,18 @@ def epsilon(u):
         l2_sigma.append(error_sigma)
         l2_div_sigma.append(error_div_sigma)
 
-    if stress_element.family().startswith("Conforming"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 2.9
-        assert min(convergence_orders(l2_div_sigma)) > 1.9
-    elif stress_element.family().startswith("Nonconforming"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 1
-        assert min(convergence_orders(l2_div_sigma)) > 1.9
-    elif stress_element.family().startswith("Johnson-Mercier"):
-        assert min(convergence_orders(l2_u)) > 1.9
-        assert min(convergence_orders(l2_sigma)) > 1
-        assert min(convergence_orders(l2_div_sigma)) > 0.9
+    sdegree = V[0].ufl_element().degree()
+    vdegree = V[1].ufl_element().degree()
+    if stress_element.family().startswith("Nonconforming"):
+        expected_rates = (vdegree+1, sdegree-1, sdegree)
+    elif stress_element.family().startswith("Conforming"):
+        expected_rates = (vdegree+1, sdegree, sdegree-1)
     else:
-        raise ValueError("Don't know what the convergence should be")
+        expected_rates = (vdegree+1, sdegree+1, sdegree)
+
+    assert min(convergence_orders(l2_u[1:])) > expected_rates[0] * 0.9
+    assert min(convergence_orders(l2_sigma[1:])) > expected_rates[1] * 0.9
+    assert min(convergence_orders(l2_div_sigma[1:])) > expected_rates[2] * 0.9
 
 
 def test_mass_conditioning(stress_element, mesh_hierarchy):

tests/regression/test_tensor_elements.py

@@ -0,0 +1,51 @@
+from firedrake import *
+import pytest
+
+
+rg = RandomGenerator(PCG64(seed=1))
+
+
+@pytest.fixture(params=("square", "cube"))
+def mesh(request):
+    if request.param == "square":
+        return UnitSquareMesh(2, 2)
+    elif request.param == "cube":
+        return UnitCubeMesh(2, 2, 2)
+    else:
+        raise ValueError("Unrecognized mesh type")
+
+
+@pytest.mark.parametrize("family, degree", [
+                         ("Regge", 0),
+                         ("Regge", 1),
+                         ("Regge", 2),
+                         ("HHJ", 0),
+                         ("HHJ", 1),
+                         ("HHJ", 2),
+                         ("GLS", 1),
+                         ("GLS", 2),
+                         ])
+def test_tensor_continuity(mesh, family, degree):
+    V = FunctionSpace(mesh, family, degree)
+    u = rg.beta(V, 1.0, 2.0)
+
+    n = FacetNormal(mesh)
+
+    space = V.ufl_element().sobolev_space
+    if space == HDivDiv:
+        utrace = dot(n, dot(u, n))
+    elif space == HEin:
+        if mesh.topological_dimension() == 2:
+            t = perp(n)
+        else:
+            t = as_matrix([[0, n[2], -n[1]], [-n[2], 0, n[0]], [n[1], -n[0], 0]])
+        utrace = dot(t, dot(u, t))
+    else:
+        if mesh.topological_dimension() == 2:
+            t = perp(n)
+            utrace = dot(t, dot(u, n))
+        else:
+            utrace = cross(n, dot(u, n))
+
+    ujump = utrace('+') - utrace('-')
+    assert assemble(inner(ujump, ujump)*dS) < 1E-10

tests/vertexonly/test_interpolation_from_parent.py

@@ -102,7 +102,7 @@ def vfs(request, parentmesh):
 
 @pytest.fixture(params=[("CG", 2, TensorFunctionSpace),
                         ("BDM", 2, VectorFunctionSpace),
-                        ("Regge", 2, FunctionSpace)],
+                        ("Regge", 2, lambda *args: FunctionSpace(*args, variant="point"))],
                 ids=lambda x: f"{x[2].__name__}({x[0]}{x[1]})")
 def tfs(request, parentmesh):
     family = request.param[0]