Merge pull request #609 from RemDelaporteMathurin/time_as_Constant

Dirichlet BC attempt 3

festim/__init__.py

@@ -18,6 +18,9 @@
 
 from .helpers import as_fenics_constant
 
+from .boundary_conditions.dirichlet_bc import DirichletBC
+from .boundary_conditions.sieverts_bc import SievertsBC
+
 from .material import Material
 
 from .mesh.mesh import Mesh

festim/boundary_conditions/dirichlet_bc.py

@@ -0,0 +1,151 @@
+import festim as F
+import ufl
+from dolfinx import fem
+import numpy as np
+
+
+class DirichletBC:
+    """
+    Dirichlet boundary condition class
+    c = value
+
+    Args:
+        subdomain (festim.Subdomain): the surface 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
+
+    Attributes:
+        subdomain (festim.Subdomain): the surface 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
+        value_fenics (fem.Function or fem.Constant): the value of the boundary condition in
+            fenics format
+        bc_expr (fem.Expression): the expression of the boundary condition that is used to
+            update the value_fenics
+
+    Usage:
+        >>> from festim import DirichletBC
+        >>> DirichletBC(subdomain=my_subdomain, value=1, species="H")
+        >>> DirichletBC(subdomain=my_subdomain, value=lambda x: 1 + x[0], species="H")
+        >>> DirichletBC(subdomain=my_subdomain, value=lambda t: 1 + t, species="H")
+        >>> DirichletBC(subdomain=my_subdomain, value=lambda T: 1 + T, species="H")
+        >>> DirichletBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + t, species="H")
+    """
+
+    def __init__(self, subdomain, value, species) -> None:
+        self.subdomain = subdomain
+        self.value = value
+        self.species = species
+
+        self.value_fenics = None
+        self.bc_expr = None
+
+    @property
+    def value_fenics(self):
+        return self._value_fenics
+
+    @value_fenics.setter
+    def value_fenics(self, value):
+        if value is None:
+            self._value_fenics = value
+            return
+        if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)):
+            raise TypeError(
+                f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not {type(value)}"
+            )
+        self._value_fenics = value
+
+    def define_surface_subdomain_dofs(self, facet_meshtags, mesh, function_space):
+        """Defines the facets and the degrees of freedom of the boundary
+        condition
+
+        Args:
+            facet_meshtags (ddolfinx.mesh.MeshTags): the mesh tags of the
+                surface facets
+            mesh (dolfinx.mesh.Mesh): the mesh
+            function_space (dolfinx.fem.FunctionSpace): the function space
+        """
+        bc_facets = facet_meshtags.find(self.subdomain.id)
+        bc_dofs = fem.locate_dofs_topological(function_space, mesh.fdim, bc_facets)
+        return bc_dofs
+
+    def create_value(
+        self, mesh, function_space: fem.FunctionSpace, temperature, t: fem.Constant
+    ):
+        """Creates the value of the boundary condition as a fenics object and sets it to
+        self.value_fenics.
+        If the value is a constant, it is converted to a fenics.Constant.
+        If the value is a function of t, it is converted to a fenics.Constant.
+        Otherwise, it is converted to a fenics.Function and the
+        expression of the function is stored in self.bc_expr.
+
+        Args:
+            mesh (dolfinx.mesh.Mesh) : the mesh
+            function_space (dolfinx.fem.FunctionSpace): the function space
+            temperature (float): the temperature
+            t (dolfinx.fem.Constant): the time
+        """
+        x = ufl.SpatialCoordinate(mesh)
+
+        if isinstance(self.value, (int, float)):
+            self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.value)
+
+        elif callable(self.value):
+            arguments = self.value.__code__.co_varnames
+
+            if "t" in arguments and "x" not in arguments and "T" not in arguments:
+                # only t is an argument
+                self.value_fenics = F.as_fenics_constant(
+                    mesh=mesh, value=self.value(t=t.value)
+                )
+            else:
+                self.value_fenics = fem.Function(function_space)
+                kwargs = {}
+                if "t" in arguments:
+                    kwargs["t"] = t
+                if "x" in arguments:
+                    kwargs["x"] = x
+                if "T" in arguments:
+                    kwargs["T"] = temperature
+
+                # store the expression of the boundary condition
+                # to update the value_fenics later
+                self.bc_expr = fem.Expression(
+                    self.value(**kwargs),
+                    function_space.element.interpolation_points(),
+                )
+                self.value_fenics.interpolate(self.bc_expr)
+
+    def create_formulation(self, dofs, function_space):
+        """Applies the boundary condition
+        Args:
+            dofs (numpy.ndarray): the degrees of freedom of surface facets
+            function_space (dolfinx.fem.FunctionSpace): the function space
+        """
+        if isinstance(self.value_fenics, fem.Function):
+            form = fem.dirichletbc(
+                value=self.value_fenics,
+                dofs=dofs,
+            )
+        else:
+            form = fem.dirichletbc(
+                value=self.value_fenics,
+                dofs=dofs,
+                V=function_space,
+            )
+        return form
+
+    def update(self, t):
+        """Updates the boundary condition value
+
+        Args:
+            t (float): the time
+        """
+        if callable(self.value):
+            arguments = self.value.__code__.co_varnames
+            if isinstance(self.value_fenics, fem.Constant) and "t" in arguments:
+                self.value_fenics.value = self.value(t=t)
+            else:
+                self.value_fenics.interpolate(self.bc_expr)

festim/boundary_conditions/sieverts_bc.py

@@ -0,0 +1,101 @@
+import festim as F
+import ufl
+
+
+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
+
+
+class SievertsBC(F.DirichletBC):
+    """
+    Sieverts boundary condition class
+
+    c = S * sqrt(pressure)
+    S = S_0 * exp(-E_S / k_B / T)
+
+    Args:
+        subdomain (festim.Subdomain): the subdomain where the boundary
+            condition is applied
+        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)
+
+    Usage:
+        >>> from festim import SievertsBC
+        >>> SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5, species="H")
+        >>> SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x: 1e5 + x[0], species="H")
+        >>> SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda t: 1e5 + t, species="H")
+        >>> SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda T: 1e5 + T, species="H")
+        >>> SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=lambda x, t: 1e5 + x[0] + t, species="H")
+    """
+
+    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
+
+        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):
+            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)

festim/hydrogen_transport_problem.py

@@ -41,8 +41,8 @@ class HydrogenTransportProblem:
         ds (dolfinx.fem.ds): the surface measure of the model
         function_space (dolfinx.fem.FunctionSpace): the function space of the
             model
-        facet_meshtags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model
-        volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the
+        facet_meshtags (dolfinx.mesh.MeshTags): the facet tags of the model
+        volume_meshtags (dolfinx.mesh.MeshTags): the volume tags of the
             model
         formulation (ufl.form.Form): the formulation of the model
         solver (dolfinx.nls.newton.NewtonSolver): the solver of the model
@@ -96,6 +96,7 @@ def __init__(
         self.volume_meshtags = None
         self.formulation = None
         self.volume_subdomains = []
+        self.bc_forms = []
 
     @property
     def temperature(self):
@@ -112,7 +113,12 @@ def initialise(self):
         self.define_function_space()
         self.define_markers_and_measures()
         self.assign_functions_to_species()
+
+        self.t = fem.Constant(self.mesh.mesh, 0.0)
+
+        self.define_boundary_conditions()
         self.create_formulation()
+        self.create_solver()
         self.defing_export_writers()
 
     def defing_export_writers(self):
@@ -177,6 +183,21 @@ def define_markers_and_measures(self):
             "dx", domain=self.mesh.mesh, subdomain_data=self.volume_meshtags
         )
 
+    def define_boundary_conditions(self):
+        """Defines the dirichlet boundary conditions of the model"""
+        for bc in self.boundary_conditions:
+            if isinstance(bc, F.DirichletBC):
+                bc_dofs = bc.define_surface_subdomain_dofs(
+                    self.facet_meshtags, self.mesh, self.function_space
+                )
+                bc.create_value(
+                    self.mesh.mesh, self.function_space, self.temperature, self.t
+                )
+                form = bc.create_formulation(
+                    dofs=bc_dofs, function_space=self.function_space
+                )
+                self.bc_forms.append(form)
+
     def assign_functions_to_species(self):
         """Creates for each species the solution, prev solution and test function"""
         if len(self.species) > 1:
@@ -233,7 +254,9 @@ def create_formulation(self):
     def create_solver(self):
         """Creates the solver of the model"""
         problem = fem.petsc.NonlinearProblem(
-            self.formulation, self.species[0].solution, bcs=self.boundary_conditions
+            self.formulation,
+            self.species[0].solution,
+            bcs=self.bc_forms,
         )
         solver = NewtonSolver(MPI.COMM_WORLD, problem)
         self.solver = solver
@@ -248,7 +271,6 @@ def run(self, final_time: float):
             list of float: the times of the simulation
             list of float: the fluxes of the simulation
         """
-        t = 0
         times, flux_values = [], []
 
         mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w")
@@ -261,19 +283,23 @@ def run(self, final_time: float):
         progress = tqdm.autonotebook.tqdm(
             desc="Solving H transport problem", total=final_time
         )
-        while t < final_time:
-            progress.update(float(self.dt))
-            t += float(self.dt)
+        while self.t.value < final_time:
+            progress.update(self.dt.value)
+            self.t.value += self.dt.value
+
+            # update boundary conditions
+            for bc in self.boundary_conditions:
+                bc.update(float(self.t))
 
             self.solver.solve(self.u)
 
-            mobile_xdmf.write_function(self.u, t)
+            mobile_xdmf.write_function(self.u, float(self.t))
 
             surface_flux = form(D * dot(grad(cm), n) * self.ds(2))
 
             flux = assemble_scalar(surface_flux)
             flux_values.append(flux)
-            times.append(t)
+            times.append(float(self.t))
 
             # update previous solution
             self.u_n.x.array[:] = self.u.x.array[:]