diff --git a/examples/multi_isotope_trapping_example.py b/examples/multi_isotope_trapping_example.py index a3528629e..9312e62c9 100644 --- a/examples/multi_isotope_trapping_example.py +++ b/examples/multi_isotope_trapping_example.py @@ -62,8 +62,7 @@ E_k=0.39, p_0=1e13, E_p=1.2, - reactant1=mobile_H, - reactant2=empty_trap, + reactant=[mobile_H, empty_trap], product=trapped_H1, ), F.Reaction( @@ -71,8 +70,7 @@ E_k=0.39, p_0=1e13, E_p=1.0, - reactant1=mobile_H, - reactant2=trapped_H1, + reactant=[mobile_H, trapped_H1], product=trapped_H2, ), F.Reaction( @@ -80,8 +78,7 @@ E_k=0.39, p_0=1e13, E_p=1.2, - reactant1=mobile_D, - reactant2=empty_trap, + reactant=[mobile_D, empty_trap], product=trapped_D1, ), F.Reaction( @@ -89,8 +86,7 @@ E_k=0.39, p_0=1e13, E_p=1.2, - reactant1=mobile_D, - reactant2=trapped_D1, + reactant=[mobile_D, trapped_D1], product=trapped_D2, ), F.Reaction( @@ -98,8 +94,7 @@ E_k=0.39, p_0=1e13, E_p=1.0, - reactant1=mobile_H, - reactant2=trapped_D1, + reactant=[mobile_H, trapped_D1], product=trapped_HD, ), F.Reaction( @@ -107,8 +102,7 @@ E_k=0.39, p_0=1e13, E_p=1.0, - reactant1=mobile_D, - reactant2=trapped_H1, + reactant=[mobile_D, trapped_H1], product=trapped_HD, ), ] diff --git a/examples/tds_example.py b/examples/tds_example.py index 0d76eed08..d0cf87013 100644 --- a/examples/tds_example.py +++ b/examples/tds_example.py @@ -52,8 +52,7 @@ E_k=0.39, p_0=1e13, E_p=0.87, - reactant1=mobile_H, - reactant2=empty_trap1, + reactant=[mobile_H, empty_trap1], product=trapped_H1, volume=my_subdomain, ), @@ -62,8 +61,7 @@ E_k=0.39, p_0=1e13, E_p=1.0, - reactant1=mobile_H, - reactant2=empty_trap2, + reactant=[mobile_H, empty_trap2], product=trapped_H2, volume=my_subdomain, ), diff --git a/examples/trapping_example.py b/examples/trapping_example.py index 8d36fc6d3..2ed82a7d7 100644 --- a/examples/trapping_example.py +++ b/examples/trapping_example.py @@ -30,8 +30,7 @@ E_p=1.2, k_0=3.8e-17, E_k=0.2, - reactant1=mobile_H, - reactant2=empty_trap, + reactant=[mobile_H, empty_trap], product=trapped_H, ) ] diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index af0efd000..5603261f6 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -655,20 +655,14 @@ def create_formulation(self): self.formulation += ((u - u_n) / self.dt) * v * self.dx(vol.id) for reaction in self.reactions: - # reactant 1 - if isinstance(reaction.reactant1, F.Species): - self.formulation += ( - reaction.reaction_term(self.temperature_fenics) - * reaction.reactant1.test_function - * self.dx(reaction.volume.id) - ) - # reactant 2 - if isinstance(reaction.reactant2, F.Species): - self.formulation += ( - reaction.reaction_term(self.temperature_fenics) - * reaction.reactant2.test_function - * self.dx(reaction.volume.id) - ) + for reactant in reaction.reactant: + if isinstance(reactant, F.Species): + self.formulation += ( + reaction.reaction_term(self.temperature_fenics) + * reactant.test_function + * self.dx(reaction.volume.id) + ) + # product if isinstance(reaction.product, list): products = reaction.product diff --git a/festim/reaction.py b/festim/reaction.py index 73abb1102..08b25ce8e 100644 --- a/festim/reaction.py +++ b/festim/reaction.py @@ -1,5 +1,5 @@ import festim as F -from typing import Union +from typing import Union, Optional, List from ufl import exp @@ -8,9 +8,8 @@ class Reaction: """A reaction between two species, with a forward and backward rate. Arguments: - reactant1 (Union[F.Species, F.ImplicitSpecies]): The first reactant. - reactant2 (Union[F.Species, F.ImplicitSpecies]): The second reactant. - product (F.Species): The product. + reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant. + product (Optional[Union[F.Species, List[F.Species]]]): The product. k_0 (float): The forward rate constant pre-exponential factor. E_k (float): The forward rate constant activation energy. p_0 (float): The backward rate constant pre-exponential factor. @@ -18,9 +17,8 @@ class Reaction: volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place. Attributes: - reactant1 (Union[F.Species, F.ImplicitSpecies]): The first reactant. - reactant2 (Union[F.Species, F.ImplicitSpecies]): The second reactant. - product (F.Species): The product. + reactant (Union[F.Species, F.ImplicitSpecies], List[Union[F.Species, F.ImplicitSpecies]]): The reactant. + product (Optional[Union[F.Species, List[F.Species]]]): The product. k_0 (float): The forward rate constant pre-exponential factor. E_k (float): The forward rate constant activation energy. p_0 (float): The backward rate constant pre-exponential factor. @@ -29,14 +27,13 @@ class Reaction: Usage: >>> # create two species - >>> species1 = F.Species("A") - >>> species2 = F.Species("B") + >>> reactant = [F.Species("A"), F.Species("B")] >>> # create a product species >>> product = F.Species("C") >>> # create a reaction between the two species - >>> reaction = Reaction(species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3) + >>> reaction = Reaction(reactant, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3) >>> print(reaction) A + B <--> C @@ -48,61 +45,59 @@ class Reaction: def __init__( self, - reactant1: Union[F.Species, F.ImplicitSpecies], - reactant2: Union[F.Species, F.ImplicitSpecies], - product: F.Species, + reactant: Union[ + F.Species, F.ImplicitSpecies, List[Union[F.Species, F.ImplicitSpecies]] + ], k_0: float, E_k: float, - p_0: float, - E_p: float, volume: F.VolumeSubdomain1D, + product: Optional[Union[F.Species, List[F.Species]]] = [], + p_0: float = None, + E_p: float = None, ) -> None: - self.reactant1 = reactant1 - self.reactant2 = reactant2 - self.product = product self.k_0 = k_0 self.E_k = E_k self.p_0 = p_0 self.E_p = E_p + self.reactant = reactant + self.product = product self.volume = volume @property - def reactant1(self): - return self._reactant1 - - @reactant1.setter - def reactant1(self, value): - if not isinstance(value, (F.Species, F.ImplicitSpecies)): - raise TypeError( - f"reactant1 must be an F.Species or F.ImplicitSpecies, not {type(value)}" + def reactant(self): + return self._reactant + + @reactant.setter + def reactant(self, value): + if not isinstance(value, list): + value = [value] + if len(value) == 0: + raise ValueError( + f"reactant must be an entry of one or more species objects, not an empty list." ) - self._reactant1 = value - - @property - def reactant2(self): - return self._reactant2 - - @reactant2.setter - def reactant2(self, value): - if not isinstance(value, (F.Species, F.ImplicitSpecies)): - raise TypeError( - f"reactant2 must be an F.Species or F.ImplicitSpecies, not {type(value)}" - ) - self._reactant2 = value + for i in value: + if not isinstance(i, (F.Species, F.ImplicitSpecies)): + raise TypeError( + f"reactant must be an F.Species or F.ImplicitSpecies, not {type(i)}" + ) + self._reactant = value def __repr__(self) -> str: + reactants = " + ".join([str(reactant) for reactant in self.reactant]) + if isinstance(self.product, list): products = " + ".join([str(product) for product in self.product]) else: products = self.product - return f"Reaction({self.reactant1} + {self.reactant2} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" + return f"Reaction({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" def __str__(self) -> str: + reactants = " + ".join([str(reactant) for reactant in self.reactant]) if isinstance(self.product, list): products = " + ".join([str(product) for product in self.product]) else: products = self.product - return f"{self.reactant1} + {self.reactant2} <--> {products}" + return f"{reactants} <--> {products}" def reaction_term(self, temperature): """Compute the reaction term at a given temperature. @@ -110,19 +105,49 @@ def reaction_term(self, temperature): Arguments: temperature (): The temperature at which the reaction term is computed. """ + + if self.product == []: + if self.p_0 is not None: + raise ValueError( + f"p_0 must be None, not {self.p_0} when no products are present." + ) + if self.E_p is not None: + raise ValueError( + f"E_p must be None, not {self.E_p} when no products are present." + ) + else: + if self.p_0 == None: + raise ValueError( + f"p_0 cannot be None when reaction products are present." + ) + elif self.E_p == None: + raise ValueError( + f"E_p cannot be None when reaction products are present." + ) + k = self.k_0 * exp(-self.E_k / (F.k_B * temperature)) - p = self.p_0 * exp(-self.E_p / (F.k_B * temperature)) - c_A = self.reactant1.concentration - c_B = self.reactant2.concentration + if self.p_0 and self.E_p: + p = self.p_0 * exp(-self.E_p / (F.k_B * temperature)) + elif self.p_0: + p = self.p_0 + else: + p = 0 + + reactants = self.reactant + product_of_reactants = reactants[0].concentration + for reactant in reactants[1:]: + product_of_reactants *= reactant.concentration if isinstance(self.product, list): products = self.product else: products = [self.product] - products_of_product = products[0].solution - for product in products[1:]: - products_of_product *= product.solution - - return k * c_A * c_B - p * products_of_product + if len(products) > 0: + product_of_products = products[0].solution + for product in products[1:]: + product_of_products *= product.solution + else: + product_of_products = 0 + return k * product_of_reactants - (p * product_of_products) diff --git a/festim/species.py b/festim/species.py index bc06467b5..04ac3f25f 100644 --- a/festim/species.py +++ b/festim/species.py @@ -117,8 +117,7 @@ class Trap(Species): ct = F.Species("trapped") trap_sites = F.ImplicitSpecies(n=1, others=[ct]) trap_reaction = F.Reaction( - reactant1=cm, - reactant2=trap_sites, + reactant=[cm, trap_sites], product=ct, k_0=1, E_k=1, @@ -153,8 +152,7 @@ def create_species_and_reaction(self): trap_site = F.ImplicitSpecies(n=self.n, others=[self.trapped_concentration]) self.reaction = F.Reaction( - reactant1=self.mobile_species, - reactant2=trap_site, + reactant=[self.mobile_species, trap_site], product=self.trapped_concentration, k_0=self.k_0, E_k=self.E_k, diff --git a/test/system_tests/test_1_mobile_1_trap.py b/test/system_tests/test_1_mobile_1_trap.py index 58aeb073c..b08f51f7b 100644 --- a/test/system_tests/test_1_mobile_1_trap.py +++ b/test/system_tests/test_1_mobile_1_trap.py @@ -76,8 +76,7 @@ def v_exact(mod): my_model.reactions = [ F.Reaction( - reactant1=mobile, - reactant2=traps, + reactant=[mobile, traps], product=trapped, k_0=k_0, E_k=E_k, @@ -273,8 +272,7 @@ def v_exact(mod): my_model.reactions = [ F.Reaction( - reactant1=mobile, - reactant2=traps, + reactant=[mobile, traps], product=trapped, k_0=k_0, E_k=E_k, @@ -403,8 +401,7 @@ def v_exact(mod): my_model.reactions = [ F.Reaction( - reactant1=mobile, - reactant2=traps, + reactant=[mobile, traps], product=trapped, k_0=k_0, E_k=E_k, diff --git a/test/system_tests/test_reaction_not_in_every_volume.py b/test/system_tests/test_reaction_not_in_every_volume.py index 5b7a5c3e6..ea315458a 100644 --- a/test/system_tests/test_reaction_not_in_every_volume.py +++ b/test/system_tests/test_reaction_not_in_every_volume.py @@ -19,8 +19,7 @@ def test_sim_reaction_not_in_every_volume(): trap = F.ImplicitSpecies(n=8.19e25, others=[ct]) my_model.reactions = [ F.Reaction( - reactant1=cm, - reactant2=trap, + reactant=[cm, trap], product=ct, k_0=8.9e-17, E_k=0.39, diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py index 6013955f9..dfaaab6a7 100644 --- a/test/test_complex_reaction.py +++ b/test/test_complex_reaction.py @@ -93,8 +93,7 @@ def model_test_reaction(stepsize=1, k=350e-4, p=120e-4, c_A_0=1): E_k=0, p_0=p, E_p=0, - reactant1=species_A, - reactant2=species_B, + reactant=[species_A, species_B], product=[species_C, species_D], volume=my_subdomain, ), diff --git a/test/test_reaction.py b/test/test_reaction.py index 741790efc..e1dc47617 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -19,12 +19,17 @@ def test_reaction_init(): # create a reaction between the two species reaction = F.Reaction( - species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, volume=my_vol + reactant=[species1, species2], + product=product, + k_0=1.0, + E_k=0.2, + p_0=0.1, + E_p=0.3, + volume=my_vol, ) # check that the attributes are set correctly - assert reaction.reactant1 == species1 - assert reaction.reactant2 == species2 + assert reaction.reactant == [species1, species2] assert reaction.product == product assert reaction.k_0 == 1.0 assert reaction.E_k == 0.2 @@ -44,7 +49,13 @@ def test_reaction_repr(): # create a reaction between the two species reaction = F.Reaction( - species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, volume=my_vol + reactant=[species1, species2], + product=product, + k_0=1.0, + E_k=0.2, + p_0=0.1, + E_p=0.3, + volume=my_vol, ) # check that the __repr__ method returns the expected string @@ -65,9 +76,8 @@ def test_reaction_repr_2_products(): # create a reaction between the two species reaction = F.Reaction( - species1, - species2, - [product1, product2], + reactant=[species1, species2], + product=[product1, product2], k_0=1.0, E_k=0.2, p_0=0.1, @@ -80,6 +90,25 @@ def test_reaction_repr_2_products(): assert repr(reaction) == expected_repr +def test_reaction_repr_0_products(): + """Test that the Reaction __repr__ method returns the expected string""" + + # create two species + species1 = F.Species("A") + + # create a reaction between the two species + reaction = F.Reaction( + reactant=species1, + k_0=1.0, + E_k=0.2, + volume=my_vol, + ) + + # check that the __repr__ method returns the expected string + expected_repr = "Reaction(A <--> , 1.0, 0.2, None, None)" + assert repr(reaction) == expected_repr + + def test_reaction_str(): """Test that the Reaction __str__ method returns the expected string""" @@ -92,7 +121,13 @@ def test_reaction_str(): # create a reaction between the two species reaction = F.Reaction( - species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, volume=my_vol + reactant=[species1, species2], + product=product, + k_0=1.0, + E_k=0.2, + p_0=0.1, + E_p=0.3, + volume=my_vol, ) # check that the __str__ method returns the expected string @@ -113,9 +148,8 @@ def test_reaction_str_2_products(): # create a reaction between the two species reaction = F.Reaction( - species1, - species2, - [product1, product2], + reactant=[species1, species2], + product=[product1, product2], k_0=1.0, E_k=0.2, p_0=0.1, @@ -128,6 +162,25 @@ def test_reaction_str_2_products(): assert str(reaction) == expected_str +def test_reaction_str_no_products(): + """Test that the Reaction __str__ method returns the expected string when there are 2 products""" + + # create two species + species1 = F.Species("A") + + # create a reaction between the two species + reaction = F.Reaction( + reactant=species1, + k_0=1.0, + E_k=0.2, + volume=my_vol, + ) + + # check that the __str__ method returns the expected string + expected_str = "A <--> " + assert str(reaction) == expected_str + + @pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0]) def test_reaction_reaction_term(temperature): """Test that the Reaction.reaction_term method returns the expected reaction term""" @@ -147,7 +200,13 @@ def test_reaction_reaction_term(temperature): # create a reaction between the two species reaction = F.Reaction( - species1, species2, product, k_0=1.0, E_k=0.2, p_0=0.1, E_p=0.3, volume=my_vol + reactant=[species1, species2], + product=product, + k_0=1.0, + E_k=0.2, + p_0=0.1, + E_p=0.3, + volume=my_vol, ) # test the reaction term at a given temperature @@ -158,9 +217,39 @@ def arrhenius(pre, act, T): p = arrhenius(reaction.p_0, reaction.E_p, temperature) expected_reaction_term = ( - k * species1.solution * species2.solution - p * product.solution + k * (species1.solution * species2.solution) - p * product.solution + ) + + assert reaction.reaction_term(temperature) == expected_reaction_term + + +@pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0]) +def test_reaction_reaction_term_no_products(temperature): + mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10) + V = functionspace(mesh, ("Lagrange", 1)) + + # create two species + species1 = F.Species("A") + species2 = F.Species("B") + species1.solution = Function(V) + species2.solution = Function(V) + + # create a reaction between the two species + reaction = F.Reaction( + reactant=[species1, species2], + k_0=1.0, + E_k=0.2, + volume=my_vol, ) + # test the reaction term at a given temperature + def arrhenius(pre, act, T): + return pre * exp(-act / (F.k_B * T)) + + k = arrhenius(reaction.k_0, reaction.E_k, temperature) + + expected_reaction_term = k * (species1.solution * species2.solution) + assert reaction.reaction_term(temperature) == expected_reaction_term @@ -185,9 +274,8 @@ def test_reaction_reaction_term_2_products(temperature): # create a reaction between the two species reaction = F.Reaction( - species1, - species2, - [product1, product2], + reactant=[species1, species2], + product=[product1, product2], k_0=1.0, E_k=0.2, p_0=0.1, @@ -204,21 +292,19 @@ def arrhenius(pre, act, T): product_of_products = product1.solution * product2.solution expected_reaction_term = ( - k * species1.solution * species2.solution - p * product_of_products + k * (species1.solution * species2.solution) - p * product_of_products ) assert reaction.reaction_term(temperature) == expected_reaction_term -def test_reactant1_setter_raises_error_with_wrong_type(): - """Test a type error is raised when the reactant1 is given a wrong type.""" +def test_reactant_setter_raises_error_with_zero_length_list(): + """Test a value error is raised when the first reactant is given a wrong type.""" with pytest.raises( - TypeError, - match="reactant1 must be an F.Species or F.ImplicitSpecies, not ", + ValueError, + match="reactant must be an entry of one or more species objects, not an empty list.", ): F.Reaction( - reactant1="A", - reactant2=F.Species("B"), - product=F.Species("C"), + reactant=[], k_0=1, E_k=0.1, p_0=2, @@ -227,15 +313,14 @@ def test_reactant1_setter_raises_error_with_wrong_type(): ) -def test_reactant2_setter_raises_error_with_wrong_type(): - """Test a type error is raised when the reactant2 is given a wrong type.""" +def test_reactant_setter_raises_error_with_wrong_type(): + """Test a type error is raised when the first reactant is given a wrong type.""" with pytest.raises( TypeError, - match="reactant2 must be an F.Species or F.ImplicitSpecies, not ", + match="reactant must be an F.Species or F.ImplicitSpecies, not ", ): F.Reaction( - reactant1=F.Species("A"), - reactant2="B", + reactant=["A", F.Species("B")], product=F.Species("C"), k_0=1, E_k=0.1, @@ -243,3 +328,65 @@ def test_reactant2_setter_raises_error_with_wrong_type(): E_p=0.2, volume=my_vol, ) + + +def test_product_setter_raise_error_p_0_no_product(): + with pytest.raises( + ValueError, + match="p_0 must be None, not 2 when no products are present.", + ): + reaction = F.Reaction( + reactant=[F.Species("A")], + k_0=1, + E_k=0.1, + p_0=2, + volume=my_vol, + ) + reaction.reaction_term(temperature=500) + + +def test_no_E_p_with_product(): + with pytest.raises( + ValueError, + match="E_p cannot be None when reaction products are present.", + ): + reaction = F.Reaction( + reactant=[F.Species("A")], + product=[F.Species("C")], + k_0=1, + E_k=0.1, + p_0=0.1, + volume=my_vol, + ) + reaction.reaction_term(temperature=500) + + +def test_no_p_0_with_product(): + with pytest.raises( + ValueError, + match="p_0 cannot be None when reaction products are present.", + ): + reaction = F.Reaction( + reactant=[F.Species("A")], + product=[F.Species("C")], + k_0=1, + E_k=0.1, + E_p=1, + volume=my_vol, + ) + reaction.reaction_term(temperature=500) + + +def test_product_setter_raise_error_E_p_no_product(): + with pytest.raises( + ValueError, + match="E_p must be None, not 2 when no products are present.", + ): + reaction = F.Reaction( + reactant=[F.Species("A")], + k_0=1, + E_k=0.1, + E_p=2, + volume=my_vol, + ) + reaction.reaction_term(temperature=500)