From df074f4166b069933fc3d1c2faf60d7dbb810b42 Mon Sep 17 00:00:00 2001 From: jhdark Date: Mon, 2 Mar 2026 15:32:31 -0500 Subject: [PATCH 01/12] 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 Date: Mon, 2 Mar 2026 15:32:40 -0500 Subject: [PATCH 02/12] 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 Date: Mon, 2 Mar 2026 15:33:15 -0500 Subject: [PATCH 03/12] 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 Date: Mon, 2 Mar 2026 15:34:09 -0500 Subject: [PATCH 04/12] 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 Date: Tue, 3 Mar 2026 08:43:30 -0500 Subject: [PATCH 05/12] 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 Date: Tue, 3 Mar 2026 09:22:17 -0500 Subject: [PATCH 06/12] 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 Date: Tue, 3 Mar 2026 09:22:48 -0500 Subject: [PATCH 07/12] 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 Date: Tue, 3 Mar 2026 13:04:21 -0500 Subject: [PATCH 08/12] 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 Date: Tue, 3 Mar 2026 16:33:44 -0500 Subject: [PATCH 09/12] 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 Date: Tue, 3 Mar 2026 16:50:49 -0500 Subject: [PATCH 10/12] 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 Date: Tue, 3 Mar 2026 16:50:56 -0500 Subject: [PATCH 11/12] 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 Date: Tue, 3 Mar 2026 16:51:03 -0500 Subject: [PATCH 12/12] 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)