From df074f4166b069933fc3d1c2faf60d7dbb810b42 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Mon, 2 Mar 2026 15:32:31 -0500
Subject: [PATCH 01/19] new outlfow bc

---
 src/festim/boundary_conditions/outflow_bc.py | 91 ++++++++++++++++++++
 1 file changed, 91 insertions(+)
 create mode 100644 src/festim/boundary_conditions/outflow_bc.py

diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py
new file mode 100644
index 000000000..ae6efd70f
--- /dev/null
+++ b/src/festim/boundary_conditions/outflow_bc.py
@@ -0,0 +1,91 @@
+import ufl
+from dolfinx import fem
+
+from festim import subdomain as _subdomain
+from festim.advection import VelocityField
+from festim.species import Species
+from festim.subdomain import VolumeSubdomain
+
+
+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_field = 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, VolumeSubdomain):
+            self._subdomain = value
+        else:
+            raise TypeError(
+                f"Subdomain must be a festim.Subdomain object, not {type(value)}"
+            )

From 4924f586a595b0988c594d9581c5488bc681c70a Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Mon, 2 Mar 2026 15:32:40 -0500
Subject: [PATCH 02/19] accesible

---
 src/festim/__init__.py                     | 1 +
 src/festim/boundary_conditions/__init__.py | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index d6df5c269..8b2b74385 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -18,6 +18,7 @@
 )
 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 (
diff --git a/src/festim/boundary_conditions/__init__.py b/src/festim/boundary_conditions/__init__.py
index 1df3d9c7a..43e465669 100644
--- a/src/festim/boundary_conditions/__init__.py
+++ b/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",

From b8027f7a2ffdcdfe4f487d26646bac33e627a9e7 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Mon, 2 Mar 2026 15:33:15 -0500
Subject: [PATCH 03/19] new advection formulation + outflow bc

---
 src/festim/hydrogen_transport_problem.py | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 805511d4a..24a5d99c3 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -869,11 +869,25 @@ def create_formulation(self):
                 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
 
+            for bc in self.boundary_conditions:
+                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
+
         # check if each species is defined in all volumes
         if not self.settings.transient:
             for spe in self.species:

From b24ac852a702597cbce35425587a00492c646b3d Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Mon, 2 Mar 2026 15:34:09 -0500
Subject: [PATCH 04/19] pull out of loop, added docs

---
 src/festim/hydrogen_transport_problem.py | 30 ++++++++++++------------
 1 file changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 24a5d99c3..289f10732 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -861,9 +861,8 @@ def create_formulation(self):
                         * self.ds(flux_bc.subdomain.id)
                     )
 
+        # add advection terms
         for adv_term in self.advection_terms:
-            # create vector functionspace based on the elements in the mesh
-
             for species in adv_term.species:
                 conc = species.solution
                 v = species.test_function
@@ -874,19 +873,20 @@ def create_formulation(self):
                 )
                 self.formulation += advection_term
 
-            for bc in self.boundary_conditions:
-                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 outflow term for outflow bcs in advection diffusion cases
+        for bc in self.boundary_conditions:
+            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
 
         # check if each species is defined in all volumes
         if not self.settings.transient:

From be69dbf49e1e021189233fa32131c7966467a66d Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 08:43:30 -0500
Subject: [PATCH 05/19] move to be with other flux terms

---
 src/festim/hydrogen_transport_problem.py | 27 ++++++++++++------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 289f10732..2979ba654 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -861,20 +861,7 @@ def create_formulation(self):
                         * self.ds(flux_bc.subdomain.id)
                     )
 
-        # 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(vel * conc, ufl.grad(v)) * self.dx(
-                    adv_term.subdomain.id
-                )
-                self.formulation += advection_term
-
-        # add outflow term for outflow bcs in advection diffusion cases
-        for bc in self.boundary_conditions:
+            # add outflow term for outflow bcs in advection diffusion cases
             if isinstance(bc, boundary_conditions.OutflowBC):
                 for species in bc.species:
                     conc = species.solution
@@ -888,6 +875,18 @@ def create_formulation(self):
                     )
                     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(vel * conc, ufl.grad(v)) * self.dx(
+                    adv_term.subdomain.id
+                )
+                self.formulation += advection_term
+
         # check if each species is defined in all volumes
         if not self.settings.transient:
             for spe in self.species:

From bf3d037a625057a8426b6a59fb465e1b252a9415 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 09:22:17 -0500
Subject: [PATCH 06/19] convert value to fenics object

---
 src/festim/hydrogen_transport_problem.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 2979ba654..aff3fbb93 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -733,6 +733,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:

From 7dad4eded56db0035730c4b9c1a3360ac2138045 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 09:22:48 -0500
Subject: [PATCH 07/19] fix subdomain typos, fix attribute name, time
 dependence

---
 src/festim/boundary_conditions/outflow_bc.py | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py
index ae6efd70f..f2eef129d 100644
--- a/src/festim/boundary_conditions/outflow_bc.py
+++ b/src/festim/boundary_conditions/outflow_bc.py
@@ -4,7 +4,6 @@
 from festim import subdomain as _subdomain
 from festim.advection import VelocityField
 from festim.species import Species
-from festim.subdomain import VolumeSubdomain
 
 
 class OutflowBC:
@@ -34,7 +33,7 @@ def __init__(
         subdomain: _subdomain.SurfaceSubdomain,
     ):
         self.subdomain = subdomain
-        self.velocity_field = velocity
+        self.velocity = velocity
         self.species = species
 
     @property
@@ -83,9 +82,21 @@ def subdomain(self):
     def subdomain(self, value):
         if value is None:
             self._subdomain = value
-        elif isinstance(value, VolumeSubdomain):
+        elif isinstance(value, _subdomain.SurfaceSubdomain):
             self._subdomain = value
         else:
             raise TypeError(
-                f"Subdomain must be a festim.Subdomain object, not {type(value)}"
+                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

From a11f25e473d528bebcf9455e38bf69aa9c5b2889 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 13:04:21 -0500
Subject: [PATCH 08/19] only interpolate if the meshes are different

---
 src/festim/advection.py | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/festim/advection.py b/src/festim/advection.py
index 2ecf944f3..6535eaccf 100644
--- a/src/festim/advection.py
+++ b/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

From 703cfa1c1c414c209158068ec90d17d6c08b4511 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 16:33:44 -0500
Subject: [PATCH 09/19] advective flux export

---
 src/festim/__init__.py               |  1 +
 src/festim/exports/__init__.py       |  2 +
 src/festim/exports/advection_flux.py | 61 ++++++++++++++++++++++++++++
 3 files changed, 64 insertions(+)
 create mode 100644 src/festim/exports/advection_flux.py

diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index 8b2b74385..5ac8afd88 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -24,6 +24,7 @@
 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
diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py
index 31a287244..cbfb5a1af 100644
--- a/src/festim/exports/__init__.py
+++ b/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",
diff --git a/src/festim/exports/advection_flux.py b/src/festim/exports/advection_flux.py
new file mode 100644
index 000000000..36da8d962
--- /dev/null
+++ b/src/festim/exports/advection_flux.py
@@ -0,0 +1,61 @@
+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, surface, filename, velocity_field):
+        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)

From 55c9082633142ea346731e1e4510143cc6fbb2ff Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 16:50:49 -0500
Subject: [PATCH 10/19] type hinting, function

---
 src/festim/exports/advection_flux.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/festim/exports/advection_flux.py b/src/festim/exports/advection_flux.py
index 36da8d962..aa2bb5742 100644
--- a/src/festim/exports/advection_flux.py
+++ b/src/festim/exports/advection_flux.py
@@ -25,7 +25,13 @@ class Advection_flux(SurfaceFlux):
     surface: SurfaceSubdomain
     velocity_field: fem.Function
 
-    def __init__(self, field, surface, filename, velocity_field):
+    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
@@ -41,7 +47,7 @@ def compute(self, u, ds: ufl.Measure, entity_maps=None):
             ds (ufl.Measure): surface measure of the model
         """
 
-        mesh = ds.ufl_domain
+        mesh = ds.ufl_domain()
         n = ufl.FacetNormal(mesh)
 
         surface_flux = assemble_scalar(

From 7f7580fe876872be4f8b1c50ec2a3fcb48d39e48 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 16:50:56 -0500
Subject: [PATCH 11/19] can remove warning now

---
 src/festim/hydrogen_transport_problem.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index aff3fbb93..fb23c2c6c 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1009,11 +1009,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()

From c2efd1856633185c51def9a94afd1c702d1ef2b3 Mon Sep 17 00:00:00 2001
From: jhdark <james.dark.95@gmail.com>
Date: Tue, 3 Mar 2026 16:51:03 -0500
Subject: [PATCH 12/19] system test

---
 .../test_advection_term_integration.py        | 79 +++++++++++++++++++
 1 file changed, 79 insertions(+)

diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index 3d3e728a0..c902bf319 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/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)

From 2686fbc0f82d8d9c5440c52d2301aba90ffb5767 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Fri, 27 Mar 2026 16:09:29 -0400
Subject: [PATCH 13/19] fix doc strings

---
 src/festim/boundary_conditions/outflow_bc.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py
index f2eef129d..8695c22a9 100644
--- a/src/festim/boundary_conditions/outflow_bc.py
+++ b/src/festim/boundary_conditions/outflow_bc.py
@@ -8,22 +8,23 @@
 
 class OutflowBC:
     """
-    Dirichlet boundary condition class
-    u = value
+    Outflow boundary condition class
+    
 
     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
-        value: The value of the boundary condition
 
     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
-        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`
 
     """
 
+    velocity: VelocityField
+    species: list[Species]
     subdomain: _subdomain.SurfaceSubdomain
 
     def __init__(

From 708a3a782ca2a49f386bf4d2c3d18ca47fc61172 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Fri, 27 Mar 2026 16:16:48 -0400
Subject: [PATCH 14/19] upadte doc strings

---
 src/festim/boundary_conditions/outflow_bc.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py
index 8695c22a9..4c24f0e7d 100644
--- a/src/festim/boundary_conditions/outflow_bc.py
+++ b/src/festim/boundary_conditions/outflow_bc.py
@@ -9,16 +9,23 @@
 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
+        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
+        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
 
     """

From 3d26353fa87e4aa0d9146b0c8627c9b617f106f4 Mon Sep 17 00:00:00 2001
From: James Dark <65899899+jhdark@users.noreply.github.com>
Date: Fri, 27 Mar 2026 16:17:39 -0400
Subject: [PATCH 15/19] Update src/festim/exports/advection_flux.py
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com>
---
 src/festim/exports/advection_flux.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/festim/exports/advection_flux.py b/src/festim/exports/advection_flux.py
index aa2bb5742..f739255cf 100644
--- a/src/festim/exports/advection_flux.py
+++ b/src/festim/exports/advection_flux.py
@@ -7,7 +7,7 @@
 from festim.subdomain.surface_subdomain import SurfaceSubdomain
 
 
-class Advection_flux(SurfaceFlux):
+class AdvectionFlux(SurfaceFlux):
     """Computes the advective flux of a field on a given surface
 
     Args:

From 64b4debb7e61277552d4162a325015b6cce6d6e5 Mon Sep 17 00:00:00 2001
From: James Dark <65899899+jhdark@users.noreply.github.com>
Date: Fri, 27 Mar 2026 16:19:10 -0400
Subject: [PATCH 16/19] Apply suggestions from code review
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com>
---
 src/festim/exports/advection_flux.py                 | 2 +-
 test/system_tests/test_advection_term_integration.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/festim/exports/advection_flux.py b/src/festim/exports/advection_flux.py
index f739255cf..e39ce2d61 100644
--- a/src/festim/exports/advection_flux.py
+++ b/src/festim/exports/advection_flux.py
@@ -13,9 +13,9 @@ class AdvectionFlux(SurfaceFlux):
     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
-        velocity_field: velocity field for which the advective flux is computed
 
     Attributes:
         see `festim.SurfaceFlux`
diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index c902bf319..b28e5912e 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/test/system_tests/test_advection_term_integration.py
@@ -288,7 +288,7 @@ def test_flux_conservation_advection():
 
     my_model = F.HydrogenTransportProblem()
 
-    test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 200, 200)
+    test_mesh_2d = create_unit_square(MPI.COMM_WORLD, 15, 15)
     my_model.mesh = F.Mesh(test_mesh_2d)
 
     # common festim objects
@@ -357,4 +357,4 @@ def test_flux_conservation_advection():
         + advection_flux_bottom.data[-1]
     )
 
-    assert np.isclose(total_source, total_outflux)
+    assert np.isclose(total_source, total_outflux, rtol=1e-3)

From 0224269576cbc6726e57afcf6db277c4dea3182d Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Fri, 27 Mar 2026 16:25:15 -0400
Subject: [PATCH 17/19] fix name of class

---
 src/festim/__init__.py                               | 2 +-
 src/festim/exports/__init__.py                       | 4 ++--
 test/system_tests/test_advection_term_integration.py | 8 ++++----
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index 5ac8afd88..8c41b4cd6 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -24,7 +24,7 @@
 from .coupled_heat_hydrogen_problem import (
     CoupledTransientHeatTransferHydrogenTransport,
 )
-from .exports.advection_flux import Advection_flux
+from .exports.advection_flux import AdvectionFlux
 from .exports.average_surface import AverageSurface
 from .exports.average_volume import AverageVolume
 from .exports.maximum_surface import MaximumSurface
diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py
index cbfb5a1af..82f6923e1 100644
--- a/src/festim/exports/__init__.py
+++ b/src/festim/exports/__init__.py
@@ -1,4 +1,4 @@
-from .advection_flux import Advection_flux
+from .advection_flux import AdvectionFlux
 from .average_surface import AverageSurface
 from .average_volume import AverageVolume
 from .maximum_surface import MaximumSurface
@@ -15,7 +15,7 @@
 from .xdmf import XDMFExport
 
 __all__ = [
-    "Advection_flux",
+    "AdvectionFlux",
     "AverageSurface",
     "AverageVolume",
     "ExportBaseClass",
diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index b28e5912e..e16161165 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/test/system_tests/test_advection_term_integration.py
@@ -324,16 +324,16 @@ def test_flux_conservation_advection():
 
     my_model.settings = F.Settings(transient=False, atol=1e-10, rtol=1e-10)
 
-    advection_flux_top = F.Advection_flux(
+    advection_flux_top = F.Advectionflux(
         velocity_field=velocity_field, field=test_H, surface=top
     )
-    advection_flux_right = F.Advection_flux(
+    advection_flux_right = F.Advectionflux(
         velocity_field=velocity_field, field=test_H, surface=right
     )
-    advection_flux_left = F.Advection_flux(
+    advection_flux_left = F.Advectionflux(
         velocity_field=velocity_field, field=test_H, surface=left
     )
-    advection_flux_bottom = F.Advection_flux(
+    advection_flux_bottom = F.Advectionflux(
         velocity_field=velocity_field, field=test_H, surface=bottom
     )
     my_model.exports = [

From fdffddf2865fa0a307b0884aa1d1cdd320ea883d Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Fri, 27 Mar 2026 16:27:17 -0400
Subject: [PATCH 18/19] fix test

---
 test/system_tests/test_advection_term_integration.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py
index e16161165..eea350fcd 100644
--- a/test/system_tests/test_advection_term_integration.py
+++ b/test/system_tests/test_advection_term_integration.py
@@ -324,16 +324,16 @@ def test_flux_conservation_advection():
 
     my_model.settings = F.Settings(transient=False, atol=1e-10, rtol=1e-10)
 
-    advection_flux_top = F.Advectionflux(
+    advection_flux_top = F.AdvectionFlux(
         velocity_field=velocity_field, field=test_H, surface=top
     )
-    advection_flux_right = F.Advectionflux(
+    advection_flux_right = F.AdvectionFlux(
         velocity_field=velocity_field, field=test_H, surface=right
     )
-    advection_flux_left = F.Advectionflux(
+    advection_flux_left = F.AdvectionFlux(
         velocity_field=velocity_field, field=test_H, surface=left
     )
-    advection_flux_bottom = F.Advectionflux(
+    advection_flux_bottom = F.AdvectionFlux(
         velocity_field=velocity_field, field=test_H, surface=bottom
     )
     my_model.exports = [

From 3cd058e2bb2fda1183a92769ea4b42f63470288d Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Sun, 29 Mar 2026 16:26:50 -0400
Subject: [PATCH 19/19] added details to doc strings

---
 src/festim/exports/advection_flux.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/festim/exports/advection_flux.py b/src/festim/exports/advection_flux.py
index e39ce2d61..c2dd945d8 100644
--- a/src/festim/exports/advection_flux.py
+++ b/src/festim/exports/advection_flux.py
@@ -10,6 +10,11 @@
 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
@@ -18,7 +23,11 @@ class AdvectionFlux(SurfaceFlux):
             exported
 
     Attributes:
-        see `festim.SurfaceFlux`
+        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