From a93482f5dc9a5603214a1bfeda94145452b6e039 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 11:53:50 -0400
Subject: [PATCH 01/10] initial Reaction class

---
 festim/__init__.py |  2 ++
 festim/reaction.py | 77 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 79 insertions(+)
 create mode 100644 festim/reaction.py

diff --git a/festim/__init__.py b/festim/__init__.py
index 2af946b92..e34146ecb 100644
--- a/festim/__init__.py
+++ b/festim/__init__.py
@@ -35,3 +35,5 @@
 
 from .exports.vtx import VTXExport
 from .exports.xdmf import XDMFExport
+
+from .reaction import Reaction
diff --git a/festim/reaction.py b/festim/reaction.py
new file mode 100644
index 000000000..aa2211145
--- /dev/null
+++ b/festim/reaction.py
@@ -0,0 +1,77 @@
+import festim as F
+from typing import Union
+
+from ufl import exp
+
+
+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.
+        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.
+
+    Attributes:
+        reactant1 (Union[F.Species, F.ImplicitSpecies]): The first reactant.
+        reactant2 (Union[F.Species, F.ImplicitSpecies]): The second reactant.
+        product (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.
+
+    Usage:
+        >>> # create two species
+        >>> species1 = F.Species("A")
+        >>> species2 = 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)
+        >>> print(reaction)
+        A + B <--> C
+
+        >>> # compute the reaction term at a given temperature
+        >>> temperature = 300.0
+        >>> reaction_term = reaction.reaction_term(temperature)
+
+    """
+
+    def __init__(
+        self,
+        reactant1: Union[F.Species, F.ImplicitSpecies],
+        reactant2: Union[F.Species, F.ImplicitSpecies],
+        product: F.Species,
+        k_0: float,
+        E_k: float,
+        p_0: float,
+        E_p: float,
+    ) -> 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
+
+    def __repr__(self) -> str:
+        return f"Reaction({self.reactant1} + {self.reactant2} <--> {self.product}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})"
+
+    def __str__(self) -> str:
+        return f"{self.reactant1} + {self.reactant2} <--> {self.product}"
+
+    def reaction_term(self, temperature):
+        k = self.k_0 * exp(-self.E_k / (F.k_B * temperature))
+        p = self.p_0 * exp(-self.E_p / (F.k_B * temperature))
+        A = self.reactant1.solution
+        B = self.reactant2.solution
+        C = self.product.solution
+        return k * A * B - p * C

From 9c79abb0263751b89e7f2600b2bf405291b3293c Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:33:38 -0400
Subject: [PATCH 02/10] added ImplicitSpecies + more functionalities

---
 festim/__init__.py |  2 +-
 festim/reaction.py | 15 +++++++++++----
 festim/species.py  | 48 ++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 60 insertions(+), 5 deletions(-)

diff --git a/festim/__init__.py b/festim/__init__.py
index e34146ecb..853092014 100644
--- a/festim/__init__.py
+++ b/festim/__init__.py
@@ -28,7 +28,7 @@
 
 from .hydrogen_transport_problem import HydrogenTransportProblem
 
-from .species import Species, Trap
+from .species import Species, Trap, ImplicitSpecies
 
 from .subdomain.surface_subdomain import SurfaceSubdomain1D
 from .subdomain.volume_subdomain import VolumeSubdomain1D
diff --git a/festim/reaction.py b/festim/reaction.py
index aa2211145..104cebbb7 100644
--- a/festim/reaction.py
+++ b/festim/reaction.py
@@ -69,9 +69,16 @@ def __str__(self) -> str:
         return f"{self.reactant1} + {self.reactant2} <--> {self.product}"
 
     def reaction_term(self, temperature):
+        """Compute the reaction term at a given temperature.
+
+        Arguments:
+            temperature (): The temperature at which the reaction term is computed.
+        """
         k = self.k_0 * exp(-self.E_k / (F.k_B * temperature))
         p = self.p_0 * exp(-self.E_p / (F.k_B * temperature))
-        A = self.reactant1.solution
-        B = self.reactant2.solution
-        C = self.product.solution
-        return k * A * B - p * C
+
+        c_A = self.reactant1.concentration
+        c_B = self.reactant2.concentration
+
+        c_C = self.product.solution
+        return k * c_A * c_B - p * c_C
diff --git a/festim/species.py b/festim/species.py
index eda76f2e1..8b23a6cdd 100644
--- a/festim/species.py
+++ b/festim/species.py
@@ -1,3 +1,6 @@
+from typing import List, Union
+
+
 class Species:
     """
     Hydrogen species class for H transport simulation.
@@ -10,6 +13,7 @@ class Species:
         solution (dolfinx.fem.Function or ...): the solution for the current timestep
         prev_solution (dolfinx.fem.Function or ...): the solution for the previous timestep
         test_function (ufl.Argument or ...): the testfunction associated with this species
+        concentration (dolfinx.fem.Function): the concentration of the species
 
     Usage:
         >>> from festim import Species, HTransportProblem
@@ -28,6 +32,10 @@ def __init__(self, name: str = None) -> None:
         self.prev_solution = None
         self.test_function = None
 
+    @property
+    def concentration(self):
+        return self.solution
+
 
 class Trap(Species):
     """Trap species class for H transport simulation.
@@ -51,3 +59,43 @@ class Trap(Species):
 
     def __init__(self, name: str = None) -> None:
         super().__init__(name)
+
+
+class ImplicitSpecies:
+    """Implicit species class for H transport simulation.
+    c = n - others
+
+    Args:
+        n (float, function): the total concentration of the species
+        others (List[Species]): the list of species from which the implicit
+            species concentration is computed (c = n - others)
+        name (str, optional): a name given to the species. Defaults to None.
+
+    Attributes:
+        name (str): a name given to the species.
+        n (float): the total concentration of the species
+        others (List[Species]): the list of species from which the implicit
+            species concentration is computed (c = n - others)
+        concentration (form): the concentration of the species
+
+    """
+
+    def __init__(
+        self,
+        n: Union[float, function],
+        others: List[Species] = None,
+        name: str = None,
+    ) -> None:
+        self.name = name
+        self.n = n
+        self.others = others
+
+    @property
+    def concentration(self):
+        if len(self.others) > 0:
+            for other in self.others:
+                if other.solution is None:
+                    raise ValueError(
+                        f"Cannot compute concentration of {self.name} because {other.name} has no solution"
+                    )
+        return self.n - sum([other.solution for other in self.others])

From bed2fba177ae85e108983adbfaf694c7072c1aa7 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:38:38 -0400
Subject: [PATCH 03/10] added __repr__ and __str__ methods + right type

---
 festim/species.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/festim/species.py b/festim/species.py
index 8b23a6cdd..2f5d36287 100644
--- a/festim/species.py
+++ b/festim/species.py
@@ -1,4 +1,4 @@
-from typing import List, Union
+from typing import List, Union, Callable
 
 
 class Species:
@@ -32,6 +32,12 @@ def __init__(self, name: str = None) -> None:
         self.prev_solution = None
         self.test_function = None
 
+    def __repr__(self) -> str:
+        return f"Species({self.name})"
+
+    def __str__(self) -> str:
+        return f"{self.name}"
+
     @property
     def concentration(self):
         return self.solution
@@ -66,7 +72,7 @@ class ImplicitSpecies:
     c = n - others
 
     Args:
-        n (float, function): the total concentration of the species
+        n (float, Callable): the total concentration of the species
         others (List[Species]): the list of species from which the implicit
             species concentration is computed (c = n - others)
         name (str, optional): a name given to the species. Defaults to None.
@@ -82,7 +88,7 @@ class ImplicitSpecies:
 
     def __init__(
         self,
-        n: Union[float, function],
+        n: Union[float, Callable],
         others: List[Species] = None,
         name: str = None,
     ) -> None:
@@ -90,6 +96,12 @@ def __init__(
         self.n = n
         self.others = others
 
+    def __repr__(self) -> str:
+        return f"ImplicitSpecies({self.name}, {self.n}, {self.others})"
+
+    def __str__(self) -> str:
+        return f"{self.name}"
+
     @property
     def concentration(self):
         if len(self.others) > 0:

From 1ee7c8c855a7466632fb718465925709254e28fe Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:41:49 -0400
Subject: [PATCH 04/10] added tests for reaction

---
 test/test_reaction.py | 96 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 96 insertions(+)
 create mode 100644 test/test_reaction.py

diff --git a/test/test_reaction.py b/test/test_reaction.py
new file mode 100644
index 000000000..efc850943
--- /dev/null
+++ b/test/test_reaction.py
@@ -0,0 +1,96 @@
+import pytest
+import festim as F
+import numpy as np
+
+
+def test_reaction_init():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # 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
+    )
+
+    # check that the attributes are set correctly
+    assert reaction.reactant1 == species1
+    assert reaction.reactant2 == species2
+    assert reaction.product == product
+    assert reaction.k_0 == 1.0
+    assert reaction.E_k == 0.2
+    assert reaction.p_0 == 0.1
+    assert reaction.E_p == 0.3
+
+
+def test_reaction_repr():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # 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
+    )
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = "Reaction(A + B <--> C, 1.0, 0.2, 0.1, 0.3)"
+    assert repr(reaction) == expected_repr
+
+
+def test_reaction_str():
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # 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
+    )
+
+    # check that the __str__ method returns the expected string
+    expected_str = "A + B <--> C"
+    assert str(reaction) == expected_str
+
+
+@pytest.mark.parametrize("temperature", [300.0, 350, 370, 500.0])
+def test_reaction_reaction_term(temperature):
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create a product species
+    product = F.Species("C")
+
+    # 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
+    )
+
+    # set the concentrations of the species
+    species1.solution = 1.0
+    species2.solution = 2.0
+    product.solution = 3.0
+
+    # test the reaction term at a given temperature
+    def arrhenius(pre, act, T):
+        return pre * np.exp(-act / (F.k_B * T))
+
+    k = arrhenius(reaction.k_0, reaction.E_k, temperature)
+    p = arrhenius(reaction.p_0, reaction.E_p, temperature)
+    expected_reaction_term = (
+        k * species1.solution * species2.solution - p * product.solution
+    )
+    assert (
+        pytest.approx(reaction.reaction_term(temperature), rel=1e-9)
+        == expected_reaction_term
+    )

From 0d25d9e68dbe3038ecef54b16edf3ff97d77540e Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:45:39 -0400
Subject: [PATCH 05/10] better test

---
 test/test_reaction.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/test/test_reaction.py b/test/test_reaction.py
index efc850943..3d93493ad 100644
--- a/test/test_reaction.py
+++ b/test/test_reaction.py
@@ -1,6 +1,9 @@
 import pytest
 import festim as F
 import numpy as np
+from dolfinx.fem import FunctionSpace, Function
+from dolfinx.mesh import create_unit_cube
+from mpi4py import MPI
 
 
 def test_reaction_init():
@@ -77,9 +80,11 @@ def test_reaction_reaction_term(temperature):
     )
 
     # set the concentrations of the species
-    species1.solution = 1.0
-    species2.solution = 2.0
-    product.solution = 3.0
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+    species2.solution = Function(V)
+    product.solution = Function(V)
 
     # test the reaction term at a given temperature
     def arrhenius(pre, act, T):
@@ -90,7 +95,4 @@ def arrhenius(pre, act, T):
     expected_reaction_term = (
         k * species1.solution * species2.solution - p * product.solution
     )
-    assert (
-        pytest.approx(reaction.reaction_term(temperature), rel=1e-9)
-        == expected_reaction_term
-    )
+    assert reaction.reaction_term(temperature) == expected_reaction_term

From bbd764ce779c1593739ea09d5ada8dc333a6d405 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:55:50 -0400
Subject: [PATCH 06/10] added tests for ImplicitSpecies

---
 test/test_species.py | 50 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/test/test_species.py b/test/test_species.py
index d4f4f3b6b..89b0049d1 100644
--- a/test/test_species.py
+++ b/test/test_species.py
@@ -2,6 +2,10 @@
 import dolfinx
 import ufl
 import numpy as np
+import pytest
+from dolfinx.fem import FunctionSpace, Function
+from dolfinx.mesh import create_unit_cube
+from mpi4py import MPI
 
 
 def test_assign_functions_to_species():
@@ -29,3 +33,49 @@ def test_assign_functions_to_species():
         assert isinstance(spe.solution, dolfinx.fem.Function)
         assert isinstance(spe.prev_solution, dolfinx.fem.Function)
         assert isinstance(spe.test_function, ufl.Argument)
+
+
+def test_implicit_species_concentration():
+    """Test that the concentration of an implicit species is computed
+    correctly.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # set the solutions of the two species
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+    species2.solution = Function(V)
+
+    # test the concentration of the implicit species
+    expected_concentration = implicit_species.n - (
+        species1.solution + species2.solution
+    )
+    assert implicit_species.concentration == expected_concentration
+
+
+def test_implicit_species_concentration_with_no_solution():
+    """Test that a ValueError is raised when on of the 'others' species
+    has no solution and the concentration of the implicit species is
+    requested.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # set the solution of the first species
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    species1.solution = Function(V)
+
+    # test that a ValueError is raised when the second species has no solution
+    with pytest.raises(ValueError):
+        implicit_species.concentration

From 8325bf68353279ec26c00d6dc66824bf8110846a Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 12:55:59 -0400
Subject: [PATCH 07/10] removed callable

---
 festim/species.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/festim/species.py b/festim/species.py
index 2f5d36287..a3f45eacd 100644
--- a/festim/species.py
+++ b/festim/species.py
@@ -1,4 +1,4 @@
-from typing import List, Union, Callable
+from typing import List
 
 
 class Species:
@@ -72,7 +72,7 @@ class ImplicitSpecies:
     c = n - others
 
     Args:
-        n (float, Callable): the total concentration of the species
+        n (float): the total concentration of the species
         others (List[Species]): the list of species from which the implicit
             species concentration is computed (c = n - others)
         name (str, optional): a name given to the species. Defaults to None.
@@ -88,7 +88,7 @@ class ImplicitSpecies:
 
     def __init__(
         self,
-        n: Union[float, Callable],
+        n: float,
         others: List[Species] = None,
         name: str = None,
     ) -> None:

From 6c7496916cc31e37c608887278cc11f2216094ad Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 13:04:14 -0400
Subject: [PATCH 08/10] print statements

---
 test/test_reaction.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/test_reaction.py b/test/test_reaction.py
index 3d93493ad..f05f9584b 100644
--- a/test/test_reaction.py
+++ b/test/test_reaction.py
@@ -95,4 +95,9 @@ def arrhenius(pre, act, T):
     expected_reaction_term = (
         k * species1.solution * species2.solution - p * product.solution
     )
+    print(expected_reaction_term)
+    print(reaction.reaction_term(temperature))
+    from dolfinx import __version__
+
+    print(__version__)
     assert reaction.reaction_term(temperature) == expected_reaction_term

From fd3b0efa7052b1ebe52b622dbd4055918c956fc0 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 13:05:21 -0400
Subject: [PATCH 09/10] ufl exp instead of np

---
 test/test_reaction.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/test/test_reaction.py b/test/test_reaction.py
index f05f9584b..58fcfbcb7 100644
--- a/test/test_reaction.py
+++ b/test/test_reaction.py
@@ -1,9 +1,9 @@
 import pytest
 import festim as F
-import numpy as np
 from dolfinx.fem import FunctionSpace, Function
 from dolfinx.mesh import create_unit_cube
 from mpi4py import MPI
+from ufl import exp
 
 
 def test_reaction_init():
@@ -88,16 +88,12 @@ def test_reaction_reaction_term(temperature):
 
     # test the reaction term at a given temperature
     def arrhenius(pre, act, T):
-        return pre * np.exp(-act / (F.k_B * T))
+        return pre * exp(-act / (F.k_B * T))
 
     k = arrhenius(reaction.k_0, reaction.E_k, temperature)
     p = arrhenius(reaction.p_0, reaction.E_p, temperature)
     expected_reaction_term = (
         k * species1.solution * species2.solution - p * product.solution
     )
-    print(expected_reaction_term)
-    print(reaction.reaction_term(temperature))
-    from dolfinx import __version__
 
-    print(__version__)
     assert reaction.reaction_term(temperature) == expected_reaction_term

From ec2d15ea9bd7764a0aef22a979ed9b54cf39f89c Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Thu, 19 Oct 2023 13:16:50 -0400
Subject: [PATCH 10/10] increased coverage

---
 test/test_species.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/test/test_species.py b/test/test_species.py
index 89b0049d1..2eafd9844 100644
--- a/test/test_species.py
+++ b/test/test_species.py
@@ -35,6 +35,42 @@ def test_assign_functions_to_species():
         assert isinstance(spe.test_function, ufl.Argument)
 
 
+def test_species_repr_and_str():
+    """Test that the __repr__ and __str__ methods of the Species class returns the
+    expected string.
+    """
+    # create a species
+    species = F.Species("A")
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = "Species(A)"
+    assert repr(species) == expected_repr
+
+    # check that the __str__ method returns the expected string
+    expected_str = "A"
+    assert str(species) == expected_str
+
+
+def test_implicit_species_repr_and_str():
+    """Test that the __repr__ and __str__ methods of the ImplicitSpecies class
+    returns the expected string.
+    """
+    # create two species
+    species1 = F.Species("A")
+    species2 = F.Species("B")
+
+    # create an implicit species that depends on the two species
+    implicit_species = F.ImplicitSpecies(3.0, [species1, species2], name="C")
+
+    # check that the __repr__ method returns the expected string
+    expected_repr = f"ImplicitSpecies(C, 3.0, {[species1, species2]})"
+    assert repr(implicit_species) == expected_repr
+
+    # check that the __str__ method returns the expected string
+    expected_str = "C"
+    assert str(implicit_species) == expected_str
+
+
 def test_implicit_species_concentration():
     """Test that the concentration of an implicit species is computed
     correctly.