From 36ed2a64d992ce10f6dfddb7f4b4fd1afb485e6d Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 09:53:17 -0400 Subject: [PATCH 01/27] accept temperature as callable --- festim/hydrogen_transport_problem.py | 80 ++++++++++++++++++++++++++-- 1 file changed, 76 insertions(+), 4 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ff171bca3..e18da565d 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -99,6 +99,7 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] + self.temperature_time_dependent = False @property def temperature(self): @@ -108,8 +109,31 @@ def temperature(self): def temperature(self, value): if value is None: self._temperature = value - else: + return + if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): + raise TypeError( + f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" + ) + self._temperature = value + + @property + def temperature(self): + return self._temperature + + @temperature.setter + def temperature(self, value): + if value is None: + self._temperature = value + elif isinstance(value, (float, int)): self._temperature = F.as_fenics_constant(value, self.mesh.mesh) + elif isinstance(value, fem.Constant): + self._temperature = value + elif callable(value): + self._temperature = value + else: + raise TypeError( + "Temperature must be a float, int, fem.Constant or callable" + ) @property def multispecies(self): @@ -123,11 +147,47 @@ def initialise(self): self.t = fem.Constant(self.mesh.mesh, 0.0) self.dt = self.settings.stepsize.get_dt(self.mesh.mesh) + if isinstance(self.temperature, fem.Constant): + pass + else: + self.define_temperature() + self.define_boundary_conditions() self.create_formulation() self.create_solver() self.defing_export_writers() + def define_temperature(self): + """If temperature value is given only dependent on t, create a + fem.Constant else create an expression to be updated and a function""" + self.temperature_value = self.temperature + arguments = self.temperature.__code__.co_varnames + if "t" in arguments and "x" not in arguments and "T" not in arguments: + # only t is an argument + self.temperature = F.as_fenics_constant( + mesh=self.mesh.mesh, value=self.temperature_value(t=float(self.t)) + ) + self.temperature_time_dependent = True + else: + x = ufl.SpatialCoordinate(self.mesh.mesh) + arguments = self.temperature.__code__.co_varnames + self.temperature_function = fem.Function(self.function_space) + kwargs = {} + if "t" in arguments: + kwargs["t"] = self.t + self.temperature_time_dependent = True + if "x" in arguments: + kwargs["x"] = x + + # store the expression of the boundary condition + # to update the value_fenics later + self.temperature_expr = fem.Expression( + self.temperature_value(**kwargs), + self.function_space.element.interpolation_points(), + ) + self.temperature_function.interpolate(self.temperature_expr) + self.temperature = self.temperature_function + def defing_export_writers(self): """Defines the export writers of the model, if field is given as a string, find species object in self.species""" @@ -375,9 +435,7 @@ def run(self): progress.update(self.dt.value) self.t.value += self.dt.value - # update boundary conditions - for bc in self.boundary_conditions: - bc.update(float(self.t)) + self.update_time_dependent_values(float(self.t)) self.solver.solve(self.u) @@ -424,3 +482,17 @@ def run(self): flux_values = [flux_values_1, flux_values_2] return times, flux_values + + def update_time_dependent_values(self, t): + """Updates the time dependent values of the model""" + + # update temperature if time dependent + if self.temperature_time_dependent: + if isinstance(self.temperature, fem.Constant): + self.temperature.value = self.temperature_value(t=t) + else: + self.temperature.interpolate(self.temperature_expr) + + # update boundary conditions + for bc in self.boundary_conditions: + bc.update(t) From 7064952e7781a764da16d8a1d468f45fcc0a3abb Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 09:53:32 -0400 Subject: [PATCH 02/27] temperature needs to be given --- test/test_dirichlet_bc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index f7f5b0e05..98d9cd731 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -297,6 +297,7 @@ def test_species_predefined(): my_model.species = [F.Species("H")] my_bc = F.DirichletBC(subdomain, 1.0, "J") my_model.boundary_conditions = [my_bc] + my_model.temperature = 1 my_model.settings = F.Settings(atol=1, rtol=0.1) my_model.settings.stepsize = 1 From 39a09cca384bdf5a68fd5981734e4a420bbcd6f0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 09:53:44 -0400 Subject: [PATCH 03/27] added tests --- test/test_temperature.py | 101 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 99 insertions(+), 2 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index 827a6960c..f5ae06e5f 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -3,13 +3,15 @@ import numpy as np import pytest import ufl +from ufl.conditional import Conditional test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) +dummy_mat = F.Material(D_0=1, E_D=1, name="dummy_mat") @pytest.mark.parametrize( - "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]] + "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", lambda x: 2 * x[0]] ) def test_temperature_type_and_processing(value): """Test that the temperature type is correctly set""" @@ -17,8 +19,103 @@ def test_temperature_type_and_processing(value): my_model.mesh = test_mesh if not isinstance(value, (fem.Constant, int, float)): - with pytest.raises(TypeError): + if callable(value): my_model.temperature = value + else: + with pytest.raises(TypeError): + my_model.temperature = value else: my_model.temperature = value assert isinstance(my_model.temperature, fem.Constant) + + +@pytest.mark.parametrize( + "value", + [ + 1.0, + lambda t: t, + lambda t: 1.0 + t, + lambda x: 1.0 + x[0], + lambda x, t: 1.0 + x[0] + t, + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ], +) +def test_time_dependent_temperature_attribute(value): + """Test that the temperature_time_dependent attribute is correctly set""" + subdomain = F.SurfaceSubdomain1D(1, x=0) + vol_subdomain = F.VolumeSubdomain1D( + 1, borders=[0, 1], material=F.Material(1, 1, "my_mat") + ) + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(vertices=np.linspace(0.0, 1.0)), + subdomains=[vol_subdomain, subdomain], + ) + my_model.species = [F.Species("H")] + my_bc = F.DirichletBC(subdomain, 0, "H") + my_model.boundary_conditions = [my_bc] + + my_model.temperature = value + + my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=2) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + my_model.initialise() + + if isinstance(value, (float, int)): + assert not my_model.temperature_time_dependent + else: + arguments = value.__code__.co_varnames + if "t" in arguments: + assert my_model.temperature_time_dependent + else: + assert not my_model.temperature_time_dependent + + +@pytest.mark.parametrize( + "value", + [ + lambda t: t, + lambda t: 1.0 + t, + lambda x, t: 1.0 + x[0] + t, + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ], +) +def test_temperature_value_updates_with_HTransportProblem(value): + """Test that different time dependent callable functions can be applied to + the temperature value, asserting in each case they match an expected value""" + subdomain = F.SurfaceSubdomain1D(1, x=4) + vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat) + + my_model = F.HydrogenTransportProblem( + mesh=test_mesh, + subdomains=[vol_subdomain, subdomain], + ) + my_model.species = [F.Species("H")] + my_bc = F.DirichletBC(subdomain, 1.0, "H") + my_model.boundary_conditions = [my_bc] + + my_model.temperature = value + + my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=3) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + # RUN + my_model.initialise() + my_model.run() + + # TEST + if callable(value): + arguments = value.__code__.co_varnames + if "x" in arguments and "t" in arguments: + expected_value = value(x=np.array([subdomain.x]), t=3.0) + computed_value = my_model.temperature.vector.array[-1] + elif "x" in arguments: + expected_value = value(x=np.array([subdomain.x])) + computed_value = my_model.temperature.vector.array[-1] + elif "t" in arguments: + expected_value = value(t=3.0) + computed_value = float(my_model.temperature) + + if isinstance(expected_value, Conditional): + expected_value = float(expected_value) + assert np.isclose(computed_value, expected_value) From 83d87fc436c5513c877c8c9df9f510cd4631bd8a Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 09:56:04 -0400 Subject: [PATCH 04/27] typo --- festim/hydrogen_transport_problem.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e18da565d..5f0215c11 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -105,21 +105,6 @@ def __init__( def temperature(self): return self._temperature - @temperature.setter - def temperature(self, value): - if value is None: - self._temperature = value - return - if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): - raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" - ) - self._temperature = value - - @property - def temperature(self): - return self._temperature - @temperature.setter def temperature(self, value): if value is None: From c686eb110024e29baf759dbc4a503a557526f466 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 11:43:02 -0400 Subject: [PATCH 05/27] new attribute fenics_value --- festim/hydrogen_transport_problem.py | 100 +++++++++++++-------------- 1 file changed, 50 insertions(+), 50 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 5f0215c11..472ded390 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -99,26 +99,23 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] + self.temperature_fenics_value = None self.temperature_time_dependent = False @property - def temperature(self): - return self._temperature + def temperature_fenics_value(self): + return self._temperature_fenics_value - @temperature.setter - def temperature(self, value): + @temperature_fenics_value.setter + def temperature_fenics_value(self, value): if value is None: - self._temperature = value - elif isinstance(value, (float, int)): - self._temperature = F.as_fenics_constant(value, self.mesh.mesh) - elif isinstance(value, fem.Constant): - self._temperature = value - elif callable(value): - self._temperature = value - else: + self._temperature_fenics_value = value + return + if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): raise TypeError( - "Temperature must be a float, int, fem.Constant or callable" + f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" ) + self._temperature_fenics_value = value @property def multispecies(self): @@ -132,11 +129,7 @@ def initialise(self): self.t = fem.Constant(self.mesh.mesh, 0.0) self.dt = self.settings.stepsize.get_dt(self.mesh.mesh) - if isinstance(self.temperature, fem.Constant): - pass - else: - self.define_temperature() - + self.define_temperature() self.define_boundary_conditions() self.create_formulation() self.create_solver() @@ -145,33 +138,40 @@ def initialise(self): def define_temperature(self): """If temperature value is given only dependent on t, create a fem.Constant else create an expression to be updated and a function""" - self.temperature_value = self.temperature - arguments = self.temperature.__code__.co_varnames - if "t" in arguments and "x" not in arguments and "T" not in arguments: - # only t is an argument - self.temperature = F.as_fenics_constant( - mesh=self.mesh.mesh, value=self.temperature_value(t=float(self.t)) + # if temperature is a float or int, create a fem.Constant + if isinstance(self.temperature, (float, int)): + self._temperature_fenics_value = F.as_fenics_constant( + self.temperature, self.mesh.mesh ) - self.temperature_time_dependent = True - else: - x = ufl.SpatialCoordinate(self.mesh.mesh) + elif isinstance(self.temperature, fem.Constant): + self._temperature_fenics_value = self.temperature + + # if temperature is callable, process accordingly + elif callable(self.temperature): arguments = self.temperature.__code__.co_varnames - self.temperature_function = fem.Function(self.function_space) - kwargs = {} - if "t" in arguments: - kwargs["t"] = self.t + if "t" in arguments and "x" not in arguments and "T" not in arguments: + # only t is an argument + self.temperature_fenics_value = F.as_fenics_constant( + mesh=self.mesh.mesh, value=self.temperature(t=float(self.t)) + ) self.temperature_time_dependent = True - if "x" in arguments: - kwargs["x"] = x - - # store the expression of the boundary condition - # to update the value_fenics later - self.temperature_expr = fem.Expression( - self.temperature_value(**kwargs), - self.function_space.element.interpolation_points(), - ) - self.temperature_function.interpolate(self.temperature_expr) - self.temperature = self.temperature_function + else: + x = ufl.SpatialCoordinate(self.mesh.mesh) + self.temperature_fenics_value = fem.Function(self.function_space) + kwargs = {} + if "t" in arguments: + kwargs["t"] = self.t + self.temperature_time_dependent = True + if "x" in arguments: + kwargs["x"] = x + + # store the expression of the boundary condition + # to update the value_fenics later + self.temperature_expr = fem.Expression( + self.temperature(**kwargs), + self.function_space.element.interpolation_points(), + ) + self.temperature_fenics_value.interpolate(self.temperature_expr) def defing_export_writers(self): """Defines the export writers of the model, if field is given as @@ -321,7 +321,7 @@ def create_dirichletbc_form(self, bc): bc.create_value( mesh=self.mesh.mesh, - temperature=self.temperature, + temperature=self.temperature_fenics_value, function_space=function_space_value, t=self.t, ) @@ -370,7 +370,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, spe + self.mesh.mesh, self.temperature_fenics_value, spe ) self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) @@ -428,7 +428,7 @@ def run(self): # TODO remove this if not self.multispecies: D_D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] + self.mesh.mesh, self.temperature_fenics_value, self.species[0] ) cm = self.u self.species[0].post_processing_solution = self.u @@ -443,10 +443,10 @@ def run(self): cm_1, cm_2 = self.u.split() D_1 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] + self.mesh.mesh, self.temperature_fenics_value, self.species[0] ) D_2 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[1] + self.mesh.mesh, self.temperature_fenics_value, self.species[1] ) surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) @@ -473,10 +473,10 @@ def update_time_dependent_values(self, t): # update temperature if time dependent if self.temperature_time_dependent: - if isinstance(self.temperature, fem.Constant): - self.temperature.value = self.temperature_value(t=t) + if isinstance(self.temperature_fenics_value, fem.Constant): + self.temperature_fenics_value.value = self.temperature(t=t) else: - self.temperature.interpolate(self.temperature_expr) + self.temperature_fenics_value.interpolate(self.temperature_expr) # update boundary conditions for bc in self.boundary_conditions: From 64fe46c1f3c3e135860b50d773f066295d1e3428 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 11:43:08 -0400 Subject: [PATCH 06/27] fix tests --- test/test_temperature.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index f5ae06e5f..b55fc5fb6 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -15,18 +15,26 @@ ) def test_temperature_type_and_processing(value): """Test that the temperature type is correctly set""" - my_model = F.HydrogenTransportProblem() - my_model.mesh = test_mesh + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(vertices=np.linspace(0.0, 1.0)), + subdomains=[ + F.VolumeSubdomain1D(1, borders=[0, 1], material=F.Material(1, 1, "my_mat")) + ], + species=[F.Species("H")], + ) + my_model.temperature = value + my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=2) + my_model.settings.stepsize = F.Stepsize(initial_value=1) if not isinstance(value, (fem.Constant, int, float)): if callable(value): - my_model.temperature = value + my_model.initialise() else: with pytest.raises(TypeError): - my_model.temperature = value + my_model.initialise() else: - my_model.temperature = value - assert isinstance(my_model.temperature, fem.Constant) + my_model.initialise() + assert isinstance(my_model.temperature_fenics_value, fem.Constant) @pytest.mark.parametrize( @@ -108,13 +116,13 @@ def test_temperature_value_updates_with_HTransportProblem(value): arguments = value.__code__.co_varnames if "x" in arguments and "t" in arguments: expected_value = value(x=np.array([subdomain.x]), t=3.0) - computed_value = my_model.temperature.vector.array[-1] + computed_value = my_model.temperature_fenics_value.vector.array[-1] elif "x" in arguments: expected_value = value(x=np.array([subdomain.x])) - computed_value = my_model.temperature.vector.array[-1] + computed_value = my_model.temperature_fenics_value.vector.array[-1] elif "t" in arguments: expected_value = value(t=3.0) - computed_value = float(my_model.temperature) + computed_value = float(my_model.temperature_fenics_value) if isinstance(expected_value, Conditional): expected_value = float(expected_value) From 10bdcc0112fef7b8a4e7175ad163ea07787bb933 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:15:51 -0400 Subject: [PATCH 07/27] updates from comments --- festim/hydrogen_transport_problem.py | 117 +++++++++++++++++++-------- 1 file changed, 83 insertions(+), 34 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 472ded390..dfa256e73 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -22,7 +22,8 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float or fem.Constant): the temperature of the model + temperature (float, int, fem.Constant or callable): the temperature of + the model sources (list of festim.Source): the hydrogen sources of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model @@ -33,7 +34,8 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (fem.Constant): the temperature of the model + temperature (float, int, fem.Constant or callable): the temperature of + the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model @@ -47,7 +49,13 @@ class HydrogenTransportProblem: model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model - multispecies (bool): True if the model has more than one species + multispecies (bool): True if the model has more than one species. + temperature_fenics (fem.Constant or fem.Function): the + temperature of the model. Defaults to None + temperature_expr (fem.Expression): the expression of the temperature + that is used to update the temperature_fenics + temperature_time_dependent (bool): True if the temperature is time + Usage: >>> import festim as F @@ -99,23 +107,47 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] - self.temperature_fenics_value = None - self.temperature_time_dependent = False + self.temperature_fenics = None + + @property + def temperature(self): + return self._temperature + + @temperature.setter + def temperature(self, value): + if value is None: + self._temperature = value + elif isinstance(value, (float, int, fem.Constant)): + self._temperature = value + elif callable(value): + self._temperature = value + else: + raise TypeError(f"Value must be a float, int, fem.Constant or callable") @property - def temperature_fenics_value(self): - return self._temperature_fenics_value + def temperature_fenics(self): + return self._temperature_fenics - @temperature_fenics_value.setter - def temperature_fenics_value(self, value): + @temperature_fenics.setter + def temperature_fenics(self, value): if value is None: - self._temperature_fenics_value = value + self._temperature_fenics = value return if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): raise TypeError( f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" ) - self._temperature_fenics_value = value + self._temperature_fenics = value + + @property + def temperature_time_dependent(self): + if self.temperature is None: + return False + if callable(self.temperature): + arguments = self.temperature.__code__.co_varnames + return "t" in arguments + else: + return False @property def multispecies(self): @@ -136,42 +168,58 @@ def initialise(self): self.defing_export_writers() def define_temperature(self): - """If temperature value is given only dependent on t, create a - fem.Constant else create an expression to be updated and a function""" + """Sets the value of temperature_fenics_value. The type depends on + self.temperature. If self.temperature is a function on t only, create + a fem.Constant. Else, create an dolfinx.fem.Expression (stored in + self.temperature_expr) to be updated, a dolfinx.fem.Function object + is created from the Expression (stored in self.temperature_fenics_value) + """ + # check if temperature is None + if self.temperature is None: + raise TypeError("Temperature needs to be defined") + # if temperature is a float or int, create a fem.Constant - if isinstance(self.temperature, (float, int)): - self._temperature_fenics_value = F.as_fenics_constant( + elif isinstance(self.temperature, (float, int)): + self.temperature_fenics = F.as_fenics_constant( self.temperature, self.mesh.mesh ) - elif isinstance(self.temperature, fem.Constant): - self._temperature_fenics_value = self.temperature + # if temperature is a fem.Constant or function, pass it to temperature_fenics + elif isinstance(self.temperature, (fem.Constant, fem.Function)): + self.temperature_fenics = self.temperature # if temperature is callable, process accordingly elif callable(self.temperature): arguments = self.temperature.__code__.co_varnames if "t" in arguments and "x" not in arguments and "T" not in arguments: # only t is an argument - self.temperature_fenics_value = F.as_fenics_constant( + self.temperature_fenics = F.as_fenics_constant( mesh=self.mesh.mesh, value=self.temperature(t=float(self.t)) ) - self.temperature_time_dependent = True else: x = ufl.SpatialCoordinate(self.mesh.mesh) - self.temperature_fenics_value = fem.Function(self.function_space) + CG_temperature = basix.ufl.element( + basix.ElementFamily.P, + self.mesh.mesh.basix_cell(), + 1, + basix.LagrangeVariant.equispaced, + ) + function_space_temperature = fem.FunctionSpace( + self.mesh.mesh, CG_temperature + ) + self.temperature_fenics = fem.Function(function_space_temperature) kwargs = {} if "t" in arguments: kwargs["t"] = self.t - self.temperature_time_dependent = True if "x" in arguments: kwargs["x"] = x - # store the expression of the boundary condition - # to update the value_fenics later + # store the expression of the temperature + # to update the temperature_fenics later self.temperature_expr = fem.Expression( self.temperature(**kwargs), self.function_space.element.interpolation_points(), ) - self.temperature_fenics_value.interpolate(self.temperature_expr) + self.temperature_fenics.interpolate(self.temperature_expr) def defing_export_writers(self): """Defines the export writers of the model, if field is given as @@ -321,7 +369,7 @@ def create_dirichletbc_form(self, bc): bc.create_value( mesh=self.mesh.mesh, - temperature=self.temperature_fenics_value, + temperature=self.temperature_fenics, function_space=function_space_value, t=self.t, ) @@ -370,7 +418,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics_value, spe + self.mesh.mesh, self.temperature_fenics, spe ) self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) @@ -428,7 +476,7 @@ def run(self): # TODO remove this if not self.multispecies: D_D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics_value, self.species[0] + self.mesh.mesh, self.temperature_fenics, self.species[0] ) cm = self.u self.species[0].post_processing_solution = self.u @@ -443,10 +491,10 @@ def run(self): cm_1, cm_2 = self.u.split() D_1 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics_value, self.species[0] + self.mesh.mesh, self.temperature_fenics, self.species[0] ) D_2 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics_value, self.species[1] + self.mesh.mesh, self.temperature_fenics, self.species[1] ) surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) @@ -469,14 +517,15 @@ def run(self): return times, flux_values def update_time_dependent_values(self, t): - """Updates the time dependent values of the model""" - + """Updates the time dependent values of the model + liketemperature, boundary conditions, sources, etc. + """ # update temperature if time dependent if self.temperature_time_dependent: - if isinstance(self.temperature_fenics_value, fem.Constant): - self.temperature_fenics_value.value = self.temperature(t=t) + if isinstance(self.temperature_fenics, fem.Constant): + self.temperature_fenics.value = self.temperature(t=t) else: - self.temperature_fenics_value.interpolate(self.temperature_expr) + self.temperature_fenics.interpolate(self.temperature_expr) # update boundary conditions for bc in self.boundary_conditions: From c427b84ff02204b7a48cc7185ffade7474657af0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:32:14 -0400 Subject: [PATCH 08/27] constant not callable --- festim/hydrogen_transport_problem.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index dfa256e73..9266a64ff 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -143,6 +143,8 @@ def temperature_fenics(self, value): def temperature_time_dependent(self): if self.temperature is None: return False + if isinstance(self.temperature, fem.Constant): + return False if callable(self.temperature): arguments = self.temperature.__code__.co_varnames return "t" in arguments From 1a43f4a740e225076739ad3ccde985d6d7aabaf3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:32:35 -0400 Subject: [PATCH 09/27] modified tests from comments --- test/test_temperature.py | 59 ++++++++++------------------------------ 1 file changed, 15 insertions(+), 44 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index b55fc5fb6..31d4af50b 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -13,34 +13,25 @@ @pytest.mark.parametrize( "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", lambda x: 2 * x[0]] ) -def test_temperature_type_and_processing(value): +def test_temperature_type(value): """Test that the temperature type is correctly set""" my_model = F.HydrogenTransportProblem( - mesh=F.Mesh1D(vertices=np.linspace(0.0, 1.0)), - subdomains=[ - F.VolumeSubdomain1D(1, borders=[0, 1], material=F.Material(1, 1, "my_mat")) - ], - species=[F.Species("H")], + mesh=test_mesh, ) - my_model.temperature = value - my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=2) - my_model.settings.stepsize = F.Stepsize(initial_value=1) if not isinstance(value, (fem.Constant, int, float)): if callable(value): - my_model.initialise() + my_model.temperature = value else: with pytest.raises(TypeError): - my_model.initialise() - else: - my_model.initialise() - assert isinstance(my_model.temperature_fenics_value, fem.Constant) + my_model.temperature = value @pytest.mark.parametrize( "value", [ 1.0, + 1, lambda t: t, lambda t: 1.0 + t, lambda x: 1.0 + x[0], @@ -50,33 +41,16 @@ def test_temperature_type_and_processing(value): ) def test_time_dependent_temperature_attribute(value): """Test that the temperature_time_dependent attribute is correctly set""" - subdomain = F.SurfaceSubdomain1D(1, x=0) - vol_subdomain = F.VolumeSubdomain1D( - 1, borders=[0, 1], material=F.Material(1, 1, "my_mat") - ) - my_model = F.HydrogenTransportProblem( - mesh=F.Mesh1D(vertices=np.linspace(0.0, 1.0)), - subdomains=[vol_subdomain, subdomain], - ) - my_model.species = [F.Species("H")] - my_bc = F.DirichletBC(subdomain, 0, "H") - my_model.boundary_conditions = [my_bc] + my_model = F.HydrogenTransportProblem() my_model.temperature = value - my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=2) - my_model.settings.stepsize = F.Stepsize(initial_value=1) - - my_model.initialise() - - if isinstance(value, (float, int)): - assert not my_model.temperature_time_dependent - else: + if callable(value): arguments = value.__code__.co_varnames if "t" in arguments: assert my_model.temperature_time_dependent - else: - assert not my_model.temperature_time_dependent + else: + assert not my_model.temperature_time_dependent @pytest.mark.parametrize( @@ -97,16 +71,13 @@ def test_temperature_value_updates_with_HTransportProblem(value): my_model = F.HydrogenTransportProblem( mesh=test_mesh, subdomains=[vol_subdomain, subdomain], + species=[F.Species("H")], + settings=F.Settings(atol=1, rtol=0.1, final_time=3), ) - my_model.species = [F.Species("H")] - my_bc = F.DirichletBC(subdomain, 1.0, "H") - my_model.boundary_conditions = [my_bc] + my_model.settings.stepsize = F.Stepsize(initial_value=1) my_model.temperature = value - my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=3) - my_model.settings.stepsize = F.Stepsize(initial_value=1) - # RUN my_model.initialise() my_model.run() @@ -116,13 +87,13 @@ def test_temperature_value_updates_with_HTransportProblem(value): arguments = value.__code__.co_varnames if "x" in arguments and "t" in arguments: expected_value = value(x=np.array([subdomain.x]), t=3.0) - computed_value = my_model.temperature_fenics_value.vector.array[-1] + computed_value = my_model.temperature_fenics.vector.array[-1] elif "x" in arguments: expected_value = value(x=np.array([subdomain.x])) - computed_value = my_model.temperature_fenics_value.vector.array[-1] + computed_value = my_model.temperature_fenics.vector.array[-1] elif "t" in arguments: expected_value = value(t=3.0) - computed_value = float(my_model.temperature_fenics_value) + computed_value = float(my_model.temperature_fenics) if isinstance(expected_value, Conditional): expected_value = float(expected_value) From 27505376d8179811515a9e11ff12655c7ed5fb95 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:36:03 -0400 Subject: [PATCH 10/27] typo --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 9266a64ff..731def860 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -520,7 +520,7 @@ def run(self): def update_time_dependent_values(self, t): """Updates the time dependent values of the model - liketemperature, boundary conditions, sources, etc. + like temperature, boundary conditions, sources, etc. """ # update temperature if time dependent if self.temperature_time_dependent: From fd5ee959d04784cd664755a36414e0af9382446c Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:43:46 -0400 Subject: [PATCH 11/27] improve coverage --- festim/hydrogen_transport_problem.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 731def860..d87013330 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -117,12 +117,14 @@ def temperature(self): def temperature(self, value): if value is None: self._temperature = value - elif isinstance(value, (float, int, fem.Constant)): + elif isinstance(value, (float, int, fem.Constant, fem.Function)): self._temperature = value elif callable(value): self._temperature = value else: - raise TypeError(f"Value must be a float, int, fem.Constant or callable") + raise TypeError( + f"Value must be a float, int, fem.Constant, fem.Function, or callable" + ) @property def temperature_fenics(self): @@ -133,10 +135,6 @@ def temperature_fenics(self, value): if value is None: self._temperature_fenics = value return - if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): - raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}" - ) self._temperature_fenics = value @property From 9a547e56a45919b281f1abbbbb0f0d5057b60c27 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:49:42 -0400 Subject: [PATCH 12/27] update doc strings --- festim/hydrogen_transport_problem.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index d87013330..5ae4d1a71 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -22,8 +22,8 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float, int, fem.Constant or callable): the temperature of - the model + temperature (float, int, fem.Constant, fem.Function or callable): the + temperature of the model sources (list of festim.Source): the hydrogen sources of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model @@ -34,8 +34,8 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float, int, fem.Constant or callable): the temperature of - the model + temperature (float, int, fem.Constant, fem.Function or callable): the + temperature of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model @@ -51,7 +51,8 @@ class HydrogenTransportProblem: solver (dolfinx.nls.newton.NewtonSolver): the solver of the model multispecies (bool): True if the model has more than one species. temperature_fenics (fem.Constant or fem.Function): the - temperature of the model. Defaults to None + temperature of the model as a fenics object (fem.Constant or + fem.Function). temperature_expr (fem.Expression): the expression of the temperature that is used to update the temperature_fenics temperature_time_dependent (bool): True if the temperature is time From bead7e3c45f355e8207f89bb81a45d3f53a2c3ea Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:51:23 -0400 Subject: [PATCH 13/27] add test for temperture as None --- test/test_temperature.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_temperature.py b/test/test_temperature.py index 31d4af50b..b0382847c 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -32,6 +32,7 @@ def test_temperature_type(value): [ 1.0, 1, + None, lambda t: t, lambda t: 1.0 + t, lambda x: 1.0 + x[0], From 22de6a6b0290eac8bdc16f7f2a67ea3623144181 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 30 Oct 2023 16:58:50 -0400 Subject: [PATCH 14/27] additional test for coverage --- test/test_temperature.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_temperature.py b/test/test_temperature.py index b0382847c..9fbb5a8f8 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -99,3 +99,18 @@ def test_temperature_value_updates_with_HTransportProblem(value): if isinstance(expected_value, Conditional): expected_value = float(expected_value) assert np.isclose(computed_value, expected_value) + + +def test_TypeError_raised_when_temperature_not_defined(): + """Test that a type error when a model is initialised without + defining a temperature""" + my_model = F.HydrogenTransportProblem( + mesh=test_mesh, + subdomains=[F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat)], + species=[F.Species("H")], + settings=F.Settings(atol=1, rtol=0.1, final_time=3), + ) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + with pytest.raises(TypeError, match="Temperature needs to be defined"): + my_model.initialise() From dc3290e1c18426f2639ed4e067b1601f0212fcfe Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 11:14:29 -0400 Subject: [PATCH 15/27] add update time dependent values --- festim/hydrogen_transport_problem.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 7521d309a..666b17c61 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -468,7 +468,6 @@ def run(self): while self.t.value < self.settings.final_time: self.iterate() - if self.multispecies: self.flux_values = [self.flux_values_1, self.flux_values_2] @@ -481,9 +480,7 @@ def iterate( self.progress.update(self.dt.value) self.t.value += self.dt.value - # update boundary conditions - for bc in self.boundary_conditions: - bc.update(float(self.t)) + self.update_time_dependent_values(t=float(self.t)) self.solver.solve(self.u) @@ -527,3 +524,12 @@ def iterate( # update previous solution self.u_n.x.array[:] = self.u.x.array[:] + def update_time_dependent_values(self, t): + if self.temperature_time_dependent: + if isinstance(self.temperature_fenics, fem.Constant): + self.temperature_fenics.value = self.temperature(t=t) + else: + self.temperature_fenics.interpolate(self.temperature_expr) + + for bc in self.boundary_conditions: + bc.update(t) From a8d0e6fd5863eac033c49c0efe70874aa087e292 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 11:14:45 -0400 Subject: [PATCH 16/27] model integration test simpler --- test/test_temperature.py | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index 9fbb5a8f8..ac3c050de 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -4,6 +4,7 @@ import pytest import ufl from ufl.conditional import Conditional +import tqdm.autonotebook test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) @@ -60,7 +61,7 @@ def test_time_dependent_temperature_attribute(value): lambda t: t, lambda t: 1.0 + t, lambda x, t: 1.0 + x[0] + t, - lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), ], ) def test_temperature_value_updates_with_HTransportProblem(value): @@ -81,24 +82,30 @@ def test_temperature_value_updates_with_HTransportProblem(value): # RUN my_model.initialise() - my_model.run() + my_model.progress = tqdm.autonotebook.tqdm( + desc="Solving H transport problem", + total=my_model.settings.final_time, + unit_scale=True, + ) + # TODO get rid of these when post processing is implemented + my_model.flux_values, my_model.times = [], [] # TEST - if callable(value): + expected_values = [] + for i in range(3): arguments = value.__code__.co_varnames if "x" in arguments and "t" in arguments: - expected_value = value(x=np.array([subdomain.x]), t=3.0) - computed_value = my_model.temperature_fenics.vector.array[-1] - elif "x" in arguments: - expected_value = value(x=np.array([subdomain.x])) - computed_value = my_model.temperature_fenics.vector.array[-1] - elif "t" in arguments: - expected_value = value(t=3.0) - computed_value = float(my_model.temperature_fenics) + expected_values.append(float(value(x=np.array([subdomain.x]), t=i + 1))) + else: + expected_values.append(float(value(t=i + 1))) - if isinstance(expected_value, Conditional): - expected_value = float(expected_value) - assert np.isclose(computed_value, expected_value) + for i in range(3): + my_model.iterate() + if isinstance(my_model.temperature_fenics, fem.Constant): + computed_value = float(my_model.temperature_fenics) + else: + computed_value = my_model.temperature_fenics.vector.array[-1] + assert np.isclose(computed_value, expected_values[i]) def test_TypeError_raised_when_temperature_not_defined(): From 72271509e3f028c7aa5768034e1613e7cd2479dd Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 11:29:08 -0400 Subject: [PATCH 17/27] expected values as arguement --- test/test_temperature.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index ac3c050de..435238258 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -56,15 +56,18 @@ def test_time_dependent_temperature_attribute(value): @pytest.mark.parametrize( - "value", + "T_function, expected_values", [ - lambda t: t, - lambda t: 1.0 + t, - lambda x, t: 1.0 + x[0] + t, - lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), + (lambda t: t, [1.0, 2.0, 3.0]), + (lambda t: 1.0 + t, [2.0, 3.0, 4.0]), + (lambda x, t: 1.0 + x[0] + t, [6.0, 7.0, 8.0]), + ( + lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), + [104.0, 0.0, 0.0], + ), ], ) -def test_temperature_value_updates_with_HTransportProblem(value): +def test_temperature_value_updates_with_HTransportProblem(T_function, expected_values): """Test that different time dependent callable functions can be applied to the temperature value, asserting in each case they match an expected value""" subdomain = F.SurfaceSubdomain1D(1, x=4) @@ -78,7 +81,7 @@ def test_temperature_value_updates_with_HTransportProblem(value): ) my_model.settings.stepsize = F.Stepsize(initial_value=1) - my_model.temperature = value + my_model.temperature = T_function # RUN my_model.initialise() @@ -91,14 +94,6 @@ def test_temperature_value_updates_with_HTransportProblem(value): my_model.flux_values, my_model.times = [], [] # TEST - expected_values = [] - for i in range(3): - arguments = value.__code__.co_varnames - if "x" in arguments and "t" in arguments: - expected_values.append(float(value(x=np.array([subdomain.x]), t=i + 1))) - else: - expected_values.append(float(value(t=i + 1))) - for i in range(3): my_model.iterate() if isinstance(my_model.temperature_fenics, fem.Constant): From 068efb57b3c2c929da5d452bd0267be6ce41e877 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 12:09:30 -0400 Subject: [PATCH 18/27] use temperature function space in temperature expr --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 666b17c61..bd48a81fc 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -218,7 +218,7 @@ def define_temperature(self): # to update the temperature_fenics later self.temperature_expr = fem.Expression( self.temperature(**kwargs), - self.function_space.element.interpolation_points(), + function_space_temperature.element.interpolation_points(), ) self.temperature_fenics.interpolate(self.temperature_expr) From 67d30cc3ee39bf7ad817cadef6ff28e2fdb90d51 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Tue, 31 Oct 2023 13:14:49 -0400 Subject: [PATCH 19/27] 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> --- festim/hydrogen_transport_problem.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index bd48a81fc..535cd4478 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -23,7 +23,7 @@ class HydrogenTransportProblem: subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model temperature (float, int, fem.Constant, fem.Function or callable): the - temperature of the model + temperature of the model (K) sources (list of festim.Source): the hydrogen sources of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model @@ -35,7 +35,7 @@ class HydrogenTransportProblem: subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model temperature (float, int, fem.Constant, fem.Function or callable): the - temperature of the model + temperature of the model (K) boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model @@ -56,6 +56,7 @@ class HydrogenTransportProblem: temperature_expr (fem.Expression): the expression of the temperature that is used to update the temperature_fenics temperature_time_dependent (bool): True if the temperature is time + dependent Usage: @@ -177,7 +178,7 @@ def define_temperature(self): """ # check if temperature is None if self.temperature is None: - raise TypeError("Temperature needs to be defined") + raise ValueError("the temperature attribute needs to be defined") # if temperature is a float or int, create a fem.Constant elif isinstance(self.temperature, (float, int)): From f69ec67e14b01be971f6763d281616449acc5bbe Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 15:36:55 -0400 Subject: [PATCH 20/27] moved temperature tests to h-transport-problem tests --- test/test_h_transport_problem.py | 155 ++++++++++++++++++++++++++++--- test/test_temperature.py | 118 ----------------------- 2 files changed, 144 insertions(+), 129 deletions(-) delete mode 100644 test/test_temperature.py diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 8d3652541..7f6bda9e4 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1,11 +1,106 @@ import festim as F import tqdm.autonotebook import mpi4py.MPI as MPI -import dolfinx +import dolfinx.mesh +from dolfinx import fem, nls import ufl import numpy as np +import pytest + +test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) +x = ufl.SpatialCoordinate(test_mesh.mesh) +dummy_mat = F.Material(D_0=1, E_D=1, name="dummy_mat") + # TODO test all the methods in the class +@pytest.mark.parametrize( + "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", lambda x: 2 * x[0]] +) +def test_temperature_setter_type(value): + """Test that the temperature type is correctly set""" + my_model = F.HydrogenTransportProblem( + mesh=test_mesh, + ) + + if not isinstance(value, (fem.Constant, int, float)): + if callable(value): + my_model.temperature = value + else: + with pytest.raises(TypeError): + my_model.temperature = value + + +@pytest.mark.parametrize( + "value", + [ + 1.0, + 1, + None, + lambda t: t, + lambda t: 1.0 + t, + lambda x: 1.0 + x[0], + lambda x, t: 1.0 + x[0] + t, + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ], +) +def test_time_dependent_temperature_attribute(value): + """Test that the temperature_time_dependent attribute is correctly set""" + + my_model = F.HydrogenTransportProblem() + my_model.temperature = value + + if callable(value): + arguments = value.__code__.co_varnames + if "t" in arguments: + assert my_model.temperature_time_dependent + else: + assert not my_model.temperature_time_dependent + + +@pytest.mark.parametrize( + "value", + [ + 1.0, + 1, + None, + fem.Constant(test_mesh.mesh, 1.0), + lambda t: t, + lambda t: 1.0 + t, + lambda x: 1.0 + x[0], + lambda x, t: 1.0 + x[0] + t, + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ], +) +def test_define_temperature(value): + """Test that the define_temperature method correctly sets the + temperature_fenics attribute to either a fem.Constant or a + fem.Function and raise a ValueError temperature is None""" + + # BUILD + my_model = F.HydrogenTransportProblem(mesh=test_mesh) + my_model.t = fem.Constant(test_mesh.mesh, 0.0) + + my_model.temperature = value + + # TEST + if value is None: + with pytest.raises( + ValueError, match="the temperature attribute needs to be defined" + ): + my_model.define_temperature() + else: + # RUN + my_model.define_temperature() + + # TEST + if isinstance(value, (fem.Constant, int, float)): + assert isinstance(my_model.temperature_fenics, fem.Constant) + elif callable(value): + arguments = value.__code__.co_varnames + if "x" in arguments: + assert isinstance(my_model.temperature_fenics, fem.Function) + else: + assert isinstance(my_model.temperature_fenics, fem.Constant) def test_iterate(): @@ -21,14 +116,12 @@ def test_iterate(): unit_scale=True, ) - my_model.boundary_conditions = [] - mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) - V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) - my_model.u = dolfinx.fem.Function(V) - my_model.u_n = dolfinx.fem.Function(V) - my_model.dt = dolfinx.fem.Constant(mesh, 2.0) + V = fem.FunctionSpace(mesh, ("CG", 1)) + my_model.u = fem.Function(V) + my_model.u_n = fem.Function(V) + my_model.dt = fem.Constant(mesh, 2.0) v = ufl.TestFunction(V) source_value = 2.0 @@ -36,10 +129,10 @@ def test_iterate(): my_model.u - my_model.u_n ) / my_model.dt * v * ufl.dx - source_value * v * ufl.dx - problem = dolfinx.fem.petsc.NonlinearProblem(form, my_model.u, bcs=[]) - my_model.solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem) + problem = fem.petsc.NonlinearProblem(form, my_model.u, bcs=[]) + my_model.solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem) - my_model.t = dolfinx.fem.Constant(mesh, 0.0) + my_model.t = fem.Constant(mesh, 0.0) for i in range(10): # RUN @@ -55,4 +148,44 @@ def test_iterate(): expected_u_value = (i + 1) * float(my_model.dt) * source_value assert np.all(np.isclose(my_model.u.x.array, expected_u_value)) - assert np.all(np.isclose(my_model.u_n.x.array, expected_u_value)) + +@pytest.mark.parametrize( + "T_function, expected_values", + [ + (lambda t: t, [1.0, 2.0, 3.0]), + (lambda t: 1.0 + t, [2.0, 3.0, 4.0]), + (lambda x, t: 1.0 + x[0] + t, [6.0, 7.0, 8.0]), + ( + lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), + [104.0, 0.0, 0.0], + ), + ], +) +def test_update_time_dependent_values_temperature(T_function, expected_values): + """Test that different time-dependent callable functions for the + temperature are updated at each time step and match an expected value""" + + # BUILD + my_model = F.HydrogenTransportProblem( + mesh=test_mesh, + ) + my_model.t = fem.Constant(my_model.mesh.mesh, 0.0) + dt = fem.Constant(test_mesh.mesh, 1.0) + + my_model.temperature = T_function + + my_model.define_temperature() + + for i in range(3): + # RUN + my_model.t.value += dt.value + my_model.update_time_dependent_values() + + # TEST + if isinstance(my_model.temperature_fenics, fem.Constant): + computed_value = float(my_model.temperature_fenics) + print(computed_value) + else: + computed_value = my_model.temperature_fenics.vector.array[-1] + print(computed_value) + assert np.isclose(computed_value, expected_values[i]) diff --git a/test/test_temperature.py b/test/test_temperature.py deleted file mode 100644 index 435238258..000000000 --- a/test/test_temperature.py +++ /dev/null @@ -1,118 +0,0 @@ -import festim as F -from dolfinx import fem -import numpy as np -import pytest -import ufl -from ufl.conditional import Conditional -import tqdm.autonotebook - -test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) -x = ufl.SpatialCoordinate(test_mesh.mesh) -dummy_mat = F.Material(D_0=1, E_D=1, name="dummy_mat") - - -@pytest.mark.parametrize( - "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", lambda x: 2 * x[0]] -) -def test_temperature_type(value): - """Test that the temperature type is correctly set""" - my_model = F.HydrogenTransportProblem( - mesh=test_mesh, - ) - - if not isinstance(value, (fem.Constant, int, float)): - if callable(value): - my_model.temperature = value - else: - with pytest.raises(TypeError): - my_model.temperature = value - - -@pytest.mark.parametrize( - "value", - [ - 1.0, - 1, - None, - lambda t: t, - lambda t: 1.0 + t, - lambda x: 1.0 + x[0], - lambda x, t: 1.0 + x[0] + t, - lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), - ], -) -def test_time_dependent_temperature_attribute(value): - """Test that the temperature_time_dependent attribute is correctly set""" - - my_model = F.HydrogenTransportProblem() - my_model.temperature = value - - if callable(value): - arguments = value.__code__.co_varnames - if "t" in arguments: - assert my_model.temperature_time_dependent - else: - assert not my_model.temperature_time_dependent - - -@pytest.mark.parametrize( - "T_function, expected_values", - [ - (lambda t: t, [1.0, 2.0, 3.0]), - (lambda t: 1.0 + t, [2.0, 3.0, 4.0]), - (lambda x, t: 1.0 + x[0] + t, [6.0, 7.0, 8.0]), - ( - lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), - [104.0, 0.0, 0.0], - ), - ], -) -def test_temperature_value_updates_with_HTransportProblem(T_function, expected_values): - """Test that different time dependent callable functions can be applied to - the temperature value, asserting in each case they match an expected value""" - subdomain = F.SurfaceSubdomain1D(1, x=4) - vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat) - - my_model = F.HydrogenTransportProblem( - mesh=test_mesh, - subdomains=[vol_subdomain, subdomain], - species=[F.Species("H")], - settings=F.Settings(atol=1, rtol=0.1, final_time=3), - ) - my_model.settings.stepsize = F.Stepsize(initial_value=1) - - my_model.temperature = T_function - - # RUN - my_model.initialise() - my_model.progress = tqdm.autonotebook.tqdm( - desc="Solving H transport problem", - total=my_model.settings.final_time, - unit_scale=True, - ) - # TODO get rid of these when post processing is implemented - my_model.flux_values, my_model.times = [], [] - - # TEST - for i in range(3): - my_model.iterate() - if isinstance(my_model.temperature_fenics, fem.Constant): - computed_value = float(my_model.temperature_fenics) - else: - computed_value = my_model.temperature_fenics.vector.array[-1] - assert np.isclose(computed_value, expected_values[i]) - - -def test_TypeError_raised_when_temperature_not_defined(): - """Test that a type error when a model is initialised without - defining a temperature""" - my_model = F.HydrogenTransportProblem( - mesh=test_mesh, - subdomains=[F.VolumeSubdomain1D(1, borders=[0, 4], material=dummy_mat)], - species=[F.Species("H")], - settings=F.Settings(atol=1, rtol=0.1, final_time=3), - ) - my_model.settings.stepsize = F.Stepsize(initial_value=1) - - with pytest.raises(TypeError, match="Temperature needs to be defined"): - my_model.initialise() From 0f63f92e1ba36ffdb9bd0f5ef00aabc7ae2c7164 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 15:37:11 -0400 Subject: [PATCH 21/27] use float(self.t) --- festim/hydrogen_transport_problem.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 535cd4478..ddfc7ae28 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -481,7 +481,7 @@ def iterate( self.progress.update(self.dt.value) self.t.value += self.dt.value - self.update_time_dependent_values(t=float(self.t)) + self.update_time_dependent_values() self.solver.solve(self.u) @@ -525,12 +525,12 @@ def iterate( # update previous solution self.u_n.x.array[:] = self.u.x.array[:] - def update_time_dependent_values(self, t): + def update_time_dependent_values(self): if self.temperature_time_dependent: if isinstance(self.temperature_fenics, fem.Constant): - self.temperature_fenics.value = self.temperature(t=t) + self.temperature_fenics.value = self.temperature(t=float(self.t)) else: self.temperature_fenics.interpolate(self.temperature_expr) for bc in self.boundary_conditions: - bc.update(t) + bc.update(float(self.t)) From 2a254b8b50305060d64aded812e385376bef4234 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 15:40:15 -0400 Subject: [PATCH 22/27] updated doc strings --- festim/hydrogen_transport_problem.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ddfc7ae28..3f7186edd 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -174,7 +174,8 @@ def define_temperature(self): self.temperature. If self.temperature is a function on t only, create a fem.Constant. Else, create an dolfinx.fem.Expression (stored in self.temperature_expr) to be updated, a dolfinx.fem.Function object - is created from the Expression (stored in self.temperature_fenics_value) + is created from the Expression (stored in self.temperature_fenics_value). + Raise a ValueError if temperature is None. """ # check if temperature is None if self.temperature is None: From 535e56fb54170710e72720c3e3c134a76db2dac7 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 17:24:09 -0400 Subject: [PATCH 23/27] remove logic from tests --- test/test_h_transport_problem.py | 93 ++++++++++++++++---------------- 1 file changed, 45 insertions(+), 48 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 7f6bda9e4..caa712fdd 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -31,76 +31,73 @@ def test_temperature_setter_type(value): @pytest.mark.parametrize( - "value", + "input, expected_value", [ - 1.0, - 1, - None, - lambda t: t, - lambda t: 1.0 + t, - lambda x: 1.0 + x[0], - lambda x, t: 1.0 + x[0] + t, - lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + (1.0, False), + (1, False), + (lambda t: t, True), + (lambda t: 1.0 + t, True), + (lambda x: 1.0 + x[0], False), + (lambda x, t: 1.0 + x[0] + t, True), + (lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), True), ], ) -def test_time_dependent_temperature_attribute(value): +def test_time_dependent_temperature_attribute(input, expected_value): """Test that the temperature_time_dependent attribute is correctly set""" my_model = F.HydrogenTransportProblem() - my_model.temperature = value + my_model.temperature = input - if callable(value): - arguments = value.__code__.co_varnames - if "t" in arguments: - assert my_model.temperature_time_dependent - else: - assert not my_model.temperature_time_dependent + assert my_model.temperature_time_dependent == expected_value + + +def test_define_temperature_value_error_raised(): + """Test that the define_temperature method correctly sets the + temperature_fenics attribute to either a fem.Constant or a + fem.Function and raise a ValueError temperature is None""" + + # BUILD + my_model = F.HydrogenTransportProblem(mesh=test_mesh) + my_model.t = fem.Constant(test_mesh.mesh, 0.0) + + my_model.temperature = None + + # TEST + with pytest.raises( + ValueError, match="the temperature attribute needs to be defined" + ): + my_model.define_temperature() @pytest.mark.parametrize( - "value", + "input, expected_type", [ - 1.0, - 1, - None, - fem.Constant(test_mesh.mesh, 1.0), - lambda t: t, - lambda t: 1.0 + t, - lambda x: 1.0 + x[0], - lambda x, t: 1.0 + x[0] + t, - lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + (1.0, fem.Constant), + (1, fem.Constant), + (fem.Constant(test_mesh.mesh, 1.0), fem.Constant), + (lambda t: t, fem.Constant), + (lambda t: 1.0 + t, fem.Constant), + (lambda x: 1.0 + x[0], fem.Function), + (lambda x, t: 1.0 + x[0] + t, fem.Function), + (lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), fem.Function), ], ) -def test_define_temperature(value): +def test_define_temperature(input, expected_type): """Test that the define_temperature method correctly sets the temperature_fenics attribute to either a fem.Constant or a - fem.Function and raise a ValueError temperature is None""" + fem.Function""" # BUILD my_model = F.HydrogenTransportProblem(mesh=test_mesh) my_model.t = fem.Constant(test_mesh.mesh, 0.0) - my_model.temperature = value + my_model.temperature = input - # TEST - if value is None: - with pytest.raises( - ValueError, match="the temperature attribute needs to be defined" - ): - my_model.define_temperature() - else: - # RUN - my_model.define_temperature() + # RUN + my_model.define_temperature() - # TEST - if isinstance(value, (fem.Constant, int, float)): - assert isinstance(my_model.temperature_fenics, fem.Constant) - elif callable(value): - arguments = value.__code__.co_varnames - if "x" in arguments: - assert isinstance(my_model.temperature_fenics, fem.Function) - else: - assert isinstance(my_model.temperature_fenics, fem.Constant) + # TEST + assert isinstance(my_model.temperature_fenics, expected_type) def test_iterate(): From ce9a411e4d32dcfbf627ba9dbb8d832c48d8d660 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 17:32:25 -0400 Subject: [PATCH 24/27] added raise value error, degree argument --- festim/hydrogen_transport_problem.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3f7186edd..d21c53874 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -137,6 +137,8 @@ def temperature_fenics(self, value): if value is None: self._temperature_fenics = value return + elif not isinstance(value, (fem.Constant, fem.Function)): + raise TypeError(f"Value must be a fem.Constant or fem.Function") self._temperature_fenics = value @property @@ -200,14 +202,15 @@ def define_temperature(self): ) else: x = ufl.SpatialCoordinate(self.mesh.mesh) - CG_temperature = basix.ufl.element( + degree = 1 + element_temperature = basix.ufl.element( basix.ElementFamily.P, self.mesh.mesh.basix_cell(), - 1, + degree, basix.LagrangeVariant.equispaced, ) function_space_temperature = fem.FunctionSpace( - self.mesh.mesh, CG_temperature + self.mesh.mesh, element_temperature ) self.temperature_fenics = fem.Function(function_space_temperature) kwargs = {} @@ -242,10 +245,11 @@ def define_function_space(self): """Creates the function space of the model, creates a mixed element if model is multispecies. Creates the main solution and previous solution function u and u_n.""" + degree = 1 element_CG = basix.ufl.element( basix.ElementFamily.P, self.mesh.mesh.basix_cell(), - 1, + degree, basix.LagrangeVariant.equispaced, ) if not self.multispecies: @@ -527,11 +531,12 @@ def iterate( self.u_n.x.array[:] = self.u.x.array[:] def update_time_dependent_values(self): + t = float(self.t) if self.temperature_time_dependent: if isinstance(self.temperature_fenics, fem.Constant): - self.temperature_fenics.value = self.temperature(t=float(self.t)) + self.temperature_fenics.value = self.temperature(t=t) else: self.temperature_fenics.interpolate(self.temperature_expr) for bc in self.boundary_conditions: - bc.update(float(self.t)) + bc.update(t) From 9f39d5197c98bba25be239bb9e3603002cbc40c9 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 17:40:52 -0400 Subject: [PATCH 25/27] update docs --- test/test_h_transport_problem.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index caa712fdd..85003b167 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -52,13 +52,10 @@ def test_time_dependent_temperature_attribute(input, expected_value): def test_define_temperature_value_error_raised(): - """Test that the define_temperature method correctly sets the - temperature_fenics attribute to either a fem.Constant or a - fem.Function and raise a ValueError temperature is None""" + """Test that a ValueError is rasied when the temperature is None""" # BUILD my_model = F.HydrogenTransportProblem(mesh=test_mesh) - my_model.t = fem.Constant(test_mesh.mesh, 0.0) my_model.temperature = None @@ -85,7 +82,7 @@ def test_define_temperature_value_error_raised(): def test_define_temperature(input, expected_type): """Test that the define_temperature method correctly sets the temperature_fenics attribute to either a fem.Constant or a - fem.Function""" + fem.Function depending on the type of input""" # BUILD my_model = F.HydrogenTransportProblem(mesh=test_mesh) From 363caf4b040a9775d5d9b6444c919078b5fd053e Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Wed, 1 Nov 2023 16:18:36 -0400 Subject: [PATCH 26/27] typo 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> --- test/test_h_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 85003b167..6045394f8 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -52,7 +52,7 @@ def test_time_dependent_temperature_attribute(input, expected_value): def test_define_temperature_value_error_raised(): - """Test that a ValueError is rasied when the temperature is None""" + """Test that a ValueError is raised when the temperature is None""" # BUILD my_model = F.HydrogenTransportProblem(mesh=test_mesh) From b92b8275a198a118defb48790e385934a82106b0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:26:20 -0400 Subject: [PATCH 27/27] is, not == in test --- test/test_h_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 6045394f8..c5a42d3f0 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -48,7 +48,7 @@ def test_time_dependent_temperature_attribute(input, expected_value): my_model = F.HydrogenTransportProblem() my_model.temperature = input - assert my_model.temperature_time_dependent == expected_value + assert my_model.temperature_time_dependent is expected_value def test_define_temperature_value_error_raised():