DO NOT SQUASH merge release into main

.github/workflows/docker.yml

@@ -75,12 +75,17 @@ jobs:
     if: ${{ ! inputs.build_dev }}
     needs: docker_merge_vanilla
     uses: ./.github/workflows/docker_build.yml
-    # Only build the 'firedrake' container for 'linux/amd64' because
-    # netgen-mesher does not have ARM wheels so many Firedrake apps cannot
-    # be installed.
+    strategy:
+      matrix:
+        os: [Linux, macOS]
+        include:
+        - os: Linux
+          platform: linux/amd64
+        - os: macOS
+          platform: linux/arm64
     with:
-      os: Linux
-      platform: linux/amd64
+      os: ${{ matrix.os }}
+      platform: ${{ matrix.platform }}
       target: firedrake
       tag: ${{ inputs.tag }}
       dockerfile: docker/Dockerfile.firedrake

AUTHORS.rst

@@ -104,7 +104,7 @@ Robert C. Kirby...............<https://www.baylor.edu/math/index.php?id=90540>
 
 Stephan C. Kramer.............<https://www.imperial.ac.uk/people/s.kramer>
 
-Tuomas Kärnä..................<https://www.tuomaskarna.com>
+Tuomas Kärnä
 
 Michael Lange.................<https://www.linkedin.com/in/michael-lange-56675994/>
 

demos/demo_references.bib

@@ -331,9 +331,21 @@ @article{Krishna:1997
 
 @article{BaierReinio:2025,
 	title={High-order finite element methods for three-dimensional multicomponent convection-diffusion},
-  author={Baier-Reinio, Aaron and Farrell, Patrick E},
+  author={Baier-Reinio, Aaron and Farrell, Patrick E.},
   journal={arXiv preprint arXiv:2408.17390},
-  year={2025}
+  year={2025},
+  doi={10.48550/arXiv.2408.17390}
+}
+
+@article{VanBrunt:2025,
+  title={Finite element methods for multicomponent convection-diffusion},
+  author={Aznaran, Francis R. A. and Farrell, Patrick E. and Monroe, Charles W. and Van-Brunt, Alexander J.},
+  journal={IMA J. Numer. Anal.},
+  volume={45},
+  number={1},
+  pages={188--222},
+  year={2025},
+  doi={10.1093/imanum/drae001}
 }
 
 @article{Brubeck2022,

demos/multicomponent/multicomponent.py.rst

@@ -78,6 +78,8 @@ The domain and mesh are visualised below.
 
 To model the mixture we employ the Stokes--Onsager--Stefan--Maxwell (SOSM) 
 partial differential equations and discretise with the method of :cite:`BaierReinio:2025`.
+We use the material property values specified in :cite:`VanBrunt:2025`,
+wherein a similar benzene-cyclohexane simulation is considered.
 In what follows species 1 refers to benzene and species 2 to cyclohexane.
 We shall discretise the following unknowns:
 

docker/Dockerfile.firedrake

@@ -1,27 +1,35 @@
-# Dockerfile for Firedrake with a full set of capabilities and applications installed.
+# Dockerfile for Firedrake with a full set of capabilities installed.
 
 FROM firedrakeproject/firedrake-vanilla-default:latest
 
 # Install optional dependencies
 RUN pip install --verbose \
         --extra-index-url https://download.pytorch.org/whl/cpu \
-        jax ngsPETSc torch vtk
+        jax torch
 
-# Set PYTHONPATH so netgen can be found
-# (see https://github.com/NGSolve/netgen/issues/213)
-ENV PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH
-
-# Install Firedrake apps
-RUN pip install --verbose --src /opt \
-        -e git+https://github.com/firedrakeproject/asQ.git#egg=asQ \
-        -e git+https://bitbucket.org/pefarrell/defcon.git#egg=defcon \
-        -e git+https://bitbucket.org/pefarrell/fascd.git#egg=fascd \
-        -e git+https://github.com/FEMlium/FEMlium.git#egg=FEMlium \
-        -e git+https://github.com/g-adopt/g-adopt.git#egg=gadopt \
-        -e git+https://github.com/firedrakeproject/gusto.git#egg=gusto \
-        -e git+https://github.com/firedrakeproject/Irksome.git#egg=Irksome \
-        -e git+https://github.com/icepack/icepack.git#egg=icepack \
-        -e git+https://github.com/thetisproject/thetis.git#egg=thetis
+# Install ngsPETSc and dependencies. Netgen do not package ARM wheels so
+# we have to build it from source in that case. If on x86 then we have to
+# set PYTHONPATH so netgen can be found (see https://github.com/NGSolve/netgen/issues/213).
+RUN \
+        if [ "$(dpkg --print-architecture)" == "arm64" ]; then \
+                apt-get update \
+                && apt-get -y install python3-tk libxmu-dev tk-dev tcl-dev cmake g++ libglu1-mesa-dev liblapacke-dev libocct-data-exchange-dev libocct-draw-dev occt-misc libtbb-dev libxi-dev \
+                && rm -rf /var/lib/apt/lists/* \
+                && mkdir /opt/ngsuite \
+                && git clone --branch=v6.2.2505 --single-branch https://github.com/NGSolve/netgen.git /opt/ngsuite/netgen-src \
+                && git -C /opt/ngsuite/netgen-src submodule update --init --recursive \
+                && mkdir /opt/ngsuite/netgen-build \
+                && mkdir /opt/ngsuite/netgen-install \
+                && cmake -DCMAKE="-cxx-flags=-flax-vector-conversions" -DCMAKE_INSTALL_PREFIX=/opt/ngsuite/netgen-install /opt/ngsuite/netgen-src \
+                && make \
+                && make install \
+                && export PATH=/opt/ngsuite/netgen-install/bin:$PATH \
+                && export PYTHONPATH=/opt/ngsuite/netgen-install/lib/python3.12/site-packages:$PYTHONPATH \
+                && pip install ngsPETSc --no-deps; \
+        else \
+                export PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH \
+                && pip install ngsPETSc; \
+        fi
 
 # Install some other niceties
 RUN apt-get update \

docs/source/conf.py

@@ -164,15 +164,12 @@
     r'http://www.cs.virginia.edu/stream/',
     r'https://www.sciencedirect.com',
     r'https://.*\.baylor\.edu.*',
-    r'https://www.tuomaskarna.com',
     r'https://www.crosscountrytrains.co.uk/',
     r'https://www.siam.org/',
     r'https://aims.ac.rw',
     r'https://mpecdt.ac.uk',
     r'https://www.hilton.com/en/hotels/leehnhn-hilton-leeds-city/',
-    r'https://www.radissonhotels.com/en-us/hotels/park-plaza-leeds',
-    r'https://www.radissonhotels.com/en-us/hotels/radisson-blu-leeds'
-    r'https://www.radissonhotels.com/en-us/hotels/radisson-blu-leeds',
+    r'https://www.radissonhotels.com/*',
     r'https://all.accor.com/hotel/*',
 ]
 linkcheck_timeout = 30

docs/source/team.ini

@@ -95,7 +95,7 @@ Christian T. Jacobs:
 Darko Janeković:
 Nick Johnson:
 Anna Kalogirou:
-Tuomas Kärnä: https://www.tuomaskarna.com
+Tuomas Kärnä:
 Stephan C. Kramer: https://www.imperial.ac.uk/people/s.kramer
 Nicolas Loriant:
 Fabio Luporini:

docs/source/visualisation.rst

@@ -332,5 +332,5 @@ matplotlib.
 .. _matplotlib: http://matplotlib.org
 .. _Arbitrary: https://www.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
 __ Arbitrary_
-.. _Tessellate: https://www.paraview.org/paraview-docs/nightly/python/paraview.simple.Tessellate.html
+.. _Tessellate: https://www.paraview.org/paraview-docs/nightly/python/paraview.simple.__init__.Tessellate.html
 .. _Tessellation: https://ieeexplore.ieee.org/document/1634311/

firedrake/assemble.py

@@ -386,6 +386,16 @@ def allocation_integral_types(self):
         else:
             return self._allocation_integral_types
 
+    @staticmethod
+    def _as_pyop2_type(tensor, indices=None):
+        if isinstance(tensor, (firedrake.Cofunction, firedrake.Function)):
+            return OneFormAssembler._as_pyop2_type(tensor, indices=indices)
+        elif isinstance(tensor, ufl.Matrix):
+            return ExplicitMatrixAssembler._as_pyop2_type(tensor, indices=indices)
+        else:
+            assert indices is None
+            return tensor
+
     def assemble(self, tensor=None, current_state=None):
         """Assemble the form.
 
@@ -410,21 +420,22 @@ def assemble(self, tensor=None, current_state=None):
         """
         def visitor(e, *operands):
             t = tensor if e is self._form else None
-            return self.base_form_assembly_visitor(e, t, *operands)
+            # Deal with 2-form bcs inside the visitor
+            bcs = self._bcs if isinstance(e, ufl.BaseForm) and len(e.arguments()) == 2 else ()
+            return self.base_form_assembly_visitor(e, t, bcs, *operands)
 
         # DAG assembly: traverse the DAG in a post-order fashion and evaluate the node on the fly.
         visited = {}
         result = BaseFormAssembler.base_form_postorder_traversal(self._form, visitor, visited)
 
-        # Apply BCs after assembly
+        # Deal with 1-form bcs outside the visitor
         rank = len(self._form.arguments())
         if rank == 1 and not isinstance(result, ufl.ZeroBaseForm):
             for bc in self._bcs:
                 OneFormAssembler._apply_bc(self, result, bc, u=current_state)
-
         return result
 
-    def base_form_assembly_visitor(self, expr, tensor, *args):
+    def base_form_assembly_visitor(self, expr, tensor, bcs, *args):
         r"""Assemble a :class:`~ufl.classes.BaseForm` object given its assembled operands.
 
             This functions contains the assembly handlers corresponding to the different nodes that
@@ -445,7 +456,7 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
                 assembler = OneFormAssembler(form, form_compiler_parameters=self._form_compiler_params,
                                              zero_bc_nodes=self._zero_bc_nodes, diagonal=self._diagonal, weight=self._weight)
             elif rank == 2:
-                assembler = TwoFormAssembler(form, bcs=self._bcs, form_compiler_parameters=self._form_compiler_params,
+                assembler = TwoFormAssembler(form, bcs=bcs, form_compiler_parameters=self._form_compiler_params,
                                              mat_type=self._mat_type, sub_mat_type=self._sub_mat_type,
                                              options_prefix=self._options_prefix, appctx=self._appctx, weight=self._weight,
                                              allocation_integral_types=self.allocation_integral_types)
@@ -456,13 +467,12 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
             if len(args) != 1:
                 raise TypeError("Not enough operands for Adjoint")
             mat, = args
-            res = tensor.petscmat if tensor else PETSc.Mat()
-            petsc_mat = mat.petscmat
+            result = tensor.petscmat if tensor else PETSc.Mat()
             # Out-of-place Hermitian transpose
-            petsc_mat.hermitianTranspose(out=res)
-            (row, col) = mat.arguments()
-            return matrix.AssembledMatrix((col, row), self._bcs, res,
-                                          options_prefix=self._options_prefix)
+            mat.petscmat.hermitianTranspose(out=result)
+            if tensor is None:
+                tensor = self.assembled_matrix(expr, bcs, result)
+            return tensor
         elif isinstance(expr, ufl.Action):
             if len(args) != 2:
                 raise TypeError("Not enough operands for Action")
@@ -480,7 +490,7 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
                     result = tensor.petscmat if tensor else PETSc.Mat()
                     lhs.petscmat.matMult(rhs.petscmat, result=result)
                     if tensor is None:
-                        tensor = self.assembled_matrix(expr, result)
+                        tensor = self.assembled_matrix(expr, bcs, result)
                     return tensor
                 else:
                     raise TypeError("Incompatible RHS for Action.")
@@ -499,9 +509,6 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
                 raise TypeError("Mismatching weights and operands in FormSum")
             if len(args) == 0:
                 raise TypeError("Empty FormSum")
-            if tensor:
-                tensor.zero()
-
             # Assemble weights
             weights = []
             for w in expr.weights():
@@ -519,27 +526,54 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
                     raise ValueError("Expecting a scalar weight expression")
                 weights.append(w)
 
+            # Scalar FormSum
             if all(isinstance(op, numbers.Complex) for op in args):
-                result = sum(weight * arg for weight, arg in zip(weights, args))
-                return tensor.assign(result) if tensor else result
-            elif (all(isinstance(op, firedrake.Cofunction) for op in args)
+                result = numpy.dot(weights, args)
+                return tensor.assign(result) if tensor else result.item()
+
+            # Accumulate coefficients in a dictionary for each unique Dat/Mat
+            terms = defaultdict(PETSc.ScalarType)
+            for arg, weight in zip(args, weights):
+                t = self._as_pyop2_type(arg)
+                terms[t] += weight
+
+            # Zero the output tensor, or rescale it if it appears in the sum
+            tensor_scale = terms.pop(self._as_pyop2_type(tensor), 0)
+            if tensor is None or tensor_scale == 1:
+                pass
+            elif tensor_scale == 0:
+                tensor.zero()
+            elif isinstance(tensor, (firedrake.Cofunction, firedrake.Function)):
+                with tensor.dat.vec as v:
+                    v.scale(tensor_scale)
+            elif isinstance(tensor, ufl.Matrix):
+                tensor.petscmat.scale(tensor_scale)
+            else:
+                raise ValueError("Expecting tensor to be None, Function, Cofunction, or Matrix")
+
+            # Compute the linear combination
+            if (all(isinstance(op, firedrake.Cofunction) for op in args)
                     or all(isinstance(op, firedrake.Function) for op in args)):
+                # Vector FormSum
                 V, = set(a.function_space() for a in args)
                 result = tensor if tensor else firedrake.Function(V)
-                result.dat.maxpy(weights, [a.dat for a in args])
+                weights = terms.values()
+                dats = terms.keys()
+                result.dat.maxpy(weights, dats)
                 return result
             elif all(isinstance(op, ufl.Matrix) for op in args):
+                # Matrix FormSum
                 result = tensor.petscmat if tensor else PETSc.Mat()
-                for (op, w) in zip(args, weights):
+                for (op, w) in terms.items():
                     if result:
                         # If result is not void, then accumulate on it
-                        result.axpy(w, op.petscmat)
+                        result.axpy(w, op.handle)
                     else:
                         # If result is void, then allocate it with first term
-                        op.petscmat.copy(result=result)
+                        op.handle.copy(result=result)
                         result.scale(w)
                 if tensor is None:
-                    tensor = self.assembled_matrix(expr, result)
+                    tensor = self.assembled_matrix(expr, bcs, result)
                 return tensor
             else:
                 raise TypeError("Mismatching FormSum shapes")
@@ -571,9 +605,8 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
                 # Occur in situations such as Interpolate composition
                 operand = assembled_operand[0]
 
-            reconstruct_interp = expr._ufl_expr_reconstruct_
             if (v, operand) != expr.argument_slots():
-                expr = reconstruct_interp(operand, v=v)
+                expr = expr._ufl_expr_reconstruct_(operand, v=v)
 
             rank = len(expr.arguments())
             if rank > 2:
@@ -586,7 +619,7 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
             default_missing_val = interp_data.pop('default_missing_val', None)
             if rank == 1 and isinstance(tensor, firedrake.Function):
                 V = tensor
-            interpolator = firedrake.Interpolator(expr, V, **interp_data)
+            interpolator = firedrake.Interpolator(expr, V, bcs=bcs, **interp_data)
             # Assembly
             return interpolator.assemble(tensor=tensor, default_missing_val=default_missing_val)
         elif tensor and isinstance(expr, (firedrake.Function, firedrake.Cofunction, firedrake.MatrixBase)):
@@ -598,8 +631,8 @@ def base_form_assembly_visitor(self, expr, tensor, *args):
         else:
             raise TypeError(f"Unrecognised BaseForm instance: {expr}")
 
-    def assembled_matrix(self, expr, petscmat):
-        return matrix.AssembledMatrix(expr.arguments(), self._bcs, petscmat,
+    def assembled_matrix(self, expr, bcs, petscmat):
+        return matrix.AssembledMatrix(expr.arguments(), bcs, petscmat,
                                       options_prefix=self._options_prefix)
 
     @staticmethod
@@ -1448,10 +1481,11 @@ def _apply_bc(self, tensor, bc, u=None):
         index = 0 if V.index is None else V.index
         space = V if V.parent is None else V.parent
         if isinstance(bc, DirichletBC):
-            if space != spaces[0]:
-                raise TypeError("bc space does not match the test function space")
-            elif space != spaces[1]:
-                raise TypeError("bc space does not match the trial function space")
+            if not any(space == fs for fs in spaces):
+                raise TypeError("bc space does not match the test or trial function space")
+            if spaces[0] != spaces[1]:
+                # Not on a diagonal block, we cannot set diagonal entries
+                return
 
             # Set diagonal entries on bc nodes to 1 if the current
             # block is on the matrix diagonal and its index matches the

firedrake/interpolation.py

@@ -1703,7 +1703,12 @@ def __init__(self, expr, V, bcs=None, **kwargs):
         super(MixedInterpolator, self).__init__(expr, V, bcs=bcs, **kwargs)
         expr = self.ufl_interpolate
         self.arguments = expr.arguments()
-        rank = len(self.arguments)
+        # Get the primal spaces
+        spaces = tuple(a.function_space().dual() if isinstance(a, Coargument) else a.function_space()
+                       for a in self.arguments)
+        # TODO consider a stricter equality test for indexed MixedFunctionSpace
+        # See https://github.com/firedrakeproject/firedrake/issues/4668
+        space_equals = lambda V1, V2: V1 == V2 and V1.parent == V2.parent and V1.index == V2.index
 
         # We need a Coargument in order to split the Interpolate
         needs_action = len([a for a in self.arguments if isinstance(a, Coargument)]) == 0
@@ -1722,12 +1727,10 @@ def __init__(self, expr, V, bcs=None, **kwargs):
                 continue
             vi, _ = form.argument_slots()
             Vtarget = vi.function_space().dual()
-            if bcs and rank != 0:
-                args = form.arguments()
-                Vsource = args[1-vi.number()].function_space()
-                sub_bcs = [bc for bc in bcs if bc.function_space() in {Vsource, Vtarget}]
-            else:
-                sub_bcs = None
+            sub_bcs = []
+            for space, index in zip(spaces, indices):
+                subspace = space.sub(index)
+                sub_bcs.extend(bc for bc in bcs if space_equals(bc.function_space(), subspace))
             if needs_action:
                 # Take the action of each sub-cofunction against each block
                 form = action(form, dual_split[indices[-1:]])