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 Advection_flux
 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,102 @@
+import ufl
+from dolfinx import fem
+
+from festim import subdomain as _subdomain
+from festim.advection import VelocityField
+from festim.species import Species
+
+
+class OutflowBC:
+    """
+    Dirichlet boundary condition class
+    u = value
+
+    Args:
+        subdomain: The surface subdomain where the boundary condition is applied
+        value: The value of the boundary condition
+
+    Attributes:
+        subdomain: The surface subdomain where the boundary condition is applied
+        value: The value of the boundary condition
+        value_fenics: The value of the boundary condition in fenics format
+        bc_expr: The expression of the boundary condition that is used to
+            update the `value_fenics`
+
+    """
+
+    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 Advection_flux
 from .average_surface import AverageSurface
 from .average_volume import AverageVolume
 from .maximum_surface import MaximumSurface
@@ -14,6 +15,7 @@
 from .xdmf import XDMFExport
 
 __all__ = [
+    "Advection_flux",
     "AverageSurface",
     "AverageVolume",
     "ExportBaseClass",

src/festim/exports/advection_flux.py

@@ -0,0 +1,67 @@
+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 Advection_flux(SurfaceFlux):
+    """Computes the advective flux of a field on a given surface
+
+    Args:
+        field: species for which the advective flux is computed
+        surface: surface subdomain
+        filename: name of the file to which the advective flux is
+            exported
+        velocity_field: velocity field for which the advective flux is computed
+
+    Attributes:
+        see `festim.SurfaceFlux`
+    """
+
+    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

@@ -754,6 +754,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:
@@ -882,15 +888,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
@@ -1011,11 +1030,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()

test/system_tests/test_advection_term_integration.py

@@ -279,3 +279,82 @@ def velocity_func(x):
 
     my_model.initialise()
     my_model.run()
+
+
+def test_flux_conservation_advection():
+    """MMS coupled heat and hydrogen test with 1 mobile species and 1 trap in a 1s
+    transient, the values of the temperature, mobile and trapped solutions at the last
+    time step is compared to an analytical solution"""
+
+    my_model = F.HydrogenTransportProblem()
+
+    test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 200, 200)
+    my_model.mesh = F.Mesh(test_mesh_2d)
+
+    # common festim objects
+    test_mat = F.Material(D_0=1.2, E_D=0.1)
+    test_vol_sub = F.VolumeSubdomain(id=1, material=test_mat)
+    left = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0))
+    right = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))
+    bottom = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[1], 0))
+    top = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[1], 1))
+    my_model.subdomains = [test_vol_sub, left, right, bottom, top]
+
+    test_H = F.Species("H", mobile=True)
+    my_model.species = [test_H]
+
+    V = dolfinx.fem.functionspace(test_mesh_2d, ("Lagrange", 1, (2,)))
+    velocity_field = dolfinx.fem.Function(V)
+    factor = 100
+    velocity_field.interpolate(lambda x: (factor * np.sin(x[0]), factor * np.cos(x[1])))
+    my_model.advection_terms = [
+        F.AdvectionTerm(velocity=velocity_field, subdomain=test_vol_sub, species=test_H)
+    ]
+
+    my_model.temperature = 100
+
+    my_model.boundary_conditions = [
+        F.OutflowBC(velocity=velocity_field, species=test_H, subdomain=right),
+        F.OutflowBC(velocity=velocity_field, species=test_H, subdomain=top),
+        F.OutflowBC(velocity=velocity_field, species=test_H, subdomain=left),
+        F.FixedConcentrationBC(species=test_H, subdomain=bottom, value=0),
+    ]
+
+    my_model.sources = [F.ParticleSource(species=test_H, value=1, volume=test_vol_sub)]
+
+    my_model.settings = F.Settings(transient=False, atol=1e-10, rtol=1e-10)
+
+    advection_flux_top = F.Advection_flux(
+        velocity_field=velocity_field, field=test_H, surface=top
+    )
+    advection_flux_right = F.Advection_flux(
+        velocity_field=velocity_field, field=test_H, surface=right
+    )
+    advection_flux_left = F.Advection_flux(
+        velocity_field=velocity_field, field=test_H, surface=left
+    )
+    advection_flux_bottom = F.Advection_flux(
+        velocity_field=velocity_field, field=test_H, surface=bottom
+    )
+    my_model.exports = [
+        advection_flux_top,
+        advection_flux_right,
+        advection_flux_left,
+        advection_flux_bottom,
+    ]
+
+    my_model.initialise()
+    my_model.run()
+
+    total_source = dolfinx.fem.assemble_scalar(
+        dolfinx.fem.form(my_model.sources[0].value.fenics_object * ufl.dx)
+    )
+
+    total_outflux = (
+        advection_flux_top.data[-1]
+        + advection_flux_right.data[-1]
+        + advection_flux_left.data[-1]
+        + advection_flux_bottom.data[-1]
+    )
+
+    assert np.isclose(total_source, total_outflux)