Merge remote-tracking branch 'upstream/fenicsx' into change-of-var-me…

…thod

docs/source/devguide/index.rst

@@ -1,3 +1,5 @@
+.. _developers_guide:
+
 =================
 Developer's Guide
 =================

docs/source/userguide/beginners_guide.rst

@@ -10,7 +10,7 @@ To make the most out of FESTIM, a basic knowledge of `Python <https://www.learnp
 
 Depending on how you install FESTIM, new users should familiarise with either `Docker <https://www.docker.com/>`_ or `Conda <https://anaconda.org/>`_.
 
-FESTIM is under version control with `git <https://git-scm.com/>`_. Even though users don't need git to install FESTIM, it is convenient to have a basic understanding of how it works. You can fin `git turorials <https://git-scm.com/doc>`_ to help you getting started. The `FESTIM source code <https://github.com/RemDelaporteMathurin/FESTIM>`_ is hosted on GitHub. It is recommended to sign up for free to GitHub in order to receive the latest updates, follow the development and submit `bug reports <https://github.com/RemDelaporteMathurin/FESTIM/issues/new/choose>`_.
+FESTIM is under version control with `git <https://git-scm.com/>`_. Even though users don't need git to install FESTIM, it is convenient to have a basic understanding of how it works. You can find `git tutorials <https://git-scm.com/doc>`_ to help you getting started. The `FESTIM source code <https://github.com/festim-dev/FESTIM>`_ is hosted on GitHub. Git knowledge and a GitHub account is required to open issues and :ref:`contribute to the code <developers_guide>`.
 
 ----------------------
 Where can I find help?

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,
+)

src/festim/__init__.py

@@ -18,6 +18,9 @@
 )
 from .boundary_conditions.flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC
 from .boundary_conditions.henrys_bc import HenrysBC
+from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC
+from .boundary_conditions.surface_reaction import SurfaceReactionBC
+
 from .boundary_conditions.sieverts_bc import SievertsBC
 from .exports.average_surface import AverageSurface
 from .exports.average_volume import AverageVolume
@@ -37,6 +40,7 @@
 from .hydrogen_transport_problem import (
     HTransportProblemDiscontinuous,
     HydrogenTransportProblem,
+    HTransportProblemPenalty,
 )
 from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 from .initial_condition import InitialCondition, InitialTemperature

src/festim/boundary_conditions/__init__.py

@@ -1,4 +1,11 @@
 from .dirichlet_bc import FixedConcentrationBC, FixedTemperatureBC
 from .flux_bc import HeatFluxBC, ParticleFluxBC
+from .surface_reaction import SurfaceReactionBC
 
-__all__ = ["FixedConcentrationBC", "FixedTemperatureBC", "ParticleFluxBC", "HeatFluxBC"]
+__all__ = [
+    "FixedConcentrationBC",
+    "FixedTemperatureBC",
+    "ParticleFluxBC",
+    "SurfaceReactionBC",
+    "HeatFluxBC",
+]

src/festim/boundary_conditions/surface_reaction.py

@@ -0,0 +1,127 @@
+from festim.boundary_conditions import ParticleFluxBC
+from festim import k_B
+from dolfinx import fem
+import ufl
+
+
+class SurfaceReactionBCpartial(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 / (k_B * temperature))
+        kd = self.k_d0 * ufl.exp(-self.E_kd / (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