diff --git a/festim/__init__.py b/festim/__init__.py index be0baf70c..502168be4 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -30,7 +30,7 @@ from .settings import Settings -from .species import Species, Trap, ImplicitSpecies +from .species import Species, Trap, ImplicitSpecies, find_species_from_name from .subdomain.surface_subdomain import SurfaceSubdomain1D from .subdomain.volume_subdomain import VolumeSubdomain1D diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index efb4d4dba..6b31ab22e 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -69,6 +69,7 @@ def define_surface_subdomain_dofs(self, facet_meshtags, mesh, 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( @@ -98,7 +99,7 @@ def create_value( 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) + mesh=mesh, value=self.value(t=float(t)) ) else: self.value_fenics = fem.Function(function_space) @@ -118,25 +119,6 @@ def create_value( ) 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 diff --git a/festim/exports/vtx.py b/festim/exports/vtx.py index 76d2d73a7..af0165435 100644 --- a/festim/exports/vtx.py +++ b/festim/exports/vtx.py @@ -48,16 +48,16 @@ def field(self): @field.setter def field(self, value): # check that field is festim.Species or list of festim.Species - if not isinstance(value, F.Species) and not isinstance(value, list): + if not isinstance(value, (F.Species, str)) and not isinstance(value, list): raise TypeError( - "field must be of type festim.Species or list of festim.Species" + "field must be of type festim.Species or str or a list of festim.Species or str" ) # check that all elements of list are festim.Species if isinstance(value, list): for element in value: - if not isinstance(element, F.Species): + if not isinstance(element, (F.Species, str)): raise TypeError( - "field must be of type festim.Species or list of festim.Species" + "field must be of type festim.Species or str or a list of festim.Species or str" ) # if field is festim.Species, convert to list if not isinstance(value, list): @@ -72,7 +72,10 @@ def define_writer(self, comm: mpi4py.MPI.Intracomm) -> None: comm (mpi4py.MPI.Intracomm): the MPI communicator """ self.writer = VTXWriter( - comm, self.filename, [field.solution for field in self.field], "BP4" + comm, + self.filename, + [field.post_processing_solution for field in self.field], + "BP4", ) def write(self, t: float): diff --git a/festim/exports/xdmf.py b/festim/exports/xdmf.py index 8318a1d78..83c8f6dce 100644 --- a/festim/exports/xdmf.py +++ b/festim/exports/xdmf.py @@ -39,16 +39,16 @@ def field(self): @field.setter def field(self, value): # check that field is festim.Species or list of festim.Species - if not isinstance(value, F.Species) and not isinstance(value, list): + if not isinstance(value, (F.Species, str)) and not isinstance(value, list): raise TypeError( - "field must be of type festim.Species or list of festim.Species" + "field must be of type festim.Species or str or a list of festim.Species or str" ) # check that all elements of list are festim.Species if isinstance(value, list): for element in value: - if not isinstance(element, F.Species): + if not isinstance(element, (F.Species, str)): raise TypeError( - "field must be of type festim.Species or list of festim.Species" + "field must be of type festim.Species or str or a list of festim.Species or str" ) # if field is festim.Species, convert to list if not isinstance(value, list): @@ -71,4 +71,4 @@ def write(self, t: float): t (float): the time of export """ for field in self.field: - self.writer.write_function(field.solution, t) + self.writer.write_function(field.post_processing_solution, t) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 63b3cf37c..ff171bca3 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -47,6 +47,7 @@ class HydrogenTransportProblem: model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model + multispecies (bool): True if the model has more than one species Usage: >>> import festim as F @@ -110,6 +111,10 @@ def temperature(self, value): else: self._temperature = F.as_fenics_constant(value, self.mesh.mesh) + @property + def multispecies(self): + return len(self.species) > 1 + def initialise(self): self.define_function_space() self.define_markers_and_measures() @@ -124,10 +129,13 @@ def initialise(self): self.defing_export_writers() def defing_export_writers(self): - """Defines the export writers of the model""" + """Defines the export writers of the model, if field is given as + a string, find species object in self.species""" for export in self.exports: - # TODO implement when export.field is an int or str - # then find solution from index of species + # if name of species is given then replace with species object + for idx, field in enumerate(export.field): + if isinstance(field, str): + export.field[idx] = F.find_species_from_name(field, self.species) if isinstance(export, (F.VTXExport, F.XDMFExport)): export.define_writer(MPI.COMM_WORLD) @@ -135,18 +143,62 @@ def defing_export_writers(self): export.writer.write_mesh(self.mesh.mesh) def define_function_space(self): - elements = basix.ufl.element( + """Creates the function space of the model, creates a mixed element if + model is multispecies. Creates the main solution and previous solution + function u and u_n.""" + element_CG = basix.ufl.element( basix.ElementFamily.P, self.mesh.mesh.basix_cell(), 1, basix.LagrangeVariant.equispaced, ) - self.function_space = fem.FunctionSpace(self.mesh.mesh, elements) + if not self.multispecies: + element = element_CG + else: + elements = [] + for spe in self.species: + if isinstance(spe, F.Species): + # TODO check if mobile or immobile for traps + elements.append(element_CG) + element = ufl.MixedElement(elements) + + self.function_space = fem.FunctionSpace(self.mesh.mesh, element) + + self.u = Function(self.function_space) + self.u_n = Function(self.function_space) + + def assign_functions_to_species(self): + """Creates the solution, prev solution, test function and + post-processing solution for each species, if model is multispecies, + created a collapsed function space for each species""" + + if not self.multispecies: + sub_solutions = [self.u] + sub_prev_solution = [self.u_n] + sub_test_functions = [ufl.TestFunction(self.function_space)] + self.species[0].sub_function_space = self.function_space + self.species[0].post_processing_solution = fem.Function(self.function_space) + else: + sub_solutions = list(ufl.split(self.u)) + sub_prev_solution = list(ufl.split(self.u_n)) + sub_test_functions = list(ufl.TestFunctions(self.function_space)) + + for idx, spe in enumerate(self.species): + spe.sub_function_space = self.function_space.sub(idx) + spe.post_processing_solution = self.u.sub(idx).collapse() + spe.collapsed_function_space, _ = self.function_space.sub( + idx + ).collapse() + + for idx, spe in enumerate(self.species): + spe.solution = sub_solutions[idx] + spe.prev_solution = sub_prev_solution[idx] + spe.test_function = sub_test_functions[idx] def define_markers_and_measures(self): """Defines the markers and measures of the model""" - dofs_facets, tags_facets = [], [] + facet_indices, tags_facets = [], [] # find all cells in domain and mark them as 0 num_cells = self.mesh.mesh.topology.index_map(self.mesh.vdim).size_local @@ -155,8 +207,10 @@ def define_markers_and_measures(self): for sub_dom in self.subdomains: if isinstance(sub_dom, F.SurfaceSubdomain1D): - dof = sub_dom.locate_dof(self.function_space) - dofs_facets.append(dof) + facet_index = sub_dom.locate_boundary_facet_indices( + self.mesh.mesh, self.mesh.fdim + ) + facet_indices.append(facet_index) tags_facets.append(sub_dom.id) if isinstance(sub_dom, F.VolumeSubdomain1D): # find all cells in subdomain and mark them as sub_dom.id @@ -171,12 +225,12 @@ def define_markers_and_measures(self): self.mesh.check_borders(self.volume_subdomains) # dofs and tags need to be in np.in32 format for meshtags - dofs_facets = np.array(dofs_facets, dtype=np.int32) + facet_indices = np.array(facet_indices, dtype=np.int32) tags_facets = np.array(tags_facets, dtype=np.int32) # define mesh tags self.facet_meshtags = meshtags( - self.mesh.mesh, self.mesh.fdim, dofs_facets, tags_facets + self.mesh.mesh, self.mesh.fdim, facet_indices, tags_facets ) self.volume_meshtags = meshtags( self.mesh.mesh, self.mesh.vdim, mesh_cell_indices, tags_volumes @@ -193,37 +247,74 @@ def define_markers_and_measures(self): def define_boundary_conditions(self): """Defines the dirichlet boundary conditions of the model""" for bc in self.boundary_conditions: + if isinstance(bc.species, str): + # if name of species is given then replace with species object + bc.species = F.find_species_from_name(bc.species, self.species) 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 - ) + form = self.create_dirichletbc_form(bc) 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: - raise NotImplementedError("Multiple species not implemented yet") - 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) + def create_dirichletbc_form(self, bc): + """Creates a dirichlet boundary condition form + + Args: + bc (festim.DirichletBC): the boundary condition + + Returns: + dolfinx.fem.bcs.DirichletBC: A representation of + the boundary condition for modifying linear systems. + """ + # create value_fenics + function_space_value = None + + if callable(bc.value): + # if bc.value is a callable then need to provide a functionspace + if not self.multispecies: + function_space_value = bc.species.sub_function_space + else: + function_space_value = bc.species.collapsed_function_space + + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=function_space_value, + t=self.t, + ) - # TODO remove this - self.u = self.species[0].solution - self.u_n = self.species[0].prev_solution + # get dofs + if self.multispecies and isinstance(bc.value_fenics, (fem.Function)): + function_space_dofs = ( + bc.species.sub_function_space, + bc.species.collapsed_function_space, + ) + else: + function_space_dofs = bc.species.sub_function_space + + bc_dofs = bc.define_surface_subdomain_dofs( + facet_meshtags=self.facet_meshtags, + mesh=self.mesh, + function_space=function_space_dofs, + ) + + # create form + if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)): + # no need to pass the functionspace since value_fenics is already a Function + function_space_form = None + else: + function_space_form = bc.species.sub_function_space + + form = fem.dirichletbc( + value=bc.value_fenics, + dofs=bc_dofs, + V=function_space_form, + ) + + return form def create_formulation(self): """Creates the formulation of the model""" if len(self.sources) > 1: raise NotImplementedError("Sources not implemented yet") - if len(self.species) > 1: - raise NotImplementedError("Multiple species not implemented yet") self.formulation = 0 @@ -234,7 +325,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature + self.mesh.mesh, self.temperature, spe ) self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) @@ -257,7 +348,7 @@ def create_solver(self): """Creates the solver of the model""" problem = fem.petsc.NonlinearProblem( self.formulation, - self.species[0].solution, + self.u, bcs=self.bc_forms, ) self.solver = NewtonSolver(MPI.COMM_WORLD, problem) @@ -273,12 +364,8 @@ def run(self): list of float: the fluxes of the simulation """ times, flux_values = [], [] + flux_values_1, flux_values_2 = [], [] - n = self.mesh.n - D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature - ) - cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( desc="Solving H transport problem", total=self.settings.final_time, @@ -295,12 +382,36 @@ def run(self): self.solver.solve(self.u) # post processing - - surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) - - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(float(self.t)) + # TODO remove this + if not self.multispecies: + D_D = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + cm = self.u + self.species[0].post_processing_solution = self.u + + surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(float(self.t)) + else: + for idx, spe in enumerate(self.species): + spe.post_processing_solution = self.u.sub(idx) + + cm_1, cm_2 = self.u.split() + D_1 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + D_2 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[1] + ) + surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) + surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) + flux_1 = assemble_scalar(surface_flux_1) + flux_2 = assemble_scalar(surface_flux_2) + flux_values_1.append(flux_1) + flux_values_2.append(flux_2) + times.append(float(self.t)) for export in self.exports: if isinstance(export, (F.VTXExport, F.XDMFExport)): @@ -309,4 +420,7 @@ def run(self): # update previous solution self.u_n.x.array[:] = self.u.x.array[:] + if self.multispecies: + flux_values = [flux_values_1, flux_values_2] + return times, flux_values diff --git a/festim/material.py b/festim/material.py index 494181e99..ceeced3ff 100644 --- a/festim/material.py +++ b/festim/material.py @@ -7,18 +7,27 @@ class Material: Material class Args: - D_0 (float or fem.Constant): the pre-exponential factor of the + D_0 (float or dict): the pre-exponential factor of the diffusion coefficient (m2/s) - E_D (float or fem.Constant): the activation energy of the diffusion + E_D (float or dict): 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 + D_0 (float or dict): the pre-exponential factor of the diffusion coefficient (m2/s) - E_D (float or fem.Constant): the activation energy of the diffusion + E_D (float or dict): the activation energy of the diffusion coeficient (eV) name (str): the name of the material + + Usage: + >>> my_mat = Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + or if several species: + >>> my_mat = Material( + D_0={"Species_1": 1.9e-7, "Species_2": 2.0e-7}, + E_D={"Species_1": 0.2, "Species_2": 0.3}, + name="my_mat" + ) """ def __init__(self, D_0, E_D, name=None) -> None: @@ -26,17 +35,46 @@ def __init__(self, D_0, E_D, name=None) -> None: self.E_D = E_D self.name = name - def get_diffusion_coefficient(self, mesh, temperature): + def get_diffusion_coefficient(self, mesh, temperature, species=None): """Defines the diffusion coefficient Args: + mesh (dolfinx.mesh.Mesh): the domain mesh temperature (dolfinx.fem.Constant): the temperature + species (festim.Species, optional): the species we want the diffusion + coefficient of. Only needed if D_0 and E_D are dicts. 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) + if isinstance(self.D_0, (float, int)) and isinstance(self.E_D, (float, int)): + 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) + + elif isinstance(self.D_0, dict) and isinstance(self.E_D, dict): + # check D_0 and E_D have the same keys + # this check should go in a setter + if list(self.D_0.keys()) != list(self.E_D.keys()): + raise ValueError("D_0 and E_D have different keys") + + if species is None: + raise ValueError("species must be provided if D_0 and E_D are dicts") + + if species in self.D_0: + D_0 = F.as_fenics_constant(self.D_0[species], mesh) + elif species.name in self.D_0: + D_0 = F.as_fenics_constant(self.D_0[species.name], mesh) + else: + raise ValueError(f"{species} is not in D_0 keys") + + if species in self.E_D: + E_D = F.as_fenics_constant(self.E_D[species], mesh) + elif species.name in self.E_D: + E_D = F.as_fenics_constant(self.E_D[species.name], mesh) + + return D_0 * ufl.exp(-E_D / F.k_B / temperature) - return D_0 * ufl.exp(-E_D / F.k_B / temperature) + else: + raise ValueError("D_0 and E_D must be either floats or dicts") diff --git a/festim/species.py b/festim/species.py index a3f45eacd..663436172 100644 --- a/festim/species.py +++ b/festim/species.py @@ -10,9 +10,18 @@ class Species: 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 + solution (dolfinx.fem.Function): the solution for the current timestep + prev_solution (dolfinx.fem.Function): the solution for the previous + timestep + test_function (ufl.Argument): the testfunction associated with this + species + sub_function_space (dolfinx.fem.function.FunctionSpaceBase): the + subspace of the function space + collapsed_function_space (dolfinx.fem.function.FunctionSpaceBase): the + collapsed function space for a species in the function space. In + case single species case, this is None. + post_processing_solution (dolfinx.fem.Function): the solution for post + processing concentration (dolfinx.fem.Function): the concentration of the species Usage: @@ -31,6 +40,9 @@ def __init__(self, name: str = None) -> None: self.solution = None self.prev_solution = None self.test_function = None + self.sub_function_space = None + self.post_processing_solution = None + self.collapsed_function_space = None def __repr__(self) -> str: return f"Species({self.name})" @@ -111,3 +123,24 @@ def concentration(self): f"Cannot compute concentration of {self.name} because {other.name} has no solution" ) return self.n - sum([other.solution for other in self.others]) + + +def find_species_from_name(name: str, species: list): + """Returns the correct species object from a list of species + based on a string + + Args: + name (str): the name of the species + species (list): the list of species + + Returns: + species (festim.Species): the species object with the correct name + + Raises: + ValueError: if the species name is not found in the list of species + + """ + for spe in species: + if spe.name == name: + return spe + raise ValueError(f"Species {name} not found in list of species") diff --git a/festim/subdomain/surface_subdomain.py b/festim/subdomain/surface_subdomain.py index 409ca6671..80fd082aa 100644 --- a/festim/subdomain/surface_subdomain.py +++ b/festim/subdomain/surface_subdomain.py @@ -1,4 +1,4 @@ -from dolfinx import fem +import dolfinx.mesh import numpy as np @@ -22,18 +22,19 @@ def __init__(self, id, x) -> None: self.id = id self.x = x - def locate_dof(self, function_space): + def locate_boundary_facet_indices(self, mesh, fdim): """Locates the dof of the surface subdomain within the function space + and return the index of the dof Args: - function_space (dolfinx.fem.FunctionSpace): the function space of - the model + mesh (dolfinx.mesh.Mesh): the mesh of the simulation + fdim (int): the dimension of the model facets Returns: - dof (np.array): the first value in the list of dofs of the surface - subdomain + index (np.array): the first value in the list of surface facet + indices of the subdomain """ - dofs = fem.locate_dofs_geometrical( - function_space, lambda x: np.isclose(x[0], self.x) + indices = dolfinx.mesh.locate_entities_boundary( + mesh, fdim, lambda x: np.isclose(x[0], self.x) ) - return dofs[0] + return indices[0] diff --git a/test/benchmark.py b/test/benchmark.py index e86360126..949fe066a 100644 --- a/test/benchmark.py +++ b/test/benchmark.py @@ -170,11 +170,12 @@ def siverts_law(T, S_0, E_S, pressure): analytical_flux = np.abs(analytical_flux) flux_values = np.array(np.abs(flux_values)) + indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) + analytical_flux = analytical_flux[indices] + flux_values = flux_values[indices] + 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_dirichlet_bc.py b/test/test_dirichlet_bc.py index 48e7ac36f..f7f5b0e05 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -57,13 +57,12 @@ def test_callable_for_value(): subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) value = lambda x, t: 1.0 + x[0] + t - species = "test" + species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) my_model = F.HydrogenTransportProblem( - mesh=F.Mesh(mesh), - subdomains=[subdomain, vol_subdomain], + mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species] ) my_model.define_function_space() @@ -96,13 +95,12 @@ def test_value_callable_x_t_T(): subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) value = lambda x, t, T: 1.0 + x[0] + t + T - species = "test" + species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) my_model = F.HydrogenTransportProblem( - mesh=F.Mesh(mesh), - subdomains=[subdomain, vol_subdomain], + mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species] ) my_model.define_function_space() @@ -132,19 +130,20 @@ def test_value_callable_x_t_T(): assert np.isclose(computed_value, expected_value) -def test_callable_t_only(): +@pytest.mark.parametrize("value", [lambda t: t, lambda t: 1.0 + t]) +def test_callable_t_only(value): """Test that the value attribute can be a callable function of t only""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) - value = lambda t: 1.0 + t - species = "test" + species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) my_model = F.HydrogenTransportProblem( mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], + species=[species], ) my_model.define_function_space() @@ -180,13 +179,14 @@ def test_callable_x_only(): subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) value = lambda x: 1.0 + x[0] - species = "test" + species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) my_model = F.HydrogenTransportProblem( mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], + species=[species], ) my_model.define_function_space() @@ -220,46 +220,7 @@ def test_callable_x_only(): "value", [ 1.0, - lambda x: 1.0 + x[0], - lambda x, t: 1.0 + x[0] + t, - lambda x, t, T: 1.0 + x[0] + t + T, - ], -) -def test_create_formulation(value): - """A test that checks that the method create_formulation can be called when value is either a callable or a float""" - # BUILD - subdomain = F.SurfaceSubdomain1D(1, x=1) - vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) - species = "test" - - bc = F.DirichletBC(subdomain, value, species) - - my_model = F.HydrogenTransportProblem( - mesh=F.Mesh(mesh), - subdomains=[subdomain, vol_subdomain], - ) - - my_model.define_function_space() - my_model.define_markers_and_measures() - - T = fem.Constant(my_model.mesh.mesh, 550.0) - t = fem.Constant(my_model.mesh.mesh, 0.0) - - dofs = bc.define_surface_subdomain_dofs( - my_model.facet_meshtags, my_model.mesh, my_model.function_space - ) - bc.create_value(my_model.mesh.mesh, my_model.function_space, T, t) - - # TEST - formulation = bc.create_formulation(dofs, my_model.function_space) - - assert isinstance(formulation, fem.DirichletBC) - - -@pytest.mark.parametrize( - "value", - [ - 1.0, + lambda t: t, lambda t: 1.0 + t, lambda x: 1.0 + x[0], lambda x, t: 1.0 + x[0] + t, @@ -268,6 +229,8 @@ def test_create_formulation(value): ], ) def test_integration_with_HTransportProblem(value): + """test that different callable functions can be applied to a dirichlet + boundary condition, asserting in each case they match an expected value""" subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -276,7 +239,7 @@ def test_integration_with_HTransportProblem(value): subdomains=[vol_subdomain, subdomain], ) my_model.species = [F.Species("H")] - my_bc = F.DirichletBC(subdomain, value, my_model.species[0]) + my_bc = F.DirichletBC(subdomain, value, "H") my_model.boundary_conditions = [my_bc] my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0) @@ -318,3 +281,98 @@ def test_integration_with_HTransportProblem(value): if isinstance(expected_value, Conditional): expected_value = float(expected_value) assert np.isclose(computed_value, expected_value) + + +def test_species_predefined(): + """Test a ValueError is raised when the species defined in the boundary + condition is not predefined in the model""" + + subdomain = F.SurfaceSubdomain1D(1, x=1) + vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh(mesh), + subdomains=[vol_subdomain, subdomain], + ) + my_model.species = [F.Species("H")] + my_bc = F.DirichletBC(subdomain, 1.0, "J") + my_model.boundary_conditions = [my_bc] + my_model.settings = F.Settings(atol=1, rtol=0.1) + my_model.settings.stepsize = 1 + + with pytest.raises(ValueError): + my_model.initialise() + + +@pytest.mark.parametrize( + "value_A, value_B", + [ + (1.0, 1.0), + (1.0, lambda t: t), + (1.0, lambda t: 1.0 + t), + (1.0, lambda x: 1.0 + x[0]), + (1.0, lambda x, t: 1.0 + x[0] + t), + (1.0, lambda x, t, T: 1.0 + x[0] + t + T), + (1.0, lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0)), + ], +) +def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B): + """test that a mixture of callable functions can be applied to dirichlet + boundary conditions in a multispecies case, asserting in each case they + match an expected value""" + subdomain_A = F.SurfaceSubdomain1D(1, x=0) + subdomain_B = F.SurfaceSubdomain1D(2, x=1) + vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh(mesh), + subdomains=[vol_subdomain, subdomain_A, subdomain_B], + ) + my_model.species = [F.Species("A"), F.Species("B")] + my_bc_A = F.DirichletBC(subdomain_A, value_A, "A") + my_bc_B = F.DirichletBC(subdomain_B, value_B, "B") + my_model.boundary_conditions = [my_bc_A, my_bc_B] + + my_model.temperature = fem.Constant(my_model.mesh.mesh, 550.0) + + my_model.settings = F.Settings(atol=1, rtol=0.1, final_time=2) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + # RUN + + my_model.initialise() + + for bc in [my_bc_A, my_bc_B]: + assert bc.value_fenics is not None + + my_model.run() + + # TEST + + expected_value = value_A + computed_value = float(my_bc_A.value_fenics) + + if isinstance(value_B, float): + expected_value = value_B + computed_value = float(my_bc_B.value_fenics) + elif callable(value_B): + arguments = value_B.__code__.co_varnames + if "x" in arguments and "t" in arguments and "T" in arguments: + expected_value = value_B(x=np.array([subdomain_B.x]), t=2.0, T=550.0) + computed_value = my_bc_B.value_fenics.vector.array[-1] + elif "x" in arguments and "t" in arguments: + expected_value = value_B(x=np.array([subdomain_B.x]), t=2.0) + computed_value = my_bc_B.value_fenics.vector.array[-1] + elif "x" in arguments: + expected_value = value_B(x=np.array([subdomain_B.x])) + computed_value = my_bc_B.value_fenics.vector.array[-1] + elif "t" in arguments: + expected_value = value_B(t=2.0) + computed_value = float(my_bc_B.value_fenics) + else: + # test fails if lambda function is not recognised + raise ValueError("value function not recognised") + + if isinstance(expected_value, Conditional): + expected_value = float(expected_value) + assert np.isclose(computed_value, expected_value) diff --git a/test/test_material.py b/test/test_material.py index e706cfe88..265bd9afb 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -1,5 +1,6 @@ import festim as F import numpy as np +import pytest test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) @@ -7,10 +8,123 @@ def test_define_diffusion_coefficient(): """Test that the diffusion coefficient is correctly defined""" T, D_0, E_D = 10, 1.2, 0.5 - + dum_spe = F.Species("dummy") my_mat = F.Material(D_0=D_0, E_D=E_D) - D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T) + D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=dum_spe) D_analytical = D_0 * np.exp(-E_D / F.k_B / T) assert np.isclose(float(D), D_analytical) + + +def test_multispecies_dict_strings(): + """Test that the diffusion coefficient is correctly defined when keys are + strings""" + T = 500 + D_0_A, D_0_B = 1, 2 + E_D_A, E_D_B = 0.1, 0.2 + A, B = F.Species("A"), F.Species("B") + + my_mat = F.Material(D_0={"A": D_0_A, "B": D_0_B}, E_D={"A": E_D_A, "B": E_D_B}) + D_A = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=A) + D_B = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=B) + + D = [float(D_A), float(D_B)] + + D_A_analytical = D_0_A * np.exp(-E_D_A / F.k_B / T) + D_B_analytical = D_0_B * np.exp(-E_D_B / F.k_B / T) + + D_analytical = [D_A_analytical, D_B_analytical] + + assert np.isclose(D, D_analytical).all() + + +def test_multispecies_dict_objects(): + """Test that the diffusion coefficient is correctly defined when keys are + festim.Species objects""" + T = 500 + D_0_A, D_0_B = 1, 2 + E_D_A, E_D_B = 0.1, 0.2 + + A = F.Species("A") + B = F.Species("B") + + my_mat = F.Material(D_0={A: D_0_A, B: D_0_B}, E_D={A: E_D_A, B: E_D_B}) + D_A = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=A) + D_B = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=B) + D = [float(D_A), float(D_B)] + + D_A_analytical = D_0_A * np.exp(-E_D_A / F.k_B / T) + D_B_analytical = D_0_B * np.exp(-E_D_B / F.k_B / T) + + D_analytical = [D_A_analytical, D_B_analytical] + + assert np.isclose(D, D_analytical).all() + + +def test_multispecies_dict_objects_and_strings(): + """Test that the diffusion coefficient is correctly defined when keys + are a mix of festim.Species objects and strings""" + T = 500 + D_0_A, D_0_B = 1, 2 + E_D_A, E_D_B = 0.1, 0.2 + + A = F.Species("A") + B = F.Species("B") + + my_mat = F.Material(D_0={A: D_0_A, "B": D_0_B}, E_D={A: E_D_A, "B": E_D_B}) + D_A = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=A) + D_B = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species=B) + D = [float(D_A), float(D_B)] + + D_A_analytical = D_0_A * np.exp(-E_D_A / F.k_B / T) + D_B_analytical = D_0_B * np.exp(-E_D_B / F.k_B / T) + + D_analytical = [D_A_analytical, D_B_analytical] + + assert np.isclose(D, D_analytical).all() + + +def test_multispecies_dict_different_keys(): + """Test that a value error is raised when the keys of the D_0 and E_D + are not the same""" + A = F.Species("A") + my_mat = F.Material(D_0={"A": 1, "B": 2}, E_D={"A": 0.1, "B": 0.2, "C": 0.3}) + + with pytest.raises(ValueError, match="D_0 and E_D have different keys"): + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) + + +def test_D_0_type_raises_error(): + """Test that a value error is raised in the get_diffusion_coefficient + function""" + A = F.Species("A") + my_mat = F.Material(D_0=[1, 1], E_D=0.1) + + with pytest.raises(ValueError, match="D_0 and E_D must be either floats or dicts"): + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) + + +def test_error_raised_when_species_not_given_with_dict(): + """Test that a value error is raised when a species has not been given in + the get_diffusion_coefficient function when using a dict for properties""" + A = F.Species("A") + B = F.Species("B") + my_mat = F.Material(D_0={A: 1, B: 1}, E_D={A: 0.1, B: 0.1}) + + with pytest.raises( + ValueError, match="species must be provided if D_0 and E_D are dicts" + ): + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500) + + +def test_error_raised_when_species_not_not_in_D_0_dict(): + """Test that a value error is raised when a species has not been given but + has no value in the dict""" + A = F.Species("A") + B = F.Species("B") + J = F.Species("J") + my_mat = F.Material(D_0={A: 1, B: 1}, E_D={A: 0.1, B: 0.1}) + + with pytest.raises(ValueError, match="J is not in D_0 keys"): + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=J) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py new file mode 100644 index 000000000..ebae0dcf5 --- /dev/null +++ b/test/test_multispecies_problem.py @@ -0,0 +1,135 @@ +import numpy as np +import festim as F +from petsc4py import PETSc +from dolfinx.fem import Constant +from ufl import exp + + +def relative_error_computed_to_analytical( + D, permeability, computed_flux, L, times, P_up +): + 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) + + # post processing + analytical_flux = np.abs(analytical_flux) + indices = np.where(analytical_flux > 0.1 * np.max(analytical_flux)) + analytical_flux = analytical_flux[indices] + computed_flux = computed_flux[indices] + + # evaulate relative error compared to analytical solution + relative_error = np.abs((computed_flux - analytical_flux) / analytical_flux) + error = relative_error.mean() + + return error + + +def test_multispecies_permeation_problem(): + """Test running a problem with 2 mobile species permeating through a 1D + domain, with different diffusion coefficients, checks that the computed + permeation flux matches the analytical solution""" + + # festim model + L = 3e-04 + my_mesh = F.Mesh1D(np.linspace(0, L, num=1001)) + my_model = F.HydrogenTransportProblem() + my_model.mesh = my_mesh + + my_mat = F.Material( + D_0={"spe_1": 1.9e-7, "spe_2": 3.8e-7}, + E_D={"spe_1": 0.2, "spe_2": 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, + ] + + spe_1 = F.Species("spe_1") + spe_2 = F.Species("spe_2") + my_model.species = [spe_1, spe_2] + + temperature = Constant(my_mesh.mesh, 500.0) + my_model.temperature = temperature + pressure = 100 + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=right_surface, value=0, species="spe_1"), + F.DirichletBC(subdomain=right_surface, value=0, species="spe_2"), + F.SievertsBC( + subdomain=left_surface, + S_0=4.02e21, + E_S=1.04, + pressure=pressure, + species="spe_1", + ), + F.SievertsBC( + subdomain=left_surface, + S_0=5.0e21, + E_S=1.2, + pressure=pressure, + species="spe_2", + ), + ] + my_model.settings = F.Settings( + atol=1e10, + rtol=1e-10, + max_iterations=30, + final_time=10, + ) + + my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) + my_model.initialise() + + my_model.solver.convergence_criterion = "incremental" + 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() + + times, flux_values = my_model.run() + + # ---------------------- analytical solutions ----------------------------- + + # common values + times = np.array(times) + P_up = float(my_model.boundary_conditions[-1].pressure) + flux_values_spe_1, flux_values_spe_2 = flux_values[0], flux_values[1] + + # ##### compute analyical solution for species 1 ##### # + D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) + S_0_spe_1 = float(my_model.boundary_conditions[-2].S_0) + E_S_spe_1 = float(my_model.boundary_conditions[-2].E_S) + S_spe_1 = S_0_spe_1 * exp(-E_S_spe_1 / F.k_B / float(temperature)) + permeability_spe_1 = float(D_spe_1) * S_spe_1 + flux_values_spe_1 = np.array(np.abs(flux_values_spe_1)) + + error_spe_1 = relative_error_computed_to_analytical( + D_spe_1, permeability_spe_1, flux_values_spe_1, L, times, P_up + ) + + # ##### compute analyical solution for species 2 ##### # + D_spe_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_2) + S_0_spe_2 = float(my_model.boundary_conditions[-1].S_0) + E_S_spe_2 = float(my_model.boundary_conditions[-1].E_S) + S_spe_2 = S_0_spe_2 * exp(-E_S_spe_2 / F.k_B / float(temperature)) + permeability_spe_2 = float(D_spe_2) * S_spe_2 + flux_values_spe_2 = np.array(np.abs(flux_values_spe_2)) + + error_spe_2 = relative_error_computed_to_analytical( + D_spe_2, permeability_spe_2, flux_values_spe_2, L, times, P_up + ) + + for err in [error_spe_1, error_spe_2]: + assert err < 0.01 diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 2911cf753..8e2eb014a 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -2,12 +2,39 @@ from dolfinx.fem import Constant from ufl import exp import numpy as np +import festim as F +import os -import festim as F +def relative_error_computed_to_analytical( + D, permeability, computed_flux, L, times, P_up +): + 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) + + # post processing + analytical_flux = np.abs(analytical_flux) + indices = np.where(analytical_flux > 0.1 * np.max(analytical_flux)) + analytical_flux = analytical_flux[indices] + computed_flux = computed_flux[indices] + + # evaulate relative error compared to analytical solution + relative_error = np.abs((computed_flux - analytical_flux) / analytical_flux) + error = relative_error.mean() + + return error def test_permeation_problem(mesh_size=1001): + """Test running a problem with a mobile species permeating through a 1D + domain, checks that the computed permeation flux matches the analytical + solution""" + + # festim model L = 3e-04 vertices = np.linspace(0, L, num=mesh_size) @@ -67,29 +94,19 @@ def test_permeation_problem(mesh_size=1001): S = S_0 * exp(-E_S / F.k_B / float(my_model.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)) - indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - analytical_flux = analytical_flux[indices] - flux_values = flux_values[indices] - - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - - error = relative_error.mean() + error = relative_error_computed_to_analytical( + D, permeability, flux_values, L, times, P_up + ) assert error < 0.01 -def test_permeation_problem_multi_volume(): +def test_permeation_problem_multi_volume(tmp_path): + """Same permeation problem as above but with 4 volume subdomains instead + of 1""" + L = 3e-04 vertices = np.linspace(0, L, num=1001) @@ -128,7 +145,9 @@ def test_permeation_problem_multi_volume(): subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H" ), ] - my_model.exports = [F.VTXExport("test.bp", field=mobile_H)] + my_model.exports = [ + F.VTXExport(os.path.join(tmp_path, "mobile_concentration_h.bp"), field=mobile_H) + ] my_model.settings = F.Settings( atol=1e10, @@ -152,7 +171,7 @@ def test_permeation_problem_multi_volume(): times, flux_values = my_model.run() - # -------------------------- analytical solution ------------------------------------- + # ---------------------- analytical solution ----------------------------- D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) S_0 = float(my_model.boundary_conditions[-1].S_0) @@ -161,23 +180,10 @@ def test_permeation_problem_multi_volume(): 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)) - indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - analytical_flux = analytical_flux[indices] - flux_values = flux_values[indices] - - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - - error = relative_error.mean() + error = relative_error_computed_to_analytical( + D, permeability, flux_values, L, times, P_up + ) assert error < 0.01 diff --git a/test/test_subdomains.py b/test/test_subdomains.py index da681bc03..e8d5f004b 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -6,6 +6,7 @@ def test_different_surface_ids(): """Checks that different surface ids are correctly set""" my_test_model = F.HydrogenTransportProblem() + my_test_model.species = [F.Species("H")] L = 1e-04 my_test_model.mesh = F.Mesh1D(np.linspace(0, L, num=3)) diff --git a/test/test_vtx.py b/test/test_vtx.py index 15c7ded6e..1f0823740 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -13,7 +13,7 @@ def test_vtx_export_one_function(tmpdir): """Test can add one function to a vtx export""" u = dolfinx.fem.Function(V) sp = F.Species("H") - sp.solution = u + sp.post_processing_solution = u filename = str(tmpdir.join("my_export.bp")) my_export = F.VTXExport(filename, field=sp) my_export.define_writer(mesh.comm) @@ -31,8 +31,8 @@ def test_vtx_export_two_functions(tmpdir): sp1 = F.Species("1") sp2 = F.Species("2") - sp1.solution = u - sp2.solution = v + sp1.post_processing_solution = u + sp2.post_processing_solution = v filename = str(tmpdir.join("my_export.bp")) my_export = F.VTXExport(filename, field=[sp1, sp2]) @@ -78,7 +78,7 @@ def test_field_attribute_is_always_list(): assert isinstance(my_export.field, list) -@pytest.mark.parametrize("field", ["H", 1, [F.Species("H"), 1]]) +@pytest.mark.parametrize("field", [["H", 1], 1, [F.Species("H"), 1]]) def test_field_attribute_raises_error_when_invalid_type(field): """Test that the field attribute raises an error if the type is not festim.Species or list""" with pytest.raises(TypeError): @@ -95,3 +95,23 @@ def test_filename_raises_error_when_wrong_type(): """Test that the filename attribute raises an error if the extension is not .bp""" with pytest.raises(TypeError): F.VTXExport(1, field=[F.Species("H")]) + + +def test_vtx_field_as_string_found_in_species(tmpdir): + """Test that the field attribute can be a string and is found in the species list""" + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) + my_mat = F.Material(D_0=1, E_D=0, name="mat") + my_model.subdomains = [ + F.VolumeSubdomain1D(1, borders=[0.0, 4.0], material=my_mat), + ] + my_model.species = [F.Species("H")] + my_model.temperature = 500 + + filename = str(tmpdir.join("my_export.bp")) + my_export = F.VTXExport(filename, field="H") + my_model.exports = [my_export] + my_model.settings = F.Settings(atol=1, rtol=0.1) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + my_model.initialise() diff --git a/test/test_xdmf.py b/test/test_xdmf.py index cc1e57a25..fd1dabc3f 100644 --- a/test/test_xdmf.py +++ b/test/test_xdmf.py @@ -27,7 +27,7 @@ def test_write(tmp_path): V = fem.FunctionSpace(domain, ("Lagrange", 1)) u = fem.Function(V) - species.solution = u + species.post_processing_solution = u for t in [0, 1, 2, 3]: my_export.write(t=t) @@ -65,7 +65,7 @@ def test_field_attribute_is_always_list(): assert isinstance(my_export.field, list) -@pytest.mark.parametrize("field", ["H", 1, [F.Species("H"), 1]]) +@pytest.mark.parametrize("field", [["H", 2], 1, [F.Species("H"), 1]]) def test_field_attribute_raises_error_when_invalid_type(field): """Test that the field attribute raises an error if the type is not festim.Species or list""" with pytest.raises(TypeError): @@ -82,3 +82,23 @@ def test_filename_raises_error_when_wrong_type(): """Test that the filename attribute raises an error if the file is not str""" with pytest.raises(TypeError): F.XDMFExport(1, field=[F.Species("H")]) + + +def test_xdmf_field_as_string_found_in_species(tmpdir): + """Test that the field attribute can be a string and is found in the species list""" + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) + my_mat = F.Material(D_0=1, E_D=0, name="mat") + my_model.subdomains = [ + F.VolumeSubdomain1D(1, borders=[0.0, 4.0], material=my_mat), + ] + my_model.species = [F.Species("H")] + my_model.temperature = 500 + + filename = str(tmpdir.join("my_export.xdmf")) + my_export = F.XDMFExport(filename, field="H") + my_model.exports = [my_export] + my_model.settings = F.Settings(atol=1, rtol=0.1) + my_model.settings.stepsize = F.Stepsize(initial_value=1) + + my_model.initialise()