diff --git a/src/festim/__init__.py b/src/festim/__init__.py index d6df5c269..5ac8afd88 100644 --- a/src/festim/__init__.py +++ b/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 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 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", diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py new file mode 100644 index 000000000..f2eef129d --- /dev/null +++ b/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 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..aa2bb5742 --- /dev/null +++ b/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) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 805511d4a..fb23c2c6c 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: @@ -861,15 +867,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 @@ -990,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() 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)