Surface reactions

examples/surface_reaction_example.py

@@ -0,0 +1,169 @@
+import festim as F
+import numpy as np
+
+import dolfinx.fem as fem
+import ufl
+
+
+class FluxFromSurfaceReaction(F.SurfaceFlux):
+    def __init__(self, reaction: F.SurfaceReactionBC):
+        super().__init__(
+            F.Species(),  # just a dummy species here
+            reaction.subdomain,
+        )
+        self.reaction = reaction.flux_bcs[0]
+
+    def compute(self, ds):
+        self.value = fem.assemble_scalar(
+            fem.form(self.reaction.value_fenics * ds(self.surface.id))
+        )
+        self.data.append(self.value)
+
+
+my_model = F.HydrogenTransportProblem()
+my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000))
+my_mat = F.Material(name="mat", D_0=1, E_D=0)
+vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)
+left = F.SurfaceSubdomain1D(id=1, x=0)
+right = F.SurfaceSubdomain1D(id=2, x=1)
+
+my_model.subdomains = [vol, left, right]
+
+H = F.Species("H")
+D = F.Species("D")
+my_model.species = [H, D]
+
+my_model.temperature = 500
+
+surface_reaction_hd = F.SurfaceReactionBC(
+    reactant=[H, D],
+    gas_pressure=0,
+    k_r0=0.01,
+    E_kr=0,
+    k_d0=0,
+    E_kd=0,
+    subdomain=right,
+)
+
+surface_reaction_hh = F.SurfaceReactionBC(
+    reactant=[H, H],
+    gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0),
+    k_r0=0.02,
+    E_kr=0,
+    k_d0=0.03,
+    E_kd=0,
+    subdomain=right,
+)
+
+surface_reaction_dd = F.SurfaceReactionBC(
+    reactant=[D, D],
+    gas_pressure=0,
+    k_r0=0.01,
+    E_kr=0,
+    k_d0=0,
+    E_kd=0,
+    subdomain=right,
+)
+
+my_model.boundary_conditions = [
+    F.DirichletBC(subdomain=left, value=2, species=H),
+    F.DirichletBC(subdomain=left, value=2, species=D),
+    surface_reaction_hd,
+    surface_reaction_hh,
+    surface_reaction_dd,
+]
+
+H_flux_right = F.SurfaceFlux(H, right)
+H_flux_left = F.SurfaceFlux(H, left)
+D_flux_right = F.SurfaceFlux(D, right)
+D_flux_left = F.SurfaceFlux(D, left)
+HD_flux = FluxFromSurfaceReaction(surface_reaction_hd)
+HH_flux = FluxFromSurfaceReaction(surface_reaction_hh)
+DD_flux = FluxFromSurfaceReaction(surface_reaction_dd)
+my_model.exports = [
+    F.XDMFExport("test.xdmf", H),
+    H_flux_left,
+    H_flux_right,
+    D_flux_left,
+    D_flux_right,
+    HD_flux,
+    HH_flux,
+    DD_flux,
+]
+
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True)
+
+my_model.settings.stepsize = 0.1
+
+my_model.initialise()
+my_model.run()
+
+import matplotlib.pyplot as plt
+
+plt.stackplot(
+    H_flux_left.t,
+    np.abs(H_flux_left.data),
+    np.abs(D_flux_left.data),
+    labels=["H_in", "D_in"],
+)
+plt.stackplot(
+    H_flux_right.t,
+    -np.abs(H_flux_right.data),
+    -np.abs(D_flux_right.data),
+    labels=["H_out", "D_out"],
+)
+plt.legend()
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (atom/m^2/s)")
+plt.figure()
+plt.stackplot(
+    HD_flux.t,
+    np.abs(HH_flux.data),
+    np.abs(HD_flux.data),
+    np.abs(DD_flux.data),
+    labels=["HH", "HD", "DD"],
+)
+plt.legend(reverse=True)
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (molecule/m^2/s)")
+
+
+plt.figure()
+plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)")
+plt.plot(
+    H_flux_right.t,
+    2 * np.array(HH_flux.data) + np.array(HD_flux.data),
+    linestyle="--",
+    label="from reaction rates (H)",
+)
+
+plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)")
+plt.plot(
+    D_flux_right.t,
+    2 * np.array(DD_flux.data) + np.array(HD_flux.data),
+    linestyle="--",
+    label="from reaction rates (D)",
+)
+plt.xlabel("Time (s)")
+plt.ylabel("Flux (atom/m^2/s)")
+plt.legend()
+plt.show()
+
+# check that H_flux_right == 2*HH_flux + HD_flux
+H_flux_from_gradient = -np.array(H_flux_right.data)
+H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data)
+assert np.allclose(
+    H_flux_from_gradient,
+    H_flux_from_reac,
+    rtol=0.5e-2,
+    atol=0.005,
+)
+# check that D_flux_right == 2*DD_flux + HD_flux
+D_flux_from_gradient = -np.array(D_flux_right.data)
+D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data)
+assert np.allclose(
+    D_flux_from_gradient,
+    D_flux_from_reac,
+    rtol=0.5e-2,
+    atol=0.005,
+)

festim/__init__.py

@@ -27,6 +27,7 @@
 from .boundary_conditions.sieverts_bc import SievertsBC
 from .boundary_conditions.henrys_bc import HenrysBC
 from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC
+from .boundary_conditions.surface_reaction import SurfaceReactionBC
 
 from .material import Material
 

festim/boundary_conditions/surface_reaction.py

@@ -0,0 +1,127 @@
+import festim as F
+from dolfinx import fem
+import ufl
+
+
+class SurfaceReactionBCpartial(F.ParticleFluxBC):
+    """Boundary condition representing a surface reaction
+    A + B <-> C
+    where A, B are the reactants and C is the product
+    the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
+    and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
+    The reaction rate is:
+    K = K_r * C_A * C_B - K_d * P_C
+    with C_A, C_B the concentration of species A and B,
+    P_C the partial pressure of species C at the surface.
+
+    This class is used to create the flux of a single species entering the surface
+    Example: The flux of species A entering the surface is K.
+
+
+    Args:
+        reactant (list): list of F.Species objects representing the reactants
+        gas_pressure (float, callable): the partial pressure of the product species
+        k_r0 (float): the pre-exponential factor of the forward reaction rate
+        E_kr (float): the activation energy of the forward reaction rate (eV)
+        k_d0 (float): the pre-exponential factor of the backward reaction rate
+        E_kd (float): the activation energy of the backward reaction rate (eV)
+        subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs
+        species (F.Species): the species to which the flux is applied
+    """
+
+    def __init__(
+        self,
+        reactant,
+        gas_pressure,
+        k_r0,
+        E_kr,
+        k_d0,
+        E_kd,
+        subdomain,
+        species,
+    ):
+        self.reactant = reactant
+        self.gas_pressure = gas_pressure
+        self.k_r0 = k_r0
+        self.E_kr = E_kr
+        self.k_d0 = k_d0
+        self.E_kd = E_kd
+        super().__init__(subdomain=subdomain, value=None, species=species)
+
+    def create_value_fenics(self, mesh, temperature, t: fem.Constant):
+        kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * temperature))
+        kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * temperature))
+        if callable(self.gas_pressure):
+            gas_pressure = self.gas_pressure(t=t)
+        else:
+            gas_pressure = self.gas_pressure
+
+        product_of_reactants = self.reactant[0].concentration
+        for reactant in self.reactant[1:]:
+            product_of_reactants *= reactant.concentration
+
+        self.value_fenics = kd * gas_pressure - kr * product_of_reactants
+
+
+class SurfaceReactionBC:
+    """Boundary condition representing a surface reaction
+    A + B <-> C
+    where A, B are the reactants and C is the product
+    the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T))
+    and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T))
+    The reaction rate is:
+    K = K_r * C_A * C_B - K_d * P_C
+    with C_A, C_B the concentration of species A and B,
+    P_C the partial pressure of species C at the surface.
+
+    The flux of species A entering the surface is K.
+    In the special case where A=B, then the flux of particle entering the surface is 2*K
+
+
+    Args:
+        reactant (list): list of F.Species objects representing the reactants
+        gas_pressure (float, callable): the partial pressure of the product species
+        k_r0 (float): the pre-exponential factor of the forward reaction rate
+        E_kr (float): the activation energy of the forward reaction rate (eV)
+        k_d0 (float): the pre-exponential factor of the backward reaction rate
+        E_kd (float): the activation energy of the backward reaction rate (eV)
+        subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs
+    """
+
+    def __init__(
+        self,
+        reactant,
+        gas_pressure,
+        k_r0,
+        E_kr,
+        k_d0,
+        E_kd,
+        subdomain,
+    ):
+
+        self.reactant = reactant
+        self.gas_pressure = gas_pressure
+        self.k_r0 = k_r0
+        self.E_kr = E_kr
+        self.k_d0 = k_d0
+        self.E_kd = E_kd
+        self.subdomain = subdomain
+
+        # create the flux boundary condition for each reactant
+        self.flux_bcs = [
+            SurfaceReactionBCpartial(
+                reactant=self.reactant,
+                gas_pressure=self.gas_pressure,
+                k_r0=self.k_r0,
+                E_kr=self.E_kr,
+                k_d0=self.k_d0,
+                E_kd=self.E_kd,
+                subdomain=self.subdomain,
+                species=species,
+            )
+            for species in self.reactant
+        ]
+
+    @property
+    def time_dependent(self):
+        return False  # no need to update if only using ufl.conditional objects

festim/hydrogen_transport_problem.py

@@ -248,7 +248,6 @@ def initialise(self):
         self.define_temperature()
         self.define_boundary_conditions()
         self.create_source_values_fenics()
-        self.create_flux_values_fenics()
         self.create_initial_conditions()
         self.create_formulation()
         self.create_solver()
@@ -512,14 +511,30 @@ def define_meshtags_and_measures(self):
         )
 
     def define_boundary_conditions(self):
-        """Defines the dirichlet boundary conditions of the model"""
+        """Create forms for DirichletBC and value_fenics for ParticleFluxBC"""
+        # @jhdark this all_bcs could be a property
+        # I just don't want to modify self.boundary_conditions
+
+        # create all_bcs which includes all flux bcs from SurfaceReactionBC
+        all_bcs = self.boundary_conditions.copy()
         for bc in self.boundary_conditions:
+            if isinstance(bc, F.SurfaceReactionBC):
+                all_bcs += bc.flux_bcs
+                all_bcs.remove(bc)
+
+        for bc in all_bcs:
             if isinstance(bc.species, str):
                 # if name of species is given then replace with species object
                 bc.species = F.find_species_from_name(bc.species, self.species)
             if isinstance(bc, F.DirichletBC):
                 form = self.create_dirichletbc_form(bc)
                 self.bc_forms.append(form)
+            if isinstance(bc, F.ParticleFluxBC):
+                bc.create_value_fenics(
+                    mesh=self.mesh.mesh,
+                    temperature=self.temperature_fenics,
+                    t=self.t,
+                )
 
     def create_dirichletbc_form(self, bc):
         """Creates a dirichlet boundary condition form
@@ -590,18 +605,6 @@ def create_source_values_fenics(self):
                     t=self.t,
                 )
 
-    def create_flux_values_fenics(self):
-        """For each particle flux create the value_fenics"""
-        for bc in self.boundary_conditions:
-            # create value_fenics for all F.ParticleFluxBC objects
-            if isinstance(bc, F.ParticleFluxBC):
-
-                bc.create_value_fenics(
-                    mesh=self.mesh.mesh,
-                    temperature=self.temperature_fenics,
-                    t=self.t,
-                )
-
     def create_initial_conditions(self):
         """For each initial condition, create the value_fenics and assign it to
         the previous solution of the condition's species"""
@@ -695,6 +698,13 @@ def create_formulation(self):
                     * bc.species.test_function
                     * self.ds(bc.subdomain.id)
                 )
+            if isinstance(bc, F.SurfaceReactionBC):
+                for flux_bc in bc.flux_bcs:
+                    self.formulation -= (
+                        flux_bc.value_fenics
+                        * flux_bc.species.test_function
+                        * self.ds(flux_bc.subdomain.id)
+                    )
 
         # check if each species is defined in all volumes
         if not self.settings.transient: