diff --git a/festim/__init__.py b/festim/__init__.py index de4a2652c..0edcfb9ea 100644 --- a/festim/__init__.py +++ b/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 diff --git a/festim/helpers.py b/festim/helpers.py new file mode 100644 index 000000000..9b17fb92b --- /dev/null +++ b/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)}" + ) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e5dfd56b2..9898f4ba2 100644 --- a/festim/hydrogen_transport_problem.py +++ b/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,9 +18,10 @@ 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 @@ -28,15 +29,18 @@ 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 - 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""" diff --git a/festim/material.py b/festim/material.py new file mode 100644 index 000000000..494181e99 --- /dev/null +++ b/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) diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py index 90f7746bc..5f6be3dac 100644 --- a/festim/mesh/mesh.py +++ b/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: diff --git a/festim/mesh/mesh_1d.py b/festim/mesh/mesh_1d.py index 6a29e4c54..f7de96bee 100644 --- a/festim/mesh/mesh_1d.py +++ b/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() diff --git a/festim/species.py b/festim/species.py index 18ed7683d..eda76f2e1 100644 --- a/festim/species.py +++ b/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 diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py index a3398993c..dd3d1bf24 100644 --- a/festim/subdomain/volume_subdomain.py +++ b/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 diff --git a/test/test_helpers.py b/test/test_helpers.py new file mode 100644 index 000000000..14b269bda --- /dev/null +++ b/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) diff --git a/test/test_material.py b/test/test_material.py new file mode 100644 index 000000000..e706cfe88 --- /dev/null +++ b/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) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 87492dec4..b8534107e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -8,8 +8,7 @@ form, assemble_scalar, ) -from ufl import dot, grad, exp, FacetNormal, ds, Measure -from dolfinx import log +from ufl import dot, grad, exp, FacetNormal import numpy as np import tqdm.autonotebook @@ -26,7 +25,8 @@ def test_permeation_problem(): my_model = F.HydrogenTransportProblem() my_model.mesh = my_mesh - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=None) + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) left_surface = F.SurfaceSubdomain1D(id=1, x=0) right_surface = F.SurfaceSubdomain1D(id=2, x=L) my_model.subdomains = [my_subdomain, left_surface, right_surface] @@ -39,6 +39,8 @@ def test_permeation_problem(): my_model.initialise() + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) + V = my_model.function_space u = mobile_H.solution @@ -98,7 +100,7 @@ def siverts_law(T, S_0, E_S, pressure): mobile_xdmf.write_function(u, t) - surface_flux = form(my_model.D * dot(grad(u), n) * my_model.ds(2)) + surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) @@ -109,13 +111,12 @@ def siverts_law(T, S_0, E_S, pressure): # analytical solution S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability = float(my_model.D) * S + permeability = float(D) * S times = np.array(times) n_array = np.arange(1, 10000)[:, np.newaxis] summation = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(my_model.D) / L**2 * times), + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), axis=0, ) analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) @@ -133,5 +134,131 @@ def siverts_law(T, S_0, E_S, pressure): assert error < 0.01 -if __name__ == "__main__": - test_permeation_problem() +def test_permeation_problem_multi_volume(): + L = 3e-04 + vertices = np.linspace(0, L, num=1001) + + my_mesh = F.Mesh1D(vertices) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = my_mesh + + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, L / 4], material=my_mat) + my_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[L / 4, L / 2], material=my_mat) + my_subdomain_3 = F.VolumeSubdomain1D( + id=3, borders=[L / 2, 3 * L / 4], material=my_mat + ) + my_subdomain_4 = F.VolumeSubdomain1D(id=4, borders=[3 * L / 4, L], material=my_mat) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=L) + my_model.subdomains = [ + my_subdomain_1, + my_subdomain_2, + my_subdomain_3, + my_subdomain_4, + left_surface, + right_surface, + ] + + mobile_H = F.Species("H") + my_model.species = [mobile_H] + + temperature = Constant(my_mesh.mesh, 500.0) + my_model.temperature = temperature + + my_model.initialise() + + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) + + V = my_model.function_space + u = mobile_H.solution + + # 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 / F.k_B / T) + return S * pressure**0.5 + + fdim = my_mesh.mesh.topology.dim - 1 + left_facets = my_model.facet_meshtags.find(1) + left_dofs = locate_dofs_topological(V, fdim, left_facets) + right_facets = my_model.facet_meshtags.find(2) + right_dofs = locate_dofs_topological(V, fdim, right_facets) + + S_0 = 4.02e21 + E_S = 1.04 + P_up = 100 + surface_conc = siverts_law(T=temperature, S_0=S_0, E_S=E_S, pressure=P_up) + 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) + my_model.boundary_conditions = [bc_sieverts, bc_outgas] + my_model.create_solver() + + my_model.solver.convergence_criterion = "incremental" + my_model.solver.rtol = 1e-10 + my_model.solver.atol = 1e10 + + my_model.solver.report = True + ksp = my_model.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() + + mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") + mobile_xdmf.write_mesh(my_model.mesh.mesh) + + final_time = 50 + + flux_values = [] + times = [] + t = 0 + progress = tqdm.autonotebook.tqdm( + desc="Solving H transport problem", total=final_time + ) + while t < final_time: + progress.update(float(my_model.dt)) + t += float(my_model.dt) + + my_model.solver.solve(u) + + mobile_xdmf.write_function(u, t) + + surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(t) + + mobile_H.prev_solution.x.array[:] = u.x.array[:] + + mobile_xdmf.close() + + # analytical solution + S = S_0 * exp(-E_S / F.k_B / float(temperature)) + permeability = float(D) * S + times = np.array(times) + + n_array = np.arange(1, 10000)[:, np.newaxis] + summation = np.sum( + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), + axis=0, + ) + analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) + + analytical_flux = np.abs(analytical_flux) + flux_values = np.array(np.abs(flux_values)) + + relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) + + relative_error = relative_error[ + np.where(analytical_flux > 0.01 * np.max(analytical_flux)) + ] + error = relative_error.mean() + + assert error < 0.01 diff --git a/test/test_temperature.py b/test/test_temperature.py new file mode 100644 index 000000000..827a6960c --- /dev/null +++ b/test/test_temperature.py @@ -0,0 +1,24 @@ +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""" + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh + + if not isinstance(value, (fem.Constant, int, float)): + with pytest.raises(TypeError): + my_model.temperature = value + else: + my_model.temperature = value + assert isinstance(my_model.temperature, fem.Constant)