new OutflowBC

src/festim/__init__.py

@@ -18,11 +18,13 @@
 )
 from .boundary_conditions.flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC
 from .boundary_conditions.henrys_bc import HenrysBC
+from .boundary_conditions.outflow_bc import OutflowBC
 from .boundary_conditions.sieverts_bc import SievertsBC
 from .boundary_conditions.surface_reaction import SurfaceReactionBC
 from .coupled_heat_hydrogen_problem import (
     CoupledTransientHeatTransferHydrogenTransport,
 )
+from .exports.advection_flux import AdvectionFlux
 from .exports.average_surface import AverageSurface
 from .exports.average_volume import AverageVolume
 from .exports.maximum_surface import MaximumSurface

src/festim/advection.py

@@ -151,18 +151,22 @@ def convert_input_value(
                     f", not a {type(vel)}"
                 )
 
-        # create vector function space and function
-        v_cg = basix.ufl.element(
-            "Lagrange",
-            function_space.mesh.topology.cell_name(),
-            1,
-            shape=(function_space.mesh.topology.dim,),
-        )
-        self.vector_function_space = fem.functionspace(function_space.mesh, v_cg)
-        self.fenics_object = fem.Function(self.vector_function_space)
-
-        # interpolate input value into fenics object function
-        nmm_interpolate(self.fenics_object, vel)
+        # check if mesh in value is same as mesh in function space, if not raise ValueError
+        if vel.function_space.mesh == function_space.mesh:
+            self.fenics_object = vel
+        else:
+            # create vector function space and function
+            v_cg = basix.ufl.element(
+                "Lagrange",
+                function_space.mesh.topology.cell_name(),
+                1,
+                shape=(function_space.mesh.topology.dim,),
+            )
+            self.vector_function_space = fem.functionspace(function_space.mesh, v_cg)
+            self.fenics_object = fem.Function(self.vector_function_space)
+
+            # interpolate input value into fenics object function
+            nmm_interpolate(self.fenics_object, vel)
 
     def update(self, t: fem.Constant):
         """Updates the velocity field

src/festim/boundary_conditions/__init__.py

@@ -1,6 +1,7 @@
 from .dirichlet_bc import DirichletBCBase, FixedConcentrationBC, FixedTemperatureBC
 from .flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC
 from .henrys_bc import HenrysBC
+from .outflow_bc import OutflowBC
 from .sieverts_bc import SievertsBC
 from .surface_reaction import SurfaceReactionBC, SurfaceReactionBCpartial
 
@@ -11,6 +12,7 @@
     "FluxBCBase",
     "HeatFluxBC",
     "HenrysBC",
+    "OutflowBC",
     "ParticleFluxBC",
     "SievertsBC",
     "SurfaceReactionBC",

src/festim/boundary_conditions/outflow_bc.py

@@ -0,0 +1,110 @@
+import ufl
+from dolfinx import fem
+
+from festim import subdomain as _subdomain
+from festim.advection import VelocityField
+from festim.species import Species
+
+
+class OutflowBC:
+    """
+    Outflow boundary condition class
+
+    dot(velocity, n) * u * v * ds
+
+    This term allows the species to flow across the boundary due to advection. It arises
+    from integrating the conservative form of the advection term by parts in the weak
+    formulation.
+
+    Args:
+        velocity: The velocity field at the outflow boundary condition
+        species: The species for which the outflow boundary condition is applied, can be
+            a list of species or a single species
+        subdomain: The surface subdomain where the boundary condition is applied
+
+    Attributes:
+        velocity: The velocity field at the outflow boundary condition
+        species: The species for which the outflow boundary condition is applied, can be
+            a list of species or a single species
+        subdomain: The surface subdomain where the boundary condition is applied
+
+    """
+
+    velocity: VelocityField
+    species: list[Species]
+    subdomain: _subdomain.SurfaceSubdomain
+
+    def __init__(
+        self,
+        velocity: fem.Function,
+        species: Species | list[Species],
+        subdomain: _subdomain.SurfaceSubdomain,
+    ):
+        self.subdomain = subdomain
+        self.velocity = velocity
+        self.species = species
+
+    @property
+    def velocity(self):
+        return self._velocity
+
+    @velocity.setter
+    def velocity(self, value):
+        err_message = f"velocity must be a fem.Function, or callable not {type(value)}"
+        if value is None:
+            self._velocity = VelocityField(value)
+        elif isinstance(
+            value,
+            fem.Function,
+        ):
+            self._velocity = VelocityField(value)
+        elif isinstance(value, fem.Constant | fem.Expression | ufl.core.expr.Expr):
+            raise TypeError(err_message)
+        elif callable(value):
+            self._velocity = VelocityField(value)
+        else:
+            raise TypeError(err_message)
+
+    @property
+    def species(self) -> list[Species]:
+        return self._species
+
+    @species.setter
+    def species(self, value):
+        if not isinstance(value, list):
+            value = [value]
+        # check that all species are of type festim.Species
+        for spe in value:
+            if not isinstance(spe, Species):
+                raise TypeError(
+                    f"elements of species must be of type festim.Species not "
+                    f"{type(spe)}"
+                )
+        self._species = value
+
+    @property
+    def subdomain(self):
+        return self._subdomain
+
+    @subdomain.setter
+    def subdomain(self, value):
+        if value is None:
+            self._subdomain = value
+        elif isinstance(value, _subdomain.SurfaceSubdomain):
+            self._subdomain = value
+        else:
+            raise TypeError(
+                f"Subdomain must be a festim.SurfaceSubdomain object, not {type(value)}"
+            )
+
+    @property
+    def time_dependent(self):
+        if self.velocity is None:
+            raise TypeError("Value must be given to determine if its time dependent")
+        if isinstance(self.velocity, fem.Constant):
+            return False
+        if callable(self.velocity):
+            arguments = self.velocity.__code__.co_varnames
+            return "t" in arguments
+        else:
+            return False

src/festim/exports/__init__.py

@@ -1,3 +1,4 @@
+from .advection_flux import AdvectionFlux
 from .average_surface import AverageSurface
 from .average_volume import AverageVolume
 from .maximum_surface import MaximumSurface
@@ -14,6 +15,7 @@
 from .xdmf import XDMFExport
 
 __all__ = [
+    "AdvectionFlux",
     "AverageSurface",
     "AverageVolume",
     "ExportBaseClass",

src/festim/exports/advection_flux.py

@@ -0,0 +1,76 @@
+import ufl
+from dolfinx import fem
+from scifem import assemble_scalar
+
+from festim.exports.surface_flux import SurfaceFlux
+from festim.species import Species
+from festim.subdomain.surface_subdomain import SurfaceSubdomain
+
+
+class AdvectionFlux(SurfaceFlux):
+    """Computes the advective flux of a field on a given surface
+
+    Advective flux is computed as the sum of the diffusive flux and the advective flux
+    at the surface:
+
+    J = ∫-D (∇u · n) ds + ∫(u · n) ds
+
+    Args:
+        field: species for which the advective flux is computed
+        surface: surface subdomain
+        velocity_field: velocity field for which the advective flux is computed
+        filename: name of the file to which the advective flux is
+            exported
+
+    Attributes:
+        field: species for which the advective flux is computed
+        surface: surface subdomain
+        velocity_field: velocity field for which the advective flux is computed
+        filename: name of the file to which the advective flux is
+            exported
+    """
+
+    field: Species
+    surface: SurfaceSubdomain
+    velocity_field: fem.Function
+
+    def __init__(
+        self,
+        field: Species,
+        surface: SurfaceSubdomain,
+        velocity_field: fem.Function,
+        filename: str | None = None,
+    ):
+        super().__init__(field=field, surface=surface, filename=filename)
+
+        self.velocity_field = velocity_field
+
+    @property
+    def title(self):
+        return f"{self.field.name} advective flux surface {self.surface.id}"
+
+    def compute(self, u, ds: ufl.Measure, entity_maps=None):
+        """Computes the value of the flux at the surface
+
+        Args:
+            ds (ufl.Measure): surface measure of the model
+        """
+
+        mesh = ds.ufl_domain()
+        n = ufl.FacetNormal(mesh)
+
+        surface_flux = assemble_scalar(
+            fem.form(
+                -self.D * ufl.dot(ufl.grad(u), n) * ds(self.surface.id),
+                entity_maps=entity_maps,
+            )
+        )
+        advective_flux = assemble_scalar(
+            fem.form(
+                u * ufl.inner(self.velocity_field, n) * ds(self.surface.id),
+                entity_maps=entity_maps,
+            )
+        )
+
+        self.value = surface_flux + advective_flux
+        self.data.append(self.value)

src/festim/hydrogen_transport_problem.py

@@ -751,6 +751,12 @@ def convert_advection_term_to_fenics_objects(self):
                 function_space=self.function_space, t=self.t
             )
 
+        for bc in self.boundary_conditions:
+            if isinstance(bc, boundary_conditions.OutflowBC):
+                bc.velocity.convert_input_value(
+                    function_space=self.function_space, t=self.t
+                )
+
     def create_flux_values_fenics(self):
         """For each particle flux create the value_fenics"""
         for bc in self.boundary_conditions:
@@ -879,15 +885,28 @@ def create_formulation(self):
                         * self.ds(flux_bc.subdomain.id)
                     )
 
-        for adv_term in self.advection_terms:
-            # create vector functionspace based on the elements in the mesh
+            # add outflow term for outflow bcs in advection diffusion cases
+            if isinstance(bc, boundary_conditions.OutflowBC):
+                for species in bc.species:
+                    conc = species.solution
+                    v = species.test_function
+                    vel = bc.velocity.fenics_object
+
+                    outflow_term = (
+                        ufl.inner(vel * conc, self.mesh.n)
+                        * v
+                        * self.ds(bc.subdomain.id)
+                    )
+                    self.formulation += outflow_term
 
+        # add advection terms
+        for adv_term in self.advection_terms:
             for species in adv_term.species:
                 conc = species.solution
                 v = species.test_function
                 vel = adv_term.velocity.fenics_object
 
-                advection_term = ufl.inner(ufl.dot(ufl.grad(conc), vel), v) * self.dx(
+                advection_term = -ufl.inner(vel * conc, ufl.grad(v)) * self.dx(
                     adv_term.subdomain.id
                 )
                 self.formulation += advection_term
@@ -1008,11 +1027,6 @@ def post_processing(self):
                     export,
                     exports.SurfaceFlux | exports.TotalSurface | exports.AverageSurface,
                 ):
-                    if len(self.advection_terms) > 0:
-                        warnings.warn(
-                            "Advection terms are not currently accounted for in the "
-                            "evaluation of surface flux values"
-                        )
                     export.compute(export.field.solution, self.ds)
                 else:
                     export.compute()