implemented SievertsBC

festim/boundary_conditions/sieverts_bc.py

@@ -1,6 +1,5 @@
 import festim as F
 import ufl
-from dolfinx.fem import Expression, Function, Constant
 
 
 def sieverts_law(T, S_0, E_S, pressure):
@@ -16,33 +15,76 @@ class SievertsBC(F.DirichletBC):
     Args:
         subdomain (festim.Subdomain): the subdomain where the boundary
             condition is applied
-        value (float or fem.Constant): the value of the boundary condition
         species (str): the name of the species
+        S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5)
+        E_S (float or fem.Constant): the Sieverts constant activation energy (eV)
+        pressure (float or callable): the pressure at the boundary (Pa)
 
     Attributes:
         subdomain (festim.Subdomain): the subdomain where the boundary
             condition is applied
         value (float or fem.Constant): the value of the boundary condition
         species (festim.Species or str): the name of the species
+        S_0 (float or fem.Constant): the Sieverts constant pre-exponential factor (H/m3/Pa0.5)
+        E_S (float or fem.Constant): the Sieverts constant activation energy (eV)
+        pressure (float or callable): the pressure at the boundary (Pa)
     """
 
     def __init__(self, subdomain, S_0, E_S, pressure, species) -> None:
+        # TODO find a way to have S_0 and E_S as fem.Constant
+        # maybe in create_value()
         self.S_0 = S_0
         self.E_S = E_S
         self.pressure = pressure
 
-        # construct value callable based on args of pressure
-        args_value_fun = ["T"]
+        value = self.create_new_value_function()
+
+        super().__init__(value=value, species=species, subdomain=subdomain)
+
+    def create_new_value_function(self):
+        """Creates a new value function based on the pressure attribute
+
+        Raises:
+            ValueError: if the pressure function is not supported
+
+        Returns:
+            callable: the value function
+        """
         if callable(self.pressure):
-            if "t" in self.pressure.__code__.co_varnames:
-                args_value_fun.append("t")
-            if "x" in self.pressure.__code__.co_varnames:
-                args_value_fun.append("x")
-
-        # FIXME
-        def value_fun(*args):
-            return sieverts_law(
-                T=args[0], S_0=self.S_0, E_S=self.E_S, pressure=pressure(*args)
-            )
-
-        super().__init__(value=value_fun, species=species, subdomain=subdomain)
+            arg_combinations = {
+                ("x",): lambda T, x=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(x=x)
+                ),
+                ("t",): lambda T, t=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(t=t)
+                ),
+                ("T",): lambda T: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(T=T)
+                ),
+                ("t", "x"): lambda T, x=None, t=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(x=x, t=t)
+                ),
+                ("T", "x"): lambda T, x=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(x=x, T=T)
+                ),
+                ("T", "t"): lambda T, t=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(t=t, T=T)
+                ),
+                ("T", "t", "x"): lambda T, x=None, t=None: sieverts_law(
+                    T, self.S_0, self.E_S, self.pressure(x=x, t=t, T=T)
+                ),
+            }
+
+            # get the arguments of the pressure function
+            args = self.pressure.__code__.co_varnames
+            key = tuple(sorted(args))
+
+            # get the lambda function based on the argument combination
+            if key not in arg_combinations:
+                raise ValueError("pressure function not supported")
+
+            func = arg_combinations[key]
+
+            return func
+        else:
+            return lambda T: sieverts_law(T, self.S_0, self.E_S, self.pressure)

test/test_sievertsbc.py

@@ -0,0 +1,102 @@
+import festim as F
+import ufl
+import pytest
+import numpy as np
+from dolfinx import fem
+import dolfinx.mesh
+from mpi4py import MPI
+
+
+def sieverts_law(T, S_0, E_S, pressure):
+    """Applies the Sieverts law to compute the concentration at the boundary"""
+    S = S_0 * ufl.exp(-E_S / F.k_B / T)
+    return S * pressure**0.5
+
+
+@pytest.mark.parametrize(
+    "pressure",
+    [
+        100,
+        lambda x: x,
+        lambda t: t,
+        lambda T: T,
+        lambda x, t: x + 2 * t,
+        lambda x, T: x + 2 * T,
+        lambda t, T: 2 * t + 3 * T,
+        lambda x, t, T: x + 2 * t + 3 * T,
+        lambda t: 1 if t < 0.5 else 2,
+    ],
+)
+def test_create_new_value_function(pressure):
+    my_BC = F.SievertsBC(
+        subdomain=None, S_0=1.0, E_S=1.0, pressure=pressure, species="H"
+    )
+    assert my_BC.value is not None
+    assert callable(my_BC.value)
+
+    pressure_kwargs, value_kwargs = {}, {}
+    if callable(pressure):
+        if "x" in pressure.__code__.co_varnames:
+            pressure_kwargs["x"] = 2.0
+            value_kwargs["x"] = 2.0
+        if "t" in pressure.__code__.co_varnames:
+            pressure_kwargs["t"] = 3.0
+            value_kwargs["t"] = 3.0
+        if "T" in pressure.__code__.co_varnames:
+            pressure_kwargs["T"] = 500.0
+    value_kwargs["T"] = 500.0
+
+    computed_value = my_BC.value(**value_kwargs)
+    if callable(pressure):
+        expected_value = sieverts_law(
+            T=500.0, S_0=1.0, E_S=1.0, pressure=pressure(**pressure_kwargs)
+        )
+    else:
+        expected_value = sieverts_law(T=500.0, S_0=1.0, E_S=1.0, pressure=pressure)
+    assert np.isclose(computed_value, expected_value)
+
+
+@pytest.mark.parametrize(
+    "pressure",
+    [
+        100,
+        lambda x: x[0],
+        lambda t: t,
+        lambda T: T,
+        lambda x, t: x[0] + 2 * t,
+        lambda x, T: x[0] + 2 * T,
+        lambda t, T: 2 * t + 3 * T,
+        lambda x, t, T: x[0] + 2 * t + 3 * T,
+        lambda t: ufl.conditional(ufl.lt(t, 1.0), 100.0, 0.0),
+        lambda t, x: ufl.conditional(ufl.lt(t, 1.0), 100.0 * x[0], 0.0),
+    ],
+)
+def test_integration_with_HTransportProblem(pressure):
+    subdomain = F.SurfaceSubdomain1D(1, x=1)
+    dummy_mat = F.Material(1.0, 1.0, "dummy")
+    vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat)
+
+    mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
+    my_model = F.HydrogenTransportProblem(
+        mesh=F.Mesh(mesh),
+        subdomains=[vol_subdomain, subdomain],
+    )
+    my_model.species = [F.Species("H")]
+    my_bc = F.SievertsBC(
+        subdomain=subdomain,
+        S_0=2.0,
+        E_S=0.5,
+        pressure=pressure,
+        species=my_model.species[0],
+    )
+    my_model.boundary_conditions = [my_bc]
+
+    my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0)
+
+    # RUN
+
+    my_model.initialise()
+
+    assert my_bc.value_fenics is not None
+
+    my_model.run(final_time=2)