Species class

.github/workflows/ci.yml

@@ -21,10 +21,10 @@ jobs:
       run: |
         mamba install -c conda-forge fenics-dolfinx
 
-    - name: Install dependencies
+    - name: Install local package and dependencies
       shell: bash -l {0}
       run: |
-        pip install pytest pytest-cov
+        pip install .[test]
 
     - name: Run tests
       shell: bash -l {0}

festim/__init__.py

@@ -13,12 +13,12 @@
     __version__ = "unknown"
 
 
-import sympy as sp
-
 R = 8.314462618  # Gas constant J.mol-1.K-1
 k_B = 8.6173303e-5  # Boltzmann constant eV.K-1
 
 from .mesh.mesh import Mesh
 from .mesh.mesh_1d import Mesh1D
 
+from .species import Species, Trap
+
 from .hydrogen_transport_problem import HydrogenTransportProblem

festim/hydrogen_transport_problem.py

@@ -1,7 +1,9 @@
-import numpy as np
 from dolfinx import fem
 import ufl
 
+from dolfinx.fem import Function
+from ufl import TestFunction
+
 
 class HydrogenTransportProblem:
     """
@@ -10,9 +12,39 @@ class HydrogenTransportProblem:
     Args:
         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
 
     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 (dolfinx.Function): 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
+        facet_tags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model
+        volume_tags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model
+
+
+    Usage:
+        >>> import festim as F
+        >>> my_model = F.HydrogenTransportProblem()
+        >>> my_model.mesh = F.Mesh(...)
+        >>> my_model.subdomains = [F.Subdomain(...)]
+        >>> my_model.species = [F.Species(name="H"), F.Species(name="Trap")]
+        >>> my_model.initialise()
+
+        or
+
+        >>> my_model = F.HydrogenTransportProblem(
+        ...     mesh=F.Mesh(...),
+        ...     subdomains=[F.Subdomain(...)],
+        ...     species=[F.Species(name="H"), F.Species(name="Trap")],
+        ... )
+        >>> my_model.initialise()
 
     """
 
@@ -52,7 +84,15 @@ def initialise(self):
             self.dx,
             self.ds,
         ) = self.mesh.create_measures_and_tags(self.function_space)
+        self.assign_functions_to_species()
 
     def define_function_space(self):
         elements = ufl.FiniteElement("CG", self.mesh.mesh.ufl_cell(), 1)
         self.function_space = fem.FunctionSpace(self.mesh.mesh, elements)
+
+    def assign_functions_to_species(self):
+        """Creates for each species the solution, prev solution and test function"""
+        for spe in self.species:
+            spe.solution = Function(self.function_space)
+            spe.prev_solution = Function(self.function_space)
+            spe.test_function = TestFunction(self.function_space)

festim/species.py

@@ -0,0 +1,58 @@
+class Species:
+    """
+    Hydrogen species class for H transport simulation.
+
+    Args:
+        name (str, optional): a name given to the species. Defaults to None.
+
+    Attributes:
+        name (str): a name given to the species.
+        solution (dolfinx.fem.Function or ...): the solution for the current timestep
+        prev_solution (dolfinx.fem.Function or ...): the solution for the previous timestep
+        test_function (ufl.Argument or ...): the testfunction associated with this species
+
+    Usage:
+        >>> from festim import Species, HTransportProblem
+        >>> species = Species(name="H")
+        >>> species.name
+        'H'
+        >>> my_model = HTransportProblem()
+        >>> my_model.species.append(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
+        self.prev_solution = None
+        self.test_function = None
+
+
+class Trap(Species):
+    """Trap species class for H transport simulation.
+
+    Args:
+        name (str, optional): a name given to the trap. Defaults to None.
+
+    Attributes:
+        name (str): a name given to the trap.
+        attributes of Species class
+
+    Usage:
+        >>> from festim import Trap, HTransportProblem
+        >>> trap = Trap(name="Trap")
+        >>> trap.name
+        'Trap'
+        >>> my_model = HTransportProblem()
+        >>> my_model.species.append(trap)
+
+    """
+
+    def __init__(self, name: str = None) -> None:
+        super().__init__(name)

setup.cfg

@@ -27,7 +27,10 @@ project_urls =
 [options]
 packages = find:
 python_requires= >=3.6
+install_requires =
+    tqdm
 
 [options.extras_require]
-tests = 
+test = 
     pytest >= 5.4.3
+    pytest-cov

test/test_permeation_problem.py

@@ -0,0 +1,136 @@
+from mpi4py import MPI
+from petsc4py import PETSc
+from dolfinx.io import XDMFFile
+from dolfinx.fem import (
+    Constant,
+    dirichletbc,
+    locate_dofs_topological,
+    form,
+    assemble_scalar,
+)
+from dolfinx.fem.petsc import (
+    NonlinearProblem,
+)
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import (
+    dot,
+    grad,
+    exp,
+    FacetNormal,
+    dx,
+    ds,
+)
+from dolfinx import log
+import numpy as np
+import tqdm.autonotebook
+
+
+import festim as F
+
+
+def test_permeation_problem():
+    # mesh nodes
+    vertices = np.linspace(0, 3e-4, num=1001)
+
+    my_mesh = F.Mesh1D(vertices)
+
+    my_model = F.HydrogenTransportProblem()
+    my_model.mesh = my_mesh
+
+    mobile_H = F.Species("H")
+    my_model.species = [mobile_H]
+
+    my_model.initialise()
+
+    V = my_model.function_space
+    u = mobile_H.solution
+    u_n = mobile_H.prev_solution
+    v = mobile_H.test_function
+
+    temperature = Constant(my_mesh.mesh, 500.0)
+    k_B = F.k_B
+
+    # TODO this should be a property of Mesh
+    n = FacetNormal(my_mesh.mesh)
+
+    def siverts_law(T, S_0, E_S, pressure):
+        S = S_0 * exp(-E_S / k_B / T)
+        return S * pressure**0.5
+
+    fdim = my_mesh.mesh.topology.dim - 1
+    left_facets = my_model.facet_tags.find(1)
+    left_dofs = locate_dofs_topological(V, fdim, left_facets)
+    right_facets = my_model.facet_tags.find(2)
+    right_dofs = locate_dofs_topological(V, fdim, right_facets)
+
+    surface_conc = siverts_law(T=temperature, S_0=4.02e21, E_S=1.04, pressure=100)
+    bc_sieverts = dirichletbc(
+        Constant(my_mesh.mesh, PETSc.ScalarType(surface_conc)), left_dofs, V
+    )
+    bc_outgas = dirichletbc(Constant(my_mesh.mesh, PETSc.ScalarType(0)), right_dofs, V)
+    bcs = [bc_sieverts, bc_outgas]
+
+    D_0 = Constant(my_mesh.mesh, 1.9e-7)
+    E_D = Constant(my_mesh.mesh, 0.2)
+
+    D = D_0 * exp(-E_D / k_B / temperature)
+
+    dt = Constant(my_mesh.mesh, 1 / 20)
+    final_time = 50
+
+    # f = Constant(my_mesh.mesh, (PETSc.ScalarType(0)))
+    variational_form = dot(D * grad(u), grad(v)) * dx
+    variational_form += ((u - u_n) / dt) * v * dx
+
+    problem = NonlinearProblem(variational_form, u, bcs=bcs)
+    solver = NewtonSolver(MPI.COMM_WORLD, problem)
+
+    solver.convergence_criterion = "incremental"
+    solver.rtol = 1e-10
+    solver.atol = 1e10
+
+    solver.report = True
+    ksp = solver.krylov_solver
+    opts = PETSc.Options()
+    option_prefix = ksp.getOptionsPrefix()
+    opts[f"{option_prefix}ksp_type"] = "cg"
+    opts[f"{option_prefix}pc_type"] = "gamg"
+    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
+    ksp.setFromOptions()
+    # log.set_log_level(log.LogLevel.INFO)
+
+    mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w")
+    mobile_xdmf.write_mesh(my_mesh.mesh)
+
+    flux_values = []
+    times = []
+    t = 0
+    progress = tqdm.autonotebook.tqdm(
+        desc="Solving H transport problem", total=final_time
+    )
+    while t < final_time:
+        progress.update(float(dt))
+        t += float(dt)
+
+        solver.solve(u)
+
+        # post process
+        surface_flux = form(D * dot(grad(u), n) * ds(2))
+        flux = assemble_scalar(surface_flux)
+        flux_values.append(flux)
+        times.append(t)
+
+        # export
+        np.savetxt("outgassing_flux.txt", np.array(flux_values))
+        np.savetxt("times.txt", np.array(times))
+
+        mobile_xdmf.write_function(u, t)
+
+        # update previous solution
+        u_n.x.array[:] = u.x.array[:]
+
+    mobile_xdmf.close()
+
+
+if __name__ == "__main__":
+    test_permeation_problem()

test/test_species.py

@@ -0,0 +1,28 @@
+import festim as F
+import dolfinx
+import ufl
+import numpy as np
+
+
+def test_assign_functions_to_species():
+    """Test that checks if the function assign_functions_to_species
+    creates the correct attributes for each species
+    """
+
+    mesh = F.Mesh1D(
+        vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
+    )
+    model = F.HydrogenTransportProblem(
+        mesh=mesh,
+        species=[F.Species(name="H"), F.Species(name="Trap")],
+    )
+    model.define_function_space()
+    model.assign_functions_to_species()
+
+    for spe in model.species:
+        assert spe.solution is not None
+        assert spe.prev_solution is not None
+        assert spe.test_function is not None
+        assert isinstance(spe.solution, dolfinx.fem.Function)
+        assert isinstance(spe.prev_solution, dolfinx.fem.Function)
+        assert isinstance(spe.test_function, ufl.Argument)