From 34ba5bf97828b076f13631ef06969c01890b16e3 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 2 Oct 2024 17:16:32 -0400 Subject: [PATCH 1/5] density is space and time dependent + test --- festim/hydrogen_transport_problem.py | 16 ++++- festim/species.py | 63 ++++++++++++++-- .../test_non_homogeneous_density.py | 71 +++++++++++++++++++ 3 files changed, 143 insertions(+), 7 deletions(-) create mode 100644 test/system_tests/test_non_homogeneous_density.py diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 7f5d6a981..ad6356f78 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -217,6 +217,14 @@ def initialise(self): self.assign_functions_to_species() self.t = fem.Constant(self.mesh.mesh, 0.0) + + for reaction in self.reactions: + for reactant in reaction.reactant: + if isinstance(reactant, F.ImplicitSpecies): + reactant.create_value_fenics( + mesh=self.mesh.mesh, + t=self.t, + ) if self.settings.transient: # TODO should raise error if no stepsize is provided # TODO Should this be an attribute of festim.Stepsize? @@ -648,12 +656,16 @@ def create_formulation(self): def update_time_dependent_values(self): super().update_time_dependent_values() + t = float(self.t) + + for reaction in self.reactions: + for reactant in reaction.reactant: + if isinstance(reactant, F.ImplicitSpecies): + reactant.update_density(t=t) if not self.temperature_time_dependent: return - t = float(self.t) - if isinstance(self.temperature_fenics, fem.Constant): self.temperature_fenics.value = self.temperature(t=t) elif isinstance(self.temperature_fenics, fem.Function): diff --git a/festim/species.py b/festim/species.py index 04ac3f25f..468a77d4c 100644 --- a/festim/species.py +++ b/festim/species.py @@ -1,4 +1,6 @@ -from typing import List +from typing import List, Union +import ufl +from dolfinx import fem import festim as F @@ -167,7 +169,7 @@ class ImplicitSpecies: c = n - others Args: - n (float): the total concentration of the species + n (Union[float, callable]): the total concentration of the species others (List[Species]): the list of species from which the implicit species concentration is computed (c = n - others) name (str, optional): a name given to the species. Defaults to None. @@ -178,12 +180,12 @@ class ImplicitSpecies: others (List[Species]): the list of species from which the implicit species concentration is computed (c = n - others) concentration (form): the concentration of the species - + value_fenics (): the total concentration as a fenics object """ def __init__( self, - n: float, + n: Union[float, callable], others: List[Species] = None, name: str = None, ) -> None: @@ -205,7 +207,58 @@ def concentration(self): raise ValueError( f"Cannot compute concentration of {self.name} because {other.name} has no solution" ) - return self.n - sum([other.solution for other in self.others]) + return self.value_fenics - sum([other.solution for other in self.others]) + + def create_value_fenics(self, mesh, t: fem.Constant): + """Creates the value of the density as a fenics object and sets it to + self.value_fenics. + If the value is a constant, it is converted to a fenics.Constant. + If the value is a function of t, it is converted to a fenics.Constant. + Otherwise, it is converted to a ufl Expression + + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + t (dolfinx.fem.Constant): the time + """ + x = ufl.SpatialCoordinate(mesh) + + if isinstance(self.n, (int, float)): + self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.n) + + elif isinstance(self.n, (fem.Function, ufl.core.expr.Expr)): + self.value_fenics = self.n + + elif callable(self.n): + arguments = self.n.__code__.co_varnames + + if "t" in arguments and "x" not in arguments and "T" not in arguments: + # only t is an argument + if not isinstance(self.n(t=float(t)), (float, int)): + raise ValueError( + f"self.value should return a float or an int, not {type(self.n(t=float(t)))} " + ) + self.value_fenics = F.as_fenics_constant( + mesh=mesh, value=self.n(t=float(t)) + ) + else: + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + kwargs["x"] = x + + self.value_fenics = self.n(**kwargs) + + def update_density(self, t): + """Updates the density value (only if the density is a function of time only) + + Args: + t (float): the time + """ + if callable(self.n): + arguments = self.n.__code__.co_varnames + if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: + self.value_fenics.value = self.n(t=t) def find_species_from_name(name: str, species: list): diff --git a/test/system_tests/test_non_homogeneous_density.py b/test/system_tests/test_non_homogeneous_density.py new file mode 100644 index 000000000..1a2c61f5a --- /dev/null +++ b/test/system_tests/test_non_homogeneous_density.py @@ -0,0 +1,71 @@ +import festim as F +import ufl +import numpy as np +import pytest + +step_function_space = lambda x: ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0) +step_function_time_non_homogeneous = lambda x, t: ufl.conditional( + ufl.gt(t, 50), 10 - 10 * x[0], 0.0 +) +step_function_time_homogeneous = lambda t: 10.0 if t > 50 else 2.0 +simple_space = lambda x: 10 - 10 * x[0] + + +@pytest.mark.parametrize( + "density_func, expected_value", + [ + (step_function_space, 5.0), + (simple_space, 5.0), + (step_function_time_non_homogeneous, 5.0), + (step_function_time_homogeneous, 10.0), + ], +) +def test_non_homogeneous_density(density_func, expected_value): + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) + my_mat = F.Material(name="mat", D_0=1, E_D=0) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=1, x=0) + right = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + trapped_H = F.Species("trapped_H", mobile=False) + + empty = F.ImplicitSpecies(n=density_func, name="empty", others=[trapped_H]) + my_model.species = [H, trapped_H] + + my_model.reactions = [ + F.Reaction( + reactant=[H, empty], + product=trapped_H, + k_0=1, + E_k=0, + p_0=0, + E_p=0, + volume=vol, + ) + ] + + my_model.temperature = 600 + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=1, species=H), + F.DirichletBC(subdomain=right, value=1, species=H), + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=100) + my_model.settings.stepsize = 1 + + total_trapped = F.TotalVolume(field=trapped_H, volume=vol) + my_model.exports = [ + F.VTXExport(filename="trapped_c.bp", field=trapped_H), + F.VTXExport(filename="c.bp", field=H), + total_trapped, + ] + + my_model.initialise() + my_model.run() + + assert np.isclose(total_trapped.value, expected_value) From 86caab0f951ba479c8b2aae0d8b48f892a198b2a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 3 Oct 2024 08:33:20 -0400 Subject: [PATCH 2/5] fixed test + refactoring + better integration with Trap class --- festim/hydrogen_transport_problem.py | 19 ++++++++++++------- festim/species.py | 7 +++++-- test/test_species.py | 2 ++ 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ad6356f78..b766f297f 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -218,13 +218,8 @@ def initialise(self): self.t = fem.Constant(self.mesh.mesh, 0.0) - for reaction in self.reactions: - for reactant in reaction.reactant: - if isinstance(reactant, F.ImplicitSpecies): - reactant.create_value_fenics( - mesh=self.mesh.mesh, - t=self.t, - ) + self.create_implicit_species_value_fenics() + if self.settings.transient: # TODO should raise error if no stepsize is provided # TODO Should this be an attribute of festim.Stepsize? @@ -241,6 +236,16 @@ def initialise(self): self.create_solver() self.initialise_exports() + def create_implicit_species_value_fenics(self): + """For each implicit species, create the value_fenics""" + for reaction in self.reactions: + for reactant in reaction.reactant: + if isinstance(reactant, F.ImplicitSpecies): + reactant.create_value_fenics( + mesh=self.mesh.mesh, + t=self.t, + ) + def create_species_from_traps(self): """Generate a species and reaction per trap defined in self.traps""" diff --git a/festim/species.py b/festim/species.py index 468a77d4c..5c9172305 100644 --- a/festim/species.py +++ b/festim/species.py @@ -86,6 +86,7 @@ class Trap(Species): volume (F.VolumeSubdomain1D): The volume subdomain where the trap is. trapped_concentration (F.Species): The immobile trapped concentration trap_reaction (F.Reaction): The reaction for trapping the mobile conc. + empty_trap_sites (F.ImplicitSpecies): The implicit species for the empty trap sites Usage: >>> import festim as F @@ -151,10 +152,12 @@ def __init__( def create_species_and_reaction(self): """create the immobile trapped species object and the reaction for trapping""" self.trapped_concentration = F.Species(name=self.name, mobile=False) - trap_site = F.ImplicitSpecies(n=self.n, others=[self.trapped_concentration]) + self.empty_trap_sites = F.ImplicitSpecies( + n=self.n, others=[self.trapped_concentration] + ) self.reaction = F.Reaction( - reactant=[self.mobile_species, trap_site], + reactant=[self.mobile_species, self.empty_trap_sites], product=self.trapped_concentration, k_0=self.k_0, E_k=self.E_k, diff --git a/test/test_species.py b/test/test_species.py index d2e672923..2f740f11e 100644 --- a/test/test_species.py +++ b/test/test_species.py @@ -88,6 +88,8 @@ def test_implicit_species_concentration(): species1.solution = Function(V) species2.solution = Function(V) + implicit_species.create_value_fenics(mesh, t=0.0) + # test the concentration of the implicit species expected_concentration = implicit_species.n - ( species1.solution + species2.solution From 3259bd26bbd62d962bc7221dd7383a66d6e45c68 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 3 Oct 2024 08:41:11 -0400 Subject: [PATCH 3/5] value fenics instead of n --- test/test_species.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_species.py b/test/test_species.py index 2f740f11e..8e09ce532 100644 --- a/test/test_species.py +++ b/test/test_species.py @@ -91,7 +91,7 @@ def test_implicit_species_concentration(): implicit_species.create_value_fenics(mesh, t=0.0) # test the concentration of the implicit species - expected_concentration = implicit_species.n - ( + expected_concentration = implicit_species.value_fenics - ( species1.solution + species2.solution ) assert implicit_species.concentration == expected_concentration From 0b9442522ab93f20e5794f445ed15c11b5156c6c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 3 Oct 2024 09:06:49 -0400 Subject: [PATCH 4/5] increased coverage --- festim/species.py | 3 + .../test_non_homogeneous_density.py | 66 +++++++++++++++++++ test/test_species.py | 15 +++++ 3 files changed, 84 insertions(+) diff --git a/festim/species.py b/festim/species.py index 5c9172305..4a9e509ee 100644 --- a/festim/species.py +++ b/festim/species.py @@ -258,6 +258,9 @@ def update_density(self, t): Args: t (float): the time """ + if isinstance(self.n, (fem.Function, ufl.core.expr.Expr)): + return + if callable(self.n): arguments = self.n.__code__.co_varnames if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: diff --git a/test/system_tests/test_non_homogeneous_density.py b/test/system_tests/test_non_homogeneous_density.py index 1a2c61f5a..1cf36ba2e 100644 --- a/test/system_tests/test_non_homogeneous_density.py +++ b/test/system_tests/test_non_homogeneous_density.py @@ -2,6 +2,8 @@ import ufl import numpy as np import pytest +import dolfinx.fem as fem +import basix step_function_space = lambda x: ufl.conditional(ufl.gt(x[0], 0.5), 10, 0.0) step_function_time_non_homogeneous = lambda x, t: ufl.conditional( @@ -69,3 +71,67 @@ def test_non_homogeneous_density(density_func, expected_value): my_model.run() assert np.isclose(total_trapped.value, expected_value) + + +def test_density_as_function(): + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) + my_mat = F.Material(name="mat", D_0=1, E_D=0) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=1, x=0) + right = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + trapped_H = F.Species("trapped_H", mobile=False) + + degree = 1 + element_CG = basix.ufl.element( + basix.ElementFamily.P, + my_model.mesh.mesh.basix_cell(), + degree, + basix.LagrangeVariant.equispaced, + ) + + density = fem.Function(fem.functionspace(my_model.mesh.mesh, element_CG)) + density.interpolate(lambda x: 10 - 10 * x[0]) + + empty = F.ImplicitSpecies(n=density, name="empty", others=[trapped_H]) + my_model.species = [H, trapped_H] + + my_model.reactions = [ + F.Reaction( + reactant=[H, empty], + product=trapped_H, + k_0=1, + E_k=0, + p_0=0, + E_p=0, + volume=vol, + ) + ] + + my_model.temperature = 600 + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=1, species=H), + F.DirichletBC(subdomain=right, value=1, species=H), + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=100) + my_model.settings.stepsize = 1 + + total_trapped = F.TotalVolume(field=trapped_H, volume=vol) + my_model.exports = [ + F.VTXExport(filename="trapped_c.bp", field=trapped_H), + F.VTXExport(filename="c.bp", field=H), + total_trapped, + ] + + my_model.initialise() + my_model.run() + + expected_value = 5.0 + + assert np.isclose(total_trapped.value, expected_value) diff --git a/test/test_species.py b/test/test_species.py index 8e09ce532..598f1d8b1 100644 --- a/test/test_species.py +++ b/test/test_species.py @@ -143,3 +143,18 @@ def test_create_species_and_reaction(): # TEST assert isinstance(my_trap.trapped_concentration, F.Species) assert isinstance(my_trap.reaction, F.Reaction) + + +def test_implicit_species_wrong_type(): + """Test that a ValueError is raised when the value of n in an ImplicitSpecies is + a function of time only but doesn't return a float or an int. + """ + mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10) + + density = lambda t: "coucou" + species = F.ImplicitSpecies(n=density) + + with pytest.raises( + ValueError, match="self.value should return a float or an int, not " + ): + species.create_value_fenics(mesh=mesh, t=0.0) From 235e75869191f33ed2f3be607540e167f61a5747 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 25 Oct 2024 09:29:55 -0400 Subject: [PATCH 5/5] update dolfinx version --- .github/workflows/ci_conda.yml | 2 +- .github/workflows/ci_docker.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci_conda.yml b/.github/workflows/ci_conda.yml index 0af9bae05..9c3d0055d 100644 --- a/.github/workflows/ci_conda.yml +++ b/.github/workflows/ci_conda.yml @@ -20,7 +20,7 @@ jobs: - name: Create Conda environment shell: bash -l {0} run: | - conda install -c conda-forge fenics-dolfinx=0.8.0 + conda install -c conda-forge fenics-dolfinx=0.9.0 - name: Install local package and dependencies shell: bash -l {0} diff --git a/.github/workflows/ci_docker.yml b/.github/workflows/ci_docker.yml index 7074e6e8f..e104da607 100644 --- a/.github/workflows/ci_docker.yml +++ b/.github/workflows/ci_docker.yml @@ -6,7 +6,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - container_version: [v0.8.0, nightly] + container_version: [v0.9.0, nightly] container: dolfinx/dolfinx:${{ matrix.container_version }} steps: - name: Checkout code