DO NOT SQUASH: merge release into main

demos/multicomponent/multicomponent.py.rst

@@ -428,7 +428,15 @@ only apply to the normal component of H(div) functions.
 Elsewhere on the boundary we enforce :math:`J_i \cdot N = 0`. Finally, instead of specifying
 the value of the barycentric velocity :math:`v` on the inflows and outflows,
 we enforce :math:`v = \rho^{-1}(G_1 + G_2)`. Boundary conditions that couple
-unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` instead of :class:`~.DirichletBC`. ::
+unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` instead of :class:`~.DirichletBC`.
+Since the barycentric velocity is in :math:`[\textrm{CG}_k]^2` it has degrees of freedom located at the points
+where the inlets/outlet meet the walls. Even though the boundary conditions on the inlets/outlet
+and the walls are compatible at these points, :class:`~.EquationBC` enforces the boundary conditions weakly whilst
+:class:`~.DirichletBC` enforces them strongly. Hence, to avoid ambiguity, we set Dirichlet boundary
+conditions on the boundary of the boundary, i.e. the points where the inlets/outlet meet the walls.
+To enforce these Dirichlet boundary conditions, tuples with the numbers of the boundary edges coincidental
+to these points need to be constructed first. This is then passed on to a Dirichlet boundary condition
+which is passed on to :class:`~.EquationBC`.::
 
     # Reference species velocities, which we choose to symmetrize so that the molar fluxes agree
     v_ref_1 = Constant(0.4e-6)                      # Reference inflow velocity of benzene, m / s
@@ -449,7 +457,9 @@ unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` ins
     # Boundary conditions on the barycentric velocity are enforced via EquationBC
     bc_data = {inlet_1_id: rho_v_inflow_1_bc_func, inlet_2_id: rho_v_inflow_2_bc_func, outlet_id: rho_v_outflow_bc_func}
     F_bc = sum(inner(v - rho_inv * flux, u) * ds(*subdomain) for subdomain, flux in bc_data.items())
-    v_bc = EquationBC(F_bc == 0, solution, (*inlet_1_id, *inlet_2_id, *outlet_id), V=Z_h.sub(2))
+    ridges = [(i, j) for i in [*inlet_1_id, *inlet_2_id, *outlet_id] for j in walls_ids]
+    sub_bcs = [DirichletBC(Z_h.sub(2), 0, ridges)]  # bc on the boundary of the boundary
+    v_bc = EquationBC(F_bc == 0, solution, (*inlet_1_id, *inlet_2_id, *outlet_id), V=Z_h.sub(2), bcs=sub_bcs)
 
     # The boundary conditions on the fluxes and barycentric velocity
     # Note that BCs on H(div) spaces only apply to the normal component

firedrake/adjoint_utils/function.py

@@ -270,7 +270,7 @@ def _ad_dot(self, other, options=None):
             return assemble(firedrake.inner(self, other)*firedrake.dx)
         elif riesz_representation == "H1":
             return assemble((firedrake.inner(self, other)
-                            + firedrake.inner(firedrake.grad(self), other))*firedrake.dx)
+                            + firedrake.inner(firedrake.grad(self), firedrake.grad(other)))*firedrake.dx)
         else:
             raise NotImplementedError(
                 "Unknown Riesz representation %s" % riesz_representation)

firedrake/assemble.py

@@ -1255,6 +1255,9 @@ def _apply_bc(self, tensor, bc, u=None):
                 bc.apply(r, u=u)
         elif isinstance(bc, EquationBCSplit):
             bc.zero(tensor)
+            if isinstance(bc.f, ufl.ZeroBaseForm) or bc.f.empty():
+                # form is empty, do nothing
+                return
             OneFormAssembler(bc.f, bcs=bc.bcs,
                              form_compiler_parameters=self._form_compiler_params,
                              needs_zeroing=False,

firedrake/bcs.py

@@ -639,9 +639,6 @@ def reconstruct(self, field=None, V=None, subu=None, u=None, row_field=None, col
             rank = len(self.f.arguments())
             splitter = ExtractSubBlock()
             form = splitter.split(self.f, argument_indices=(row_field, col_field)[:rank])
-            if isinstance(form, ufl.ZeroBaseForm) or form.empty():
-                # form is empty, do nothing
-                return
             if u is not None:
                 form = firedrake.replace(form, {self.u: u})
         if action_x is not None:

firedrake/fml/form_manipulation_language.py

@@ -93,7 +93,7 @@ class Term(object):
 
     __slots__ = ["form", "labels"]
 
-    def __init__(self, form: ufl.Form, label_dict: Mapping = None):
+    def __init__(self, form: ufl.BaseForm, label_dict: Mapping = None):
         """
 
         Parameters
@@ -216,6 +216,9 @@ def __mul__(
 
     __rmul__ = __mul__
 
+    def __neg__(self):
+        return -1 * self
+
     def __truediv__(
         self,
         other: Union[float, Constant, ufl.algebra.Product]
@@ -274,7 +277,7 @@ def __init__(self, *terms: Sequence[Term]):
 
     def __add__(
         self,
-        other: Union[ufl.Form, Term, "LabelledForm"]
+        other: Union[ufl.BaseForm, Term, "LabelledForm"]
     ) -> "LabelledForm":
         """Add a form, term or labelled form to this labelled form.
 
@@ -289,12 +292,12 @@ def __add__(
             A labelled form containing the terms.
 
         """
-        if isinstance(other, ufl.Form):
-            return LabelledForm(*self, Term(other))
-        elif type(other) is Term:
+        if type(other) is Term:
             return LabelledForm(*self, other)
         elif type(other) is LabelledForm:
             return LabelledForm(*self, *other)
+        elif isinstance(other, ufl.BaseForm):
+            return LabelledForm(*self, Term(other))
         elif other is None:
             return self
         else:
@@ -304,7 +307,7 @@ def __add__(
 
     def __sub__(
         self,
-        other: Union[ufl.Form, Term, "LabelledForm"]
+        other: Union[ufl.BaseForm, Term, "LabelledForm"]
     ) -> "LabelledForm":
         """Subtract a form, term or labelled form from this labelled form.
 
@@ -320,15 +323,18 @@ def __sub__(
 
         """
         if type(other) is Term:
-            return LabelledForm(*self, Constant(-1.)*other)
+            return LabelledForm(*self, -other)
         elif type(other) is LabelledForm:
-            return LabelledForm(*self, *[Constant(-1.)*t for t in other])
+            return LabelledForm(*self, *[-t for t in other])
         elif other is None:
             return self
         else:
             # Make new Term for other and subtract it
             return LabelledForm(*self, Term(Constant(-1.)*other))
 
+    def __rsub__(self, other):
+        return other + (-self)
+
     def __mul__(
         self,
         other: Union[float, Constant, ufl.algebra.Product]
@@ -371,6 +377,9 @@ def __truediv__(
 
     __rmul__ = __mul__
 
+    def __neg__(self):
+        return -1 * self
+
     def __iter__(self) -> Sequence:
         """Iterable of the terms in the labelled form."""
         return iter(self.terms)
@@ -427,7 +436,7 @@ def label_map(
         return new_labelled_form
 
     @property
-    def form(self) -> ufl.Form:
+    def form(self) -> ufl.BaseForm:
         """Provide the whole form from the labelled form.
 
         Raises
@@ -437,7 +446,7 @@ def form(self) -> ufl.Form:
 
         Returns
         -------
-        ufl.Form
+        ufl.BaseForm
             The whole form corresponding to all the terms.
 
         """
@@ -479,7 +488,7 @@ def __init__(
 
     def __call__(
         self,
-        target: Union[ufl.Form, Term, LabelledForm],
+        target: Union[ufl.BaseForm, Term, LabelledForm],
         value: Any = None
     ) -> Union[Term, LabelledForm]:
         """Apply the label to a form or term.
@@ -494,7 +503,7 @@ def __call__(
         Raises
         ------
         ValueError
-            If the `target` is not a ufl.Form, Term or
+            If the `target` is not a ufl.BaseForm, Term or
             LabelledForm.
 
         Returns
@@ -514,7 +523,7 @@ def __call__(
             self.value = self.default_value
         if isinstance(target, LabelledForm):
             return LabelledForm(*(self(t, value) for t in target.terms))
-        elif isinstance(target, ufl.Form):
+        elif isinstance(target, ufl.BaseForm):
             return LabelledForm(Term(target, {self.label: self.value}))
         elif isinstance(target, Term):
             new_labels = target.labels.copy()

firedrake/interpolation.py

@@ -441,7 +441,11 @@ def __init__(self, expr: Interpolate):
         else:
             self.access = op2.WRITE
 
-        if self.target_space.ufl_element().mapping() != "identity":
+        # TODO check V.finat_element.is_lagrange() once https://github.com/firedrakeproject/fiat/pull/200 is released
+        target_element = self.target_space.ufl_element()
+        if not ((isinstance(target_element, MixedElement)
+                 and all(sub.mapping() == "identity" for sub in target_element.sub_elements))
+                or target_element.mapping() == "identity"):
             # Identity mapping between reference cell and physical coordinates
             # implies point evaluation nodes.
             raise NotImplementedError(
@@ -510,7 +514,8 @@ def _get_symbolic_expressions(self) -> tuple[Interpolate, Interpolate]:
         elif len(shape) == 1:
             fs_type = partial(VectorFunctionSpace, dim=shape[0])
         else:
-            fs_type = partial(TensorFunctionSpace, shape=shape)
+            symmetry = self.target_space.ufl_element().symmetry()
+            fs_type = partial(TensorFunctionSpace, shape=shape, symmetry=symmetry)
 
         # Get expression for point evaluation at the dest_node_coords
         P0DG_vom = fs_type(vom, "DG", 0)

firedrake/mg/embedded.py

@@ -10,7 +10,7 @@
 __all__ = ("TransferManager", )
 
 
-native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement", "Crouzeix-Raviart"])
+native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement", "Crouzeix-Raviart", "Kong-Mulder-Veldhuizen"])
 alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced-Hsieh-Clough-Tocher", "Johnson-Mercier",
                              "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu",
                              "Guzman-Neilan", "Guzman-Neilan Bubble"])

firedrake/scripts/firedrake-zenodo

@@ -26,14 +26,12 @@ DESCRIPTIONS = {
     "ufl": "UFL: the Unified Form Language",
     "fiat": "FIAT: the Finite Element Automated Tabulator",
     "petsc": "PETSc: the Portable, Extensible Toolkit for Scientific Computation",
-    "loopy": "loopy: Transformation-Based Generation of High-Performance CPU/GPU Code",
 }
 
 PYPI_PACKAGE_NAMES = {
     "firedrake": "firedrake",
     "ufl": "fenics-ufl",
     "fiat": "firedrake-fiat",
-    "loopy": "loopy",
     "petsc": "petsc4py",
 }
 
@@ -295,11 +293,11 @@ def zenodo_records():
     result = None
     i = 1
     while True:
-        if i > 20:
+        if i > 8000//25:
             raise RuntimeError("More than 8000 uploads on zenodo?")
         response = requests.get(f"{ZENODO_URL}/records", params={"q": 'owners:19586 OR owners:19587',
                                                                  "all_versions": True,
-                                                                 "size": 400,
+                                                                 "size": 25,
                                                                  "sort": "mostrecent",
                                                                  "page": i})
         if response.status_code == 200:

firedrake/supermeshing.py

@@ -20,15 +20,18 @@
 from pyop2.compilation import load
 from pyop2.mpi import COMM_SELF
 from pyop2.utils import get_petsc_dir
+from collections import defaultdict
 
 
 __all__ = ["assemble_mixed_mass_matrix", "intersection_finder"]
 
 
+# TODO replace with KAIJ (we require petsc4py wrappers)
 class BlockMatrix(object):
-    def __init__(self, mat, dimension):
+    def __init__(self, mat, dimension, block_scale=None):
         self.mat = mat
         self.dimension = dimension
+        self.block_scale = block_scale
 
     def mult(self, mat, x, y):
         sizes = self.mat.getSizes()
@@ -41,6 +44,8 @@ def mult(self, mat, x, y):
             xi = PETSc.Vec().createWithArray(xa, size=sizes[1], comm=x.comm)
             yi = PETSc.Vec().createWithArray(ya, size=sizes[0], comm=y.comm)
             self.mat.mult(xi, yi)
+            if self.block_scale is not None:
+                yi.scale(self.block_scale[i])
             y.array[start::stride] = yi.array_r
 
     def multTranspose(self, mat, x, y):
@@ -54,6 +59,8 @@ def multTranspose(self, mat, x, y):
             xi = PETSc.Vec().createWithArray(xa, size=sizes[0], comm=x.comm)
             yi = PETSc.Vec().createWithArray(ya, size=sizes[1], comm=y.comm)
             self.mat.multTranspose(xi, yi)
+            if self.block_scale is not None:
+                yi.scale(self.block_scale[i])
             y.array[start::stride] = yi.array_r
 
 
@@ -68,14 +75,6 @@ def assemble_mixed_mass_matrix(V_A, V_B):
     if len(V_A) > 1 or len(V_B) > 1:
         raise NotImplementedError("Sorry, only implemented for non-mixed spaces")
 
-    if V_A.ufl_element().mapping() != "identity" or V_B.ufl_element().mapping() != "identity":
-        msg = """
-Sorry, only implemented for affine maps for now. To do non-affine, we'd need to
-import much more of the assembly engine of UFL/TSFC/etc to do the assembly on
-each supermesh cell.
-"""
-        raise NotImplementedError(msg)
-
     mesh_A = V_A.mesh()
     mesh_B = V_B.mesh()
 
@@ -116,15 +115,39 @@ def likely(cell_A):
                 def likely(cell_A):
                     return cell_map[cell_A]
 
-    assert V_A.value_size == V_B.value_size
-    orig_value_size = V_A.value_size
-    if V_A.value_size > 1:
+    assert V_A.block_size == V_B.block_size
+    orig_block_size = V_A.block_size
+
+    # To deal with symmetry, each block of the mass matrix must be rescaled by the multiplicity
+    if V_A.ufl_element().mapping() == "symmetries":
+        symmetry = V_A.ufl_element().symmetry()
+        assert V_B.ufl_element().mapping() == "symmetries"
+        assert V_B.ufl_element().symmetry() == symmetry
+
+        multiplicity = defaultdict(int)
+        for idx in numpy.ndindex(V_A.value_shape):
+            idx = symmetry.get(idx, idx)
+            multiplicity[idx] += 1
+
+        block_scale = tuple(scale for idx, scale in multiplicity.items())
+    else:
+        block_scale = None
+
+    if V_A.block_size > 1:
         V_A = firedrake.FunctionSpace(mesh_A, V_A.ufl_element().sub_elements[0])
-    if V_B.value_size > 1:
+    if V_B.block_size > 1:
         V_B = firedrake.FunctionSpace(mesh_B, V_B.ufl_element().sub_elements[0])
 
-    assert V_A.value_size == 1
-    assert V_B.value_size == 1
+    if V_A.ufl_element().mapping() != "identity" or V_B.ufl_element().mapping() != "identity":
+        msg = """
+Sorry, only implemented for affine maps for now. To do non-affine, we'd need to
+import much more of the assembly engine of UFL/TSFC/etc to do the assembly on
+each supermesh cell.
+"""
+        raise NotImplementedError(msg)
+
+    assert V_A.block_size == 1
+    assert V_B.block_size == 1
 
     preallocator = PETSc.Mat().create(comm=mesh_A.comm)
     preallocator.setType(PETSc.Mat.Type.PREALLOCATOR)
@@ -155,7 +178,7 @@ def likely(cell_A):
     onnz = numpy.repeat(onnz, cset.cdim)
     preallocator.destroy()
 
-    assert V_A.value_size == V_B.value_size
+    assert V_A.block_size == V_B.block_size
     rdim = V_B.dof_dset.cdim
     cdim = V_A.dof_dset.cdim
 
@@ -445,16 +468,16 @@ def likely(cell_A):
     lib.restype = ctypes.c_int
 
     ammm(V_A, V_B, likely, node_locations_A, node_locations_B, M_SS, ctypes.addressof(lib), mat)
-    if orig_value_size == 1:
+    if orig_block_size == 1:
         return mat
     else:
         (lrows, grows), (lcols, gcols) = mat.getSizes()
-        lrows *= orig_value_size
-        grows *= orig_value_size
-        lcols *= orig_value_size
-        gcols *= orig_value_size
+        lrows *= orig_block_size
+        grows *= orig_block_size
+        lcols *= orig_block_size
+        gcols *= orig_block_size
         size = ((lrows, grows), (lcols, gcols))
-        context = BlockMatrix(mat, orig_value_size)
+        context = BlockMatrix(mat, orig_block_size, block_scale=block_scale)
         blockmat = PETSc.Mat().createPython(size, context=context, comm=mat.comm)
         blockmat.setUp()
         return blockmat

scripts/firedrake-check

@@ -24,6 +24,10 @@ TESTS = {
         # near nullspace
         "tests/firedrake/regression/test_nullspace.py::test_near_nullspace",
     ),
+    2: (
+        # HDF5/checkpointing
+        "tests/firedrake/output/test_io_function.py::test_io_function_base[cell_family_degree2]",
+    ),
     3: (
         "tests/firedrake/regression/test_dg_advection.py::test_dg_advection_icosahedral_sphere[nprocs=3]",
         # vertex-only mesh

setup.py

@@ -238,6 +238,7 @@ def extensions():
     "tests/firedrake/regression/test_nullspace.py",
     "tests/firedrake/regression/test_dg_advection.py",
     "tests/firedrake/regression/test_interpolate_cross_mesh.py",
+    "tests/firedrake/output/test_io_function.py",
 )