removed SpeciesChangeVar

Author: RemDelaporteMathurin

src/festim/__init__.py

@@ -59,7 +59,7 @@
 from .reaction import Reaction
 from .settings import Settings
 from .source import HeatSource, ParticleSource, SourceBase
-from .species import ImplicitSpecies, Species, SpeciesChangeVar, find_species_from_name
+from .species import ImplicitSpecies, Species, find_species_from_name
 from .stepsize import Stepsize
 from .subdomain.interface import Interface
 from .subdomain.surface_subdomain import SurfaceSubdomain, find_surface_from_id

src/festim/hydrogen_transport_problem.py

@@ -1510,38 +1510,8 @@ def create_formulation(self):
                     self.formulation += ((c - c_n) / self.dt) * v * self.dx(vol.id)
 
         for reaction in self.reactions:
-            if isinstance(reaction.product, list):
-                products = reaction.product
-            else:
-                products = [reaction.product]
-
-            # hack enforce the concentration attribute of the species for all species
-            # to be used in reaction.reaction_term
-
-            for spe in self.species:
-                if spe.mobile:
-                    K_S = reaction.volume.material.get_solubility_coefficient(
-                        self.mesh.mesh, self.temperature_fenics, spe
-                    )
-                    assert isinstance(spe, _species.SpeciesChangeVar)
-                    spe.concentration = spe.solution * K_S
+            self.add_reaction_term(reaction)
 
-            # reactant
-            for reactant in reaction.reactant:
-                if isinstance(reactant, festim.species.Species):
-                    self.formulation += (
-                        reaction.reaction_term(self.temperature_fenics)
-                        * reactant.test_function
-                        * self.dx(reaction.volume.id)
-                    )
-
-            # product
-            for product in products:
-                self.formulation += (
-                    -reaction.reaction_term(self.temperature_fenics)
-                    * product.test_function
-                    * self.dx(reaction.volume.id)
-                )
         # add sources
         for source in self.sources:
             self.formulation -= (
@@ -1581,6 +1551,55 @@ def create_formulation(self):
                             spe.solution * spe.test_function * self.dx(vol.id)
                         )
 
+    def add_reaction_term(self, reaction: _reaction.Reaction):
+        """Adds the reaction term to the formulation"""
+
+        products = (
+            reaction.product
+            if isinstance(reaction.product, list)
+            else [reaction.product]
+        )
+
+        # we cannot use the `concentration` attribute of the mobile species and need to use u * K_S instead
+
+        def get_concentrations(species_list) -> List:
+            concentrations = []
+            for spe in species_list:
+                if isinstance(spe, _species.ImplicitSpecies):
+                    concentrations.append(None)
+                elif spe.mobile:
+                    K_S = reaction.volume.material.get_solubility_coefficient(
+                        self.mesh.mesh, self.temperature_fenics, spe
+                    )
+                    concentrations.append(spe.solution * K_S)
+                else:
+                    concentrations.append(None)
+            return concentrations
+
+        reactant_concentrations = get_concentrations(reaction.reactant)
+        product_concentrations = get_concentrations(products)
+
+        # get the reaction term from the reaction
+        reaction_term = reaction.reaction_term(
+            temperature=self.temperature_fenics,
+            reactant_concentrations=reactant_concentrations,
+            product_concentrations=product_concentrations,
+        )
+
+        # add reaction term to formulation
+        # reactant
+        for reactant in reaction.reactant:
+            if isinstance(reactant, festim.species.Species):
+                self.formulation += (
+                    reaction_term * reactant.test_function * self.dx(reaction.volume.id)
+                )
+
+        # product
+        for product in products:
+            self.formulation += (
+                -reaction_term * product.test_function * self.dx(reaction.volume.id)
+            )
+
     def initialise_exports(self):
         self.override_post_processing_solution()
         super().initialise_exports()

src/festim/reaction.py

@@ -1,6 +1,7 @@
-from typing import Optional, Union
+from typing import Optional, Union, List
 
 from ufl import exp
+from ufl.core.expr import Expr
 
 from festim import k_B as _k_B
 from festim.species import ImplicitSpecies as _ImplicitSpecies
@@ -102,11 +103,27 @@ def __str__(self) -> str:
             products = self.product
         return f"{reactants} <--> {products}"
 
-    def reaction_term(self, temperature):
+    def reaction_term(
+        self,
+        temperature,
+        reactant_concentrations: List = None,
+        product_concentrations: List = None,
+    ) -> Expr:
         """Compute the reaction term at a given temperature.
 
         Arguments:
-            temperature (): The temperature at which the reaction term is computed.
+            temperature: The temperature at which the reaction term is computed.
+            reactant_concentrations: The concentrations of the reactants. Must
+                be the same length as the reactants. If None, the ``concentration``
+                attribute of each reactant is used. If an element is None, the
+                ``concentration`` attribute of the reactant is used.
+            product_concentrations: The concentrations of the products. Must
+                be the same length as the products. If None, the ``concentration``
+                attribute of each product is used. If an element is None, the
+                ``concentration`` attribute of the product is used.
+
+        Returns:
+            The reaction term to be used in a formulation.
         """
 
         if self.product == []:
@@ -130,6 +147,9 @@ def reaction_term(self, temperature):
                     "E_p cannot be None when reaction products are present."
                 )
 
+        products = self.product if isinstance(self.product, list) else [self.product]
+
+        # reaction rates
         k = self.k_0 * exp(-self.E_k / (_k_B * temperature))
 
         if self.p_0 and self.E_p:
@@ -139,20 +159,35 @@ def reaction_term(self, temperature):
         else:
             p = 0
 
+        # if reactant_concentrations is provided, use these concentrations
         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
+        if reactant_concentrations is not None:
+            assert len(reactant_concentrations) == len(reactants)
+            for i, reactant in enumerate(reactants):
+                if reactant_concentrations[i] is None:
+                    reactant_concentrations[i] = reactant.concentration
+        else:
+            reactant_concentrations = [reactant.concentration for reactant in reactants]
+
+        # if product_concentrations is provided, use these concentrations
+        if product_concentrations is not None:
+            assert len(product_concentrations) == len(products)
+            for i, product in enumerate(products):
+                if product_concentrations[i] is None:
+                    product_concentrations[i] = product.concentration
         else:
-            products = [self.product]
+            product_concentrations = [product.concentration for product in products]
 
-        if len(products) > 0:
-            product_of_products = products[0].concentration
-            for product in products[1:]:
-                product_of_products *= product.concentration
+        # multiply all concentrations to be used in the term
+        product_of_reactants = reactant_concentrations[0]
+        for reactant_conc in reactant_concentrations[1:]:
+            product_of_reactants *= reactant_conc
+
+        if products:
+            product_of_products = product_concentrations[0]
+            for product_conc in product_concentrations[1:]:
+                product_of_products *= product_conc
         else:
             product_of_products = 0
+
         return k * product_of_reactants - (p * product_of_products)

test/system_tests/test_multi_material_change_of_var.py

@@ -103,8 +103,8 @@ def test_run():
         right_surface,
     ]
 
-    H = F.SpeciesChangeVar("H", mobile=True)
-    trapped_H = F.SpeciesChangeVar("H_trapped", mobile=False)
+    H = F.Species("H", mobile=True)
+    trapped_H = F.Species("H_trapped", mobile=False)
     empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H])
 
     my_model.species = [H, trapped_H]
@@ -193,7 +193,7 @@ def c_exact_bot_np(x):
         bottom_surface,
     ]
 
-    H = F.SpeciesChangeVar("H", mobile=True)
+    H = F.Species("H", mobile=True)
 
     my_model.species = [H]