Arbitrary reaction

examples/multi_isotope_trapping_example.py

@@ -62,53 +62,47 @@
         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(
         k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density),
         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(
         k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density),
         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(
         k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density),
         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(
         k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density),
         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(
         k_0=4.1e-7 / (1.1e-10**2 * 6 * w_atom_density),
         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,
     ),
 ]

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,
     ),

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,
     )
 ]

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

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,19 +8,17 @@ 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.
         E_p (float): The backward rate constant activation energy.
         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,81 +45,109 @@ 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.
 
         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)

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,

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,

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,

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,
         ),