SurfaceKineticsBC

festim/__init__.py

@@ -56,6 +56,7 @@
 from .boundary_conditions.fluxes.convective_flux import ConvectiveFlux
 from .boundary_conditions.fluxes.flux_custom import CustomFlux
 from .boundary_conditions.fluxes.mass_flux import MassFlux
+from .boundary_conditions.fluxes.surface_kinetics import SurfaceKinetics
 
 from .exports.exports import Exports
 from .exports.export import Export
@@ -83,6 +84,7 @@
 from .exports.derived_quantities.total_volume import TotalVolume
 from .exports.derived_quantities.average_surface import AverageSurface
 from .exports.derived_quantities.point_value import PointValue
+from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen
 
 from .exports.derived_quantities.derived_quantities import DerivedQuantities
 

festim/boundary_conditions/fluxes/surface_kinetics.py

@@ -0,0 +1,92 @@
+from festim import FluxBC
+from fenics import *
+import sympy as sp
+
+
+class SurfaceKinetics(FluxBC):
+    """ """
+
+    def __init__(
+        self, k_sb, k_bs, l_abs, N_s, N_b, J_vs, surfaces, initial_condition, **prms
+    ):
+        super().__init__(surfaces=surfaces, field=0)
+        self.k_sb = k_sb
+        self.k_bs = k_bs
+        self.J_vs = J_vs
+        self.l_abs = l_abs
+        self.N_s = N_s
+        self.N_b = N_b
+        self.J_vs = J_vs
+        self.initial_condition = initial_condition
+        self.prms = prms
+        self.convert_prms()
+
+        self.solutions = [None] * len(self.surfaces)
+        self.previous_solutions = [None] * len(self.surfaces)
+        self.test_functions = [None] * len(self.surfaces)
+        self.post_processing_solutions = [None] * len(self.surfaces)
+
+    def create_form(self, solute, solute_prev, solute_test_function, T, ds, dt):
+        """Creates the general form associated with the surface species
+        d c_s/ dt = k_bs l_abs c_m (1 - c_s/N_s) - k_sb c_s (1 - c_b/N_b) + J_vs
+
+        Args:
+            solution (fenics.Function or ufl.Indexed): mobile solution for "current"
+                timestep
+            previous_solution (fenics.Function or ufl.Indexed): mobile solution for
+                "previous" timestep
+            test_function (fenics.TestFunction or ufl.Indexed): mobile test function
+            T (festim.Temperature): the temperature of the simulation
+            ds (fenics.Measure): the ds measure of the sim
+            dt (festim.Stepsize): the step-size
+        """
+
+        l_abs = self.l_abs
+        N_s = self.N_s
+        N_b = self.N_b
+        self.F = 0
+
+        for i, surf in enumerate(self.surfaces):
+
+            J_vs = self.J_vs
+            if callable(J_vs):
+                J_vs = J_vs(self.solutions[i], T.T, **self.prms)
+            k_sb = self.k_sb
+            if callable(k_sb):
+                k_sb = k_sb(self.solutions[i], T.T, **self.prms)
+            k_bs = self.k_bs
+            if callable(k_bs):
+                k_bs = k_bs(self.solutions[i], T.T, **self.prms)
+
+            J_sb = k_sb * self.solutions[i] * (1 - solute / N_b)
+            J_bs = k_bs * (solute * l_abs) * (1 - self.solutions[i] / N_s)
+
+            if dt is not None:
+                # Surface concentration form
+                self.F += (
+                    (self.solutions[i] - self.previous_solutions[i])
+                    / dt.value
+                    * self.test_functions[i]
+                    * ds(surf)
+                )
+                # Flux to solute species
+                self.F += (
+                    -l_abs
+                    * (solute - solute_prev)
+                    / dt.value
+                    * solute_test_function
+                    * ds(surf)
+                )
+
+            self.F += -(J_vs + J_bs - J_sb) * self.test_functions[i] * ds(surf)
+            self.F += (J_bs - J_sb) * solute_test_function * ds(surf)
+
+        self.sub_expressions += [expression for expression in self.prms.values()]
+
+    def convert_prms(self):
+        # create Expressions or Constant for all parameters
+        for key, value in self.prms.items():
+            if isinstance(value, (int, float)):
+                self.prms[key] = Constant(value)
+            else:
+                self.prms[key] = Expression(sp.printing.ccode(value), t=0, degree=1)

festim/concentration/mobile.py

@@ -1,4 +1,4 @@
-from festim import Concentration, FluxBC, k_B, RadioactiveDecay
+from festim import Concentration, FluxBC, k_B, RadioactiveDecay, SurfaceKinetics
 from fenics import *
 
 
@@ -38,7 +38,7 @@ def create_form(self, materials, mesh, T, dt=None, traps=None, soret=False):
         self.F = 0
         self.create_diffusion_form(materials, mesh, T, dt=dt, traps=traps, soret=soret)
         self.create_source_form(mesh.dx)
-        self.create_fluxes_form(T, mesh.ds)
+        self.create_fluxes_form(T, mesh.ds, dt)
 
     def create_diffusion_form(
         self, materials, mesh, T, dt=None, traps=None, soret=False
@@ -171,7 +171,7 @@ def create_source_form(self, dx):
         self.F += F_source
         self.sub_expressions += expressions_source
 
-    def create_fluxes_form(self, T, ds):
+    def create_fluxes_form(self, T, ds, dt=None):
         """Modifies the formulation and adds fluxes based
         on parameters in self.boundary_conditions
         """
@@ -184,12 +184,25 @@ def create_fluxes_form(self, T, ds):
         for bc in self.boundary_conditions:
             if bc.field != "T":
                 if isinstance(bc, FluxBC):
-                    bc.create_form(T.T, solute)
+                    if isinstance(bc, SurfaceKinetics):
+                        bc.create_form(
+                            solute,
+                            self.previous_solution,
+                            self.test_function,
+                            T,
+                            ds,
+                            dt,
+                        )
+                        F += bc.F
+                        expressions_fluxes += bc.sub_expressions
+
+                    else:
+                        bc.create_form(T.T, solute)
+                        for surf in bc.surfaces:
+                            F += -self.test_function * bc.form * ds(surf)
                     # TODO : one day we will get rid of this huge expressions list
                     expressions_fluxes += bc.sub_expressions
 
-                    for surf in bc.surfaces:
-                        F += -self.test_function * bc.form * ds(surf)
         self.F_fluxes = F
         self.F += F
         self.sub_expressions += expressions_fluxes

festim/exports/derived_quantities/adsorbed_hydrogen.py

@@ -0,0 +1,39 @@
+from festim import SurfaceQuantity, k_B
+import fenics as f
+
+
+class AdsorbedHydrogen(SurfaceQuantity):
+    """
+    Computes the hydrogen surface concentration on a surface
+
+    Args:
+        surface (int): the surface id
+
+    Attribtutes
+        surface (int): the surface id
+        export_unit (str): the unit of the derived quantity in the export file
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen adsorbed field
+
+    """
+
+    def __init__(self, surface) -> None:
+        super().__init__(field="adsorbed", surface=surface)
+
+    @property
+    def export_unit(self):
+        return f"H m-2"
+
+    @property
+    def title(self):
+        quantity_title = f"Concentration of adsorbed H on surface {self.surface}"
+        if self.show_units:
+            return quantity_title + f" ({self.export_unit})"
+        else:
+            return quantity_title
+
+    def compute(self):
+        return f.assemble(self.function * self.ds(self.surface))

festim/exports/exports.py

@@ -93,7 +93,13 @@ def write(self, label_to_function, dx):
                                 label_to_function[quantity.field] = f.project(
                                     label_to_function[quantity.field], self.V_DG1
                                 )
-                        quantity.function = label_to_function[quantity.field]
+                        if isinstance(quantity, festim.AdsorbedHydrogen):
+                            for surf_funcs in label_to_function[quantity.field]:
+                                if quantity.surface in surf_funcs[1]:
+                                    ind = surf_funcs[1].index(quantity.surface)
+                                    quantity.function = surf_funcs[0][ind]
+                        else:
+                            quantity.function = label_to_function[quantity.field]
                     export.compute(self.t)
                 # export derived quantities
                 if export.is_export(self.t, self.final_time, self.nb_iterations):

festim/generic_simulation.py

@@ -485,6 +485,12 @@ def update_post_processing_solutions(self):
                 [self.mobile.post_processing_solution]
                 + [trap.post_processing_solution for trap in self.traps]
             ),
+            # add tuples (pp_solutions, surfaces) for each SurfaceKinetics boundary condition
+            "adsorbed": [
+                (bc.post_processing_solutions, bc.surfaces)
+                for bc in self.boundary_conditions
+                if isinstance(bc, festim.SurfaceKinetics)
+            ],
         }
         for trap in self.traps:
             label_to_function[trap.id] = trap.post_processing_solution

festim/h_transport_problem.py

@@ -1,6 +1,7 @@
 from fenics import *
 import festim
 import warnings
+import sympy as sp
 
 
 class HTransportProblem:
@@ -110,15 +111,32 @@ def define_function_space(self, mesh):
 
         # function space for H concentrations
         nb_traps = len(self.traps)
-        if nb_traps == 0:
+
+        # the number of surfaces where SurfaceKinetics is used
+        nb_adsorbed = sum(
+            [
+                len(bc.surfaces)
+                for bc in self.boundary_conditions
+                if isinstance(bc, festim.SurfaceKinetics)
+            ]
+        )
+
+        if nb_traps == 0 and nb_adsorbed == 0:
             V = FunctionSpace(mesh.mesh, element_solute, order_solute)
+        elif nb_traps == 0 and nb_adsorbed != 0:
+            solute = FiniteElement(element_solute, mesh.mesh.ufl_cell(), order_solute)
+            adsorbed = FiniteElement("R", mesh.mesh.ufl_cell(), 0)
+            element = [solute] + [adsorbed] * nb_adsorbed
+            V = FunctionSpace(mesh.mesh, MixedElement(element))
         else:
             solute = FiniteElement(element_solute, mesh.mesh.ufl_cell(), order_solute)
             traps = FiniteElement(
                 self.settings.traps_element_type, mesh.mesh.ufl_cell(), order_trap
             )
-            element = [solute] + [traps] * nb_traps
+            adsorbed = FiniteElement("R", mesh.mesh.ufl_cell(), 0)
+            element = [solute] + [traps] * nb_traps + [adsorbed] * nb_adsorbed
             V = FunctionSpace(mesh.mesh, MixedElement(element))
+
         self.V = V
         self.V_CG1 = FunctionSpace(mesh.mesh, "CG", 1)
         self.V_DG1 = FunctionSpace(mesh.mesh, "DG", 1)
@@ -141,11 +159,31 @@ def initialise_concentrations(self):
             self.mobile.previous_solution = self.u_n
             self.mobile.test_function = self.v
         else:
-            for i, concentration in enumerate([self.mobile, *self.traps]):
-                concentration.solution = self.u.sub(i)
-                # concentration.solution = list(split(self.u))[i]
-                concentration.previous_solution = self.u_n.sub(i)
-                concentration.test_function = list(split(self.v))[i]
+            conc_list = [self.mobile]
+            if self.traps:
+                conc_list += [*self.traps]
+            if any(
+                isinstance(bc, festim.SurfaceKinetics)
+                for bc in self.boundary_conditions
+            ):
+                conc_list += [
+                    bc
+                    for bc in self.boundary_conditions
+                    if isinstance(bc, festim.SurfaceKinetics)
+                ]
+
+            for i, concentration in enumerate(conc_list):
+                if isinstance(concentration, festim.SurfaceKinetics):
+                    # iterate through each surface of each SurfaceKinetics
+                    for j in range(len(concentration.surfaces)):
+                        concentration.solutions[j] = self.u.sub(i + j)
+                        concentration.previous_solutions[j] = self.u_n.sub(i + j)
+                        concentration.test_functions[j] = list(split(self.v))[i + j]
+                else:
+                    concentration.solution = self.u.sub(i)
+                    # concentration.solution = list(split(self.u))[i]
+                    concentration.previous_solution = self.u_n.sub(i)
+                    concentration.test_function = list(split(self.v))[i]
 
         print("Defining initial values")
         field_to_component = {
@@ -176,6 +214,23 @@ def initialise_concentrations(self):
                     functionspace, value, label=ini.label, time_step=ini.time_step
                 )
 
+        if any(
+            isinstance(bc, festim.SurfaceKinetics) for bc in self.boundary_conditions
+        ):
+            # iterate through each surface of each SurfaceKinetics
+            for i, bc in enumerate(
+                [
+                    bc
+                    for bc in self.boundary_conditions
+                    if isinstance(bc, festim.SurfaceKinetics)
+                ],
+                len(self.traps) + 1,
+            ):
+                for j in range(len(bc.previous_solutions)):
+                    functionspace = self.V.sub(i + j).collapse()
+                    comp = interpolate(Constant(bc.initial_condition), functionspace)
+                    assign(bc.previous_solutions[j], comp)
+
         # initial guess needs to be non zero if chemical pot
         if self.settings.chemical_pot:
             if self.V.num_sub_spaces() == 0:
@@ -189,9 +244,22 @@ def initialise_concentrations(self):
         # this is needed to correctly create the formulation
         # TODO: write a test for this?
         if self.V.num_sub_spaces() != 0:
-            for i, concentration in enumerate([self.mobile, *self.traps]):
+            """
+            for i, concentration in enumerate(сonc_list):
                 concentration.previous_solution = list(split(self.u_n))[i]
                 concentration.solution = list(split(self.u))[i]
+            """
+
+            for i, concentration in enumerate(conc_list):
+                if isinstance(concentration, festim.SurfaceKinetics):
+                    for j in range(len(concentration.surfaces)):
+                        concentration.solutions[j] = list(split(self.u))[i + j]
+                        concentration.previous_solutions[j] = list(split(self.u_n))[
+                            i + j
+                        ]
+                else:
+                    concentration.solution = list(split(self.u))[i]
+                    concentration.previous_solution = list(split(self.u_n))[i]
 
     def define_variational_problem(self, materials, mesh, dt=None):
         """Creates the variational problem for hydrogen transport (form,
@@ -326,6 +394,17 @@ def update_post_processing_solutions(self, exports):
         for i, trap in enumerate(self.traps, 1):
             trap.post_processing_solution = res[i]
 
+        for i, bc in enumerate(
+            [
+                bc
+                for bc in self.boundary_conditions
+                if isinstance(bc, festim.SurfaceKinetics)
+            ],
+            len(self.traps) + 1,
+        ):
+            for j in range(len(bc.post_processing_solutions)):
+                bc.post_processing_solutions[j] = res[i + j]
+
         if self.settings.chemical_pot:
             self.mobile.post_processing_solution_to_concentration()
         else: