From af9dcb1eeaa7b05f64ae2a09aab9206eda8a8156 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 9 Jul 2026 11:43:37 -0400 Subject: [PATCH 1/2] start of rework --- src/festim/__init__.py | 2 +- src/festim/reaction.py | 251 ++++++++++++++++++++++++++++------------- 2 files changed, 174 insertions(+), 79 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index b2095f861..50fbf2298 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -73,7 +73,7 @@ from .mesh.mesh_1d import Mesh1D from .mesh.mesh_from_xdmf import MeshFromXDMF from .problem import ProblemBase -from .reaction import Reaction +from .reaction import BaseReaction, Reaction from .settings import Settings from .source import HeatSource, ParticleSource, SourceBase from .species import ImplicitSpecies, Species, find_species_from_name diff --git a/src/festim/reaction.py b/src/festim/reaction.py index 588fef7ea..2a6bcbb84 100644 --- a/src/festim/reaction.py +++ b/src/festim/reaction.py @@ -9,72 +9,48 @@ from festim.subdomain.volume_subdomain import VolumeSubdomain1D as VS1D -class Reaction: - """A reaction between two species, with a forward and backward rate. +class BaseReaction: + """A generic reaction between one or more reactant species and zero or more + product species, taking place within a volume. + + The forward and backward rates are not restricted to Arrhenius laws: they can + be a float, a ufl expression, or a callable that takes the temperature and + returns either of those. This makes it possible to build other reaction types + (eg. trapping, hydride formation) on top of this class. Arguments: 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. + F.ImplicitSpecies]]): The reactant(s). + product (Optional[Union[F.Species, List[F.Species]]]): The product(s). + forward_rate: The forward reaction rate. Can be a float, a ufl expression, + or a callable of temperature returning either of those. volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place. + backward_rate: The backward reaction rate, following the same rules as + forward_rate. If None, the reaction is irreversible. Attributes: 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. + F.ImplicitSpecies]]): The reactant(s). + product (Optional[Union[F.Species, List[F.Species]]]): The product(s). + forward_rate: The forward reaction rate. + backward_rate: The backward reaction rate. volume (F.VolumeSubdomain1D): The volume subdomain where the reaction takes place. - - Examples: - - :: testsetup:: Reaction - - from festim import Reaction, Species, ImplicitSpecies - - :: testcode:: Reaction - - # create a volume subdomain - # create two species - reactant = [F.Species("A"), F.Species("B")] - - # create a product species - product = F.Species("C") - - # create a reaction between the two species - 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 - - # compute the reaction term at a given temperature - temperature = 300.0 - reaction_term = reaction.reaction_term(temperature) """ def __init__( self, reactant: _Species | _ImplicitSpecies | list[_Species | _ImplicitSpecies], - k_0: float, - E_k: float, + product: Union[_Species, list[_Species]] | None, + forward_rate, volume: VS1D, - product: Union[_Species, list[_Species]] | None = [], - p_0: float | None = None, - E_p: float | None = None, + backward_rate=None, ) -> None: - 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.forward_rate = forward_rate + self.backward_rate = backward_rate self.volume = volume @property @@ -99,12 +75,14 @@ def reactant(self, 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({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" # noqa: E501 + return ( + f"{type(self).__name__}({reactants} <--> {products}, " + f"{self.forward_rate}, {self.backward_rate})" + ) def __str__(self) -> str: reactants = " + ".join([str(reactant) for reactant in self.reactant]) @@ -114,6 +92,12 @@ def __str__(self) -> str: products = self.product return f"{reactants} <--> {products}" + @staticmethod + def _evaluate_rate(rate, temperature): + """Evaluate a rate at a given temperature, allowing rates expressed as + floats, ufl expressions, or callables of temperature.""" + return rate(temperature) if callable(rate) else rate + def reaction_term( self, temperature, @@ -155,36 +139,13 @@ def get_concentration(species): else: return species.concentration - 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 is None: - raise ValueError( - "p_0 cannot be None when reaction products are present." - ) - elif self.E_p is None: - raise ValueError( - "E_p cannot be None when reaction products are present." - ) - # reaction rates - k = self.k_0 * exp(-self.E_k / (_k_B * temperature)) - - if self.p_0 and self.E_p: - p = self.p_0 * exp(-self.E_p / (_k_B * temperature)) - elif self.p_0: - p = self.p_0 - else: - p = 0 + k = self._evaluate_rate(self.forward_rate, temperature) + p = ( + self._evaluate_rate(self.backward_rate, temperature) + if self.backward_rate is not None + else 0 + ) # if reactant_concentrations is provided, use these concentrations reactants = self.reactant @@ -222,3 +183,137 @@ def get_concentration(species): product_of_products = 0 return k * product_of_reactants - (p * product_of_products) + + +class Reaction(BaseReaction): + """A reaction between two species, with a forward and backward rate built + from Arrhenius laws. This is typically used to model trapping/detrapping. + + Arguments: + 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: + 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. + + Examples: + + :: testsetup:: Reaction + + from festim import Reaction, Species, ImplicitSpecies + + :: testcode:: Reaction + + # create a volume subdomain + # create two species + reactant = [F.Species("A"), F.Species("B")] + + # create a product species + product = F.Species("C") + + # create a reaction between the two species + 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 + + # compute the reaction term at a given temperature + temperature = 300.0 + reaction_term = reaction.reaction_term(temperature) + """ + + def __init__( + self, + reactant: _Species | _ImplicitSpecies | list[_Species | _ImplicitSpecies], + k_0: float, + E_k: float, + volume: VS1D, + product: Union[_Species, list[_Species]] | None = None, + p_0: float | None = None, + E_p: float | None = None, + ) -> None: + if product is None: + product = [] + + self.k_0 = k_0 + self.E_k = E_k + self.p_0 = p_0 + self.E_p = E_p + + def forward_rate(temperature): + return self.k_0 * exp(-self.E_k / (_k_B * temperature)) + + backward_rate = None + if self.p_0 is not None: + + def backward_rate(temperature): + if self.p_0 and self.E_p: + return self.p_0 * exp(-self.E_p / (_k_B * temperature)) + elif self.p_0: + return self.p_0 + else: + return 0 + + super().__init__( + reactant=reactant, + product=product, + forward_rate=forward_rate, + backward_rate=backward_rate, + volume=volume, + ) + + def reaction_term( + self, + temperature, + reactant_concentrations: list | None = None, + product_concentrations: list | None = None, + ) -> Expr: + # validate that p_0/E_p are consistent with the presence of products + # only at this point (not at __init__) so that Reaction objects can be + # built ahead of time even if p_0/E_p aren't set yet + 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 is None: + raise ValueError( + "p_0 cannot be None when reaction products are present." + ) + elif self.E_p is None: + raise ValueError( + "E_p cannot be None when reaction products are present." + ) + + return super().reaction_term( + temperature, reactant_concentrations, product_concentrations + ) + + 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"{type(self).__name__}({reactants} <--> {products}, {self.k_0}, {self.E_k}, {self.p_0}, {self.E_p})" # noqa: E501 From d2cca0186c8f7ec8592f3c8b94dbdf1eca7e5ae9 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jul 2026 12:07:31 -0400 Subject: [PATCH 2/2] rename to GenericReaction --- src/festim/__init__.py | 2 +- src/festim/hydrogen_transport_problem.py | 51 +++++++----------------- src/festim/reaction.py | 35 +++++++++++++++- 3 files changed, 49 insertions(+), 39 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 5375605ca..1795206f0 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -72,7 +72,7 @@ from .mesh.mesh_1d import Mesh1D from .mesh.mesh_from_xdmf import MeshFromXDMF from .problem import ProblemBase -from .reaction import BaseReaction, Reaction +from .reaction import GenericReaction, Reaction from .settings import Settings from .source import HeatSource, ParticleSource, SourceBase from .species import ImplicitSpecies, Species, find_species_from_name diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 19cc7b11f..c7bad6286 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -853,24 +853,13 @@ def create_formulation(self): self.formulation += ((u - u_n) / self.dt) * v * self.dx(vol.id) for reaction in self.reactions: - for reactant in reaction.reactant: - if isinstance(reactant, _species.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 - else: - products = [reaction.product] - for product in products: - self.formulation += ( - -reaction.reaction_term(self.temperature_fenics) - * product.test_function - * self.dx(reaction.volume.id) + # a reaction enters the formulation as a set of volumetric source + # contributions (-R for reactants, +R for products) + for species, value in reaction.source_contributions( + self.temperature_fenics + ): + self.formulation -= ( + value * species.test_function * self.dx(reaction.volume.id) ) # add sources for source in self.sources: @@ -1531,24 +1520,14 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): if reaction.volume != subdomain: continue - # reactant - for reactant in reaction.reactant: - if isinstance(reactant, _species.Species): - form += ( - reaction.reaction_term(self.temperature_fenics) - * reactant.subdomain_to_test_function[subdomain] - * self.dx(subdomain.id) - ) - - # product - if isinstance(reaction.product, list): - products = reaction.product - else: - products = [reaction.product] - for product in products: - form += ( - -reaction.reaction_term(self.temperature_fenics) - * product.subdomain_to_test_function[subdomain] + # a reaction enters the formulation as a set of volumetric source + # contributions (-R for reactants, +R for products) + for species, value in reaction.source_contributions( + self.temperature_fenics + ): + form -= ( + value + * species.subdomain_to_test_function[subdomain] * self.dx(subdomain.id) ) diff --git a/src/festim/reaction.py b/src/festim/reaction.py index 2a6bcbb84..d3e85db3e 100644 --- a/src/festim/reaction.py +++ b/src/festim/reaction.py @@ -9,7 +9,7 @@ from festim.subdomain.volume_subdomain import VolumeSubdomain1D as VS1D -class BaseReaction: +class GenericReaction: """A generic reaction between one or more reactant species and zero or more product species, taking place within a volume. @@ -184,8 +184,39 @@ def get_concentration(species): return k * product_of_reactants - (p * product_of_products) + def source_contributions(self, temperature): + """Express the reaction as a set of volumetric source contributions, one + per participating :class:`~festim.species.Species`. -class Reaction(BaseReaction): + Each contribution is returned in the same sign convention as a + :class:`~festim.source.SourceBase`, i.e. the returned ``value`` is meant + to enter the residual as ``-value * test_function * dx``. Reactants are + consumed (``value = -R``) and products are produced (``value = +R``), + where ``R`` is :meth:`reaction_term`. Implicit species (which have no + governing equation) are skipped. + + Arguments: + temperature: The temperature at which the reaction term is computed. + + Returns: + A list of ``(species, value)`` tuples. + """ + rate = self.reaction_term(temperature) + + products = self.product if isinstance(self.product, list) else [self.product] + + contributions = [ + (reactant, -rate) + for reactant in self.reactant + if isinstance(reactant, _Species) + ] + contributions += [ + (product, rate) for product in products if isinstance(product, _Species) + ] + return contributions + + +class Reaction(GenericReaction): """A reaction between two species, with a forward and backward rate built from Arrhenius laws. This is typically used to model trapping/detrapping.