Merge pull request #596 from jhdark/materials_class

Materials class and temperature processing

festim/__init__.py

@@ -16,12 +16,16 @@
 R = 8.314462618  # Gas constant J.mol-1.K-1
 k_B = 8.6173303e-5  # Boltzmann constant eV.K-1
 
+from .helpers import as_fenics_constant
+
+from .material import Material
+
 from .mesh.mesh import Mesh
 from .mesh.mesh_1d import Mesh1D
 
+from .hydrogen_transport_problem import HydrogenTransportProblem
+
 from .species import Species, Trap
 
 from .subdomain.surface_subdomain import SurfaceSubdomain1D
 from .subdomain.volume_subdomain import VolumeSubdomain1D
-
-from .hydrogen_transport_problem import HydrogenTransportProblem

festim/helpers.py

@@ -0,0 +1,24 @@
+from dolfinx import fem
+
+
+def as_fenics_constant(value, mesh):
+    """Converts a value to a dolfinx.Constant
+
+    Args:
+        value (float, int or dolfinx.Constant): the value to convert
+        mesh (dolfinx.mesh.Mesh): the mesh of the domiain
+
+    Returns:
+        dolfinx.Constant: the converted value
+
+    Raises:
+        TypeError: if the value is not a float, an int or a dolfinx.Constant
+    """
+    if isinstance(value, (float, int)):
+        return fem.Constant(mesh, float(value))
+    elif isinstance(value, fem.Constant):
+        return value
+    else:
+        raise TypeError(
+            f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}"
+        )

festim/hydrogen_transport_problem.py

@@ -3,8 +3,8 @@
 import ufl
 from mpi4py import MPI
 from dolfinx.fem import Function
-from dolfinx.mesh import meshtags, locate_entities
-from ufl import TestFunction, dot, grad, exp, Measure
+from dolfinx.mesh import meshtags
+from ufl import TestFunction, dot, grad, Measure
 import numpy as np
 
 import festim as F
@@ -18,25 +18,29 @@ class HydrogenTransportProblem:
         mesh (festim.Mesh): the mesh of the model
         subdomains (list of festim.Subdomain): the subdomains of the model
         species (list of festim.Species): the species of the model
-        temperature (float or dolfinx.Function): the temperature of the model
+        temperature (float or fem.Constant): the temperature of the model
         sources (list of festim.Source): the hydrogen sources of the model
-        boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model
+        boundary_conditions (list of festim.BoundaryCondition): the boundary
+            conditions of the model
         solver_parameters (dict): the solver parameters of the model
         exports (list of festim.Export): the exports of the model
 
     Attributes:
         mesh (festim.Mesh): the mesh of the model
         subdomains (list of festim.Subdomain): the subdomains of the model
         species (list of festim.Species): the species of the model
-        temperature (float or dolfinx.Function): the temperature of the model
-        boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model
+        temperature (fem.Constant): the temperature of the model
+        boundary_conditions (list of festim.BoundaryCondition): the boundary
+            conditions of the model
         solver_parameters (dict): the solver parameters of the model
         exports (list of festim.Export): the exports of the model
         dx (dolfinx.fem.dx): the volume measure of the model
         ds (dolfinx.fem.ds): the surface measure of the model
-        function_space (dolfinx.fem.FunctionSpace): the function space 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 model
+        volume_meshtags (dolfinx.cpp.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
 
@@ -88,12 +92,20 @@ def __init__(
         self.facet_meshtags = None
         self.volume_meshtags = None
         self.formulation = None
+        self.volume_subdomains = []
 
-    def initialise(self):
-        """Initialise the model. Creates suitable function
-        spaces, facet and volume tags...
-        """
+    @property
+    def temperature(self):
+        return self._temperature
+
+    @temperature.setter
+    def temperature(self, value):
+        if value is None:
+            self._temperature = value
+        else:
+            self._temperature = F.as_fenics_constant(value, self.mesh.mesh)
 
+    def initialise(self):
         self.define_function_space()
         self.define_markers_and_measures()
         self.assign_functions_to_species()
@@ -104,6 +116,8 @@ def define_function_space(self):
         self.function_space = fem.FunctionSpace(self.mesh.mesh, elements)
 
     def define_markers_and_measures(self):
+        """Defines the markers and measures of the model"""
+
         dofs_facets, tags_facets = [], []
 
         # TODO this should be a property of mesh
@@ -122,6 +136,7 @@ def define_markers_and_measures(self):
                 tags_facets.append(sub_dom.id)
             if isinstance(sub_dom, F.VolumeSubdomain1D):
                 # find all cells in subdomain and mark them as sub_dom.id
+                self.volume_subdomains.append(sub_dom)
                 entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim)
                 tags_volumes[entities] = sub_dom.id
 
@@ -159,39 +174,38 @@ def create_formulation(self):
         if len(self.species) > 1:
             raise NotImplementedError("Multiple species not implemented yet")
 
-        # TODO expose D_0 and E_D as parameters of a Material class
-        D_0 = fem.Constant(self.mesh.mesh, 1.9e-7)
-        E_D = fem.Constant(self.mesh.mesh, 0.2)
-
-        D = D_0 * exp(-E_D / F.k_B / self.temperature)
-
         # TODO expose dt as parameter of the model
         dt = fem.Constant(self.mesh.mesh, 1 / 20)
 
-        self.D = D  # TODO remove this
         self.dt = dt  # TODO remove this
 
+        self.formulation = 0
+
         for spe in self.species:
             u = spe.solution
             u_n = spe.prev_solution
             v = spe.test_function
 
-            formulation = dot(D * grad(u), grad(v)) * self.dx
-            formulation += ((u - u_n) / dt) * v * self.dx
-
-            # add sources
-            for source in self.sources:
-                # f = Constant(my_mesh.mesh, (PETSc.ScalarType(0)))
-                if source.species == spe:
-                    formulation += source * v * self.dx
-            # add fluxes
-            # TODO implement this
-            # for bc in self.boundary_conditions:
-            #     pass
-            #     if bc.species == spe and bc.type != "dirichlet":
-            #         formulation += bc * v * self.ds
-
-        self.formulation = formulation
+            for vol in self.volume_subdomains:
+                D = vol.material.get_diffusion_coefficient(
+                    self.mesh.mesh, self.temperature
+                )
+
+                self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id)
+                self.formulation += ((u - u_n) / dt) * v * self.dx(vol.id)
+
+                # add sources
+                # TODO implement this
+                # for source in self.sources:
+                #     # f = Constant(my_mesh.mesh, (PETSc.ScalarType(0)))
+                #     if source.species == spe:
+                #         formulation += source * v * self.dx
+                # add fluxes
+                # TODO implement this
+                # for bc in self.boundary_conditions:
+                #     pass
+                #     if bc.species == spe and bc.type != "dirichlet":
+                #         formulation += bc * v * self.ds
 
     def create_solver(self):
         """Creates the solver of the model"""

festim/material.py

@@ -0,0 +1,42 @@
+import ufl
+import festim as F
+
+
+class Material:
+    """
+    Material class
+
+    Args:
+        D_0 (float or fem.Constant): the pre-exponential factor of the
+            diffusion coefficient (m2/s)
+        E_D (float or fem.Constant): the activation energy of the diffusion
+            coeficient (eV)
+        name (str): the name of the material
+
+    Attributes:
+        D_0 (float or fem.Constant): the pre-exponential factor of the
+            diffusion coefficient (m2/s)
+        E_D (float or fem.Constant): the activation energy of the diffusion
+            coeficient (eV)
+        name (str): the name of the material
+    """
+
+    def __init__(self, D_0, E_D, name=None) -> None:
+        self.D_0 = D_0
+        self.E_D = E_D
+        self.name = name
+
+    def get_diffusion_coefficient(self, mesh, temperature):
+        """Defines the diffusion coefficient
+        Args:
+            mesh (dolfinx.mesh.Mesh): the domain mesh
+            temperature (dolfinx.fem.Constant): the temperature
+        Returns:
+            ufl.algebra.Product: the diffusion coefficient
+        """
+
+        # check type of values
+        D_0 = F.as_fenics_constant(self.D_0, mesh)
+        E_D = F.as_fenics_constant(self.E_D, mesh)
+
+        return D_0 * ufl.exp(-E_D / F.k_B / temperature)

festim/mesh/mesh.py

@@ -5,15 +5,14 @@ class Mesh:
     """
     Mesh class
 
+    Args:
+            mesh (dolfinx.mesh.Mesh, optional): the mesh. Defaults to None.
+
     Attributes:
         mesh (dolfinx.mesh.Mesh): the mesh
     """
 
     def __init__(self, mesh=None):
-        """Inits Mesh
-        Args:
-            mesh (dolfinx.mesh.Mesh, optional): the mesh. Defaults to None.
-        """
         self.mesh = mesh
 
         if self.mesh is not None:

festim/mesh/mesh_1d.py

@@ -9,17 +9,14 @@ class Mesh1D(Mesh):
     """
     1D Mesh
 
+    Args:
+            vertices (list): the mesh x-coordinates (m)
+
     Attributes:
         vertices (list): the mesh x-coordinates (m)
     """
 
     def __init__(self, vertices, **kwargs) -> None:
-        """Inits Mesh1D
-
-        Args:
-            vertices (list): the mesh x-coordinates (m)
-        """
-
         self.vertices = vertices
 
         mesh = self.generate_mesh()

festim/species.py

@@ -22,11 +22,6 @@ class Species:
     """
 
     def __init__(self, name: str = None) -> None:
-        """_summary_
-
-        Args:
-            name (str, optional): a name given to the species. Defaults to None.
-        """
         self.name = name
 
         self.solution = None

festim/subdomain/volume_subdomain.py

@@ -30,7 +30,7 @@ def locate_subdomain_entities(self, mesh, vdim):
         """Locates all cells in subdomain borders within domain
 
         Args:
-            mesh (dolfinx.cpp.mesh.Mesh): the mesh of the model
+            mesh (dolfinx.mesh.Mesh): the mesh of the model
             vdim (int): the dimension of the volumes of the mesh,
                 for 1D this is always 1
 

test/test_helpers.py

@@ -0,0 +1,21 @@
+import festim as F
+from dolfinx import fem
+import numpy as np
+import pytest
+import ufl
+
+test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
+x = ufl.SpatialCoordinate(test_mesh.mesh)
+
+
+@pytest.mark.parametrize(
+    "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]]
+)
+def test_temperature_type_and_processing(value):
+    """Test that the temperature type is correctly set"""
+
+    if not isinstance(value, (fem.Constant, int, float)):
+        with pytest.raises(TypeError):
+            F.as_fenics_constant(value, test_mesh.mesh)
+    else:
+        assert isinstance(F.as_fenics_constant(value, test_mesh.mesh), fem.Constant)

test/test_material.py

@@ -0,0 +1,16 @@
+import festim as F
+import numpy as np
+
+test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0]))
+
+
+def test_define_diffusion_coefficient():
+    """Test that the diffusion coefficient is correctly defined"""
+    T, D_0, E_D = 10, 1.2, 0.5
+
+    my_mat = F.Material(D_0=D_0, E_D=E_D)
+    D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T)
+
+    D_analytical = D_0 * np.exp(-E_D / F.k_B / T)
+
+    assert np.isclose(float(D), D_analytical)