Refactoring

src/festim/__init__.py

@@ -44,22 +44,22 @@
     Value,
     get_interpolation_points,
 )
+from .problem import ProblemBase
 from .hydrogen_transport_problem import (
     HTransportProblemDiscontinuous,
     HTransportProblemPenalty,
     HydrogenTransportProblem,
+    HydrogenTransportProblemDiscontinuousChangeVar,
 )
 from .initial_condition import InitialCondition, InitialTemperature
 from .material import Material
 from .mesh.mesh import Mesh
 from .mesh.mesh_1d import Mesh1D
 from .mesh.mesh_from_xdmf import MeshFromXDMF
-from .problem import ProblemBase
-from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 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/coupled_heat_hydrogen_problem.py

@@ -7,8 +7,8 @@
     HTransportProblemDiscontinuous,
     HTransportProblemPenalty,
     HydrogenTransportProblem,
+    HydrogenTransportProblemDiscontinuousChangeVar,
 )
-from festim.problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 
 
 class CoupledTransientHeatTransferHydrogenTransport:

src/festim/hydrogen_transport_problem.py

@@ -1,4 +1,5 @@
 from collections.abc import Callable
+from typing import List
 
 from mpi4py import MPI
 
@@ -1471,3 +1472,259 @@
             entity_maps=entity_maps,
         )
         self.J = dolfinx.fem.form(J, entity_maps=entity_maps)
+
+
+class HydrogenTransportProblemDiscontinuousChangeVar(HydrogenTransportProblem):
+    species: List[_species.Species]
+
+    def initialise(self):
+        # check if a SurfaceReactionBC is given
+        for bc in self.boundary_conditions:
+            if isinstance(bc, (boundary_conditions.SurfaceReactionBC)):
+                raise ValueError(
+                    f"{type(bc)} not implemented for HydrogenTransportProblemDiscontinuousChangeVar"
+                )
+            if isinstance(bc, boundary_conditions.ParticleFluxBC):
+                if bc.species_dependent_value:
+                    raise ValueError(
+                        f"{type(bc)} concentration-dependent not implemented for HydrogenTransportProblemDiscontinuousChangeVar"
+                    )
+
+        super().initialise()
+
+    def create_formulation(self):
+        """Creates the formulation of the model"""
+
+        self.formulation = 0
+
+        # add diffusion and time derivative for each species
+        for spe in self.species:
+            u = spe.solution
+            u_n = spe.prev_solution
+            v = spe.test_function
+
+            for vol in self.volume_subdomains:
+                D = vol.material.get_diffusion_coefficient(
+                    self.mesh.mesh, self.temperature_fenics, spe
+                )
+                if spe.mobile:
+                    K_S = vol.material.get_solubility_coefficient(
+                        self.mesh.mesh, self.temperature_fenics, spe
+                    )
+                    c = u * K_S
+                    c_n = u_n * K_S
+                else:
+                    c = u
+                    c_n = u_n
+                if spe.mobile:
+                    self.formulation += ufl.dot(D * ufl.grad(c), ufl.grad(v)) * self.dx(
+                        vol.id
+                    )
+
+                if self.settings.transient:
+                    self.formulation += ((c - c_n) / self.dt) * v * self.dx(vol.id)
+
+        for reaction in self.reactions:
+            self.add_reaction_term(reaction)
+
+        # add sources
+        for source in self.sources:
+            self.formulation -= (
+                source.value.fenics_object
+                * source.species.test_function
+                * self.dx(source.volume.id)
+            )
+
+        # add fluxes
+        for bc in self.boundary_conditions:
+            if isinstance(bc, boundary_conditions.ParticleFluxBC):
+                self.formulation -= (
+                    bc.value_fenics
+                    * bc.species.test_function
+                    * self.ds(bc.subdomain.id)
+                )
+
+        # check if each species is defined in all volumes
+        if not self.settings.transient:
+            for spe in self.species:
+                # if species mobile, already defined in diffusion term
+                if not spe.mobile:
+                    not_defined_in_volume = self.volume_subdomains.copy()
+                    for vol in self.volume_subdomains:
+                        # check reactions
+                        for reaction in self.reactions:
+                            if (
+                                spe in reaction.product
+                            ):  # TODO we probably need this in HydrogenTransportProblem too no?
+                                if vol == reaction.volume:
+                                    if vol in not_defined_in_volume:
+                                        not_defined_in_volume.remove(vol)
+
+                    # add c = 0 to formulation where needed
+                    for vol in not_defined_in_volume:
+                        self.formulation += (
+                            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()
+
+    def override_post_processing_solution(self):
+        # override the post-processing solution c = theta * K_S
+        Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
+        Q1 = fem.functionspace(self.mesh.mesh, ("DG", 1))
+
+        for spe in self.species:
+            if not spe.mobile:
+                continue
+            K_S0 = fem.Function(Q0)
+            E_KS = fem.Function(Q0)
+            for subdomain in self.volume_subdomains:
+                entities = subdomain.locate_subdomain_entities_correct(
+                    self.volume_meshtags
+                )
+                K_S0.x.array[entities] = subdomain.material.get_K_S_0(spe)
+                E_KS.x.array[entities] = subdomain.material.get_E_K_S(spe)
+
+            K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+
+            theta = spe.solution
+
+            spe.dg_expr = fem.Expression(
+                theta * K_S, get_interpolation_points(Q1.element)
+            )
+            spe.post_processing_solution = fem.Function(Q1)
+            spe.post_processing_solution.interpolate(
+                spe.dg_expr
+            )  # NOTE: do we need this line since it's in initialise?
+
+    def post_processing(self):
+        # need to compute c = theta * K_S
+        # this expression is stored in species.dg_expr
+        for spe in self.species:
+            if not spe.mobile:
+                continue
+            spe.post_processing_solution.interpolate(spe.dg_expr)
+
+        super().post_processing()
+
+    def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
+        """Creates a dirichlet boundary condition form
+
+        Args:
+            bc (festim.DirichletBC): the boundary condition
+
+        Returns:
+            dolfinx.fem.bcs.DirichletBC: A representation of
+                the boundary condition for modifying linear systems.
+        """
+        # create value_fenics
+        if not self.multispecies:
+            function_space_value = bc.species.sub_function_space
+        else:
+            function_space_value = bc.species.collapsed_function_space
+
+        # create K_S function
+        Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
+        K_S0 = fem.Function(Q0)
+        E_KS = fem.Function(Q0)
+        for subdomain in self.volume_subdomains:
+            entities = subdomain.locate_subdomain_entities_correct(self.volume_meshtags)
+            K_S0.x.array[entities] = subdomain.material.get_K_S_0(bc.species)
+            E_KS.x.array[entities] = subdomain.material.get_E_K_S(bc.species)
+
+        K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+
+        bc.create_value(
+            temperature=self.temperature_fenics,
+            function_space=function_space_value,
+            t=self.t,
+            K_S=K_S,
+        )
+
+        # get dofs
+        if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
+            function_space_dofs = (
+                bc.species.sub_function_space,
+                bc.species.collapsed_function_space,
+            )
+        else:
+            function_space_dofs = bc.species.sub_function_space
+
+        bc_dofs = bc.define_surface_subdomain_dofs(
+            facet_meshtags=self.facet_meshtags,
+            function_space=function_space_dofs,
+        )
+
+        # create form
+        if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
+            # no need to pass the functionspace since value_fenics is already a Function
+            function_space_form = None
+        else:
+            function_space_form = bc.species.sub_function_space
+
+        form = fem.dirichletbc(
+            value=bc.value_fenics,
+            dofs=bc_dofs,
+            V=function_space_form,
+        )
+
+        return form
+
+    def update_time_dependent_values(self):
+        super().update_time_dependent_values()
+
+        if self.temperature_time_dependent:
+            for bc in self.boundary_conditions:
+                if isinstance(bc, boundary_conditions.FixedConcentrationBC):
+                    bc.update(self.t)

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].solution
-            for product in products[1:]:
-                product_of_products *= product.solution
+        # 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)

src/festim/species.py

@@ -225,13 +225,3 @@ def find_species_from_name(name: str, species: list):
         if spe.name == name:
             return spe
     raise ValueError(f"Species {name} not found in list of species")
-
-
-class SpeciesChangeVar(Species):
-    @property
-    def concentration(self):
-        return self._concentration
-
-    @concentration.setter
-    def concentration(self, value):
-        self._concentration = value

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]