From 31f24e3dec0725860b9123f12bdf3509e88a8a30 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 13:33:32 -0400 Subject: [PATCH 01/77] find based on boundary not geometry --- festim/subdomain/surface_subdomain.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/festim/subdomain/surface_subdomain.py b/festim/subdomain/surface_subdomain.py index 409ca6671..f1ee32839 100644 --- a/festim/subdomain/surface_subdomain.py +++ b/festim/subdomain/surface_subdomain.py @@ -1,4 +1,5 @@ from dolfinx import fem +import dolfinx.mesh import numpy as np @@ -22,18 +23,18 @@ def __init__(self, id, x) -> None: self.id = id self.x = x - def locate_dof(self, function_space): + def locate_dof(self, mesh, fdim): """Locates the dof of the surface subdomain within the function space 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 """ - dofs = fem.locate_dofs_geometrical( - function_space, lambda x: np.isclose(x[0], self.x) + dofs = dolfinx.mesh.locate_entities_boundary( + mesh, fdim, lambda x: np.isclose(x[0], self.x) ) return dofs[0] From 08d6a58e915a29bcf40a5d8411ed07ceafbd5446 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 13:33:52 -0400 Subject: [PATCH 02/77] accept dict --- festim/material.py | 41 ++++++++++++++++++++++++++++++++--------- 1 file changed, 32 insertions(+), 9 deletions(-) diff --git a/festim/material.py b/festim/material.py index 494181e99..e8a5142c5 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 (1 species): + >>> my_mat = Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + Usage (multispecies): + >>> my_mat = Material( + D_0={"Species_1": 1.9e-7, "D": 2.0e-7}, + E_D={"H": 0.2, "D": 0.3}, + name="my_mat" + ) """ def __init__(self, D_0, E_D, name=None) -> None: @@ -26,7 +35,7 @@ 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): """Defines the diffusion coefficient Args: mesh (dolfinx.mesh.Mesh): the domain mesh @@ -34,9 +43,23 @@ def get_diffusion_coefficient(self, mesh, temperature): Returns: ufl.algebra.Product: the diffusion coefficient """ + if isinstance(self.D_0, float) and isinstance(self.E_D, float): + D_0 = F.as_fenics_constant(self.D_0, mesh) + E_D = F.as_fenics_constant(self.E_D, mesh) - # 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) - return D_0 * ufl.exp(-E_D / F.k_B / temperature) + elif isinstance(self.D_0, dict) and isinstance(self.E_D, dict): + # check lengths of dicts are the same + if len(self.D_0) != len(self.E_D): + raise ValueError( + "The number of pre-exponential factors and activation energies " + "must be the same" + ) + # check E_D keys for species or species.name + if species in self.D_0.keys(): + D_0 = F.as_fenics_constant(self.D_0[species], mesh) + elif species.name in self.D_0.keys(): + D_0 = F.as_fenics_constant(self.D_0[species.name], mesh) + else: + raise ValueError("Species not defined") From 08d26e9011407fb55035f228bd3c96b8867eafda Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 13:34:08 -0400 Subject: [PATCH 03/77] attriute sub_function_space --- festim/species.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/festim/species.py b/festim/species.py index a3f45eacd..075c40268 100644 --- a/festim/species.py +++ b/festim/species.py @@ -13,6 +13,7 @@ class 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 + sub_function_space (dolfinx.fem.FunctionSpace): the subspace of the function space concentration (dolfinx.fem.Function): the concentration of the species Usage: @@ -31,6 +32,7 @@ def __init__(self, name: str = None) -> None: self.solution = None self.prev_solution = None self.test_function = None + self.sub_function_space = None def __repr__(self) -> str: return f"Species({self.name})" From a83bd48b5b766e1f562fb363765da227bd6109f4 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:24:21 -0400 Subject: [PATCH 04/77] find species from name --- festim/__init__.py | 2 +- festim/species.py | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) 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/species.py b/festim/species.py index 075c40268..3d2817263 100644 --- a/festim/species.py +++ b/festim/species.py @@ -113,3 +113,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") From fa3ae112903aa0335d84b7da3d389630eff4a217 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:25:06 -0400 Subject: [PATCH 05/77] updated for multispecies --- festim/hydrogen_transport_problem.py | 88 ++++++++++++++++++++-------- 1 file changed, 63 insertions(+), 25 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 63b3cf37c..680110a80 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -135,13 +135,45 @@ def defing_export_writers(self): export.writer.write_mesh(self.mesh.mesh) def define_function_space(self): - elements = basix.ufl.element( + 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) + elements = [] + if len(self.species) <= 1: + mixed_element = element_CG + else: + for spe in self.species: + if isinstance(spe, F.Species): + # TODO check if mobile or immobile for traps + elements.append(element_CG) + mixed_element = ufl.MixedElement(elements) + + self.function_space = fem.FunctionSpace(self.mesh.mesh, mixed_element) + + self.u = Function(self.function_space) + self.u_n = Function(self.function_space) + + def assign_functions_to_species(self): + """Creates for each species the solution, prev solution and test + function""" + sub_solutions = list(ufl.split(self.u)) + sub_prev_solution = list(ufl.split(self.u_n)) + + if len(self.species) == 1: + sub_test_functions = [ufl.TestFunction(self.function_space)] + self.species[0].sub_function_space = self.function_space + else: + 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) + + 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""" @@ -155,7 +187,7 @@ 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) + dof = sub_dom.locate_dof(self.mesh.mesh, self.mesh.fdim) dofs_facets.append(dof) tags_facets.append(sub_dom.id) if isinstance(sub_dom, F.VolumeSubdomain1D): @@ -193,37 +225,37 @@ 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): + # for spe in self.species: + # if spe.name == bc.species: + # bc.species.assign(spe) + # else: + # raise ValueError( + # f"Species {bc.species} not found in model species" + # ) + if isinstance(bc.species, str): + # if name of species is given then replace with species object + name = bc.species + bc.species = F.find_species_from_name(name, self.species) if isinstance(bc, F.DirichletBC): bc_dofs = bc.define_surface_subdomain_dofs( - self.facet_meshtags, self.mesh, self.function_space + self.facet_meshtags, self.mesh, bc.species.sub_function_space ) bc.create_value( - self.mesh.mesh, self.function_space, self.temperature, self.t + self.mesh.mesh, + bc.species.sub_function_space, + self.temperature, + self.t, ) form = bc.create_formulation( - dofs=bc_dofs, function_space=self.function_space + dofs=bc_dofs, function_space=bc.species.sub_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: - 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) - - # TODO remove this - self.u = self.species[0].solution - self.u_n = self.species[0].prev_solution - 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 +266,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 +289,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) @@ -276,7 +308,7 @@ def run(self): n = self.mesh.n D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature + self.mesh.mesh, self.temperature, self.species[0] ) cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( @@ -294,8 +326,14 @@ def run(self): self.solver.solve(self.u) - # post processing + if len(self.species) == 1: + cm = self.u + else: + res = list(self.u.split()) + for idx, spe in enumerate(self.species): + spe.solution = res[idx] + # post processing surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) flux = assemble_scalar(surface_flux) From f64076dfa111b105eff393de12d2023a46f354f0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:25:15 -0400 Subject: [PATCH 06/77] accept dict --- festim/material.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/festim/material.py b/festim/material.py index e8a5142c5..494a361fb 100644 --- a/festim/material.py +++ b/festim/material.py @@ -43,7 +43,7 @@ def get_diffusion_coefficient(self, mesh, temperature, species): Returns: ufl.algebra.Product: the diffusion coefficient """ - if isinstance(self.D_0, float) and isinstance(self.E_D, float): + 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) @@ -56,10 +56,20 @@ def get_diffusion_coefficient(self, mesh, temperature, species): "The number of pre-exponential factors and activation energies " "must be the same" ) - # check E_D keys for species or species.name + # check D_0 keys for species or species.name if species in self.D_0.keys(): D_0 = F.as_fenics_constant(self.D_0[species], mesh) elif species.name in self.D_0.keys(): D_0 = F.as_fenics_constant(self.D_0[species.name], mesh) else: raise ValueError("Species not defined") + + # check E_D keys for species or species.name + if species in self.E_D.keys(): + E_D = F.as_fenics_constant(self.E_D[species], mesh) + elif species.name in self.E_D.keys(): + E_D = F.as_fenics_constant(self.E_D[species.name], mesh) + else: + raise ValueError("Species not definedy") + + return D_0 * ufl.exp(-E_D / F.k_B / temperature) From 9bf437380b43801756f647ed5e8c296631a969fd Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:25:34 -0400 Subject: [PATCH 07/77] updated tests for multispecies --- test/test_dirichlet_bc.py | 21 ++-- test/test_material.py | 2 +- test/test_permeation_problem.py | 188 ++++++++++++++++---------------- 3 files changed, 107 insertions(+), 104 deletions(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 48e7ac36f..d56c158a2 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() @@ -138,13 +136,14 @@ def test_callable_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() @@ -230,13 +230,14 @@ def test_create_formulation(value): # BUILD subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) - 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() @@ -276,7 +277,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) diff --git a/test/test_material.py b/test/test_material.py index e706cfe88..4e68c18a4 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -9,7 +9,7 @@ def test_define_diffusion_coefficient(): 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 = my_mat.get_diffusion_coefficient(test_mesh.mesh, T, species="dummy") D_analytical = D_0 * np.exp(-E_D / F.k_B / T) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 2911cf753..7ced0cde2 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -59,7 +59,9 @@ def test_permeation_problem(mesh_size=1001): # -------------------------- analytical solution ------------------------------------- - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature) + D = my_mat.get_diffusion_coefficient( + my_mesh.mesh, my_model.temperature, my_model.species[0] + ) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) @@ -89,95 +91,95 @@ def test_permeation_problem(mesh_size=1001): assert error < 0.01 -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.boundary_conditions = [ - F.DirichletBC(subdomain=right_surface, value=0, species="H"), - F.SievertsBC( - 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.settings = F.Settings( - atol=1e10, - rtol=1e-10, - max_iterations=30, - final_time=50, - ) - - 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 solution ------------------------------------- - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) - - S_0 = float(my_model.boundary_conditions[-1].S_0) - E_S = float(my_model.boundary_conditions[-1].E_S) - P_up = float(my_model.boundary_conditions[-1].pressure) - 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() - - assert error < 0.01 +# 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.boundary_conditions = [ +# F.DirichletBC(subdomain=right_surface, value=0, species="H"), +# F.SievertsBC( +# 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.settings = F.Settings( +# atol=1e10, +# rtol=1e-10, +# max_iterations=30, +# final_time=50, +# ) + +# 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 solution ------------------------------------- +# D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) + +# S_0 = float(my_model.boundary_conditions[-1].S_0) +# E_S = float(my_model.boundary_conditions[-1].E_S) +# P_up = float(my_model.boundary_conditions[-1].pressure) +# 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() + +# assert error < 0.01 From 441aeeaec358c2f8418cc2b00bae1b797de44366 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:48:23 -0400 Subject: [PATCH 08/77] added tests --- test/test_material.py | 18 ++++++++++++++++++ test/test_multispecies_problem.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) create mode 100644 test/test_multispecies_problem.py diff --git a/test/test_material.py b/test/test_material.py index 4e68c18a4..8c57553c6 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -14,3 +14,21 @@ def test_define_diffusion_coefficient(): D_analytical = D_0 * np.exp(-E_D / F.k_B / T) assert np.isclose(float(D), D_analytical) + + +def test_multispecies_dict(): + T = 500 + D_0_A, D_0_B = 1, 2 + E_D_A, E_D_B = 0.1, 0.2 + + 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() diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py new file mode 100644 index 000000000..cac552185 --- /dev/null +++ b/test/test_multispecies_problem.py @@ -0,0 +1,28 @@ +import numpy as np +import festim as F + + +def test_multispecies_problem_initialisation(): + """Test that the multispecies problem is correctly initialised""" + my_model = F.HydrogenTransportProblem() + L = 1e-04 + my_model.mesh = F.Mesh1D(np.linspace(0, L, num=1001)) + my_mat = F.Material( + D_0={"H": 1.9e-7, "D": 1.9e-07}, E_D={"H": 0.2, "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] + mobile_H = F.Species("H") + mobile_D = F.Species("D") + my_model.species = [mobile_H, mobile_D] + my_model.temperature = 600 + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left_surface, value=1e18, species="H"), + F.DirichletBC(subdomain=right_surface, value=1e18, species="D"), + ] + my_model.settings = F.Settings(atol=1e10, rtol=1e-10) + my_model.settings.stepsize = 0.1 + + my_model.initialise() From 8480c60b05732f568e5641bd2614fb528fe609f8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:48:53 -0400 Subject: [PATCH 09/77] remove commented snip, post_processing solution --- festim/hydrogen_transport_problem.py | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 680110a80..b309dc4eb 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -225,14 +225,6 @@ 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): - # for spe in self.species: - # if spe.name == bc.species: - # bc.species.assign(spe) - # else: - # raise ValueError( - # f"Species {bc.species} not found in model species" - # ) if isinstance(bc.species, str): # if name of species is given then replace with species object name = bc.species @@ -331,7 +323,7 @@ def run(self): else: res = list(self.u.split()) for idx, spe in enumerate(self.species): - spe.solution = res[idx] + spe.post_processing_solution = res[idx] # post processing surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) From b0e887570d6bb3f2e17ba137ae518d3375d587bf Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 14:49:06 -0400 Subject: [PATCH 10/77] post_processing solution --- festim/species.py | 1 + 1 file changed, 1 insertion(+) diff --git a/festim/species.py b/festim/species.py index 3d2817263..dc35ad35d 100644 --- a/festim/species.py +++ b/festim/species.py @@ -33,6 +33,7 @@ def __init__(self, name: str = None) -> None: self.prev_solution = None self.test_function = None self.sub_function_space = None + self.post_processing_solution = None def __repr__(self) -> str: return f"Species({self.name})" From 8ee41f62f1daf70a3eac1809deabb981aa5548a2 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 15:12:37 -0400 Subject: [PATCH 11/77] uncomment test --- test/test_permeation_problem.py | 184 ++++++++++++++++---------------- 1 file changed, 92 insertions(+), 92 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 7ced0cde2..7a541753d 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -91,95 +91,95 @@ def test_permeation_problem(mesh_size=1001): assert error < 0.01 -# 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.boundary_conditions = [ -# F.DirichletBC(subdomain=right_surface, value=0, species="H"), -# F.SievertsBC( -# 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.settings = F.Settings( -# atol=1e10, -# rtol=1e-10, -# max_iterations=30, -# final_time=50, -# ) - -# 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 solution ------------------------------------- -# D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) - -# S_0 = float(my_model.boundary_conditions[-1].S_0) -# E_S = float(my_model.boundary_conditions[-1].E_S) -# P_up = float(my_model.boundary_conditions[-1].pressure) -# 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() - -# assert error < 0.01 +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.boundary_conditions = [ + F.DirichletBC(subdomain=right_surface, value=0, species="H"), + F.SievertsBC( + 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.settings = F.Settings( + atol=1e10, + rtol=1e-10, + max_iterations=30, + final_time=50, + ) + + 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 solution ------------------------------------- + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, species="H") + + S_0 = float(my_model.boundary_conditions[-1].S_0) + E_S = float(my_model.boundary_conditions[-1].E_S) + P_up = float(my_model.boundary_conditions[-1].pressure) + 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]m + flux_values = flux_values[indices] + + relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) + + error = relative_error.mean() + + assert error < 0.01 From b4bef5f817b29a29f797c45209fbbf0f516b9f08 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 15:12:48 -0400 Subject: [PATCH 12/77] refactoring --- festim/material.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/material.py b/festim/material.py index 494a361fb..1e2f3daa9 100644 --- a/festim/material.py +++ b/festim/material.py @@ -53,7 +53,7 @@ def get_diffusion_coefficient(self, mesh, temperature, species): # check lengths of dicts are the same if len(self.D_0) != len(self.E_D): raise ValueError( - "The number of pre-exponential factors and activation energies " + "The number of pre-exponential factors and activation energies" "must be the same" ) # check D_0 keys for species or species.name From 42f6cb548889702f135a1abd1ba189a1ac4f6377 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 15:12:59 -0400 Subject: [PATCH 13/77] more tests and doc strings --- test/test_material.py | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/test/test_material.py b/test/test_material.py index 8c57553c6..a0d56fd27 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])) @@ -16,7 +17,9 @@ def test_define_diffusion_coefficient(): assert np.isclose(float(D), D_analytical) -def test_multispecies_dict(): +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 @@ -32,3 +35,35 @@ def test_multispecies_dict(): 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_different_lengths(): + """Test that a value error is rasied when the length of the D_0 and E_D + are not the same""" + 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): + D = my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species="A") From 449abf45c60e8850e4f88c75fec62c8a3c45c044 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 15:14:51 -0400 Subject: [PATCH 14/77] typo --- test/test_permeation_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 7a541753d..906919892 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -175,7 +175,7 @@ def test_permeation_problem_multi_volume(): flux_values = np.array(np.abs(flux_values)) indices = np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - analytical_flux = analytical_flux[indices]m + analytical_flux = analytical_flux[indices] flux_values = flux_values[indices] relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) From 00befddde535a0a1c8c66d70927dec48336e99f2 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 20 Oct 2023 17:31:13 -0400 Subject: [PATCH 15/77] test new helper --- test/test_dirichlet_bc.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index d56c158a2..14ca8e5c8 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -319,3 +319,24 @@ 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 rasied 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() From 69b2eccc80580aaa15770a9d39282e52e29734a8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:34:34 -0400 Subject: [PATCH 16/77] test materials more throughly --- test/test_material.py | 93 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 83 insertions(+), 10 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index a0d56fd27..326d866eb 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -8,9 +8,11 @@ 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, species="dummy") + D = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=dum_spe, model_species=[dum_spe] + ) D_analytical = D_0 * np.exp(-E_D / F.k_B / T) @@ -23,10 +25,17 @@ def test_multispecies_dict_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") + spe_list = [A, 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_A = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=A, model_species=spe_list + ) + D_B = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=B, model_species=spe_list + ) + D = [float(D_A), float(D_B)] D_A_analytical = D_0_A * np.exp(-E_D_A / F.k_B / T) @@ -46,10 +55,15 @@ def test_multispecies_dict_objects(): A = F.Species("A") B = F.Species("B") + spe_list = [A, 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_A = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=A, model_species=spe_list + ) + D_B = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=B, model_species=spe_list + ) D = [float(D_A), float(D_B)] D_A_analytical = D_0_A * np.exp(-E_D_A / F.k_B / T) @@ -60,10 +74,69 @@ def test_multispecies_dict_objects(): assert np.isclose(D, D_analytical).all() -def test_multispecies_dict_different_lengths(): - """Test that a value error is rasied when the length of the D_0 and E_D +def test_multispecies_dict_objects_and_strings(): + """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") + spe_list = [A, 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, model_species=spe_list + ) + D_B = my_mat.get_diffusion_coefficient( + test_mesh.mesh, T, species=B, model_species=spe_list + ) + 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 rasied when the keys of the D_0 and E_D are not the same""" + A = F.Species("A") + spe_list = [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): - D = my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species="A") + with pytest.raises(ValueError, match="D_0 and E_D have different keys"): + my_mat.get_diffusion_coefficient( + test_mesh.mesh, 500, species=A, model_species=spe_list + ) + + +def test_multispecies_dict_wrong_name_species_not_found(): + """Test that a value error is rasied when the length of the D_0 and E_D + are not the same""" + J = F.Species("J") + spe_list = [J] + my_mat = F.Material(D_0={"A": 1, "B": 2}, E_D={"A": 0.1, "B": 0.2}) + + with pytest.raises(ValueError, match="Species A not found in list of species"): + my_mat.get_diffusion_coefficient( + test_mesh.mesh, 500, species=J, model_species=spe_list + ) + + +def test_multispecies_dict_contains_species_not_in_species_list(): + """Test that a value error is rasied in the get_diffusion_coefficient + function""" + J = F.Species("J") + A = F.Species("A") + spe_list = [J] + my_mat = F.Material(D_0={A: 1, "B": 2}, E_D={A: 0.1, "B": 0.2}) + + with pytest.raises(ValueError, match="Species A not found in model species"): + my_mat.get_diffusion_coefficient( + test_mesh.mesh, 500, species=J, model_species=spe_list + ) From 03aaa1f18c5592f83219f390a20e0e940e7ed5a3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:48:46 -0400 Subject: [PATCH 17/77] use solution, pass species list --- festim/hydrogen_transport_problem.py | 34 +++++++++++++++++++--------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index b309dc4eb..9d607f0ca 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -159,19 +159,29 @@ def define_function_space(self): def assign_functions_to_species(self): """Creates for each species the solution, prev solution and test function""" - sub_solutions = list(ufl.split(self.u)) - sub_prev_solution = list(ufl.split(self.u_n)) if len(self.species) == 1: + 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.species[0].sub_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) + post_processing_function_space, _ = spe.sub_function_space.collapse() + spe.post_processing_solution = fem.Function( + post_processing_function_space + ) for idx, spe in enumerate(self.species): spe.solution = sub_solutions[idx] + spe.post_processing_solution = sub_solutions[idx] spe.prev_solution = sub_prev_solution[idx] spe.test_function = sub_test_functions[idx] @@ -258,7 +268,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, spe + self.mesh.mesh, self.temperature, spe, self.species ) self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) @@ -300,7 +310,7 @@ def run(self): n = self.mesh.n D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] + self.mesh.mesh, self.temperature, self.species[0], self.species ) cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( @@ -320,17 +330,19 @@ def run(self): if len(self.species) == 1: cm = self.u + self.species[0].post_processing_solution = 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)) + else: res = list(self.u.split()) for idx, spe in enumerate(self.species): - spe.post_processing_solution = res[idx] - - # post processing - surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) + spe.solution = res[idx] - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(float(self.t)) + times = flux_values = 0 for export in self.exports: if isinstance(export, (F.VTXExport, F.XDMFExport)): From 6c63813022ec3e02706aded1c384e94b0f1b7479 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:49:05 -0400 Subject: [PATCH 18/77] better material processing --- festim/material.py | 47 ++++++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/festim/material.py b/festim/material.py index 1e2f3daa9..72e61525a 100644 --- a/festim/material.py +++ b/festim/material.py @@ -35,14 +35,20 @@ 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, species): + def get_diffusion_coefficient(self, mesh, temperature, species, model_species): """Defines the diffusion coefficient Args: + mesh (dolfinx.mesh.Mesh): the domain mesh temperature (dolfinx.fem.Constant): the temperature + species (festim.Species): the species + model_species (list): the list of species in the model Returns: ufl.algebra.Product: the diffusion coefficient """ + if species not in model_species: + raise ValueError(f"Species {species} not found in model species") + 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) @@ -50,26 +56,27 @@ def get_diffusion_coefficient(self, mesh, temperature, species): return D_0 * ufl.exp(-E_D / F.k_B / temperature) elif isinstance(self.D_0, dict) and isinstance(self.E_D, dict): - # check lengths of dicts are the same - if len(self.D_0) != len(self.E_D): - raise ValueError( - "The number of pre-exponential factors and activation energies" - "must be the same" - ) - # check D_0 keys for species or species.name - if species in self.D_0.keys(): - D_0 = F.as_fenics_constant(self.D_0[species], mesh) - elif species.name in self.D_0.keys(): + # check D_0 and E_D have the same keys + if list(self.D_0.keys()) != list(self.E_D.keys()): + raise ValueError("D_0 and E_D have different keys") + + for key in self.D_0.keys(): + if isinstance(key, str): + F.find_species_from_name(key, model_species) + elif key not in model_species: + raise ValueError(f"Species {key} not found in model species") + + try: D_0 = F.as_fenics_constant(self.D_0[species.name], mesh) - else: - raise ValueError("Species not defined") + except KeyError: + D_0 = F.as_fenics_constant(self.D_0[species], mesh) - # check E_D keys for species or species.name - if species in self.E_D.keys(): - E_D = F.as_fenics_constant(self.E_D[species], mesh) - elif species.name in self.E_D.keys(): + try: E_D = F.as_fenics_constant(self.E_D[species.name], mesh) - else: - raise ValueError("Species not definedy") + except KeyError: + E_D = F.as_fenics_constant(self.E_D[species], 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") From 1cf208be5120b6c20b642c68b91b8827fe92e9ac Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:49:18 -0400 Subject: [PATCH 19/77] pass species list --- test/test_permeation_problem.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 906919892..587235e99 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -60,7 +60,7 @@ def test_permeation_problem(mesh_size=1001): # -------------------------- analytical solution ------------------------------------- D = my_mat.get_diffusion_coefficient( - my_mesh.mesh, my_model.temperature, my_model.species[0] + my_mesh.mesh, my_model.temperature, my_model.species[0], my_model.species ) S_0 = float(my_model.boundary_conditions[-1].S_0) @@ -154,8 +154,10 @@ def test_permeation_problem_multi_volume(): times, flux_values = my_model.run() - # -------------------------- analytical solution ------------------------------------- - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, species="H") + # ---------------------- analytical solution ----------------------------- + D = my_mat.get_diffusion_coefficient( + my_mesh.mesh, temperature, my_model.species[0], model_species=my_model.species + ) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) From 8f894eb97fcacea389ccffe521ce0df969e36050 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:56:49 -0400 Subject: [PATCH 20/77] more tests for material to increase coverage --- test/test_material.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/test_material.py b/test/test_material.py index 326d866eb..23a780d9d 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -140,3 +140,30 @@ def test_multispecies_dict_contains_species_not_in_species_list(): my_mat.get_diffusion_coefficient( test_mesh.mesh, 500, species=J, model_species=spe_list ) + + +def test_contains_species_not_in_species_list(): + """Test that a value error is rasied in the get_diffusion_coefficient + function with one species""" + J = F.Species("J") + A = F.Species("A") + spe_list = [J] + my_mat = F.Material(D_0=1, E_D=0) + + with pytest.raises(ValueError, match="Species A not found in model species"): + my_mat.get_diffusion_coefficient( + test_mesh.mesh, 500, species=A, model_species=spe_list + ) + + +def test_D_0_type_rasies_error(): + """Test that a value error is rasied in the get_diffusion_coefficient + function""" + A = F.Species("A") + spe_list = [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, model_species=spe_list + ) From c7eb997215d102cdd40b1fb4529ce8a9194f8807 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 10:58:16 -0400 Subject: [PATCH 21/77] line not needed --- festim/hydrogen_transport_problem.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 9d607f0ca..a885cfe17 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -342,8 +342,6 @@ def run(self): for idx, spe in enumerate(self.species): spe.solution = res[idx] - times = flux_values = 0 - for export in self.exports: if isinstance(export, (F.VTXExport, F.XDMFExport)): export.write(float(self.t)) From 5d7458014f13817b133cdfe6c682c5f178eda9e1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 21 Oct 2023 11:01:31 -0400 Subject: [PATCH 22/77] comment for now --- festim/hydrogen_transport_problem.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index a885cfe17..e3c211aab 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -337,10 +337,11 @@ def run(self): flux_values.append(flux) times.append(float(self.t)) - else: - res = list(self.u.split()) - for idx, spe in enumerate(self.species): - spe.solution = res[idx] + # TODO in multi-species export all functions + # else: + # res = list(self.u.split()) + # for idx, spe in enumerate(self.species): + # spe.solution = res[idx] for export in self.exports: if isinstance(export, (F.VTXExport, F.XDMFExport)): From 07d35ef292fc206169e174639b6d19f0dbdbb37c Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 10:57:23 -0400 Subject: [PATCH 23/77] export post_processing_solution --- festim/exports/vtx.py | 5 ++++- festim/exports/xdmf.py | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/festim/exports/vtx.py b/festim/exports/vtx.py index 76d2d73a7..62ba0fe5e 100644 --- a/festim/exports/vtx.py +++ b/festim/exports/vtx.py @@ -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..a45514956 100644 --- a/festim/exports/xdmf.py +++ b/festim/exports/xdmf.py @@ -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) From e32d49a8afc6aa9a4e1058c55d80aa33b4c5a95b Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 10:58:32 -0400 Subject: [PATCH 24/77] updated doc strings --- festim/material.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/material.py b/festim/material.py index 72e61525a..59d8deba1 100644 --- a/festim/material.py +++ b/festim/material.py @@ -24,8 +24,8 @@ class Material: >>> my_mat = Material(D_0=1.9e-7, E_D=0.2, name="my_mat") Usage (multispecies): >>> my_mat = Material( - D_0={"Species_1": 1.9e-7, "D": 2.0e-7}, - E_D={"H": 0.2, "D": 0.3}, + 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" ) """ From a2a4d87def8d12b603610e4c801752dadacf1300 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 10:58:46 -0400 Subject: [PATCH 25/77] updated tests with new export method --- test/test_multispecies_problem.py | 13 ++++++++++++- test/test_vtx.py | 6 +++--- test/test_xdmf.py | 2 +- 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index cac552185..6dc60170e 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -22,7 +22,18 @@ def test_multispecies_problem_initialisation(): F.DirichletBC(subdomain=left_surface, value=1e18, species="H"), F.DirichletBC(subdomain=right_surface, value=1e18, species="D"), ] - my_model.settings = F.Settings(atol=1e10, rtol=1e-10) + my_model.settings = F.Settings(atol=1e10, rtol=1e-10, final_time=50) my_model.settings.stepsize = 0.1 + my_model.exports = [ + F.XDMFExport("results/multispecies/test_H.xdmf", field=mobile_H), + F.XDMFExport("results/multispecies/test_D.xdmf", field=mobile_D), + F.VTXExport("results/multispecies/test_H_vts.bp", field=mobile_H), + F.VTXExport("results/multispecies/test_D_vtx.bp", field=mobile_D), + ] my_model.initialise() + my_model.run() + + +if __name__ == "__main__": + test_multispecies_problem_initialisation() diff --git a/test/test_vtx.py b/test/test_vtx.py index 15c7ded6e..1366a10e0 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]) diff --git a/test/test_xdmf.py b/test/test_xdmf.py index cc1e57a25..8a3f8bd53 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) From 5254a20f2b2b79797185f3479fcd48eee6a9bfea Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 10:59:02 -0400 Subject: [PATCH 26/77] export multispecies cases --- festim/hydrogen_transport_problem.py | 50 +++++++++++++--------------- 1 file changed, 23 insertions(+), 27 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e3c211aab..e51590cb6 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -128,7 +128,6 @@ def defing_export_writers(self): for export in self.exports: # TODO implement when export.field is an int or str # then find solution from index of species - if isinstance(export, (F.VTXExport, F.XDMFExport)): export.define_writer(MPI.COMM_WORLD) if isinstance(export, F.XDMFExport): @@ -142,16 +141,16 @@ def define_function_space(self): basix.LagrangeVariant.equispaced, ) elements = [] - if len(self.species) <= 1: - mixed_element = element_CG - else: + if len(self.species) > 1: for spe in self.species: if isinstance(spe, F.Species): # TODO check if mobile or immobile for traps elements.append(element_CG) - mixed_element = ufl.MixedElement(elements) + element_data = ufl.MixedElement(elements) + else: + element_data = element_CG - self.function_space = fem.FunctionSpace(self.mesh.mesh, mixed_element) + self.function_space = fem.FunctionSpace(self.mesh.mesh, element_data) self.u = Function(self.function_space) self.u_n = Function(self.function_space) @@ -160,28 +159,24 @@ def assign_functions_to_species(self): """Creates for each species the solution, prev solution and test function""" - if len(self.species) == 1: - 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.species[0].sub_function_space - ) - else: + if len(self.species) > 1: 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) - post_processing_function_space, _ = spe.sub_function_space.collapse() - spe.post_processing_solution = fem.Function( - post_processing_function_space - ) + spe.post_processing_solution = self.u.sub(idx).collapse() + + else: + 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) for idx, spe in enumerate(self.species): spe.solution = sub_solutions[idx] - spe.post_processing_solution = sub_solutions[idx] spe.prev_solution = sub_prev_solution[idx] spe.test_function = sub_test_functions[idx] @@ -328,7 +323,14 @@ def run(self): self.solver.solve(self.u) - if len(self.species) == 1: + if len(self.species) > 1: + # res = list(self.u.split()) + for idx, spe in enumerate(self.species): + pass + # print(self.u.sub(idx).collapse()) + # quit() + # spe.post_processing_solution = self.u.sub(idx).collapse() + else: cm = self.u self.species[0].post_processing_solution = self.u # post processing @@ -337,12 +339,6 @@ def run(self): flux_values.append(flux) times.append(float(self.t)) - # TODO in multi-species export all functions - # else: - # res = list(self.u.split()) - # for idx, spe in enumerate(self.species): - # spe.solution = res[idx] - for export in self.exports: if isinstance(export, (F.VTXExport, F.XDMFExport)): export.write(float(self.t)) From 6ee60b7c4c6178c050824577d061a392246dfc6d Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 11:30:30 -0400 Subject: [PATCH 27/77] vtx export field as string --- festim/exports/vtx.py | 8 ++++---- festim/hydrogen_transport_problem.py | 9 +++++++++ test/test_vtx.py | 21 ++++++++++++++++++++- 3 files changed, 33 insertions(+), 5 deletions(-) diff --git a/festim/exports/vtx.py b/festim/exports/vtx.py index 62ba0fe5e..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): diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e51590cb6..87f01b5ed 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -128,6 +128,15 @@ def defing_export_writers(self): for export in self.exports: # TODO implement when export.field is an int or str # then find solution from index of species + if isinstance(export.field, str): + export.field = F.find_species_from_name(field, self.species) + elif isinstance(export.field, list): + 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) if isinstance(export, F.XDMFExport): diff --git a/test/test_vtx.py b/test/test_vtx.py index 1366a10e0..5681bf878 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -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,22 @@ 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): + 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() From 6e32e79b80d2fc547e1ff760eecb9ccb953bf223 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 11:33:15 -0400 Subject: [PATCH 28/77] xdmf accept field as str --- festim/exports/xdmf.py | 8 ++++---- test/test_xdmf.py | 21 ++++++++++++++++++++- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/festim/exports/xdmf.py b/festim/exports/xdmf.py index a45514956..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): diff --git a/test/test_xdmf.py b/test/test_xdmf.py index 8a3f8bd53..3e86ca5db 100644 --- a/test/test_xdmf.py +++ b/test/test_xdmf.py @@ -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,22 @@ 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): + 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() From 9a4475c4f18000bd00ce452e9a39d88f10c117ce Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 11:42:31 -0400 Subject: [PATCH 29/77] field always list, refactoring --- festim/hydrogen_transport_problem.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 87f01b5ed..9ed1217e4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -128,14 +128,9 @@ def defing_export_writers(self): for export in self.exports: # TODO implement when export.field is an int or str # then find solution from index of species - if isinstance(export.field, str): - export.field = F.find_species_from_name(field, self.species) - elif isinstance(export.field, list): - for idx, field in enumerate(export.field): - if isinstance(field, str): - export.field[idx] = F.find_species_from_name( - field, self.species - ) + 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) @@ -241,8 +236,7 @@ def define_boundary_conditions(self): for bc in self.boundary_conditions: if isinstance(bc.species, str): # if name of species is given then replace with species object - name = bc.species - bc.species = F.find_species_from_name(name, self.species) + 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, bc.species.sub_function_space From 42386b5bddf587f203634b9173f24dbf21ff0aeb Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 11:59:46 -0400 Subject: [PATCH 30/77] refactoring --- festim/hydrogen_transport_problem.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 9ed1217e4..5b2d4f1b1 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -326,17 +326,10 @@ def run(self): self.solver.solve(self.u) - if len(self.species) > 1: - # res = list(self.u.split()) - for idx, spe in enumerate(self.species): - pass - # print(self.u.sub(idx).collapse()) - # quit() - # spe.post_processing_solution = self.u.sub(idx).collapse() - else: + if len(self.species) == 1: cm = self.u self.species[0].post_processing_solution = self.u - # post processing + surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) From ad0af6b4f8b67a99896a745c7110373834a4ec82 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 12:02:11 -0400 Subject: [PATCH 31/77] tempdr --- test/test_multispecies_problem.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 6dc60170e..103142e46 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -2,7 +2,7 @@ import festim as F -def test_multispecies_problem_initialisation(): +def test_multispecies_problem_initialisation(tmpdir): """Test that the multispecies problem is correctly initialised""" my_model = F.HydrogenTransportProblem() L = 1e-04 @@ -25,15 +25,11 @@ def test_multispecies_problem_initialisation(): my_model.settings = F.Settings(atol=1e10, rtol=1e-10, final_time=50) my_model.settings.stepsize = 0.1 my_model.exports = [ - F.XDMFExport("results/multispecies/test_H.xdmf", field=mobile_H), - F.XDMFExport("results/multispecies/test_D.xdmf", field=mobile_D), - F.VTXExport("results/multispecies/test_H_vts.bp", field=mobile_H), - F.VTXExport("results/multispecies/test_D_vtx.bp", field=mobile_D), + F.XDMFExport(str(tmpdir.join("my_export_1.xdmf")), field=mobile_H), + F.XDMFExport(str(tmpdir.join("my_export_2.xdmf")), field="D"), + F.VTXExport(str(tmpdir.join("my_export_3.bp")), field="H"), + F.VTXExport(str(tmpdir.join("my_export_4.bp")), field=mobile_D), ] my_model.initialise() my_model.run() - - -if __name__ == "__main__": - test_multispecies_problem_initialisation() From 89de52f514c5c2ebbb131d588831ce37924945d7 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 12:02:20 -0400 Subject: [PATCH 32/77] added doc strings --- test/test_vtx.py | 1 + test/test_xdmf.py | 1 + 2 files changed, 2 insertions(+) diff --git a/test/test_vtx.py b/test/test_vtx.py index 5681bf878..1f0823740 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -98,6 +98,7 @@ def test_filename_raises_error_when_wrong_type(): 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") diff --git a/test/test_xdmf.py b/test/test_xdmf.py index 3e86ca5db..fd1dabc3f 100644 --- a/test/test_xdmf.py +++ b/test/test_xdmf.py @@ -85,6 +85,7 @@ def test_filename_raises_error_when_wrong_type(): 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") From e3e5feb3ad52cda488b60e08a6e9f80a18829452 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 12:04:13 -0400 Subject: [PATCH 33/77] updated comment --- festim/hydrogen_transport_problem.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 5b2d4f1b1..3fac2caf1 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -126,8 +126,7 @@ def initialise(self): def defing_export_writers(self): """Defines the export writers of the model""" 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) From f0cd736e4ebfda1bd6968062c8c78d8585d8b9a3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 12:06:21 -0400 Subject: [PATCH 34/77] updated doc strings --- festim/species.py | 1 + 1 file changed, 1 insertion(+) diff --git a/festim/species.py b/festim/species.py index dc35ad35d..6cbc3f5a6 100644 --- a/festim/species.py +++ b/festim/species.py @@ -14,6 +14,7 @@ class Species: prev_solution (dolfinx.fem.Function or ...): the solution for the previous timestep test_function (ufl.Argument or ...): the testfunction associated with this species sub_function_space (dolfinx.fem.FunctionSpace): the subspace of the function space + post_processing_solution (dolfinx.fem.Function): the solution for post processing concentration (dolfinx.fem.Function): the concentration of the species Usage: From 2c12d444563d92eeec70df01c712a3cd5fad724f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 23 Oct 2023 13:36:26 -0400 Subject: [PATCH 35/77] get_diffusion_coeff doesn't require list of species --- festim/hydrogen_transport_problem.py | 4 +- festim/material.py | 32 +++++------ test/test_material.py | 81 ++++------------------------ test/test_permeation_problem.py | 6 +-- 4 files changed, 29 insertions(+), 94 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3fac2caf1..cfbffd189 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -265,7 +265,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, spe, self.species + self.mesh.mesh, self.temperature, spe ) self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) @@ -307,7 +307,7 @@ def run(self): n = self.mesh.n D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0], self.species + self.mesh.mesh, self.temperature, self.species[0] ) cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( diff --git a/festim/material.py b/festim/material.py index 59d8deba1..811766f6d 100644 --- a/festim/material.py +++ b/festim/material.py @@ -35,19 +35,17 @@ 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, species, model_species): + 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): the species - model_species (list): the list of species in the model + 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 """ - if species not in model_species: - raise ValueError(f"Species {species} not found in model species") if isinstance(self.D_0, (float, int)) and isinstance(self.E_D, (float, int)): D_0 = F.as_fenics_constant(self.D_0, mesh) @@ -57,24 +55,26 @@ def get_diffusion_coefficient(self, mesh, temperature, species, model_species): 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") - for key in self.D_0.keys(): - if isinstance(key, str): - F.find_species_from_name(key, model_species) - elif key not in model_species: - raise ValueError(f"Species {key} not found in model species") + if species is None: + raise ValueError("species must be provided if D_0 and E_D are dicts") - try: - D_0 = F.as_fenics_constant(self.D_0[species.name], mesh) - except KeyError: + 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") - try: - E_D = F.as_fenics_constant(self.E_D[species.name], mesh) - except KeyError: + 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) + else: + raise ValueError(f"{species} is not in E_D keys") return D_0 * ufl.exp(-E_D / F.k_B / temperature) diff --git a/test/test_material.py b/test/test_material.py index 23a780d9d..df8f70c08 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -10,9 +10,7 @@ def test_define_diffusion_coefficient(): 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, species=dum_spe, model_species=[dum_spe] - ) + 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) @@ -26,15 +24,10 @@ def test_multispecies_dict_strings(): 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") - spe_list = [A, 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, model_species=spe_list - ) - D_B = my_mat.get_diffusion_coefficient( - test_mesh.mesh, T, species=B, model_species=spe_list - ) + 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)] @@ -55,15 +48,10 @@ def test_multispecies_dict_objects(): A = F.Species("A") B = F.Species("B") - spe_list = [A, 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, model_species=spe_list - ) - D_B = my_mat.get_diffusion_coefficient( - test_mesh.mesh, T, species=B, model_species=spe_list - ) + 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) @@ -83,15 +71,10 @@ def test_multispecies_dict_objects_and_strings(): A = F.Species("A") B = F.Species("B") - spe_list = [A, 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, model_species=spe_list - ) - D_B = my_mat.get_diffusion_coefficient( - test_mesh.mesh, T, species=B, model_species=spe_list - ) + 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) @@ -106,54 +89,10 @@ def test_multispecies_dict_different_keys(): """Test that a value error is rasied when the keys of the D_0 and E_D are not the same""" A = F.Species("A") - spe_list = [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, model_species=spe_list - ) - - -def test_multispecies_dict_wrong_name_species_not_found(): - """Test that a value error is rasied when the length of the D_0 and E_D - are not the same""" - J = F.Species("J") - spe_list = [J] - my_mat = F.Material(D_0={"A": 1, "B": 2}, E_D={"A": 0.1, "B": 0.2}) - - with pytest.raises(ValueError, match="Species A not found in list of species"): - my_mat.get_diffusion_coefficient( - test_mesh.mesh, 500, species=J, model_species=spe_list - ) - - -def test_multispecies_dict_contains_species_not_in_species_list(): - """Test that a value error is rasied in the get_diffusion_coefficient - function""" - J = F.Species("J") - A = F.Species("A") - spe_list = [J] - my_mat = F.Material(D_0={A: 1, "B": 2}, E_D={A: 0.1, "B": 0.2}) - - with pytest.raises(ValueError, match="Species A not found in model species"): - my_mat.get_diffusion_coefficient( - test_mesh.mesh, 500, species=J, model_species=spe_list - ) - - -def test_contains_species_not_in_species_list(): - """Test that a value error is rasied in the get_diffusion_coefficient - function with one species""" - J = F.Species("J") - A = F.Species("A") - spe_list = [J] - my_mat = F.Material(D_0=1, E_D=0) - - with pytest.raises(ValueError, match="Species A not found in model species"): - my_mat.get_diffusion_coefficient( - test_mesh.mesh, 500, species=A, model_species=spe_list - ) + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) def test_D_0_type_rasies_error(): @@ -164,6 +103,4 @@ def test_D_0_type_rasies_error(): 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, model_species=spe_list - ) + my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 587235e99..d69b3dfa6 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -60,7 +60,7 @@ def test_permeation_problem(mesh_size=1001): # -------------------------- analytical solution ------------------------------------- D = my_mat.get_diffusion_coefficient( - my_mesh.mesh, my_model.temperature, my_model.species[0], my_model.species + my_mesh.mesh, my_model.temperature, my_model.species[0] ) S_0 = float(my_model.boundary_conditions[-1].S_0) @@ -155,9 +155,7 @@ def test_permeation_problem_multi_volume(): times, flux_values = my_model.run() # ---------------------- analytical solution ----------------------------- - D = my_mat.get_diffusion_coefficient( - my_mesh.mesh, temperature, my_model.species[0], model_species=my_model.species - ) + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, my_model.species[0]) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) From 8787e8e91c7dae3da891a3bfe66e4faf568f5221 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Mon, 23 Oct 2023 14:08:07 -0400 Subject: [PATCH 36/77] only use list if multispecies MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3fac2caf1..f8f02717e 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -143,8 +143,8 @@ def define_function_space(self): 1, basix.LagrangeVariant.equispaced, ) - elements = [] if len(self.species) > 1: + elements = [] for spe in self.species: if isinstance(spe, F.Species): # TODO check if mobile or immobile for traps From 00c63f8a2e619602ff0e4bae38bfda0a1c860a3b Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Mon, 23 Oct 2023 14:08:37 -0400 Subject: [PATCH 37/77] simpler naming MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- festim/hydrogen_transport_problem.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index f8f02717e..4f0ac625b 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -149,11 +149,11 @@ def define_function_space(self): if isinstance(spe, F.Species): # TODO check if mobile or immobile for traps elements.append(element_CG) - element_data = ufl.MixedElement(elements) + element = ufl.MixedElement(elements) else: - element_data = element_CG + element = element_CG - self.function_space = fem.FunctionSpace(self.mesh.mesh, element_data) + self.function_space = fem.FunctionSpace(self.mesh.mesh, element) self.u = Function(self.function_space) self.u_n = Function(self.function_space) From 3c2cba887c669f587fcf50df419f313b10a6b50e Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Mon, 23 Oct 2023 14:08:50 -0400 Subject: [PATCH 38/77] typo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- test/test_dirichlet_bc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 14ca8e5c8..eb5e26bc5 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -322,7 +322,7 @@ def test_integration_with_HTransportProblem(value): def test_species_predefined(): - """Test a ValueError is rasied when the species defined in the boundary + """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) From 68212f331c191541bb59d53c10fe8b704e4bb757 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 14:11:15 -0400 Subject: [PATCH 39/77] update test doc string --- test/test_material.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index 23a780d9d..652783922 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -75,8 +75,8 @@ def test_multispecies_dict_objects(): def test_multispecies_dict_objects_and_strings(): - """Test that the diffusion coefficient is correctly defined when keys are - festim.Species objects""" + """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 From 4fdcf70cdeadceaa4c950160728f3cb37057342f Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 23 Oct 2023 14:29:29 -0400 Subject: [PATCH 40/77] update doc strings --- festim/hydrogen_transport_problem.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 4f0ac625b..01a8185c6 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -124,7 +124,8 @@ 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: # if name of species is given then replace with species object for idx, field in enumerate(export.field): @@ -137,6 +138,8 @@ def defing_export_writers(self): export.writer.write_mesh(self.mesh.mesh) def define_function_space(self): + """Creates the function space of the model, creates a mixed element if + model is multispecies""" element_CG = basix.ufl.element( basix.ElementFamily.P, self.mesh.mesh.basix_cell(), @@ -159,8 +162,8 @@ def define_function_space(self): self.u_n = Function(self.function_space) def assign_functions_to_species(self): - """Creates for each species the solution, prev solution and test - function""" + """Creates the solution, prev solution, test function and + post-processing solution for each species""" if len(self.species) > 1: sub_solutions = list(ufl.split(self.u)) From 337a35862ace8accccb5bc6248a14b695b5baf55 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 23 Oct 2023 14:43:07 -0400 Subject: [PATCH 41/77] removed unused variable --- test/test_material.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_material.py b/test/test_material.py index df8f70c08..53cc3f645 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -99,7 +99,6 @@ def test_D_0_type_rasies_error(): """Test that a value error is rasied in the get_diffusion_coefficient function""" A = F.Species("A") - spe_list = [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"): From d91dd751cfacae8935559862bb1c0a4d29539099 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 24 Oct 2023 14:12:57 -0400 Subject: [PATCH 42/77] updated tests to improve coverage --- festim/material.py | 2 -- test/test_material.py | 25 +++++++++++++++++++++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/festim/material.py b/festim/material.py index 811766f6d..36984e97f 100644 --- a/festim/material.py +++ b/festim/material.py @@ -73,8 +73,6 @@ def get_diffusion_coefficient(self, mesh, temperature, species=None): 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) - else: - raise ValueError(f"{species} is not in E_D keys") return D_0 * ufl.exp(-E_D / F.k_B / temperature) diff --git a/test/test_material.py b/test/test_material.py index 653f4f37a..a7a9dc774 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -103,3 +103,28 @@ def test_D_0_type_rasies_error(): 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_rasied_when_species_not_given_with_dict(): + """Test that a value error is rasied 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_rasied_when_species_not_not_in_D_0_dict(): + """Test that a value error is rasied 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) From a6bef0c64aad5bbcf25fde4d2a029d1db18fa3ba Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 09:52:36 -0400 Subject: [PATCH 43/77] accept callable functions in multispecies --- festim/hydrogen_transport_problem.py | 87 +++++++++++++++----- test/test_multispecies_problem.py | 116 ++++++++++++++++++++++----- test/test_permeation_problem.py | 2 - 3 files changed, 166 insertions(+), 39 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 411d9532a..fc6a7782e 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -173,6 +173,9 @@ def assign_functions_to_species(self): 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() else: sub_solutions = [self.u] @@ -240,18 +243,43 @@ def define_boundary_conditions(self): # 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, bc.species.sub_function_space - ) - bc.create_value( - self.mesh.mesh, - bc.species.sub_function_space, - self.temperature, - self.t, - ) - form = bc.create_formulation( - dofs=bc_dofs, function_space=bc.species.sub_function_space - ) + if (len(self.species) > 1) and ( + not isinstance(bc.value, (int, float, fem.Constant)) + ): + bc_dofs = bc.define_surface_subdomain_dofs( + facet_meshtags=self.facet_meshtags, + mesh=self.mesh, + function_space=( + bc.species.sub_function_space, + bc.species.collapsed_function_space, + ), + ) + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=bc.species.collapsed_function_space, + t=self.t, + ) + else: + bc_dofs = bc.define_surface_subdomain_dofs( + facet_meshtags=self.facet_meshtags, + mesh=self.mesh, + function_space=bc.species.sub_function_space, + ) + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=bc.species.sub_function_space, + t=self.t, + ) + if (len(self.species) == 1) and ( + isinstance(bc.value_fenics, (fem.Function)) + ): + form = bc.create_formulation(dofs=bc_dofs, function_space=None) + else: + form = bc.create_formulation( + dofs=bc_dofs, function_space=bc.species.sub_function_space + ) self.bc_forms.append(form) def create_formulation(self): @@ -307,12 +335,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, self.species[0] - ) - cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( desc="Solving H transport problem", total=self.settings.final_time, @@ -328,14 +352,36 @@ def run(self): self.solver.solve(self.u) - if len(self.species) == 1: + if len(self.species) <= 1: + 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 * dot(grad(cm), n) * self.ds(2)) + 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: + # res = list(self.u.split()) + 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)): @@ -344,4 +390,7 @@ def run(self): # update previous solution self.u_n.x.array[:] = self.u.x.array[:] + if len(self.species) == 2: + flux_values = [flux_values_1, flux_values_2] + return times, flux_values diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 103142e46..1ccc870ff 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -1,35 +1,115 @@ import numpy as np import festim as F +from petsc4py import PETSc +from dolfinx.fem import Constant +from ufl import exp -def test_multispecies_problem_initialisation(tmpdir): - """Test that the multispecies problem is correctly initialised""" +def test_multispecies_permeation_problem(): + L = 3e-04 + my_mesh = F.Mesh1D(np.linspace(0, L, num=1001)) my_model = F.HydrogenTransportProblem() - L = 1e-04 - my_model.mesh = F.Mesh1D(np.linspace(0, L, num=1001)) + my_model.mesh = my_mesh + my_mat = F.Material( - D_0={"H": 1.9e-7, "D": 1.9e-07}, E_D={"H": 0.2, "D": 0.2}, name="my_mat" + D_0={"D": 1.9e-7, "T": 3.8e-7}, + E_D={"D": 0.2, "T": 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] - mobile_H = F.Species("H") + my_model.subdomains = [ + my_subdomain, + left_surface, + right_surface, + ] + mobile_D = F.Species("D") - my_model.species = [mobile_H, mobile_D] - my_model.temperature = 600 + mobile_T = F.Species("T") + my_model.species = [mobile_D, mobile_T] + + temperature = Constant(my_mesh.mesh, 500.0) + my_model.temperature = temperature + my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left_surface, value=1e18, species="H"), - F.DirichletBC(subdomain=right_surface, value=1e18, species="D"), + F.DirichletBC(subdomain=right_surface, value=0, species="D"), + F.DirichletBC(subdomain=right_surface, value=0, species="T"), + F.SievertsBC( + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="D" + ), + F.SievertsBC( + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="T" + ), ] - my_model.settings = F.Settings(atol=1e10, rtol=1e-10, final_time=50) - my_model.settings.stepsize = 0.1 my_model.exports = [ - F.XDMFExport(str(tmpdir.join("my_export_1.xdmf")), field=mobile_H), - F.XDMFExport(str(tmpdir.join("my_export_2.xdmf")), field="D"), - F.VTXExport(str(tmpdir.join("my_export_3.bp")), field="H"), - F.VTXExport(str(tmpdir.join("my_export_4.bp")), field=mobile_D), + F.XDMFExport("mobile_concentration_D.xdmf", field=mobile_D), + F.XDMFExport("mobile_concentration_T.xdmf", field=mobile_T), ] + 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.run() + + 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 solution ----------------------------- + D_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_D) + D_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_T) + S_0 = float(my_model.boundary_conditions[-1].S_0) + E_S = float(my_model.boundary_conditions[-1].E_S) + P_up = float(my_model.boundary_conditions[-1].pressure) + S = S_0 * exp(-E_S / F.k_B / float(temperature)) + permeability_1 = float(D_1) * S + permeability_2 = float(D_2) * S + times = np.array(times) + + n_array_1 = np.arange(1, 10000)[:, np.newaxis] + n_array_2 = np.arange(1, 10000)[:, np.newaxis] + summation_1 = np.sum( + (-1) ** n_array_1 + * np.exp(-((np.pi * n_array_1) ** 2) * float(D_1) / L**2 * times), + axis=0, + ) + summation_2 = np.sum( + (-1) ** n_array_2 + * np.exp(-((np.pi * n_array_2) ** 2) * float(D_2) / L**2 * times), + axis=0, + ) + analytical_flux_1 = P_up**0.5 * permeability_1 / L * (2 * summation_1 + 1) + analytical_flux_2 = P_up**0.5 * permeability_2 / L * (2 * summation_2 + 1) + analytical_flux_1 = np.abs(analytical_flux_1) + analytical_flux_2 = np.abs(analytical_flux_2) + flux_values_1 = np.array(np.abs(flux_values[0])) + flux_values_2 = np.array(np.abs(flux_values[1])) + + indices_1 = np.where(analytical_flux_1 > 0.1 * np.max(analytical_flux_1)) + indices_2 = np.where(analytical_flux_2 > 0.1 * np.max(analytical_flux_2)) + analytical_flux_1 = analytical_flux_1[indices_1] + analytical_flux_2 = analytical_flux_2[indices_2] + flux_values_1 = flux_values_1[indices_1] + flux_values_2 = flux_values_2[indices_2] + + relative_error_1 = np.abs((flux_values_1 - analytical_flux_1) / analytical_flux_1) + relative_error_2 = np.abs((flux_values_2 - analytical_flux_2) / analytical_flux_2) + + error_1 = relative_error_1.mean() + error_2 = relative_error_2.mean() + + for err in [error_1, error_2]: + assert err < 0.01 diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index d69b3dfa6..47bed233b 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -2,8 +2,6 @@ from dolfinx.fem import Constant from ufl import exp import numpy as np - - import festim as F From 4c67bd998a6cf84e9d72ef50d4999b13a061114b Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 09:53:00 -0400 Subject: [PATCH 44/77] test not needed --- festim/boundary_conditions/dirichlet_bc.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index efb4d4dba..cd334c483 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( @@ -124,17 +125,11 @@ def create_formulation(self, dofs, function_space): 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, - ) + form = fem.dirichletbc( + value=self.value_fenics, + dofs=dofs, + V=function_space, + ) return form def update(self, t): From 7ca6763bbc3d2e67553fbbdf96a9ba822839f181 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 09:53:10 -0400 Subject: [PATCH 45/77] test not needed --- test/test_dirichlet_bc.py | 41 --------------------------------------- 1 file changed, 41 deletions(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index eb5e26bc5..065a97b25 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -216,47 +216,6 @@ def test_callable_x_only(): assert np.isclose(computed_value, expected_value) -@pytest.mark.parametrize( - "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 = 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() - 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", [ From 619a0fb214a2deece7b663985ccbe13ba35ea01b Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 10:07:06 -0400 Subject: [PATCH 46/77] updated doc strings --- festim/hydrogen_transport_problem.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index fc6a7782e..c5172441e 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -163,7 +163,8 @@ def define_function_space(self): def assign_functions_to_species(self): """Creates the solution, prev solution, test function and - post-processing solution for each species""" + post-processing solution for each species, if model is multispecies, + created a collapsed function space for each species""" if len(self.species) > 1: sub_solutions = list(ufl.split(self.u)) From 0652bf1e7870d1db7120e3e72e0626a374cb3ead Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 10:15:59 -0400 Subject: [PATCH 47/77] dont export results locally --- test/test_multispecies_problem.py | 11 ++++++++--- test/test_permeation_problem.py | 15 +++++++++++---- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 1ccc870ff..8f39f2520 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -3,9 +3,10 @@ from petsc4py import PETSc from dolfinx.fem import Constant from ufl import exp +import os -def test_multispecies_permeation_problem(): +def test_multispecies_permeation_problem(tmp_path): L = 3e-04 my_mesh = F.Mesh1D(np.linspace(0, L, num=1001)) my_model = F.HydrogenTransportProblem() @@ -43,8 +44,12 @@ def test_multispecies_permeation_problem(): ), ] my_model.exports = [ - F.XDMFExport("mobile_concentration_D.xdmf", field=mobile_D), - F.XDMFExport("mobile_concentration_T.xdmf", field=mobile_T), + F.XDMFExport( + os.path.join(tmp_path, "mobile_concentration_D.xdmf"), field=mobile_D + ), + F.XDMFExport( + os.path.join(tmp_path, "mobile_concentration_T.xdmf"), field=mobile_T + ), ] my_model.settings = F.Settings( atol=1e10, diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 47bed233b..3244edb9e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -3,9 +3,10 @@ from ufl import exp import numpy as np import festim as F +import os -def test_permeation_problem(mesh_size=1001): +def test_permeation_problem(tmp_path, mesh_size=1001): L = 3e-04 vertices = np.linspace(0, L, num=mesh_size) @@ -31,7 +32,11 @@ def test_permeation_problem(mesh_size=1001): subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H" ), ] - my_model.exports = [F.XDMFExport("mobile_concentration.xdmf", field=mobile_H)] + my_model.exports = [ + F.XDMFExport( + os.path.join(tmp_path, "mobile_concentration_H.xdmf"), field=mobile_H + ) + ] my_model.settings = F.Settings( atol=1e10, @@ -89,7 +94,7 @@ def test_permeation_problem(mesh_size=1001): assert error < 0.01 -def test_permeation_problem_multi_volume(): +def test_permeation_problem_multi_volume(tmp_path): L = 3e-04 vertices = np.linspace(0, L, num=1001) @@ -128,7 +133,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, From 68fc9ddc03fc087fcf963afa0692550c4db3a81a Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 10:32:51 -0400 Subject: [PATCH 48/77] dont use tmppath in benchmark case --- test/test_permeation_problem.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 3244edb9e..55e532fd2 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -6,7 +6,7 @@ import os -def test_permeation_problem(tmp_path, mesh_size=1001): +def test_permeation_problem(mesh_size=1001): L = 3e-04 vertices = np.linspace(0, L, num=mesh_size) @@ -32,11 +32,7 @@ def test_permeation_problem(tmp_path, mesh_size=1001): subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H" ), ] - my_model.exports = [ - F.XDMFExport( - os.path.join(tmp_path, "mobile_concentration_H.xdmf"), field=mobile_H - ) - ] + my_model.exports = [F.XDMFExport("mobile_concentration.xdmf", field=mobile_H)] my_model.settings = F.Settings( atol=1e10, From 487ea24e10f1f37b1b089e49f08bddab84e0c0fd Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Wed, 25 Oct 2023 10:39:51 -0400 Subject: [PATCH 49/77] update define_function_space doc strings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index c5172441e..471e42e4d 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -139,7 +139,7 @@ def defing_export_writers(self): def define_function_space(self): """Creates the function space of the model, creates a mixed element if - model is multispecies""" + 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(), From eab51811c7c1195836d21a2c427463f361bd3979 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 10:52:42 -0400 Subject: [PATCH 50/77] facet indices not dofs --- festim/hydrogen_transport_problem.py | 15 +++++++++------ festim/subdomain/surface_subdomain.py | 12 ++++++------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 471e42e4d..a9a9447ac 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -139,7 +139,8 @@ def defing_export_writers(self): def define_function_space(self): """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.""" + 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(), @@ -193,7 +194,7 @@ def assign_functions_to_species(self): 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 @@ -202,8 +203,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.mesh.mesh, self.mesh.fdim) - 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 @@ -218,12 +221,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 diff --git a/festim/subdomain/surface_subdomain.py b/festim/subdomain/surface_subdomain.py index f1ee32839..80fd082aa 100644 --- a/festim/subdomain/surface_subdomain.py +++ b/festim/subdomain/surface_subdomain.py @@ -1,4 +1,3 @@ -from dolfinx import fem import dolfinx.mesh import numpy as np @@ -23,18 +22,19 @@ def __init__(self, id, x) -> None: self.id = id self.x = x - def locate_dof(self, mesh, fdim): + 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: 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 = dolfinx.mesh.locate_entities_boundary( + indices = dolfinx.mesh.locate_entities_boundary( mesh, fdim, lambda x: np.isclose(x[0], self.x) ) - return dofs[0] + return indices[0] From 5eb7ebf355910be618ad0afb450d1d32cd81b9a4 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 10:58:22 -0400 Subject: [PATCH 51/77] species list needs to be more than 0 --- test/test_subdomains.py | 1 + 1 file changed, 1 insertion(+) 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)) From 33fc5b63da222b73e78ef4a7fce43b8b1416d42d Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 11:08:49 -0400 Subject: [PATCH 52/77] more readable for multispecies --- festim/hydrogen_transport_problem.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index a9a9447ac..5bd5cb5f5 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -147,15 +147,15 @@ def define_function_space(self): 1, basix.LagrangeVariant.equispaced, ) - if len(self.species) > 1: + if len(self.species) == 1: + 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) - else: - element = element_CG self.function_space = fem.FunctionSpace(self.mesh.mesh, element) @@ -167,7 +167,13 @@ def assign_functions_to_species(self): post-processing solution for each species, if model is multispecies, created a collapsed function space for each species""" - if len(self.species) > 1: + if len(self.species) == 1: + 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)) @@ -179,13 +185,6 @@ def assign_functions_to_species(self): idx ).collapse() - else: - 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) - for idx, spe in enumerate(self.species): spe.solution = sub_solutions[idx] spe.prev_solution = sub_prev_solution[idx] @@ -356,7 +355,7 @@ def run(self): self.solver.solve(self.u) - if len(self.species) <= 1: + if len(self.species) == 1: D_D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature, self.species[0] ) From 5fe2b0dc91d2e4a329008f58dad121130eeffe92 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 11:12:52 -0400 Subject: [PATCH 53/77] typo --- test/test_material.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index a7a9dc774..265bd9afb 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -86,7 +86,7 @@ def test_multispecies_dict_objects_and_strings(): def test_multispecies_dict_different_keys(): - """Test that a value error is rasied when the keys of the D_0 and E_D + """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}) @@ -95,8 +95,8 @@ def test_multispecies_dict_different_keys(): my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) -def test_D_0_type_rasies_error(): - """Test that a value error is rasied in the get_diffusion_coefficient +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) @@ -105,8 +105,8 @@ def test_D_0_type_rasies_error(): my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=A) -def test_error_rasied_when_species_not_given_with_dict(): - """Test that a value error is rasied when a species has not been given in +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") @@ -118,8 +118,8 @@ def test_error_rasied_when_species_not_given_with_dict(): my_mat.get_diffusion_coefficient(test_mesh.mesh, 500) -def test_error_rasied_when_species_not_not_in_D_0_dict(): - """Test that a value error is rasied when a species has not been given but +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") From 39d50017f2378f236322da6c085797dc01a309af Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 11:13:12 -0400 Subject: [PATCH 54/77] not needed --- festim/hydrogen_transport_problem.py | 1 - 1 file changed, 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 5bd5cb5f5..a90f5eb3d 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -367,7 +367,6 @@ def run(self): flux_values.append(flux) times.append(float(self.t)) else: - # res = list(self.u.split()) for idx, spe in enumerate(self.species): spe.post_processing_solution = self.u.sub(idx) From 2515c8d7a0055673276cad9d3c63e289a70f35b9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:33:35 -0400 Subject: [PATCH 55/77] split methods --- festim/hydrogen_transport_problem.py | 75 ++++++++++++++-------------- 1 file changed, 38 insertions(+), 37 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index a90f5eb3d..cf0d91fcd 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -246,45 +246,46 @@ def define_boundary_conditions(self): # 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): - if (len(self.species) > 1) and ( - not isinstance(bc.value, (int, float, fem.Constant)) - ): - bc_dofs = bc.define_surface_subdomain_dofs( - facet_meshtags=self.facet_meshtags, - mesh=self.mesh, - function_space=( - bc.species.sub_function_space, - bc.species.collapsed_function_space, - ), - ) - bc.create_value( - mesh=self.mesh.mesh, - temperature=self.temperature, - function_space=bc.species.collapsed_function_space, - t=self.t, - ) - else: - bc_dofs = bc.define_surface_subdomain_dofs( - facet_meshtags=self.facet_meshtags, - mesh=self.mesh, - function_space=bc.species.sub_function_space, - ) - bc.create_value( - mesh=self.mesh.mesh, - temperature=self.temperature, - function_space=bc.species.sub_function_space, - t=self.t, - ) - if (len(self.species) == 1) and ( - isinstance(bc.value_fenics, (fem.Function)) - ): - form = bc.create_formulation(dofs=bc_dofs, function_space=None) - else: - form = bc.create_formulation( - dofs=bc_dofs, function_space=bc.species.sub_function_space - ) + form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) + def create_dirichletbc_form(self, bc): + if (len(self.species) > 1) and ( + not isinstance(bc.value, (int, float, fem.Constant)) + ): + bc_dofs = bc.define_surface_subdomain_dofs( + facet_meshtags=self.facet_meshtags, + mesh=self.mesh, + function_space=( + bc.species.sub_function_space, + bc.species.collapsed_function_space, + ), + ) + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=bc.species.collapsed_function_space, + t=self.t, + ) + else: + bc_dofs = bc.define_surface_subdomain_dofs( + facet_meshtags=self.facet_meshtags, + mesh=self.mesh, + function_space=bc.species.sub_function_space, + ) + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=bc.species.sub_function_space, + t=self.t, + ) + if (len(self.species) == 1) and (isinstance(bc.value_fenics, (fem.Function))): + return bc.create_formulation(dofs=bc_dofs, function_space=None) + else: + return bc.create_formulation( + dofs=bc_dofs, function_space=bc.species.sub_function_space + ) + def create_formulation(self): """Creates the formulation of the model""" if len(self.sources) > 1: From 80c7ec6634e28fcab5f9bd8222fdaa120f4c1a26 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:35:16 -0400 Subject: [PATCH 56/77] refactoring --- festim/hydrogen_transport_problem.py | 43 ++++++++++++---------------- 1 file changed, 19 insertions(+), 24 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index cf0d91fcd..85ad6f4b4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -253,32 +253,27 @@ def create_dirichletbc_form(self, bc): if (len(self.species) > 1) and ( not isinstance(bc.value, (int, float, fem.Constant)) ): - bc_dofs = bc.define_surface_subdomain_dofs( - facet_meshtags=self.facet_meshtags, - mesh=self.mesh, - function_space=( - bc.species.sub_function_space, - bc.species.collapsed_function_space, - ), - ) - bc.create_value( - mesh=self.mesh.mesh, - temperature=self.temperature, - function_space=bc.species.collapsed_function_space, - t=self.t, + function_space_dofs = ( + bc.species.sub_function_space, + bc.species.collapsed_function_space, ) + function_space_value = bc.species.collapsed_function_space else: - bc_dofs = bc.define_surface_subdomain_dofs( - facet_meshtags=self.facet_meshtags, - mesh=self.mesh, - function_space=bc.species.sub_function_space, - ) - bc.create_value( - mesh=self.mesh.mesh, - temperature=self.temperature, - function_space=bc.species.sub_function_space, - t=self.t, - ) + function_space_dofs = bc.species.sub_function_space + function_space_value = 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, + ) + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=function_space_value, + t=self.t, + ) + if (len(self.species) == 1) and (isinstance(bc.value_fenics, (fem.Function))): return bc.create_formulation(dofs=bc_dofs, function_space=None) else: From 03b5b7364301aac236ab2f2489674f79eaf648b9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:38:09 -0400 Subject: [PATCH 57/77] callable --- festim/hydrogen_transport_problem.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 85ad6f4b4..1237992a5 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -250,9 +250,7 @@ def define_boundary_conditions(self): self.bc_forms.append(form) def create_dirichletbc_form(self, bc): - if (len(self.species) > 1) and ( - not isinstance(bc.value, (int, float, fem.Constant)) - ): + if len(self.species) > 1 and callable(bc.value): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, From 5687d5bb2b640425bd8a0a270cda44d550354954 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:47:37 -0400 Subject: [PATCH 58/77] separated things a bit --- festim/hydrogen_transport_problem.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 1237992a5..aee3319dc 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -250,28 +250,38 @@ def define_boundary_conditions(self): self.bc_forms.append(form) def create_dirichletbc_form(self, bc): + # create value_fenics + if callable(bc.value): + if len(self.species) == 1: + function_space_value = bc.species.sub_function_space + else: + function_space_value = bc.species.collapsed_function_space + else: + function_space_value = None + + bc.create_value( + mesh=self.mesh.mesh, + temperature=self.temperature, + function_space=function_space_value, + t=self.t, + ) + + # get dofs if len(self.species) > 1 and callable(bc.value): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, ) - function_space_value = bc.species.collapsed_function_space else: function_space_dofs = bc.species.sub_function_space - function_space_value = 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, ) - bc.create_value( - mesh=self.mesh.mesh, - temperature=self.temperature, - function_space=function_space_value, - t=self.t, - ) + # create form if (len(self.species) == 1) and (isinstance(bc.value_fenics, (fem.Function))): return bc.create_formulation(dofs=bc_dofs, function_space=None) else: From 2a81986d945b375ae8c3f80499c94b44285bbc24 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:50:26 -0400 Subject: [PATCH 59/77] removed DirichletBC.create_form method --- festim/boundary_conditions/dirichlet_bc.py | 13 ------------- festim/hydrogen_transport_problem.py | 20 +++++++++++++++++--- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index cd334c483..aea5c2753 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -119,19 +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 - """ - 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/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index aee3319dc..00fa90bf4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -250,6 +250,14 @@ def define_boundary_conditions(self): self.bc_forms.append(form) def create_dirichletbc_form(self, bc): + """Creates a dirichlet boundary condition form + + Args: + bc (festim.DirichletBC): _description_ + + Returns: + dolfinx.fem.bcs.DirichletBC: _description_ + """ # create value_fenics if callable(bc.value): if len(self.species) == 1: @@ -283,11 +291,17 @@ def create_dirichletbc_form(self, bc): # create form if (len(self.species) == 1) and (isinstance(bc.value_fenics, (fem.Function))): - return bc.create_formulation(dofs=bc_dofs, function_space=None) + form = fem.dirichletbc( + value=bc.value_fenics, + dofs=bc_dofs, + ) else: - return bc.create_formulation( - dofs=bc_dofs, function_space=bc.species.sub_function_space + form = fem.dirichletbc( + value=bc.value_fenics, + dofs=bc_dofs, + V=bc.species.sub_function_space, ) + return form def create_formulation(self): """Creates the formulation of the model""" From d31c5a0f1704e4dca04852cd44db58d417b65a6c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 11:58:42 -0400 Subject: [PATCH 60/77] correct type for t --- festim/boundary_conditions/dirichlet_bc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index aea5c2753..6b31ab22e 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -99,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) From 69edb78671f0d447202567093fc3f9daa32b83e9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 12:01:28 -0400 Subject: [PATCH 61/77] added a test that catches the previous bug --- test/test_dirichlet_bc.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 065a97b25..c73742c15 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -130,12 +130,13 @@ 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 + # value = lambda t: 1.0 + t species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) From 93160375926a0d49482c49cc3c32d106c7c3a9ac Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 12:01:40 -0400 Subject: [PATCH 62/77] correct condition for test --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 00fa90bf4..6bff602e9 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -275,7 +275,7 @@ def create_dirichletbc_form(self, bc): ) # get dofs - if len(self.species) > 1 and callable(bc.value): + if len(self.species) > 1 and isinstance(bc.value_fenics, (fem.Function)): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, From f4d017fde50ce5a28d24ab03c274b394e46d9b1c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 25 Oct 2023 12:13:34 -0400 Subject: [PATCH 63/77] documentation --- festim/hydrogen_transport_problem.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 6bff602e9..3a45adec7 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -253,19 +253,22 @@ def create_dirichletbc_form(self, bc): """Creates a dirichlet boundary condition form Args: - bc (festim.DirichletBC): _description_ + bc (festim.DirichletBC): the boundary condition Returns: - dolfinx.fem.bcs.DirichletBC: _description_ + 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 len(self.species) == 1: function_space_value = bc.species.sub_function_space else: function_space_value = bc.species.collapsed_function_space - else: - function_space_value = None bc.create_value( mesh=self.mesh.mesh, @@ -290,10 +293,11 @@ def create_dirichletbc_form(self, bc): ) # create form - if (len(self.species) == 1) and (isinstance(bc.value_fenics, (fem.Function))): + if len(self.species) == 1 and isinstance(bc.value_fenics, (fem.Function)): form = fem.dirichletbc( value=bc.value_fenics, dofs=bc_dofs, + # no need to pass the functionspace since value_fenics is already a Function ) else: form = fem.dirichletbc( From d254c871fb8c692b6391bddff7c419ba9b73a5b0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:11:37 -0400 Subject: [PATCH 64/77] refactoring for clarity --- test/test_multispecies_problem.py | 124 +++++++++++++++++------------- 1 file changed, 71 insertions(+), 53 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 8f39f2520..e2c05ae3e 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -6,15 +6,20 @@ import os -def test_multispecies_permeation_problem(tmp_path): +def test_multispecies_permeation_problem(): + """Test running a permeation problem with two species permeating from the + same side of a 1D domain with different diffusion coefficients, asserting + that both species match their respecitive analytcial solutions""" + + # 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={"D": 1.9e-7, "T": 3.8e-7}, - E_D={"D": 0.2, "T": 0.2}, + 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) @@ -26,29 +31,21 @@ def test_multispecies_permeation_problem(tmp_path): right_surface, ] - mobile_D = F.Species("D") - mobile_T = F.Species("T") - my_model.species = [mobile_D, mobile_T] + 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 my_model.boundary_conditions = [ - F.DirichletBC(subdomain=right_surface, value=0, species="D"), - F.DirichletBC(subdomain=right_surface, value=0, species="T"), + 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=100, species="D" + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_1" ), F.SievertsBC( - subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="T" - ), - ] - my_model.exports = [ - F.XDMFExport( - os.path.join(tmp_path, "mobile_concentration_D.xdmf"), field=mobile_D - ), - F.XDMFExport( - os.path.join(tmp_path, "mobile_concentration_T.xdmf"), field=mobile_T + subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_2" ), ] my_model.settings = F.Settings( @@ -59,7 +56,6 @@ def test_multispecies_permeation_problem(tmp_path): ) my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) - my_model.initialise() my_model.solver.convergence_criterion = "incremental" @@ -73,48 +69,70 @@ def test_multispecies_permeation_problem(tmp_path): times, flux_values = my_model.run() - # ---------------------- analytical solution ----------------------------- - D_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_D) - D_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, mobile_T) + # ---------------------- analytical solutions ----------------------------- + D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) + D_spe_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_2) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) P_up = float(my_model.boundary_conditions[-1].pressure) S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability_1 = float(D_1) * S - permeability_2 = float(D_2) * S + permeability_spe_1 = float(D_spe_1) * S + permeability_spe_2 = float(D_spe_2) * S times = np.array(times) - n_array_1 = np.arange(1, 10000)[:, np.newaxis] - n_array_2 = np.arange(1, 10000)[:, np.newaxis] - summation_1 = np.sum( - (-1) ** n_array_1 - * np.exp(-((np.pi * n_array_1) ** 2) * float(D_1) / L**2 * times), + # ##### compute analyical solution for species 1 ##### # + n_array = np.arange(1, 10000)[:, np.newaxis] + + summation_spe_1 = np.sum( + (-1) ** n_array + * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_1) / L**2 * times), axis=0, ) - summation_2 = np.sum( - (-1) ** n_array_2 - * np.exp(-((np.pi * n_array_2) ** 2) * float(D_2) / L**2 * times), + analytical_flux_spe_1 = ( + P_up**0.5 * permeability_spe_1 / L * (2 * summation_spe_1 + 1) + ) + + # post processing + analytical_flux_spe_1 = np.abs(analytical_flux_spe_1) + flux_values_spe_1 = np.array(np.abs(flux_values[0])) + indices_spe_1 = np.where( + analytical_flux_spe_1 > 0.1 * np.max(analytical_flux_spe_1) + ) + + analytical_flux_spe_1 = analytical_flux_spe_1[indices_spe_1] + flux_values_spe_1 = flux_values_spe_1[indices_spe_1] + + # evaluate relative error compared to analytical solution + relative_error_spe_1 = np.abs( + (flux_values_spe_1 - analytical_flux_spe_1) / analytical_flux_spe_1 + ) + error_spe_1 = relative_error_spe_1.mean() + + # ##### compute analyical solution for species 2 ##### # + summation_spe_2 = np.sum( + (-1) ** n_array + * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_2) / L**2 * times), axis=0, ) - analytical_flux_1 = P_up**0.5 * permeability_1 / L * (2 * summation_1 + 1) - analytical_flux_2 = P_up**0.5 * permeability_2 / L * (2 * summation_2 + 1) - analytical_flux_1 = np.abs(analytical_flux_1) - analytical_flux_2 = np.abs(analytical_flux_2) - flux_values_1 = np.array(np.abs(flux_values[0])) - flux_values_2 = np.array(np.abs(flux_values[1])) - - indices_1 = np.where(analytical_flux_1 > 0.1 * np.max(analytical_flux_1)) - indices_2 = np.where(analytical_flux_2 > 0.1 * np.max(analytical_flux_2)) - analytical_flux_1 = analytical_flux_1[indices_1] - analytical_flux_2 = analytical_flux_2[indices_2] - flux_values_1 = flux_values_1[indices_1] - flux_values_2 = flux_values_2[indices_2] - - relative_error_1 = np.abs((flux_values_1 - analytical_flux_1) / analytical_flux_1) - relative_error_2 = np.abs((flux_values_2 - analytical_flux_2) / analytical_flux_2) - - error_1 = relative_error_1.mean() - error_2 = relative_error_2.mean() - - for err in [error_1, error_2]: + analytical_flux_spe_2 = ( + P_up**0.5 * permeability_spe_2 / L * (2 * summation_spe_2 + 1) + ) + + # post processing + analytical_flux_spe_2 = np.abs(analytical_flux_spe_2) + flux_values_spe_2 = np.array(np.abs(flux_values[1])) + indices_spe_2 = np.where( + analytical_flux_spe_2 > 0.1 * np.max(analytical_flux_spe_2) + ) + + analytical_flux_spe_2 = analytical_flux_spe_2[indices_spe_2] + flux_values_spe_2 = flux_values_spe_2[indices_spe_2] + + # evaluate relative error compared to analytical solution + relative_error_spe_2 = np.abs( + (flux_values_spe_2 - analytical_flux_spe_2) / analytical_flux_spe_2 + ) + error_spe_2 = relative_error_spe_2.mean() + + for err in [error_spe_1, error_spe_2]: assert err < 0.01 From a5b44483e6cd061381e272c7e793f628349e9fb8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:11:53 -0400 Subject: [PATCH 65/77] post processing same as festim test --- test/benchmark.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) 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 From 2ca1b1b91d89557baff49c1b321126d181910684 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:28:30 -0400 Subject: [PATCH 66/77] comments and doc strings in permeation tests --- test/test_multispecies_problem.py | 7 ++++--- test/test_permeation_problem.py | 15 ++++++++++++--- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index e2c05ae3e..f7153a881 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -7,9 +7,10 @@ def test_multispecies_permeation_problem(): - """Test running a permeation problem with two species permeating from the - same side of a 1D domain with different diffusion coefficients, asserting - that both species match their respecitive analytcial solutions""" + """Test running a problem with 2 mobile species permeating through a 1D + 0.3mm domain, with different diffusion coefficients, asserting that the + resulting concentration fields are less than 1% different from their + respecitive analytical solutions""" # festim model L = 3e-04 diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 55e532fd2..0898e412a 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -7,6 +7,11 @@ def test_permeation_problem(mesh_size=1001): + """Test running a problem with a mobile species permeating through a 1D 0.3mm domain + asserting that the resulting concentration field is less than 1% different from a + respecitive analytical solution""" + + # festim model L = 3e-04 vertices = np.linspace(0, L, num=mesh_size) @@ -76,21 +81,25 @@ def test_permeation_problem(mesh_size=1001): ) analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) + # post processing 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] + # evaulate relative error compared to analytical solution relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - error = relative_error.mean() assert error < 0.01 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) @@ -179,8 +188,8 @@ def test_permeation_problem_multi_volume(tmp_path): analytical_flux = analytical_flux[indices] flux_values = flux_values[indices] + # evaulate relative error compared to analytical solution relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - error = relative_error.mean() assert error < 0.01 From 03911010c3aa75eb24339c37c86d32c041e39b52 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:30:00 -0400 Subject: [PATCH 67/77] species arguement null in single species case --- test/test_permeation_problem.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 0898e412a..3ec522846 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -63,9 +63,7 @@ def test_permeation_problem(mesh_size=1001): # -------------------------- analytical solution ------------------------------------- - D = my_mat.get_diffusion_coefficient( - my_mesh.mesh, my_model.temperature, my_model.species[0] - ) + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, my_model.temperature) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) @@ -165,7 +163,7 @@ def test_permeation_problem_multi_volume(tmp_path): times, flux_values = my_model.run() # ---------------------- analytical solution ----------------------------- - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, my_model.species[0]) + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) S_0 = float(my_model.boundary_conditions[-1].S_0) E_S = float(my_model.boundary_conditions[-1].E_S) From d182007039b56e1a8a343a0c0b0849e7d391ed20 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:35:18 -0400 Subject: [PATCH 68/77] updates species doc strings --- festim/species.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/festim/species.py b/festim/species.py index 6cbc3f5a6..d397e6f89 100644 --- a/festim/species.py +++ b/festim/species.py @@ -10,11 +10,17 @@ 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 - sub_function_space (dolfinx.fem.FunctionSpace): the subspace of the function space - post_processing_solution (dolfinx.fem.Function): the solution for post processing + 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.FunctionSpace): the subspace of the + function space + collapsed_function_space (dolfinx.fem.FunctionSpace): the collapsed + function space for a species in the function space + post_processing_solution (dolfinx.fem.Function): the solution for post # + processing concentration (dolfinx.fem.Function): the concentration of the species Usage: @@ -35,6 +41,7 @@ def __init__(self, name: str = None) -> 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})" From 8277fec4ff0421657a6fbb9918d10c29a4d79326 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:44:22 -0400 Subject: [PATCH 69/77] multispecies property --- festim/hydrogen_transport_problem.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3a45adec7..eecf70d61 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -110,6 +110,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() @@ -147,7 +151,7 @@ def define_function_space(self): 1, basix.LagrangeVariant.equispaced, ) - if len(self.species) == 1: + if not self.multispecies: element = element_CG else: elements = [] @@ -167,7 +171,7 @@ def assign_functions_to_species(self): post-processing solution for each species, if model is multispecies, created a collapsed function space for each species""" - if len(self.species) == 1: + if not self.multispecies: sub_solutions = [self.u] sub_prev_solution = [self.u_n] sub_test_functions = [ufl.TestFunction(self.function_space)] @@ -265,7 +269,7 @@ def create_dirichletbc_form(self, bc): if callable(bc.value): # if bc.value is a callable then need to provide a functionspace - if len(self.species) == 1: + if not self.multispecies: function_space_value = bc.species.sub_function_space else: function_space_value = bc.species.collapsed_function_space @@ -278,7 +282,7 @@ def create_dirichletbc_form(self, bc): ) # get dofs - if len(self.species) > 1 and isinstance(bc.value_fenics, (fem.Function)): + if self.multispecies and isinstance(bc.value_fenics, (fem.Function)): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, @@ -293,7 +297,7 @@ def create_dirichletbc_form(self, bc): ) # create form - if len(self.species) == 1 and isinstance(bc.value_fenics, (fem.Function)): + if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)): form = fem.dirichletbc( value=bc.value_fenics, dofs=bc_dofs, @@ -377,7 +381,7 @@ def run(self): self.solver.solve(self.u) - if len(self.species) == 1: + if not self.multispecies: D_D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature, self.species[0] ) @@ -414,7 +418,7 @@ def run(self): # update previous solution self.u_n.x.array[:] = self.u.x.array[:] - if len(self.species) == 2: + if self.multispecies: flux_values = [flux_values_1, flux_values_2] return times, flux_values From 9eb5d655cef6f682034d0b35c4c1ecbe53ecff25 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:56:08 -0400 Subject: [PATCH 70/77] test dirichlet bc with multispecies --- test/test_dirichlet_bc.py | 72 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index c73742c15..2d0feb48b 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -221,6 +221,7 @@ def test_callable_x_only(): "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, @@ -300,3 +301,74 @@ def test_species_predefined(): 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): + 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) From 7ba4d012a1c720058052bd9e05ba2ec7f64b7b76 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 13:57:39 -0400 Subject: [PATCH 71/77] typo --- festim/species.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/species.py b/festim/species.py index d397e6f89..be5e528a6 100644 --- a/festim/species.py +++ b/festim/species.py @@ -19,7 +19,7 @@ class Species: function space collapsed_function_space (dolfinx.fem.FunctionSpace): the collapsed function space for a species in the function space - post_processing_solution (dolfinx.fem.Function): the solution for post # + post_processing_solution (dolfinx.fem.Function): the solution for post processing concentration (dolfinx.fem.Function): the concentration of the species From 335ce16c648035e9edf438ff62a91ed6a3f3f158 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 14:02:40 -0400 Subject: [PATCH 72/77] updated doc strings --- festim/material.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/material.py b/festim/material.py index 36984e97f..ceeced3ff 100644 --- a/festim/material.py +++ b/festim/material.py @@ -20,9 +20,9 @@ class Material: coeficient (eV) name (str): the name of the material - Usage (1 species): + Usage: >>> my_mat = Material(D_0=1.9e-7, E_D=0.2, name="my_mat") - Usage (multispecies): + 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}, From 87e7ab779df6406f9e09bbfef00d5e7f7876d389 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 14:12:37 -0400 Subject: [PATCH 73/77] add doc string to dirichletbc tests --- test/test_dirichlet_bc.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 2d0feb48b..a29dbcaad 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -230,6 +230,8 @@ def test_callable_x_only(): ], ) 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) @@ -316,6 +318,9 @@ def test_species_predefined(): ], ) 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) From ccec32192346f7cf257dcd1ec204f3a5ff1ebdfd Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 14:28:27 -0400 Subject: [PATCH 74/77] refactoring --- festim/hydrogen_transport_problem.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index eecf70d61..6e3a7f7b7 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -268,7 +268,6 @@ def create_dirichletbc_form(self, bc): 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: @@ -298,17 +297,16 @@ def create_dirichletbc_form(self, bc): # create form if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)): - form = fem.dirichletbc( - value=bc.value_fenics, - dofs=bc_dofs, - # no need to pass the functionspace since value_fenics is already a Function - ) + # no need to pass the functionspace since value_fenics is already a Function + function_space_form = None else: - form = fem.dirichletbc( - value=bc.value_fenics, - dofs=bc_dofs, - V=bc.species.sub_function_space, - ) + 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): From 59be9e74931d696a5b642ea0a42a89112a9adc4f Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 14:29:36 -0400 Subject: [PATCH 75/77] refactoring --- festim/hydrogen_transport_problem.py | 1 + 1 file changed, 1 insertion(+) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 6e3a7f7b7..e2aff9a16 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -301,6 +301,7 @@ def create_dirichletbc_form(self, bc): function_space_form = None else: function_space_form = bc.species.sub_function_space + form = fem.dirichletbc( value=bc.value_fenics, dofs=bc_dofs, From 2dd86f813dabca9ef80839b491fea5cd5cfe871e Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Wed, 25 Oct 2023 15:09:58 -0400 Subject: [PATCH 76/77] remove comment MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- test/test_dirichlet_bc.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index a29dbcaad..f7f5b0e05 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -136,7 +136,6 @@ def test_callable_t_only(value): subdomain = F.SurfaceSubdomain1D(1, x=1) vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) - # value = lambda t: 1.0 + t species = F.Species("test") bc = F.DirichletBC(subdomain, value, species) From d5c4feffef2b987dd3d78dec110d5a7783022ea0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 16:17:19 -0400 Subject: [PATCH 77/77] updates based on comments --- festim/hydrogen_transport_problem.py | 3 + festim/species.py | 9 +- test/test_multispecies_problem.py | 120 +++++++++++++-------------- test/test_permeation_problem.py | 68 +++++++-------- 4 files changed, 98 insertions(+), 102 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e2aff9a16..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 @@ -380,6 +381,8 @@ def run(self): self.solver.solve(self.u) + # post processing + # TODO remove this if not self.multispecies: D_D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature, self.species[0] diff --git a/festim/species.py b/festim/species.py index be5e528a6..663436172 100644 --- a/festim/species.py +++ b/festim/species.py @@ -15,10 +15,11 @@ class Species: timestep test_function (ufl.Argument): the testfunction associated with this species - sub_function_space (dolfinx.fem.FunctionSpace): the subspace of the - function space - collapsed_function_space (dolfinx.fem.FunctionSpace): the collapsed - function space for a species in the function space + 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 diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index f7153a881..ebae0dcf5 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -3,14 +3,35 @@ from petsc4py import PETSc from dolfinx.fem import Constant from ufl import exp -import os + + +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 - 0.3mm domain, with different diffusion coefficients, asserting that the - resulting concentration fields are less than 1% different from their - respecitive analytical solutions""" + domain, with different diffusion coefficients, checks that the computed + permeation flux matches the analytical solution""" # festim model L = 3e-04 @@ -38,15 +59,24 @@ def test_multispecies_permeation_problem(): 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=100, species="spe_1" + subdomain=left_surface, + S_0=4.02e21, + E_S=1.04, + pressure=pressure, + species="spe_1", ), F.SievertsBC( - subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="spe_2" + subdomain=left_surface, + S_0=5.0e21, + E_S=1.2, + pressure=pressure, + species="spe_2", ), ] my_model.settings = F.Settings( @@ -71,69 +101,35 @@ def test_multispecies_permeation_problem(): times, flux_values = my_model.run() # ---------------------- analytical solutions ----------------------------- - D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) - D_spe_2 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_2) - S_0 = float(my_model.boundary_conditions[-1].S_0) - E_S = float(my_model.boundary_conditions[-1].E_S) - P_up = float(my_model.boundary_conditions[-1].pressure) - S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability_spe_1 = float(D_spe_1) * S - permeability_spe_2 = float(D_spe_2) * S + + # 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 ##### # - n_array = np.arange(1, 10000)[:, np.newaxis] - - summation_spe_1 = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_1) / L**2 * times), - axis=0, - ) - analytical_flux_spe_1 = ( - P_up**0.5 * permeability_spe_1 / L * (2 * summation_spe_1 + 1) - ) - - # post processing - analytical_flux_spe_1 = np.abs(analytical_flux_spe_1) - flux_values_spe_1 = np.array(np.abs(flux_values[0])) - indices_spe_1 = np.where( - analytical_flux_spe_1 > 0.1 * np.max(analytical_flux_spe_1) - ) - - analytical_flux_spe_1 = analytical_flux_spe_1[indices_spe_1] - flux_values_spe_1 = flux_values_spe_1[indices_spe_1] - - # evaluate relative error compared to analytical solution - relative_error_spe_1 = np.abs( - (flux_values_spe_1 - analytical_flux_spe_1) / analytical_flux_spe_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 ) - error_spe_1 = relative_error_spe_1.mean() # ##### compute analyical solution for species 2 ##### # - summation_spe_2 = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(D_spe_2) / L**2 * times), - axis=0, - ) - analytical_flux_spe_2 = ( - P_up**0.5 * permeability_spe_2 / L * (2 * summation_spe_2 + 1) - ) - - # post processing - analytical_flux_spe_2 = np.abs(analytical_flux_spe_2) - flux_values_spe_2 = np.array(np.abs(flux_values[1])) - indices_spe_2 = np.where( - analytical_flux_spe_2 > 0.1 * np.max(analytical_flux_spe_2) - ) - - analytical_flux_spe_2 = analytical_flux_spe_2[indices_spe_2] - flux_values_spe_2 = flux_values_spe_2[indices_spe_2] - - # evaluate relative error compared to analytical solution - relative_error_spe_2 = np.abs( - (flux_values_spe_2 - analytical_flux_spe_2) / analytical_flux_spe_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 ) - error_spe_2 = relative_error_spe_2.mean() 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 3ec522846..8e2eb014a 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -6,10 +6,33 @@ import os +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 0.3mm domain - asserting that the resulting concentration field is less than 1% different from a - respecitive analytical solution""" + """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 @@ -71,25 +94,11 @@ 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) - - # post processing - 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] - - # evaulate relative error compared to analytical solution - 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 @@ -171,23 +180,10 @@ def test_permeation_problem_multi_volume(tmp_path): 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] - - # evaulate relative error compared to analytical solution - 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