Knook/multicomponent (#4322)

Multicomponent microfluidics demo

---------

Co-authored-by: Patrick E. Farrell <patrick.farrell@maths.ox.ac.uk>
Co-authored-by: Connor Ward <c.ward20@imperial.ac.uk>
Co-authored-by: Aaron Baier-Reinio <130484015+ABaierReinio@users.noreply.github.com>
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
Co-authored-by: ABaierReinio <baierreinio@maths.ox.ac.uk>

Author: 5 people

demos/boussinesq/boussinesq.py.rst

@@ -160,43 +160,46 @@ cannot always be sure that the linear solver at hand is correctly utilising the
 To directly eliminate the nullspace we introduce a class :code:`FixAtPointBC` which
 implements a boundary condition that fixes a field at a single point. ::
 
-    import firedrake.utils as firedrake_utils
-
-    class FixAtPointBC(firedrake.DirichletBC):
-        r'''A special BC object for pinning a function at a point.
-
-        :arg V: the :class:`.FunctionSpace` on which the boundary condition should be applied.
-        :arg g: the boundary condition value.
-        :arg bc_point: the point at which to pin the function.
-            The location of the finite element DOF nearest to bc_point is actually used.
-        '''
-        def __init__(self, V, g, bc_point):
-            super(FixAtPointBC, self).__init__(V, g, bc_point)
-            if isinstance(bc_point, tuple):
-                bc_point = as_vector(bc_point)
-            self.bc_point = bc_point
-
-        @firedrake_utils.cached_property
-        def nodes(self):
-            V = self.function_space()
-            x = firedrake.SpatialCoordinate(V.mesh())
-            xdist = x - self.bc_point
-
-            test = firedrake.TestFunction(V)
-            trial = firedrake.TrialFunction(V)
-            xphi = firedrake.assemble(ufl.inner(xdist * test, xdist * trial) * ufl.dx, diagonal=True)
-            phi = firedrake.assemble(ufl.inner(test, trial) * ufl.dx, diagonal=True)
-            with xphi.dat.vec as xu, phi.dat.vec as u:
-                xu.pointwiseDivide(xu, u)
-                min_index, min_value = xu.min()     # Find the index of the DOF closest to bc_point
-
-            nodes = V.dof_dset.lgmap.applyInverse([min_index])
-            nodes = nodes[nodes >= 0]
-            return nodes
-
+    import functools
+
+    class FixAtPointBC(DirichletBC):
+       r'''A special BC object for pinning a function at a point.
+
+       :arg V: the :class:`.FunctionSpace` on which the boundary condition should be applied.
+       :arg g: the boundary condition value.
+       :arg bc_point: the point at which to pin the function.
+           The location of the finite element DOF nearest to bc_point is actually used.
+       '''
+       def __init__(self, V, g, bc_point):
+           super().__init__(V, g, bc_point)
+
+       @functools.cached_property
+       def nodes(self):
+           V = self.function_space()
+           if V.mesh().ufl_coordinate_element().degree() != 1:
+               # Ensure a P1 mesh
+               coordinates = V.mesh().coordinates
+               P1 = coordinates.function_space().reconstruct(degree=1)
+               P1_mesh = Mesh(Function(P1).interpolate(coordinates))
+               V = V.reconstruct(mesh=P1_mesh)
+
+           point = [tuple(self.sub_domain)]
+           vom = VertexOnlyMesh(V.mesh(), point)
+           P0 = FunctionSpace(vom, "DG", 0)
+           Fvom = Cofunction(P0.dual()).assign(1)
+
+           # Take the basis function with the largest abs value at bc_point
+           v = TestFunction(V)
+           F = assemble(Interpolate(inner(v, v), Fvom))
+           with F.dat.vec as Fvec:
+               max_index, _ = Fvec.max()
+           nodes = V.dof_dset.lgmap.applyInverse([max_index])
+           nodes = nodes[nodes >= 0]
+           return nodes
+ 
 We use this to fix the pressure and auxiliary temperature at the origin::
 
-    aux_bcs = [FixAtPointBC(Z.sub(1), 0, as_vector([0, 0])), 
+    aux_bcs = [FixAtPointBC(Z.sub(1), 0, (0, 0)), 
                FixAtPointBC(Z.sub(2), 0, as_vector([0, 0]))]
 
 :code:`FixAtPointBC` takes three arguments: the function space to fix, the value with which it

demos/demo_references.bib

@@ -34,6 +34,14 @@ @Book{2006:SV
   year =         2006
 }
 
+@book{Perry:2007,
+	title =        {Perry’s Chemical Engineers’ Handbook},
+	edition =      {8th},
+	author =       {Green, Don W. and Perry, Robert H.},
+	year =         2007,
+	publisher =   {McGraw Hill Professional}
+}
+
 @InBook{2015:lmscup,
   author =       {O. Bokhove and A. Kalogirou},
   title =        {Lectures on the Theory of Water Waves},
@@ -309,6 +317,25 @@ @article{Ricker:1953
   doi={https://doi.org/10.1190/1.1437843}
 }
 
+@article{Krishna:1997,
+	title={The {M}axwell--{S}tefan approach to mass transfer},
+	author={Krishna, Rajamani and Wesselingh, Johannes A},
+	journal={Chemical Engineering Science},
+	volume={52},
+  number={6},
+	pages={861--911},
+	year={1997},
+  publisher={Elsevier},
+	doi={10.1016/S0009-2509(96)00458-7}
+}
+
+@article{BaierReinio:2025,
+	title={High-order finite element methods for three-dimensional multicomponent convection-diffusion},
+  author={Baier-Reinio, Aaron and Farrell, Patrick E},
+  journal={arXiv preprint arXiv:2408.17390},
+  year={2025}
+}
+
 @article{Brubeck2022,
   author        = {Brubeck, Pablo D. and Farrell, Patrick E.},
   doi           = {10.1137/21M1444187},