From 56948fb2706d80054bc979d13e2deed2a0f69220 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 25 Oct 2023 15:12:07 -0400 Subject: [PATCH 01/62] time dependent bcs updated --- festim/boundary_conditions/dirichlet_bc.py | 3 ++ festim/hydrogen_transport_problem.py | 9 +++- test/test_time_depenent_objects.py | 56 ++++++++++++++++++++++ 3 files changed, 66 insertions(+), 2 deletions(-) create mode 100644 test/test_time_depenent_objects.py diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index 6b31ab22e..3f56638ea 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -41,6 +41,7 @@ def __init__(self, subdomain, value, species) -> None: self.value_fenics = None self.bc_expr = None + self.t_dependence = False @property def value_fenics(self): @@ -101,11 +102,13 @@ def create_value( self.value_fenics = F.as_fenics_constant( mesh=mesh, value=self.value(t=float(t)) ) + self.t_dependence = True else: self.value_fenics = fem.Function(function_space) kwargs = {} if "t" in arguments: kwargs["t"] = t + self.t_dependence = True if "x" in arguments: kwargs["x"] = x if "T" in arguments: diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e2aff9a16..cb99d3571 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -98,6 +98,7 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] + self.t_dependent_objects = [] @property def temperature(self): @@ -253,6 +254,10 @@ def define_boundary_conditions(self): form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) + for bc in self.boundary_conditions: + if bc.t_dependence: + self.t_dependent_objects.append(bc) + def create_dirichletbc_form(self, bc): """Creates a dirichlet boundary condition form @@ -375,8 +380,8 @@ def run(self): self.t.value += self.dt.value # update boundary conditions - for bc in self.boundary_conditions: - bc.update(float(self.t)) + for object in self.t_dependent_objects: + object.update(float(self.t)) self.solver.solve(self.u) diff --git a/test/test_time_depenent_objects.py b/test/test_time_depenent_objects.py new file mode 100644 index 000000000..5c4cfccbb --- /dev/null +++ b/test/test_time_depenent_objects.py @@ -0,0 +1,56 @@ +from petsc4py import PETSc +from dolfinx.fem import Constant +from ufl import exp +import numpy as np +import festim as F + + +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 + my_model = F.HydrogenTransportProblem() + L = 3e-04 + my_model.mesh = F.Mesh1D(np.linspace(0, L, num=mesh_size)) + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=L) + my_model.subdomains = [my_subdomain, left_surface, right_surface] + + mobile_H = F.Species("H") + my_model.species = [mobile_H] + + my_model.temperature = 500 + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=right_surface, value=0, species="H"), + F.DirichletBC( + subdomain=left_surface, value=lambda t: 1e17 + 1e17 * t, species="H" + ), + ] + my_model.exports = [F.XDMFExport("mobile_concentration.xdmf", 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() From 751841c07fcaf0e2ace0361742dc37aecf1f139a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 26 Oct 2023 13:24:04 -0400 Subject: [PATCH 02/62] dont add to another list, update bcs --- festim/hydrogen_transport_problem.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 4539f6dd5..1382a3739 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -99,7 +99,6 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] - self.t_dependent_objects = [] @property def temperature(self): @@ -255,10 +254,6 @@ def define_boundary_conditions(self): form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) - for bc in self.boundary_conditions: - if bc.t_dependence: - self.t_dependent_objects.append(bc) - def create_dirichletbc_form(self, bc): """Creates a dirichlet boundary condition form @@ -381,8 +376,9 @@ def run(self): self.t.value += self.dt.value # update boundary conditions - for object in self.t_dependent_objects: - object.update(float(self.t)) + for bc in self.boundary_conditions: + if bc.time_dependent: + bc.update(float(self.t)) self.solver.solve(self.u) From d5cd7c0948ed70167829682d3e78600ff13e2b64 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 26 Oct 2023 13:24:13 -0400 Subject: [PATCH 03/62] renaming --- festim/boundary_conditions/dirichlet_bc.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index 3f56638ea..ad11224fe 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -41,7 +41,7 @@ def __init__(self, subdomain, value, species) -> None: self.value_fenics = None self.bc_expr = None - self.t_dependence = False + self.time_dependent = False @property def value_fenics(self): @@ -102,13 +102,13 @@ def create_value( self.value_fenics = F.as_fenics_constant( mesh=mesh, value=self.value(t=float(t)) ) - self.t_dependence = True + self.time_dependent = True else: self.value_fenics = fem.Function(function_space) kwargs = {} if "t" in arguments: kwargs["t"] = t - self.t_dependence = True + self.time_dependent = True if "x" in arguments: kwargs["x"] = x if "T" in arguments: From b87518e1f58a41d923234f8734cd374e6bc24c9e Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 26 Oct 2023 13:25:03 -0400 Subject: [PATCH 04/62] system test not needed --- test/test_time_depenent_objects.py | 56 ------------------------------ 1 file changed, 56 deletions(-) delete mode 100644 test/test_time_depenent_objects.py diff --git a/test/test_time_depenent_objects.py b/test/test_time_depenent_objects.py deleted file mode 100644 index 5c4cfccbb..000000000 --- a/test/test_time_depenent_objects.py +++ /dev/null @@ -1,56 +0,0 @@ -from petsc4py import PETSc -from dolfinx.fem import Constant -from ufl import exp -import numpy as np -import festim as F - - -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 - my_model = F.HydrogenTransportProblem() - L = 3e-04 - my_model.mesh = F.Mesh1D(np.linspace(0, L, num=mesh_size)) - my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) - left_surface = F.SurfaceSubdomain1D(id=1, x=0) - right_surface = F.SurfaceSubdomain1D(id=2, x=L) - my_model.subdomains = [my_subdomain, left_surface, right_surface] - - mobile_H = F.Species("H") - my_model.species = [mobile_H] - - my_model.temperature = 500 - - my_model.boundary_conditions = [ - F.DirichletBC(subdomain=right_surface, value=0, species="H"), - F.DirichletBC( - subdomain=left_surface, value=lambda t: 1e17 + 1e17 * t, species="H" - ), - ] - my_model.exports = [F.XDMFExport("mobile_concentration.xdmf", 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() From 86cd7e76fe2f5dc5bf71084c9783d107f6a36aec Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 26 Oct 2023 13:55:36 -0400 Subject: [PATCH 05/62] test time dependent attribute --- test/test_dirichlet_bc.py | 41 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index f7f5b0e05..d8ab301d7 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -376,3 +376,44 @@ def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B): if isinstance(expected_value, Conditional): expected_value = float(expected_value) assert np.isclose(computed_value, expected_value) + + +@pytest.mark.parametrize( + "value", + [ + 1.0, + lambda t: t, + lambda t: 1.0 + t, + lambda x: 1.0 + x[0], + lambda x, t: 1.0 + x[0] + t, + lambda x, t, T: 1.0 + x[0] + t + T, + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ], +) +def test_time_dependence_attribute(value): + subdomain = F.SurfaceSubdomain1D(1, x=0) + 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, value, "H") + my_model.boundary_conditions = [my_bc] + + 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) + + my_model.initialise() + + if isinstance(value, (float, int)): + assert not my_model.boundary_conditions[0].time_dependent + else: + arguments = value.__code__.co_varnames + if "t" in arguments: + assert my_model.boundary_conditions[0].time_dependent + else: + assert not my_model.boundary_conditions[0].time_dependent From 6874755eae022693b76d294fa1aa37ba5ce7fde1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 26 Oct 2023 13:57:56 -0400 Subject: [PATCH 06/62] update doc strings --- festim/boundary_conditions/dirichlet_bc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index ad11224fe..95a628849 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -24,6 +24,7 @@ class DirichletBC: fenics format bc_expr (fem.Expression): the expression of the boundary condition that is used to update the value_fenics + time_dependent (bool): True if the boundary condition is time dependent Usage: >>> from festim import DirichletBC From 334f190a329b6af51745cba87489499fbbf3585e Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 28 Oct 2023 09:43:08 -0400 Subject: [PATCH 07/62] post processing --- festim/hydrogen_transport_problem.py | 110 ++++++++++++++++++--------- 1 file changed, 73 insertions(+), 37 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ff171bca3..caec672c4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -99,6 +99,7 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] + self.derived_quantities = {} @property def temperature(self): @@ -126,22 +127,29 @@ def initialise(self): self.define_boundary_conditions() self.create_formulation() self.create_solver() - self.defing_export_writers() + self.initialise_exports() - def defing_export_writers(self): + def initialise_exports(self): """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): if isinstance(field, str): - export.field[idx] = F.find_species_from_name(field, self.species) + # export.field[idx] = F.find_species_from_name(field, self.species) + export.field = 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): export.writer.write_mesh(self.mesh.mesh) + if isinstance(export, F.SurfaceFlux): + self.derived_quantities["t(s)"] = [] + self.derived_quantities[ + f"Surface_flux_subdomain_{export.subdomain.id}" + ] = [] + 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 @@ -218,6 +226,7 @@ def define_markers_and_measures(self): entities = sub_dom.locate_subdomain_entities( self.mesh.mesh, self.mesh.vdim ) + print(entities) tags_volumes[entities] = sub_dom.id # check if all borders are defined @@ -382,40 +391,7 @@ 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] - ) - cm = self.u - self.species[0].post_processing_solution = self.u - - surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(float(self.t)) - else: - for idx, spe in enumerate(self.species): - spe.post_processing_solution = self.u.sub(idx) - - cm_1, cm_2 = self.u.split() - D_1 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] - ) - D_2 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[1] - ) - surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) - surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) - flux_1 = assemble_scalar(surface_flux_1) - flux_2 = assemble_scalar(surface_flux_2) - flux_values_1.append(flux_1) - flux_values_2.append(flux_2) - times.append(float(self.t)) - - for export in self.exports: - if isinstance(export, (F.VTXExport, F.XDMFExport)): - export.write(float(self.t)) + self.post_processing(times, flux_values, flux_values_1, flux_values_2) # update previous solution self.u_n.x.array[:] = self.u.x.array[:] @@ -424,3 +400,63 @@ def run(self): flux_values = [flux_values_1, flux_values_2] return times, flux_values + + def post_processing(self, times, flux_values, flux_values_1, flux_values_2): + if not self.multispecies: + D_D = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + cm = self.u + self.species[0].post_processing_solution = self.u + + surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(float(self.t)) + else: + for idx, spe in enumerate(self.species): + spe.post_processing_solution = self.u.sub(idx) + + cm_1, cm_2 = self.u.split() + D_1 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[0] + ) + D_2 = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, self.species[1] + ) + surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) + surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) + flux_1 = assemble_scalar(surface_flux_1) + flux_2 = assemble_scalar(surface_flux_2) + flux_values_1.append(flux_1) + flux_values_2.append(flux_2) + times.append(float(self.t)) + + # for export in self.exports: + # if isinstance(export, F.SurfaceQuantity): + # export.compute(self.mesh, self.ds) + # if isinstance(export, F.VolumeQuantity): + # export.compute(self.mesh, self.dx) + + for export in self.exports: + # if isinstance(export, F.SurfaceFlux): + # self.derived_quantities["t(s)"].append(float(self.t)) + # for vol in self.volume_subdomains: + # if export.subdomain.x == vol.borders[0] or vol.borders[1]: + # D_export = vol.material.get_diffusion_coefficient( + # self.mesh.mesh, self.temperature, export.field + # ) + # if not hasattr(export.subdomain, "x"): + # raise NotImplementedError("only works for 1D") + # export_value = export.compute_quantity( + # D_export, + # self.mesh, + # self.ds, + # ) + # self.derived_quantities[ + # f"Surface_flux_subdomain_{export.subdomain.id}" + # ].append(export_value) + # export.write(t=float(self.t)) + + if isinstance(export, (F.VTXExport, F.XDMFExport)): + export.write(float(self.t)) From 0146875de6daf3d8b6d18416be1688e6177f696d Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 10:32:07 -0400 Subject: [PATCH 08/62] simpler method --- festim/hydrogen_transport_problem.py | 51 ++++++++++++++-------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index caec672c4..83a1059c2 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -134,10 +134,14 @@ def initialise_exports(self): 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): - if isinstance(field, str): - # export.field[idx] = F.find_species_from_name(field, self.species) - export.field = F.find_species_from_name(field, self.species) + if 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 + ) + elif isinstance(export.field, str): + export.field = F.find_species_from_name(export.field, self.species) if isinstance(export, (F.VTXExport, F.XDMFExport)): export.define_writer(MPI.COMM_WORLD) @@ -146,9 +150,7 @@ def initialise_exports(self): if isinstance(export, F.SurfaceFlux): self.derived_quantities["t(s)"] = [] - self.derived_quantities[ - f"Surface_flux_subdomain_{export.subdomain.id}" - ] = [] + self.derived_quantities[f"{export.title}"] = [] def define_function_space(self): """Creates the function space of the model, creates a mixed element if @@ -226,7 +228,6 @@ def define_markers_and_measures(self): entities = sub_dom.locate_subdomain_entities( self.mesh.mesh, self.mesh.vdim ) - print(entities) tags_volumes[entities] = sub_dom.id # check if all borders are defined @@ -439,24 +440,22 @@ def post_processing(self, times, flux_values, flux_values_1, flux_values_2): # export.compute(self.mesh, self.dx) for export in self.exports: - # if isinstance(export, F.SurfaceFlux): - # self.derived_quantities["t(s)"].append(float(self.t)) - # for vol in self.volume_subdomains: - # if export.subdomain.x == vol.borders[0] or vol.borders[1]: - # D_export = vol.material.get_diffusion_coefficient( - # self.mesh.mesh, self.temperature, export.field - # ) - # if not hasattr(export.subdomain, "x"): - # raise NotImplementedError("only works for 1D") - # export_value = export.compute_quantity( - # D_export, - # self.mesh, - # self.ds, - # ) - # self.derived_quantities[ - # f"Surface_flux_subdomain_{export.subdomain.id}" - # ].append(export_value) - # export.write(t=float(self.t)) + if isinstance(export, F.SurfaceFlux): + # evaluate value of export + D_export = export.volume_subdomain.material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature, export.field + ) + export_value = export.compute_quantity( + D_export, + self.mesh, + self.ds, + ) + + # update derived quantities dict + self.derived_quantities["t(s)"].append(float(self.t)) + self.derived_quantities[f"{export.title}"].append(export_value) + + export.write(t=float(self.t)) if isinstance(export, (F.VTXExport, F.XDMFExport)): export.write(float(self.t)) From 176f42ddcddb78703d985f784571a3cec1ab66fc Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 31 Oct 2023 10:32:19 -0400 Subject: [PATCH 09/62] import surface quantity and flux --- festim/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/festim/__init__.py b/festim/__init__.py index 502168be4..bbe6f48da 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -37,6 +37,8 @@ from .stepsize import Stepsize +from .exports.surface_quantity import SurfaceQuantity +from .exports.surface_flux import SurfaceFlux from .exports.vtx import VTXExport from .exports.xdmf import XDMFExport From 1d031a302b1458f61d0ec4b335517da061f3f86e Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 11:39:25 -0400 Subject: [PATCH 10/62] surface flux class --- festim/exports/surface_flux.py | 88 ++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 festim/exports/surface_flux.py diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py new file mode 100644 index 000000000..bbcad76ce --- /dev/null +++ b/festim/exports/surface_flux.py @@ -0,0 +1,88 @@ +from dolfinx import fem +import festim as F +import ufl +import csv + + +class SurfaceFlux(F.SurfaceQuantity): + """Exports surface flux at a given subdomain + + Args: + + Attributes: + + Usage: + """ + + def __init__( + self, + field, + surface_subdomain: int or F.SurfaceSubdomain1D, + filename: str = None, + ) -> None: + self.field = field + self.surface_subdomain = surface_subdomain + self.filename = filename + + self.D = None + + if self.filename is None: + self.filename = f"Surface_Flux_subdomain_{surface_subdomain.id}.csv" + + self.title = "Flux surface {}: {}".format( + self.surface_subdomain.id, self.field.name + ) + + with open(self.filename, mode="w", newline="") as file: + writer = csv.writer(file) + writer.writerow(["Time", f"{self.title}"]) + + @property + def filename(self): + return self._filename + + @filename.setter + def filename(self, value): + if not isinstance(value, str): + raise TypeError("filename must be of type str") + if not value.endswith(".csv"): + raise ValueError("filename must end with .csv") + self._filename = value + + @property + def surface_subdomain(self): + return self._surface_subdomain + + @surface_subdomain.setter + def surface_subdomain(self, value): + if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance( + value, bool + ): + raise TypeError("surface should be an int or F.SurfaceSubdomain1D") + self._surface_subdomain = value + + @property + def field(self): + return self._field + + @field.setter + def field(self, value): + # check that field is festim.Species + if not isinstance(value, (F.Species, str)): + raise TypeError("field must be of type festim.Species") + + self._field = value + + def compute(self, mesh, ds): + self.value = fem.assemble_scalar( + fem.form( + self.D + * ufl.dot(ufl.grad(self.field.solution), mesh.n) + * ds(self.surface_subdomain.id) + ) + ) + + def write(self, t): + with open(self.filename, mode="a", newline="") as file: + writer = csv.writer(file) + writer.writerow([t, self.value]) From 183405474fb4ca7bfc41ba5b95e3150ed40b6a00 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 11:39:39 -0400 Subject: [PATCH 11/62] use D global for exports --- festim/hydrogen_transport_problem.py | 94 +++++++++++++++------------- 1 file changed, 51 insertions(+), 43 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 83a1059c2..b0d9a887f 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -100,6 +100,8 @@ def __init__( self.volume_subdomains = [] self.bc_forms = [] self.derived_quantities = {} + self.D_global = [] + self.D_global_expr = [] @property def temperature(self): @@ -117,7 +119,7 @@ def multispecies(self): return len(self.species) > 1 def initialise(self): - self.define_function_space() + self.define_function_spaces() self.define_markers_and_measures() self.assign_functions_to_species() @@ -149,10 +151,44 @@ def initialise_exports(self): export.writer.write_mesh(self.mesh.mesh) if isinstance(export, F.SurfaceFlux): + self.define_D_global(export) self.derived_quantities["t(s)"] = [] self.derived_quantities[f"{export.title}"] = [] - def define_function_space(self): + def define_D_global(self, export): + for spe in self.species: + D_0 = fem.Function(self.V_DG_0) + E_D = fem.Function(self.V_DG_0) + for vol in self.volume_subdomains: + entities = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) + if not self.multispecies: + if isinstance(vol.material.D_0, (float, int, fem.Constant)): + D_0.x.array[entities] = vol.material.D_0 + E_D.x.array[entities] = vol.material.E_D + else: + raise NotImplementedError( + "diffusion values as functions not supported" + ) + else: + D_0.x.array[entities] = vol.material.D_0[f"{spe}"] + E_D.x.array[entities] = vol.material.E_D[f"{spe}"] + + # create global D function + D = fem.Function(self.V_DG_1) + expr = D_0 * ufl.exp( + -E_D / F.as_fenics_constant(F.k_B, self.mesh.mesh) / self.temperature + ) + D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) + D.interpolate(D_expr) + + # add the global D to the export + export.D = D + + # add to global list for updating later + self.D_global_expr.append(D_expr) + self.D_global.append(D) + + def define_function_spaces(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.""" @@ -173,6 +209,8 @@ def define_function_space(self): element = ufl.MixedElement(elements) self.function_space = fem.FunctionSpace(self.mesh.mesh, element) + self.V_DG_0 = fem.FunctionSpace(self.mesh.mesh, ("DG", 0)) + self.V_DG_1 = fem.FunctionSpace(self.mesh.mesh, ("DG", 1)) self.u = Function(self.function_space) self.u_n = Function(self.function_space) @@ -389,10 +427,15 @@ def run(self): for bc in self.boundary_conditions: bc.update(float(self.t)) + # update global D if temperature time dependent or internal + # variables time dependent + for D, expr in zip(self.D_global, self.D_global_expr): + D.interpolate(expr) + self.solver.solve(self.u) # post processing - self.post_processing(times, flux_values, flux_values_1, flux_values_2) + self.post_processing() # update previous solution self.u_n.x.array[:] = self.u.x.array[:] @@ -402,58 +445,23 @@ def run(self): return times, flux_values - def post_processing(self, times, flux_values, flux_values_1, flux_values_2): + def post_processing(self): if not self.multispecies: - D_D = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] - ) - cm = self.u self.species[0].post_processing_solution = self.u - - surface_flux = form(D_D * dot(grad(cm), self.mesh.n) * self.ds(2)) - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(float(self.t)) else: for idx, spe in enumerate(self.species): spe.post_processing_solution = self.u.sub(idx) - cm_1, cm_2 = self.u.split() - D_1 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[0] - ) - D_2 = self.subdomains[0].material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, self.species[1] - ) - surface_flux_1 = form(D_1 * dot(grad(cm_1), self.mesh.n) * self.ds(2)) - surface_flux_2 = form(D_2 * dot(grad(cm_2), self.mesh.n) * self.ds(2)) - flux_1 = assemble_scalar(surface_flux_1) - flux_2 = assemble_scalar(surface_flux_2) - flux_values_1.append(flux_1) - flux_values_2.append(flux_2) - times.append(float(self.t)) - - # for export in self.exports: - # if isinstance(export, F.SurfaceQuantity): - # export.compute(self.mesh, self.ds) - # if isinstance(export, F.VolumeQuantity): - # export.compute(self.mesh, self.dx) - + self.derived_quantities["t(s)"].append(float(self.t)) for export in self.exports: + # TODO if export type derived quantity if isinstance(export, F.SurfaceFlux): - # evaluate value of export - D_export = export.volume_subdomain.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature, export.field - ) - export_value = export.compute_quantity( - D_export, + export.compute( self.mesh, self.ds, ) - # update derived quantities dict - self.derived_quantities["t(s)"].append(float(self.t)) - self.derived_quantities[f"{export.title}"].append(export_value) + self.derived_quantities[f"{export.title}"].append(export.value) export.write(t=float(self.t)) From 67dabb9140126044e26413088b0d408fda55f2c9 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 11:39:55 -0400 Subject: [PATCH 12/62] use derived quantites dict --- test/test_multispecies_problem.py | 22 +++++++++++-- test/test_permeation_problem.py | 22 +++++++++++-- test/test_surface_flux.py | 53 +++++++++++++++++++++++++++++++ 3 files changed, 92 insertions(+), 5 deletions(-) create mode 100644 test/test_surface_flux.py diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index ebae0dcf5..5db5892d2 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -85,8 +85,20 @@ def test_multispecies_permeation_problem(): max_iterations=30, final_time=10, ) - my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) + + my_model.exports = [ + F.SurfaceFlux( + filename="flux_spe_1.csv", + field=spe_1, + surface_subdomain=right_surface, + ), + F.SurfaceFlux( + filename="flux_spe_2.csv", + field=spe_2, + surface_subdomain=right_surface, + ), + ] my_model.initialise() my_model.solver.convergence_criterion = "incremental" @@ -98,14 +110,18 @@ def test_multispecies_permeation_problem(): opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" ksp.setFromOptions() - times, flux_values = my_model.run() + my_model.run() + + times = my_model.derived_quantities["t(s)"] + flux_values_spe_1 = my_model.derived_quantities["Flux surface 2: spe_1"] + flux_values_spe_2 = my_model.derived_quantities["Flux surface 2: spe_2"] # ---------------------- analytical solutions ----------------------------- # common values times = np.array(times) P_up = float(my_model.boundary_conditions[-1].pressure) - flux_values_spe_1, flux_values_spe_2 = flux_values[0], flux_values[1] + # flux_values_spe_1, flux_values_spe_2 = flux_values[0], flux_values[1] # ##### compute analyical solution for species 1 ##### # D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 8e2eb014a..6b074c775 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -4,6 +4,7 @@ import numpy as np import festim as F import os +import ufl def relative_error_computed_to_analytical( @@ -60,7 +61,14 @@ 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("mobile_concentration.xdmf", field=mobile_H), + F.SurfaceFlux( + filename="my_surface_flux.csv", + field=mobile_H, + surface_subdomain=right_surface, + ), + ] my_model.settings = F.Settings( atol=1e10, @@ -82,7 +90,12 @@ def test_permeation_problem(mesh_size=1001): opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" ksp.setFromOptions() - times, flux_values = my_model.run() + my_model.run() + + # print(my_model.derived_quantities.keys()) + + times = my_model.derived_quantities["t(s)"] + flux_values = my_model.derived_quantities["Flux surface 2: H"] # -------------------------- analytical solution ------------------------------------- @@ -187,3 +200,8 @@ def test_permeation_problem_multi_volume(tmp_path): ) assert error < 0.01 + + +if __name__ == "__main__": + test_permeation_problem() + # test_permeation_problem_multi_volume(tmp_path=".") diff --git a/test/test_surface_flux.py b/test/test_surface_flux.py new file mode 100644 index 000000000..a6dc7ec42 --- /dev/null +++ b/test/test_surface_flux.py @@ -0,0 +1,53 @@ +import festim as F +import numpy as np +from petsc4py import PETSc + + +def surface_flux_export(): + """Test that the field attribute can be a string and is found in the species list""" + L = 3e-04 + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=L) + mobile_H = F.Species("H") + + my_model = F.HydrogenTransportProblem() + + my_model.mesh = F.Mesh1D(vertices=np.linspace(0, L, num=1001)) + my_model.subdomains = [my_subdomain, left_surface, right_surface] + my_model.species = [mobile_H] + my_model.temperature = 500 + 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.settings = F.Settings(atol=1e10, rtol=1e-10, final_time=50) + my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) + + my_export = F.SurfaceFlux( + filename="my_surface_flux.csv", + field=mobile_H, + surface_subdomain=right_surface, + volume_subdomain=my_subdomain, + ) + my_model.exports = [my_export, F.XDMFExport("mobile_H.xdmf", field=mobile_H)] + + 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() + + my_model.run() + + +if __name__ == "__main__": + surface_flux_export() From 1db757f18b2c98d135b2d3f523e5519f513fb36f Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 11:49:46 -0400 Subject: [PATCH 13/62] work in multispecies --- 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 b0d9a887f..581c3655d 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -119,7 +119,7 @@ def multispecies(self): return len(self.species) > 1 def initialise(self): - self.define_function_spaces() + self.define_function_space() self.define_markers_and_measures() self.assign_functions_to_species() @@ -134,6 +134,7 @@ def initialise(self): def initialise_exports(self): """Defines the export writers of the model, if field is given as a string, find species object in self.species""" + self.derived_quantities["t(s)"] = [] for export in self.exports: # if name of species is given then replace with species object if isinstance(export.field, list): @@ -152,7 +153,6 @@ def initialise_exports(self): if isinstance(export, F.SurfaceFlux): self.define_D_global(export) - self.derived_quantities["t(s)"] = [] self.derived_quantities[f"{export.title}"] = [] def define_D_global(self, export): @@ -182,13 +182,14 @@ def define_D_global(self, export): D.interpolate(D_expr) # add the global D to the export - export.D = D + if export.field == spe: + export.D = D # add to global list for updating later self.D_global_expr.append(D_expr) self.D_global.append(D) - def define_function_spaces(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.""" From 202225590a18016706fe3ae14d6ae38b7184a51e Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 11:50:00 -0400 Subject: [PATCH 14/62] updated permeation test --- test/test_permeation_problem.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 6b074c775..0b2464b7f 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -159,7 +159,14 @@ def test_permeation_problem_multi_volume(tmp_path): ), ] my_model.exports = [ - F.VTXExport(os.path.join(tmp_path, "mobile_concentration_h.bp"), field=mobile_H) + F.VTXExport( + os.path.join(tmp_path, "mobile_concentration_h.bp"), field=mobile_H + ), + F.SurfaceFlux( + filename=os.path.join(tmp_path, "my_surface_flux.csv"), + field=mobile_H, + surface_subdomain=right_surface, + ), ] my_model.settings = F.Settings( @@ -182,7 +189,10 @@ def test_permeation_problem_multi_volume(tmp_path): opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" ksp.setFromOptions() - times, flux_values = my_model.run() + my_model.run() + + times = my_model.derived_quantities["t(s)"] + flux_values = my_model.derived_quantities["Flux surface 2: H"] # ---------------------- analytical solution ----------------------------- D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) From 42a93b9fb1d27604f70093f6e1d2ba831c539674 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 12:42:47 -0400 Subject: [PATCH 15/62] blank lines --- festim/hydrogen_transport_problem.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 581c3655d..4eb25e610 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -135,6 +135,7 @@ def initialise_exports(self): """Defines the export writers of the model, if field is given as a string, find species object in self.species""" self.derived_quantities["t(s)"] = [] + for export in self.exports: # if name of species is given then replace with species object if isinstance(export.field, list): @@ -151,6 +152,7 @@ def initialise_exports(self): if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) + # initialise surface flux derived quantity if isinstance(export, F.SurfaceFlux): self.define_D_global(export) self.derived_quantities[f"{export.title}"] = [] @@ -463,8 +465,8 @@ def post_processing(self): ) # update derived quantities dict self.derived_quantities[f"{export.title}"].append(export.value) - - export.write(t=float(self.t)) + if export.write_to_file: + export.write(t=float(self.t)) if isinstance(export, (F.VTXExport, F.XDMFExport)): export.write(float(self.t)) From 42e2e9085ea3737c5bc80294e47ca75339b44dbd Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:28:02 -0400 Subject: [PATCH 16/62] split methods --- festim/exports/surface_flux.py | 2 +- festim/hydrogen_transport_problem.py | 63 +++++++++++++--------------- 2 files changed, 31 insertions(+), 34 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index bbcad76ce..cfb28aece 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -4,7 +4,7 @@ import csv -class SurfaceFlux(F.SurfaceQuantity): +class SurfaceFlux: """Exports surface flux at a given subdomain Args: diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 581c3655d..fa5419369 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -152,42 +152,39 @@ def initialise_exports(self): export.writer.write_mesh(self.mesh.mesh) if isinstance(export, F.SurfaceFlux): - self.define_D_global(export) + D, D_expr = self.define_D_global(export.field) self.derived_quantities[f"{export.title}"] = [] - - def define_D_global(self, export): - for spe in self.species: - D_0 = fem.Function(self.V_DG_0) - E_D = fem.Function(self.V_DG_0) - for vol in self.volume_subdomains: - entities = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) - if not self.multispecies: - if isinstance(vol.material.D_0, (float, int, fem.Constant)): - D_0.x.array[entities] = vol.material.D_0 - E_D.x.array[entities] = vol.material.E_D - else: - raise NotImplementedError( - "diffusion values as functions not supported" - ) - else: - D_0.x.array[entities] = vol.material.D_0[f"{spe}"] - E_D.x.array[entities] = vol.material.E_D[f"{spe}"] - - # create global D function - D = fem.Function(self.V_DG_1) - expr = D_0 * ufl.exp( - -E_D / F.as_fenics_constant(F.k_B, self.mesh.mesh) / self.temperature - ) - D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) - D.interpolate(D_expr) - - # add the global D to the export - if export.field == spe: + # add the global D to the export export.D = D + self.D_global_expr.append(D_expr) + self.D_global.append(D) + + def define_D_global(self, spe): + assert isinstance(spe, F.Species) + D_0 = fem.Function(self.V_DG_0) + E_D = fem.Function(self.V_DG_0) + for vol in self.volume_subdomains: + entities = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) + if not self.multispecies: + if isinstance(vol.material.D_0, (float, int, fem.Constant)): + D_0.x.array[entities] = vol.material.D_0 + E_D.x.array[entities] = vol.material.E_D + else: + raise NotImplementedError( + "diffusion values as functions not supported" + ) + else: + D_0.x.array[entities] = vol.material.D_0[f"{spe}"] + E_D.x.array[entities] = vol.material.E_D[f"{spe}"] - # add to global list for updating later - self.D_global_expr.append(D_expr) - self.D_global.append(D) + # create global D function + D = fem.Function(self.V_DG_1) + expr = D_0 * ufl.exp( + -E_D / F.as_fenics_constant(F.k_B, self.mesh.mesh) / self.temperature + ) + D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) + D.interpolate(D_expr) + return D, D_expr def define_function_space(self): """Creates the function space of the model, creates a mixed element if From 176274f33e471eb0cd27bb1f5a7dda0a0c96bf54 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:36:51 -0400 Subject: [PATCH 17/62] commented out import --- festim/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/__init__.py b/festim/__init__.py index bbe6f48da..9f35ae4a6 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -37,7 +37,7 @@ from .stepsize import Stepsize -from .exports.surface_quantity import SurfaceQuantity +# from .exports.surface_quantity import SurfaceQuantity from .exports.surface_flux import SurfaceFlux from .exports.vtx import VTXExport from .exports.xdmf import XDMFExport From a5665c962d885a2375deaa97e68f6da732031964 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:37:05 -0400 Subject: [PATCH 18/62] no need for D_global attribute --- festim/hydrogen_transport_problem.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index fa5419369..ff9bfedfc 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -100,8 +100,6 @@ def __init__( self.volume_subdomains = [] self.bc_forms = [] self.derived_quantities = {} - self.D_global = [] - self.D_global_expr = [] @property def temperature(self): @@ -151,13 +149,18 @@ def initialise_exports(self): if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) + spe_to_D_global = {} + spe_to_D_global_expr = {} if isinstance(export, F.SurfaceFlux): - D, D_expr = self.define_D_global(export.field) + if export.field in spe_to_D_global: + D = spe_to_D_global[export.field] + D_expr = spe_to_D_global_expr[export.field] + else: + D, D_expr = self.define_D_global(export.field) self.derived_quantities[f"{export.title}"] = [] # add the global D to the export export.D = D - self.D_global_expr.append(D_expr) - self.D_global.append(D) + export.D_expr = D_expr def define_D_global(self, spe): assert isinstance(spe, F.Species) @@ -427,9 +430,13 @@ def run(self): # update global D if temperature time dependent or internal # variables time dependent - for D, expr in zip(self.D_global, self.D_global_expr): - D.interpolate(expr) + species_not_updated = self.species.copy() + for export in self.exports: + if isinstance(export, F.SurfaceFlux): + if export.field in species_not_updated: + export.D.interpolate(export.D_expr) + # solve main problem self.solver.solve(self.u) # post processing From f6e87c65c5ff0d857ff2fc06edb37fcfec4e1a21 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:39:31 -0400 Subject: [PATCH 19/62] dict initialisation outside of loop --- festim/hydrogen_transport_problem.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index ff9bfedfc..149081033 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -133,6 +133,7 @@ def initialise_exports(self): """Defines the export writers of the model, if field is given as a string, find species object in self.species""" self.derived_quantities["t(s)"] = [] + for export in self.exports: # if name of species is given then replace with species object if isinstance(export.field, list): @@ -149,8 +150,12 @@ def initialise_exports(self): if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) - spe_to_D_global = {} - spe_to_D_global_expr = {} + # compute diffusivity function for surface fluxes + # correspondance dicts + spe_to_D_global = {} + spe_to_D_global_expr = {} + + for export in self.exports: if isinstance(export, F.SurfaceFlux): if export.field in spe_to_D_global: D = spe_to_D_global[export.field] From 8482ea1c317e2de3b6147969ced6032172b7f00a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:41:58 -0400 Subject: [PATCH 20/62] add D and D_expr to dicts + comments --- festim/hydrogen_transport_problem.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 149081033..469200a30 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -158,11 +158,17 @@ def initialise_exports(self): for export in self.exports: if isinstance(export, F.SurfaceFlux): if export.field in spe_to_D_global: + # if already computed then use the same D D = spe_to_D_global[export.field] D_expr = spe_to_D_global_expr[export.field] else: + # if not computed then compute D and add it to the dict D, D_expr = self.define_D_global(export.field) + spe_to_D_global[export.field] = D + spe_to_D_global_expr[export.field] = D_expr + self.derived_quantities[f"{export.title}"] = [] + # add the global D to the export export.D = D export.D_expr = D_expr From 42e00ab6b82c98071c5b0be87e23c50b8fcc7bc9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:47:41 -0400 Subject: [PATCH 21/62] cell_indices instead of entities --- festim/hydrogen_transport_problem.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 469200a30..3b21adeb7 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -151,9 +151,9 @@ def initialise_exports(self): export.writer.write_mesh(self.mesh.mesh) # compute diffusivity function for surface fluxes - # correspondance dicts - spe_to_D_global = {} - spe_to_D_global_expr = {} + + spe_to_D_global = {} # links species to global D function + spe_to_D_global_expr = {} # links species to D expression for export in self.exports: if isinstance(export, F.SurfaceFlux): @@ -175,21 +175,23 @@ def initialise_exports(self): def define_D_global(self, spe): assert isinstance(spe, F.Species) + D_0 = fem.Function(self.V_DG_0) E_D = fem.Function(self.V_DG_0) for vol in self.volume_subdomains: - entities = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) + cell_indices = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) if not self.multispecies: if isinstance(vol.material.D_0, (float, int, fem.Constant)): - D_0.x.array[entities] = vol.material.D_0 - E_D.x.array[entities] = vol.material.E_D + # replace values of D_0 and E_D by values from the material + D_0.x.array[cell_indices] = vol.material.D_0 + E_D.x.array[cell_indices] = vol.material.E_D else: raise NotImplementedError( "diffusion values as functions not supported" ) else: - D_0.x.array[entities] = vol.material.D_0[f"{spe}"] - E_D.x.array[entities] = vol.material.E_D[f"{spe}"] + D_0.x.array[cell_indices] = vol.material.D_0[f"{spe}"] + E_D.x.array[cell_indices] = vol.material.E_D[f"{spe}"] # create global D function D = fem.Function(self.V_DG_1) From 43a587a7d45cc8546cb68257ae1bef01279c044d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 13:50:56 -0400 Subject: [PATCH 22/62] refactored define_D_global --- festim/hydrogen_transport_problem.py | 16 +++------ festim/material.py | 52 ++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 12 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3b21adeb7..62f2b5aa3 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -180,18 +180,10 @@ def define_D_global(self, spe): E_D = fem.Function(self.V_DG_0) for vol in self.volume_subdomains: cell_indices = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) - if not self.multispecies: - if isinstance(vol.material.D_0, (float, int, fem.Constant)): - # replace values of D_0 and E_D by values from the material - D_0.x.array[cell_indices] = vol.material.D_0 - E_D.x.array[cell_indices] = vol.material.E_D - else: - raise NotImplementedError( - "diffusion values as functions not supported" - ) - else: - D_0.x.array[cell_indices] = vol.material.D_0[f"{spe}"] - E_D.x.array[cell_indices] = vol.material.E_D[f"{spe}"] + + # replace values of D_0 and E_D by values from the material + D_0.x.array[cell_indices] = vol.material.get_D_0(species=spe) + E_D.x.array[cell_indices] = vol.material.get_E_D(species=spe) # create global D function D = fem.Function(self.V_DG_1) diff --git a/festim/material.py b/festim/material.py index ceeced3ff..c46ff74d1 100644 --- a/festim/material.py +++ b/festim/material.py @@ -35,6 +35,58 @@ def __init__(self, D_0, E_D, name=None) -> None: self.E_D = E_D self.name = name + def get_D_0(self, species=None): + """Returns the pre-exponential factor of the diffusion coefficient + Args: + species (festim.Species or str, optional): the species we want the pre-exponential + factor of the diffusion coefficient of. Only needed if D_0 is a dict. + Returns: + float: the pre-exponential factor of the diffusion coefficient + """ + + if isinstance(self.D_0, (float, int)): + return self.D_0 + + elif isinstance(self.D_0, dict): + if species is None: + raise ValueError("species must be provided if D_0 is a dict") + + if species in self.D_0: + return self.D_0[species] + elif species.name in self.D_0: + return self.D_0[species.name] + else: + raise ValueError(f"{species} is not in D_0 keys") + + else: + raise ValueError("D_0 must be either a float or a dict") + + def get_E_D(self, species=None): + """Returns the activation energy of the diffusion coefficient + Args: + species (festim.Species or str, optional): the species we want the activation + energy of the diffusion coefficient of. Only needed if E_D is a dict. + Returns: + float: the activation energy of the diffusion coefficient + """ + + if isinstance(self.E_D, (float, int)): + return self.E_D + + elif isinstance(self.E_D, dict): + if species is None: + raise ValueError("species must be provided if E_D is a dict") + + if species in self.E_D: + return self.E_D[species] + elif species.name in self.E_D: + return self.E_D[species.name] + else: + raise ValueError(f"{species} is not in E_D keys") + + else: + raise ValueError("E_D must be either a float or a dict") + def get_diffusion_coefficient(self, mesh, temperature, species=None): """Defines the diffusion coefficient Args: From 2709196a8316c8ad9af6b58f8fd1433e34132f59 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 14:04:01 -0400 Subject: [PATCH 23/62] added TODO --- festim/material.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/festim/material.py b/festim/material.py index c46ff74d1..057760434 100644 --- a/festim/material.py +++ b/festim/material.py @@ -98,6 +98,14 @@ def get_diffusion_coefficient(self, mesh, temperature, species=None): Returns: ufl.algebra.Product: the diffusion coefficient """ + # TODO use get_D_0 and get_E_D to refactore this method, something like: + # D_0 = self.get_D_0(species=species) + # E_D = self.get_E_D(species=species) + + # D_0 = F.as_fenics_constant(D_0, mesh) + # E_D = F.as_fenics_constant(E_D, mesh) + + # return D_0 * ufl.exp(-E_D / F.k_B / temperature) if isinstance(self.D_0, (float, int)) and isinstance(self.E_D, (float, int)): D_0 = F.as_fenics_constant(self.D_0, mesh) From 876d5e7babf70771f65724f9f0b149bd7cb0edaf Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 14:04:11 -0400 Subject: [PATCH 24/62] docstrings + species instead of spe --- festim/hydrogen_transport_problem.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 62f2b5aa3..339593027 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -173,8 +173,17 @@ def initialise_exports(self): export.D = D export.D_expr = D_expr - def define_D_global(self, spe): - assert isinstance(spe, F.Species) + def define_D_global(self, species): + """Defines the global diffusion coefficient for a given species + + Args: + species (F.Species): the species + + Returns: + dolfinx.fem.Function, dolfinx.fem.Expression: the global diffusion + coefficient and the expression of the global diffusion coefficient + """ + assert isinstance(species, F.Species) D_0 = fem.Function(self.V_DG_0) E_D = fem.Function(self.V_DG_0) @@ -182,8 +191,8 @@ def define_D_global(self, spe): cell_indices = vol.locate_subdomain_entities(self.mesh.mesh, self.mesh.vdim) # replace values of D_0 and E_D by values from the material - D_0.x.array[cell_indices] = vol.material.get_D_0(species=spe) - E_D.x.array[cell_indices] = vol.material.get_E_D(species=spe) + D_0.x.array[cell_indices] = vol.material.get_D_0(species=species) + E_D.x.array[cell_indices] = vol.material.get_E_D(species=species) # create global D function D = fem.Function(self.V_DG_1) From d7e9d0d502657230b6090f590fd99815df63cdb3 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 1 Nov 2023 14:09:50 -0400 Subject: [PATCH 25/62] pop species from species_not_updated --- festim/hydrogen_transport_problem.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 339593027..ab8b2d576 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -444,11 +444,13 @@ def run(self): # update global D if temperature time dependent or internal # variables time dependent - species_not_updated = self.species.copy() + species_not_updated = self.species.copy() # make a copy of the species for export in self.exports: if isinstance(export, F.SurfaceFlux): + # if the D of the species has not been updated yet if export.field in species_not_updated: export.D.interpolate(export.D_expr) + species_not_updated.remove(export.field) # solve main problem self.solver.solve(self.u) From 45ed02d0341bd6ea8aa88c2a564270530225089d Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:17:02 -0400 Subject: [PATCH 26/62] store data and t, move setters to parent --- festim/exports/surface_flux.py | 52 ++-------------------------------- 1 file changed, 3 insertions(+), 49 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index 5ddbe9c1a..efb8b204c 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -20,11 +20,10 @@ def __init__( surface_subdomain: int or F.SurfaceSubdomain1D, filename: str = None, ) -> None: - self.field = field - self.surface_subdomain = surface_subdomain - self.filename = filename + super().__init__(field, surface_subdomain, filename) - self.D = None + self.t = [] + self.data = [] self.title = "Flux surface {}: {}".format( self.surface_subdomain.id, self.field.name @@ -35,51 +34,6 @@ def __init__( writer = csv.writer(file) writer.writerow(["Time", f"{self.title}"]) - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is None: - self._filename = None - if not isinstance(value, str): - raise TypeError("filename must be of type str") - if not value.endswith(".csv"): - raise ValueError("filename must end with .csv") - self._filename = value - - @property - def surface_subdomain(self): - return self._surface_subdomain - - @surface_subdomain.setter - def surface_subdomain(self, value): - if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance( - value, bool - ): - raise TypeError("surface should be an int or F.SurfaceSubdomain1D") - self._surface_subdomain = value - - @property - def field(self): - return self._field - - @field.setter - def field(self, value): - # check that field is festim.Species - if not isinstance(value, (F.Species, str)): - raise TypeError("field must be of type festim.Species") - - self._field = value - - @property - def write_to_file(self): - if self.filename is None: - return False - else: - return True - def compute(self, mesh, ds): self.value = fem.assemble_scalar( fem.form( From 068efdc9531678d5872eadbb14c54d74320f8df8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:17:10 -0400 Subject: [PATCH 27/62] add setters --- festim/exports/surface_quantity.py | 52 +++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index 066b1afb7..a1cb38eec 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -9,14 +9,56 @@ class SurfaceQuantity: Attributes: """ - def __init__( - self, field, surface_subdomain, volume_subdomain, filename: str = None - ) -> None: + def __init__(self, field, surface_subdomain, filename: str = None) -> None: self.field = field - self.subdomain = surface_subdomain - self.volume_subdomain = volume_subdomain + self.surface_subdomain = surface_subdomain self.filename = filename self.ds = None self.n = None self.D = None self.S = None + + @property + def filename(self): + return self._filename + + @filename.setter + def filename(self, value): + if value is None: + self._filename = None + elif not isinstance(value, str): + raise TypeError("filename must be of type str") + elif not value.endswith(".csv"): + raise ValueError("filename must end with .csv") + self._filename = value + + @property + def surface_subdomain(self): + return self._surface_subdomain + + @surface_subdomain.setter + def surface_subdomain(self, value): + if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance( + value, bool + ): + raise TypeError("surface should be an int or F.SurfaceSubdomain1D") + self._surface_subdomain = value + + @property + def field(self): + return self._field + + @field.setter + def field(self, value): + # check that field is festim.Species + if not isinstance(value, (F.Species, str)): + raise TypeError("field must be of type festim.Species") + + self._field = value + + @property + def write_to_file(self): + if self.filename is None: + return False + else: + return True From a341ecdf8f0b77a0cc7b1bd4f45ad101fe91aea8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:17:22 -0400 Subject: [PATCH 28/62] include surface quantity --- festim/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/__init__.py b/festim/__init__.py index 9f35ae4a6..bbe6f48da 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -37,7 +37,7 @@ from .stepsize import Stepsize -# from .exports.surface_quantity import SurfaceQuantity +from .exports.surface_quantity import SurfaceQuantity from .exports.surface_flux import SurfaceFlux from .exports.vtx import VTXExport from .exports.xdmf import XDMFExport From 3461803f6fc33c2f103e6e5a7d6f13f7e781fe57 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:17:39 -0400 Subject: [PATCH 29/62] refactoring --- festim/hydrogen_transport_problem.py | 40 ++++++++++++++-------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 66b6ba07e..795decce8 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -1,16 +1,12 @@ from dolfinx import fem +from dolfinx.mesh import meshtags from dolfinx.nls.petsc import NewtonSolver -from dolfinx.io import XDMFFile import basix import ufl from mpi4py import MPI -from dolfinx.fem import Function, form, assemble_scalar from dolfinx.mesh import meshtags -from ufl import TestFunction, dot, grad, Measure, FacetNormal import numpy as np import tqdm.autonotebook - - import festim as F @@ -99,7 +95,6 @@ def __init__( self.formulation = None self.volume_subdomains = [] self.bc_forms = [] - self.derived_quantities = {} @property def temperature(self): @@ -117,7 +112,7 @@ def multispecies(self): return len(self.species) > 1 def initialise(self): - self.define_function_space() + self.define_function_spaces() self.define_markers_and_measures() self.assign_functions_to_species() @@ -132,7 +127,6 @@ def initialise(self): def initialise_exports(self): """Defines the export writers of the model, if field is given as a string, find species object in self.species""" - self.derived_quantities["t(s)"] = [] for export in self.exports: # if name of species is given then replace with species object @@ -167,8 +161,6 @@ def initialise_exports(self): spe_to_D_global[export.field] = D spe_to_D_global_expr[export.field] = D_expr - self.derived_quantities[f"{export.title}"] = [] - # add the global D to the export export.D = D export.D_expr = D_expr @@ -182,6 +174,7 @@ def define_D_global(self, species): Returns: dolfinx.fem.Function, dolfinx.fem.Expression: the global diffusion coefficient and the expression of the global diffusion coefficient + for a given species """ assert isinstance(species, F.Species) @@ -203,10 +196,12 @@ def define_D_global(self, species): D.interpolate(D_expr) return D, D_expr - def define_function_space(self): + def define_function_spaces(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.""" + function u and u_n. Create global DG function spaces of degree 0 and 1 + for the global diffusion coefficient""" + element_CG = basix.ufl.element( basix.ElementFamily.P, self.mesh.mesh.basix_cell(), @@ -224,11 +219,13 @@ def define_function_space(self): element = ufl.MixedElement(elements) self.function_space = fem.FunctionSpace(self.mesh.mesh, element) + + # create global DG function spaces of degree 0 and 1 self.V_DG_0 = fem.FunctionSpace(self.mesh.mesh, ("DG", 0)) self.V_DG_1 = fem.FunctionSpace(self.mesh.mesh, ("DG", 1)) - self.u = Function(self.function_space) - self.u_n = Function(self.function_space) + self.u = fem.Function(self.function_space) + self.u_n = fem.Function(self.function_space) def assign_functions_to_species(self): """Creates the solution, prev solution, test function and @@ -300,10 +297,10 @@ def define_markers_and_measures(self): ) # define measures - self.ds = Measure( + self.ds = ufl.Measure( "ds", domain=self.mesh.mesh, subdomain_data=self.facet_meshtags ) - self.dx = Measure( + self.dx = ufl.Measure( "dx", domain=self.mesh.mesh, subdomain_data=self.volume_meshtags ) @@ -391,7 +388,9 @@ def create_formulation(self): self.mesh.mesh, self.temperature, spe ) - self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) + self.formulation += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * self.dx( + vol.id + ) self.formulation += ((u - u_n) / self.dt) * v * self.dx(vol.id) # add sources @@ -473,7 +472,6 @@ def post_processing(self): for idx, spe in enumerate(self.species): spe.post_processing_solution = self.u.sub(idx) - self.derived_quantities["t(s)"].append(float(self.t)) for export in self.exports: # TODO if export type derived quantity if isinstance(export, F.SurfaceFlux): @@ -481,8 +479,10 @@ def post_processing(self): self.mesh, self.ds, ) - # update derived quantities dict - self.derived_quantities[f"{export.title}"].append(export.value) + # update export data + export.t.append(float(self.t)) + + # if filename given write export data to file if export.write_to_file: export.write(t=float(self.t)) From a52d76873ed052787a8a31ce929b2af0cfceebf5 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 1 Nov 2023 16:19:14 -0400 Subject: [PATCH 30/62] updated tests --- test/test_dirichlet_bc.py | 8 +++---- test/test_multispecies_problem.py | 30 +++++++++++++------------- test/test_permeation_problem.py | 35 +++++++++++++------------------ test/test_species.py | 2 +- test/test_subdomains.py | 2 +- 5 files changed, 34 insertions(+), 43 deletions(-) diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index f7f5b0e05..486950b7d 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -65,7 +65,7 @@ def test_callable_for_value(): mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species] ) - my_model.define_function_space() + my_model.define_function_spaces() my_model.define_markers_and_measures() T = fem.Constant(my_model.mesh.mesh, 550.0) @@ -103,7 +103,7 @@ def test_value_callable_x_t_T(): mesh=F.Mesh(mesh), subdomains=[subdomain, vol_subdomain], species=[species] ) - my_model.define_function_space() + my_model.define_function_spaces() my_model.define_markers_and_measures() T = fem.Constant(my_model.mesh.mesh, 550.0) @@ -146,7 +146,7 @@ def test_callable_t_only(value): species=[species], ) - my_model.define_function_space() + my_model.define_function_spaces() my_model.define_markers_and_measures() T = fem.Constant(my_model.mesh.mesh, 550.0) @@ -189,7 +189,7 @@ def test_callable_x_only(): species=[species], ) - my_model.define_function_space() + my_model.define_function_spaces() my_model.define_markers_and_measures() T = fem.Constant(my_model.mesh.mesh, 550.0) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 5db5892d2..cbe2ce110 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -87,18 +87,17 @@ def test_multispecies_permeation_problem(): ) my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) - my_model.exports = [ - F.SurfaceFlux( - filename="flux_spe_1.csv", - field=spe_1, - surface_subdomain=right_surface, - ), - F.SurfaceFlux( - filename="flux_spe_2.csv", - field=spe_2, - surface_subdomain=right_surface, - ), - ] + outgassing_flux_spe_1 = F.SurfaceFlux( + filename="flux_spe_1.csv", + field=spe_1, + surface_subdomain=right_surface, + ) + outgassing_flux_spe_2 = F.SurfaceFlux( + filename="flux_spe_2.csv", + field=spe_2, + surface_subdomain=right_surface, + ) + my_model.exports = [outgassing_flux_spe_1, outgassing_flux_spe_2] my_model.initialise() my_model.solver.convergence_criterion = "incremental" @@ -112,16 +111,15 @@ def test_multispecies_permeation_problem(): my_model.run() - times = my_model.derived_quantities["t(s)"] - flux_values_spe_1 = my_model.derived_quantities["Flux surface 2: spe_1"] - flux_values_spe_2 = my_model.derived_quantities["Flux surface 2: spe_2"] + times = outgassing_flux_spe_1.t + flux_values_spe_1 = outgassing_flux_spe_1.data + flux_values_spe_2 = outgassing_flux_spe_2.data # ---------------------- analytical solutions ----------------------------- # common values times = np.array(times) P_up = float(my_model.boundary_conditions[-1].pressure) - # flux_values_spe_1, flux_values_spe_2 = flux_values[0], flux_values[1] # ##### compute analyical solution for species 1 ##### # D_spe_1 = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature, spe_1) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 0b2464b7f..cdf8382ef 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -61,13 +61,13 @@ def test_permeation_problem(mesh_size=1001): subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H" ), ] + outgassing_flux = F.SurfaceFlux( + field=mobile_H, + surface_subdomain=right_surface, + ) my_model.exports = [ F.XDMFExport("mobile_concentration.xdmf", field=mobile_H), - F.SurfaceFlux( - filename="my_surface_flux.csv", - field=mobile_H, - surface_subdomain=right_surface, - ), + outgassing_flux, ] my_model.settings = F.Settings( @@ -92,10 +92,8 @@ def test_permeation_problem(mesh_size=1001): my_model.run() - # print(my_model.derived_quantities.keys()) - - times = my_model.derived_quantities["t(s)"] - flux_values = my_model.derived_quantities["Flux surface 2: H"] + times = outgassing_flux.t + flux_values = outgassing_flux.data # -------------------------- analytical solution ------------------------------------- @@ -158,15 +156,15 @@ def test_permeation_problem_multi_volume(tmp_path): subdomain=left_surface, S_0=4.02e21, E_S=1.04, pressure=100, species="H" ), ] + outgassing_flux = F.SurfaceFlux( + field=mobile_H, + surface_subdomain=right_surface, + ) my_model.exports = [ F.VTXExport( os.path.join(tmp_path, "mobile_concentration_h.bp"), field=mobile_H ), - F.SurfaceFlux( - filename=os.path.join(tmp_path, "my_surface_flux.csv"), - field=mobile_H, - surface_subdomain=right_surface, - ), + outgassing_flux, ] my_model.settings = F.Settings( @@ -191,8 +189,8 @@ def test_permeation_problem_multi_volume(tmp_path): my_model.run() - times = my_model.derived_quantities["t(s)"] - flux_values = my_model.derived_quantities["Flux surface 2: H"] + times = outgassing_flux.t + flux_values = outgassing_flux.data # ---------------------- analytical solution ----------------------------- D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) @@ -210,8 +208,3 @@ def test_permeation_problem_multi_volume(tmp_path): ) assert error < 0.01 - - -if __name__ == "__main__": - test_permeation_problem() - # test_permeation_problem_multi_volume(tmp_path=".") diff --git a/test/test_species.py b/test/test_species.py index 2eafd9844..d081d1ad0 100644 --- a/test/test_species.py +++ b/test/test_species.py @@ -23,7 +23,7 @@ def test_assign_functions_to_species(): # F.Species(name="Trap"), ], ) - model.define_function_space() + model.define_function_spaces() model.assign_functions_to_species() for spe in model.species: diff --git a/test/test_subdomains.py b/test/test_subdomains.py index e8d5f004b..5e8489944 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -20,7 +20,7 @@ def test_different_surface_ids(): volume_subdomain, ] - my_test_model.define_function_space() + my_test_model.define_function_spaces() my_test_model.define_markers_and_measures() for surf_id in surface_subdomains_ids: From a550c54eccb2a4c9d11d841856aed554a6dc000d Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 13:58:16 -0400 Subject: [PATCH 31/62] dont use skip_post_processing --- test/test_h_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index c5a42d3f0..a8dfc570f 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -130,7 +130,7 @@ def test_iterate(): for i in range(10): # RUN - my_model.iterate(skip_post_processing=True) + my_model.iterate() # TEST From 2356e1a17a007d188b788f353b22daf183c6f5c4 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:08:16 -0400 Subject: [PATCH 32/62] added docs, only give n --- festim/exports/surface_flux.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index efb8b204c..e189a14cb 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -8,10 +8,18 @@ class SurfaceFlux(F.SurfaceQuantity): """Exports surface flux at a given subdomain Args: + field (festim.Species): species for which the surface flux is computed + surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + filename (str, optional): name of the file to which the surface flux is exported Attributes: - - Usage: + field (festim.Species): species for which the surface flux is computed + surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + filename (str): name of the file to which the surface flux is exported + t (list): list of time values + data (list): list of surface flux values + title (str): title of the exported data + write_to_file (bool): True if the data is exported to a file """ def __init__( @@ -34,17 +42,12 @@ def __init__( writer = csv.writer(file) writer.writerow(["Time", f"{self.title}"]) - def compute(self, mesh, ds): + def compute(self, n, ds): self.value = fem.assemble_scalar( fem.form( - self.D - * ufl.dot(ufl.grad(self.field.solution), mesh.n) + -self.D + * ufl.dot(ufl.grad(self.field.solution), n) * ds(self.surface_subdomain.id) ) ) self.data.append(self.value) - - def write(self, t): - with open(self.filename, mode="a", newline="") as file: - writer = csv.writer(file) - writer.writerow([t, self.value]) From 22e3bf43f80a93db47f37382be4213a034cffc26 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:08:24 -0400 Subject: [PATCH 33/62] updated docs --- festim/exports/surface_quantity.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index a1cb38eec..9070aaa20 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -1,22 +1,26 @@ import festim as F +import csv class SurfaceQuantity: """Export SurfaceQuantity Args: + field (festim.Species): species for which the surface flux is computed + surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + filename (str, optional): name of the file to which the surface flux is exported Attributes: + field (festim.Species): species for which the surface flux is computed + surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + filename (str): name of the file to which the surface flux is exported + write_to_file (bool): True if the data is exported to a file """ def __init__(self, field, surface_subdomain, filename: str = None) -> None: self.field = field self.surface_subdomain = surface_subdomain self.filename = filename - self.ds = None - self.n = None - self.D = None - self.S = None @property def filename(self): @@ -62,3 +66,8 @@ def write_to_file(self): return False else: return True + + def write(self, t): + with open(self.filename, mode="a", newline="") as file: + writer = csv.writer(file) + writer.writerow([t, self.value]) From fbfaea865c049074ea0153388471933ba102118d Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:09:00 -0400 Subject: [PATCH 34/62] refactoring --- festim/hydrogen_transport_problem.py | 50 ++++++++++++---------------- 1 file changed, 22 insertions(+), 28 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 7584731b4..14a808eaf 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -249,13 +249,13 @@ def initialise_exports(self): spe_to_D_global_expr = {} # links species to D expression for export in self.exports: - if isinstance(export, F.SurfaceFlux): + if isinstance(export, F.SurfaceQuantity): if export.field in spe_to_D_global: # if already computed then use the same D D = spe_to_D_global[export.field] D_expr = spe_to_D_global_expr[export.field] else: - # if not computed then compute D and add it to the dict + # compute D and add it to the dict D, D_expr = self.define_D_global(export.field) spe_to_D_global[export.field] = D spe_to_D_global_expr[export.field] = D_expr @@ -541,30 +541,18 @@ def run(self): return self.times, self.flux_values - def iterate( - self, skip_post_processing=False - ): # TODO remove skip_post_processing flag, just temporary + def iterate(self): """Iterates the model for a given time step""" self.progress.update(self.dt.value) self.t.value += self.dt.value self.update_time_dependent_values() - # update global D if temperature time dependent or internal - # variables time dependent - species_not_updated = self.species.copy() # make a copy of the species - for export in self.exports: - if isinstance(export, F.SurfaceFlux): - # if the D of the species has not been updated yet - if export.field in species_not_updated: - export.D.interpolate(export.D_expr) - species_not_updated.remove(export.field) + # solve main problem + self.solver.solve(self.u) - # solve main problem - self.solver.solve(self.u) - - # post processing - self.post_processing() + # post processing + self.post_processing() # update previous solution self.u_n.x.array[:] = self.u.x.array[:] @@ -577,20 +565,27 @@ def update_time_dependent_values(self): else: self.temperature_fenics.interpolate(self.temperature_expr) - return times, flux_values + # update global D if temperature time dependent or internal + # variables time dependent + species_not_updated = self.species.copy() # make a copy of the species + for export in self.exports: + if isinstance(export, F.SurfaceFlux): + # if the D of the species has not been updated yet + if export.field in species_not_updated: + export.D.interpolate(export.D_expr) + species_not_updated.remove(export.field) + + for bc in self.boundary_conditions: + bc.update(t=t) def post_processing(self): - if not self.multispecies: - self.species[0].post_processing_solution = self.u - else: - for idx, spe in enumerate(self.species): - spe.post_processing_solution = self.u.sub(idx) + """Post processes the model""" for export in self.exports: # TODO if export type derived quantity - if isinstance(export, F.SurfaceFlux): + if isinstance(export, F.SurfaceQuantity): export.compute( - self.mesh, + self.mesh.n, self.ds, ) # update export data @@ -602,4 +597,3 @@ def post_processing(self): if isinstance(export, (F.VTXExport, F.XDMFExport)): export.write(float(self.t)) - From b5804e8d84a7ee0723334fec838ba9c5b85c80ee Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:09:26 -0400 Subject: [PATCH 35/62] rename test file --- test/test_surface_flux.py | 53 --------------- test/test_surface_quantity.py | 120 ++++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+), 53 deletions(-) delete mode 100644 test/test_surface_flux.py create mode 100644 test/test_surface_quantity.py diff --git a/test/test_surface_flux.py b/test/test_surface_flux.py deleted file mode 100644 index a6dc7ec42..000000000 --- a/test/test_surface_flux.py +++ /dev/null @@ -1,53 +0,0 @@ -import festim as F -import numpy as np -from petsc4py import PETSc - - -def surface_flux_export(): - """Test that the field attribute can be a string and is found in the species list""" - L = 3e-04 - my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) - left_surface = F.SurfaceSubdomain1D(id=1, x=0) - right_surface = F.SurfaceSubdomain1D(id=2, x=L) - mobile_H = F.Species("H") - - my_model = F.HydrogenTransportProblem() - - my_model.mesh = F.Mesh1D(vertices=np.linspace(0, L, num=1001)) - my_model.subdomains = [my_subdomain, left_surface, right_surface] - my_model.species = [mobile_H] - my_model.temperature = 500 - 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.settings = F.Settings(atol=1e10, rtol=1e-10, final_time=50) - my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) - - my_export = F.SurfaceFlux( - filename="my_surface_flux.csv", - field=mobile_H, - surface_subdomain=right_surface, - volume_subdomain=my_subdomain, - ) - my_model.exports = [my_export, F.XDMFExport("mobile_H.xdmf", field=mobile_H)] - - 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() - - my_model.run() - - -if __name__ == "__main__": - surface_flux_export() diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py new file mode 100644 index 000000000..d38a5d949 --- /dev/null +++ b/test/test_surface_quantity.py @@ -0,0 +1,120 @@ +import festim as F +import numpy as np +import ufl +from dolfinx.mesh import meshtags +from dolfinx import fem +import pytest +import os + + +def surface_flux_export_compute(): + """Test that the surface flux export computes the correct value""" + + # BUILD + L = 4.0 + D = 1.5 + my_mesh = F.Mesh1D(np.linspace(0, L, 10000)) + dummy_surface = F.SurfaceSubdomain1D(id=1, x=4) + + # define mesh ds measure + facet_indices = np.array( + dummy_surface.locate_boundary_facet_indices(my_mesh.mesh, 0), + dtype=np.int32, + ) + tags_facets = np.array( + [1], + dtype=np.int32, + ) + facet_meshtags = meshtags(my_mesh.mesh, 0, facet_indices, tags_facets) + ds = ufl.Measure("ds", domain=my_mesh.mesh, subdomain_data=facet_meshtags) + + # give function to species + V = fem.FunctionSpace(my_mesh.mesh, ("CG", 1)) + u = fem.Function(V) + u.interpolate(lambda x: 2 * x[0] ** 2 + 1) + + my_species = F.Species("H") + my_species.solution = u + + my_export = F.SurfaceFlux( + filename="my_surface_flux.csv", + field=my_species, + surface_subdomain=dummy_surface, + ) + my_export.D = D + + # RUN + my_export.compute(n=my_mesh.n, ds=ds) + + # TEST + expected_value = -D * 4 * dummy_surface.x + computed_value = my_export.value + + assert np.isclose(computed_value, expected_value, rtol=1e-2) + + +@pytest.mark.parametrize( + "input, expected_value", + [("export.csv", True), (None, False)], +) +def test_write_to_file_attribute(input, expected_value): + """Test that the write_to_file attribute is correctly set when a filename is given""" + my_export = F.SurfaceQuantity( + filename=input, + field="H", + surface_subdomain=1, + ) + + assert my_export.write_to_file is expected_value + + +def test_title_generation(tmp_path): + """Test that the title is made to be written to the header""" + my_export = F.SurfaceFlux( + filename=os.path.join(tmp_path, "my_export.csv"), + field=F.Species("TEST"), + surface_subdomain=F.SurfaceSubdomain1D(id=35, x=1), + ) + assert my_export.title == "Flux surface 35: TEST" + + +def test_filename_setter_raises_TypeError(tmp_path): + """Test that a TypeError is raised when the filename is not a string""" + + with pytest.raises(TypeError): + F.SurfaceFlux( + filename=os.path.join(tmp_path, 1), + field=F.Species("test"), + surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + ) + + +def test_filename_setter_raises_ValueError(tmp_path): + """Test that a ValueError is raised when the filename does not end with .csv""" + + with pytest.raises(ValueError): + F.SurfaceFlux( + filename=os.path.join(tmp_path, "my_export.xdmf"), + field=F.Species("test"), + surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + ) + + +def test_writer(tmp_path): + """Test that the writes values at each timestep""" + my_export = F.SurfaceFlux( + filename=os.path.join(tmp_path, "my_export.csv"), + field=F.Species("test"), + surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + ) + my_export.value = 2.0 + + for i in range(10): + my_export.write(i) + + # computed value should be range + 1 for the header + computed_value = len(np.genfromtxt(my_export.filename, delimiter=",")) + + expected_value = 11 + + assert computed_value == expected_value From 6c41bf39705fd20e168cb275ea9beec2503ebd25 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:14:35 -0400 Subject: [PATCH 36/62] remove time dependent bcs --- festim/boundary_conditions/dirichlet_bc.py | 3 -- test/test_dirichlet_bc.py | 41 ---------------------- 2 files changed, 44 deletions(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index 95a628849..6bfe2f466 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -42,7 +42,6 @@ def __init__(self, subdomain, value, species) -> None: self.value_fenics = None self.bc_expr = None - self.time_dependent = False @property def value_fenics(self): @@ -103,13 +102,11 @@ def create_value( self.value_fenics = F.as_fenics_constant( mesh=mesh, value=self.value(t=float(t)) ) - self.time_dependent = True else: self.value_fenics = fem.Function(function_space) kwargs = {} if "t" in arguments: kwargs["t"] = t - self.time_dependent = True if "x" in arguments: kwargs["x"] = x if "T" in arguments: diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 5727a4032..b34cd6a7f 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -377,44 +377,3 @@ def test_integration_with_a_multispecies_HTransportProblem(value_A, value_B): if isinstance(expected_value, Conditional): expected_value = float(expected_value) assert np.isclose(computed_value, expected_value) - - -@pytest.mark.parametrize( - "value", - [ - 1.0, - lambda t: t, - lambda t: 1.0 + t, - lambda x: 1.0 + x[0], - lambda x, t: 1.0 + x[0] + t, - lambda x, t, T: 1.0 + x[0] + t + T, - lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), - ], -) -def test_time_dependence_attribute(value): - subdomain = F.SurfaceSubdomain1D(1, x=0) - 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, value, "H") - my_model.boundary_conditions = [my_bc] - - 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) - - my_model.initialise() - - if isinstance(value, (float, int)): - assert not my_model.boundary_conditions[0].time_dependent - else: - arguments = value.__code__.co_varnames - if "t" in arguments: - assert my_model.boundary_conditions[0].time_dependent - else: - assert not my_model.boundary_conditions[0].time_dependent From 1a8ac516bc876acbc1d45efec7c7d10360006bff Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:15:50 -0400 Subject: [PATCH 37/62] remove from doc --- festim/boundary_conditions/dirichlet_bc.py | 1 - 1 file changed, 1 deletion(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index 2ab2d61b6..776535e24 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -24,7 +24,6 @@ class DirichletBC: fenics format bc_expr (fem.Expression): the expression of the boundary condition that is used to update the value_fenics - time_dependent (bool): True if the boundary condition is time dependent Usage: >>> from festim import DirichletBC From b7898a84a47b03424e9a2dc0b80c6c512ceb702a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 16:33:15 -0400 Subject: [PATCH 38/62] only export when needed, use surface quantity in tests --- test/test_multispecies_problem.py | 2 -- test/test_permeation_problem.py | 1 + test/test_surface_quantity.py | 9 ++++----- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index cbe2ce110..7e85a8fe1 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -88,12 +88,10 @@ def test_multispecies_permeation_problem(): my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20) outgassing_flux_spe_1 = F.SurfaceFlux( - filename="flux_spe_1.csv", field=spe_1, surface_subdomain=right_surface, ) outgassing_flux_spe_2 = F.SurfaceFlux( - filename="flux_spe_2.csv", field=spe_2, surface_subdomain=right_surface, ) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index cdf8382ef..a1571356f 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -157,6 +157,7 @@ def test_permeation_problem_multi_volume(tmp_path): ), ] outgassing_flux = F.SurfaceFlux( + filename=os.path.join(tmp_path, "outgassing_flux.csv"), field=mobile_H, surface_subdomain=right_surface, ) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index d38a5d949..4dc2310d7 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -82,7 +82,7 @@ def test_filename_setter_raises_TypeError(tmp_path): """Test that a TypeError is raised when the filename is not a string""" with pytest.raises(TypeError): - F.SurfaceFlux( + F.SurfaceQuantity( filename=os.path.join(tmp_path, 1), field=F.Species("test"), surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), @@ -93,7 +93,7 @@ def test_filename_setter_raises_ValueError(tmp_path): """Test that a ValueError is raised when the filename does not end with .csv""" with pytest.raises(ValueError): - F.SurfaceFlux( + F.SurfaceQuantity( filename=os.path.join(tmp_path, "my_export.xdmf"), field=F.Species("test"), surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), @@ -102,7 +102,7 @@ def test_filename_setter_raises_ValueError(tmp_path): def test_writer(tmp_path): """Test that the writes values at each timestep""" - my_export = F.SurfaceFlux( + my_export = F.SurfaceQuantity( filename=os.path.join(tmp_path, "my_export.csv"), field=F.Species("test"), surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), @@ -112,9 +112,8 @@ def test_writer(tmp_path): for i in range(10): my_export.write(i) - # computed value should be range + 1 for the header computed_value = len(np.genfromtxt(my_export.filename, delimiter=",")) - expected_value = 11 + expected_value = 10 assert computed_value == expected_value From de7440caf14e209b8f204b5f1f61a58f233681f1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:20:51 -0400 Subject: [PATCH 39/62] test species in list --- test/test_h_transport_problem.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index ab30b9e77..3a83869cd 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -207,3 +207,22 @@ def test_update_time_dependent_values_temperature(T_function, expected_values): computed_value = my_model.temperature_fenics.vector.array[-1] print(computed_value) assert np.isclose(computed_value, expected_values[i]) + + +def test_initialise_exports_find_species_with_one_field(): + """Test that a species can be found from the model species if given as a string""" + + # BUILD + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])), + species=[F.Species("H"), F.Species("D")], + temperature=1, + ) + surf = F.SurfaceSubdomain1D(id=1, x=1) + + my_model.exports = [F.SurfaceFlux(field="J", surface_subdomain=surf)] + my_model.define_function_spaces() + + # TEST + with pytest.raises(ValueError, match="Species J not found in list of species"): + my_model.initialise_exports() From c610f436a8e66ded65259bc487ddf59a79ef29b3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:21:23 -0400 Subject: [PATCH 40/62] new tests --- test/test_surface_quantity.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index 4dc2310d7..e8e5133c8 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -78,12 +78,12 @@ def test_title_generation(tmp_path): assert my_export.title == "Flux surface 35: TEST" -def test_filename_setter_raises_TypeError(tmp_path): +def test_filename_setter_raises_TypeError(): """Test that a TypeError is raised when the filename is not a string""" - with pytest.raises(TypeError): + with pytest.raises(TypeError, match="filename must be of type str"): F.SurfaceQuantity( - filename=os.path.join(tmp_path, 1), + filename=1, field=F.Species("test"), surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), ) @@ -100,6 +100,16 @@ def test_filename_setter_raises_ValueError(tmp_path): ) +def test_field_setter_raises_TypeError(): + """Test that a TypeError is raised when the field is not a F.Species""" + + with pytest.raises(TypeError): + F.SurfaceQuantity( + field=1, + surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + ) + + def test_writer(tmp_path): """Test that the writes values at each timestep""" my_export = F.SurfaceQuantity( From 41e25e2ad90ddc6745ac48bd5835cd2bde418794 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:21:43 -0400 Subject: [PATCH 41/62] only create file after species matches --- festim/hydrogen_transport_problem.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index dd50acbfb..75e055dd8 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -246,6 +246,8 @@ def initialise_exports(self): export.define_writer(MPI.COMM_WORLD) if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) + elif isinstance(export, F.SurfaceQuantity): + export.create_file() # compute diffusivity function for surface fluxes @@ -596,7 +598,7 @@ def post_processing(self): export.t.append(float(self.t)) # if filename given write export data to file - if export.write_to_file: + if export.filename is not None: export.write(t=float(self.t)) if isinstance(export, (F.VTXExport, F.XDMFExport)): From 62f653d279a9e6c1dc1ee366dad0ccca27d87afb Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:21:53 -0400 Subject: [PATCH 42/62] create file as method --- festim/exports/surface_flux.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index e189a14cb..b5e2434b1 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -33,15 +33,6 @@ def __init__( self.t = [] self.data = [] - self.title = "Flux surface {}: {}".format( - self.surface_subdomain.id, self.field.name - ) - - if self.write_to_file: - with open(self.filename, mode="w", newline="") as file: - writer = csv.writer(file) - writer.writerow(["Time", f"{self.title}"]) - def compute(self, n, ds): self.value = fem.assemble_scalar( fem.form( @@ -51,3 +42,13 @@ def compute(self, n, ds): ) ) self.data.append(self.value) + + def create_file(self): + self.title = "Flux surface {}: {}".format( + self.surface_subdomain.id, self.field.name + ) + + if self.write_to_file: + with open(self.filename, mode="w", newline="") as file: + writer = csv.writer(file) + writer.writerow(["t(s)", f"{self.title}"]) From 6233d388ee0d47c8f2c04f0783377619b09bd5a7 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:23:26 -0400 Subject: [PATCH 43/62] not needed --- festim/hydrogen_transport_problem.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 75e055dd8..fcec00df8 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -542,11 +542,6 @@ def run(self): while self.t.value < self.settings.final_time: self.iterate() - if self.multispecies: - self.flux_values = [self.flux_values_1, self.flux_values_2] - - return self.times, self.flux_values - def iterate(self): """Iterates the model for a given time step""" self.progress.update(self.dt.value) From f385a71fc0b269f56d2665e9f5631ad16fb1abf9 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:24:51 -0400 Subject: [PATCH 44/62] call create file --- test/test_surface_quantity.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index e8e5133c8..f07423c1c 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -75,6 +75,7 @@ def test_title_generation(tmp_path): field=F.Species("TEST"), surface_subdomain=F.SurfaceSubdomain1D(id=35, x=1), ) + my_export.create_file() assert my_export.title == "Flux surface 35: TEST" From 1e5c3f54b72170560c72ae7488718b22e4c878c1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:26:04 -0400 Subject: [PATCH 45/62] not accepting ints for now --- festim/exports/surface_flux.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index b5e2434b1..9d4083d4c 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -25,7 +25,7 @@ class SurfaceFlux(F.SurfaceQuantity): def __init__( self, field, - surface_subdomain: int or F.SurfaceSubdomain1D, + surface_subdomain: F.SurfaceSubdomain1D, filename: str = None, ) -> None: super().__init__(field, surface_subdomain, filename) From 5eef473df5bb48b9f03732407ae65bf82f512760 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:29:53 -0400 Subject: [PATCH 46/62] update global D in post processing --- festim/hydrogen_transport_problem.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index fcec00df8..9cc8e606f 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -566,6 +566,13 @@ def update_time_dependent_values(self): else: self.temperature_fenics.interpolate(self.temperature_expr) + for bc in self.boundary_conditions: + bc.update(t=t) + + def post_processing(self): + """Post processes the model""" + + if self.temperature_time_dependent: # update global D if temperature time dependent or internal # variables time dependent species_not_updated = self.species.copy() # make a copy of the species @@ -576,12 +583,6 @@ def update_time_dependent_values(self): export.D.interpolate(export.D_expr) species_not_updated.remove(export.field) - for bc in self.boundary_conditions: - bc.update(t=t) - - def post_processing(self): - """Post processes the model""" - for export in self.exports: # TODO if export type derived quantity if isinstance(export, F.SurfaceQuantity): From 2414d62b950fd4ae98958a66c6016a6f5ef0a09a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:39:09 -0400 Subject: [PATCH 47/62] surface instead of surface subdomain --- festim/exports/surface_flux.py | 19 ++++++++---------- festim/exports/surface_quantity.py | 26 +++++++++---------------- festim/hydrogen_transport_problem.py | 2 +- test/test_h_transport_problem.py | 2 +- test/test_multispecies_problem.py | 4 ++-- test/test_permeation_problem.py | 4 ++-- test/test_surface_quantity.py | 29 +++++++--------------------- 7 files changed, 30 insertions(+), 56 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index 9d4083d4c..1a3fdad88 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -9,26 +9,25 @@ class SurfaceFlux(F.SurfaceQuantity): Args: field (festim.Species): species for which the surface flux is computed - surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + surface (festim.SurfaceSubdomain1D): surface subdomain filename (str, optional): name of the file to which the surface flux is exported Attributes: field (festim.Species): species for which the surface flux is computed - surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + surface (festim.SurfaceSubdomain1D): surface subdomain filename (str): name of the file to which the surface flux is exported t (list): list of time values data (list): list of surface flux values title (str): title of the exported data - write_to_file (bool): True if the data is exported to a file """ def __init__( self, field, - surface_subdomain: F.SurfaceSubdomain1D, + surface: F.SurfaceSubdomain1D, filename: str = None, ) -> None: - super().__init__(field, surface_subdomain, filename) + super().__init__(field, surface, filename) self.t = [] self.data = [] @@ -38,17 +37,15 @@ def compute(self, n, ds): fem.form( -self.D * ufl.dot(ufl.grad(self.field.solution), n) - * ds(self.surface_subdomain.id) + * ds(self.surface.id) ) ) self.data.append(self.value) - def create_file(self): - self.title = "Flux surface {}: {}".format( - self.surface_subdomain.id, self.field.name - ) + def initialise_export(self): + self.title = "Flux surface {}: {}".format(self.surface.id, self.field.name) - if self.write_to_file: + if self.filename is not None: with open(self.filename, mode="w", newline="") as file: writer = csv.writer(file) writer.writerow(["t(s)", f"{self.title}"]) diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index 9070aaa20..3045d6bab 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -7,19 +7,18 @@ class SurfaceQuantity: Args: field (festim.Species): species for which the surface flux is computed - surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + surface (festim.SurfaceSubdomain1D): surface subdomain filename (str, optional): name of the file to which the surface flux is exported Attributes: field (festim.Species): species for which the surface flux is computed - surface_subdomain (festim.SurfaceSubdomain1D): surface subdomain + surface (festim.SurfaceSubdomain1D): surface subdomain filename (str): name of the file to which the surface flux is exported - write_to_file (bool): True if the data is exported to a file """ - def __init__(self, field, surface_subdomain, filename: str = None) -> None: + def __init__(self, field, surface, filename: str = None) -> None: self.field = field - self.surface_subdomain = surface_subdomain + self.surface = surface self.filename = filename @property @@ -37,16 +36,16 @@ def filename(self, value): self._filename = value @property - def surface_subdomain(self): - return self._surface_subdomain + def surface(self): + return self._surface - @surface_subdomain.setter - def surface_subdomain(self, value): + @surface.setter + def surface(self, value): if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance( value, bool ): raise TypeError("surface should be an int or F.SurfaceSubdomain1D") - self._surface_subdomain = value + self._surface = value @property def field(self): @@ -60,13 +59,6 @@ def field(self, value): self._field = value - @property - def write_to_file(self): - if self.filename is None: - return False - else: - return True - def write(self, t): with open(self.filename, mode="a", newline="") as file: writer = csv.writer(file) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 9cc8e606f..2610d43a8 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -247,7 +247,7 @@ def initialise_exports(self): if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) elif isinstance(export, F.SurfaceQuantity): - export.create_file() + export.initialise_export() # compute diffusivity function for surface fluxes diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 3a83869cd..55a5c07e2 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -220,7 +220,7 @@ def test_initialise_exports_find_species_with_one_field(): ) surf = F.SurfaceSubdomain1D(id=1, x=1) - my_model.exports = [F.SurfaceFlux(field="J", surface_subdomain=surf)] + my_model.exports = [F.SurfaceFlux(field="J", surface=surf)] my_model.define_function_spaces() # TEST diff --git a/test/test_multispecies_problem.py b/test/test_multispecies_problem.py index 7e85a8fe1..c2408a26d 100644 --- a/test/test_multispecies_problem.py +++ b/test/test_multispecies_problem.py @@ -89,11 +89,11 @@ def test_multispecies_permeation_problem(): outgassing_flux_spe_1 = F.SurfaceFlux( field=spe_1, - surface_subdomain=right_surface, + surface=right_surface, ) outgassing_flux_spe_2 = F.SurfaceFlux( field=spe_2, - surface_subdomain=right_surface, + surface=right_surface, ) my_model.exports = [outgassing_flux_spe_1, outgassing_flux_spe_2] my_model.initialise() diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index a1571356f..eb90c22c0 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -63,7 +63,7 @@ def test_permeation_problem(mesh_size=1001): ] outgassing_flux = F.SurfaceFlux( field=mobile_H, - surface_subdomain=right_surface, + surface=right_surface, ) my_model.exports = [ F.XDMFExport("mobile_concentration.xdmf", field=mobile_H), @@ -159,7 +159,7 @@ def test_permeation_problem_multi_volume(tmp_path): outgassing_flux = F.SurfaceFlux( filename=os.path.join(tmp_path, "outgassing_flux.csv"), field=mobile_H, - surface_subdomain=right_surface, + surface=right_surface, ) my_model.exports = [ F.VTXExport( diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index f07423c1c..25b8ad782 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -39,7 +39,7 @@ def surface_flux_export_compute(): my_export = F.SurfaceFlux( filename="my_surface_flux.csv", field=my_species, - surface_subdomain=dummy_surface, + surface=dummy_surface, ) my_export.D = D @@ -53,29 +53,14 @@ def surface_flux_export_compute(): assert np.isclose(computed_value, expected_value, rtol=1e-2) -@pytest.mark.parametrize( - "input, expected_value", - [("export.csv", True), (None, False)], -) -def test_write_to_file_attribute(input, expected_value): - """Test that the write_to_file attribute is correctly set when a filename is given""" - my_export = F.SurfaceQuantity( - filename=input, - field="H", - surface_subdomain=1, - ) - - assert my_export.write_to_file is expected_value - - def test_title_generation(tmp_path): """Test that the title is made to be written to the header""" my_export = F.SurfaceFlux( filename=os.path.join(tmp_path, "my_export.csv"), field=F.Species("TEST"), - surface_subdomain=F.SurfaceSubdomain1D(id=35, x=1), + surface=F.SurfaceSubdomain1D(id=35, x=1), ) - my_export.create_file() + my_export.initialise_export() assert my_export.title == "Flux surface 35: TEST" @@ -86,7 +71,7 @@ def test_filename_setter_raises_TypeError(): F.SurfaceQuantity( filename=1, field=F.Species("test"), - surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + surface=F.SurfaceSubdomain1D(id=1, x=0), ) @@ -97,7 +82,7 @@ def test_filename_setter_raises_ValueError(tmp_path): F.SurfaceQuantity( filename=os.path.join(tmp_path, "my_export.xdmf"), field=F.Species("test"), - surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + surface=F.SurfaceSubdomain1D(id=1, x=0), ) @@ -107,7 +92,7 @@ def test_field_setter_raises_TypeError(): with pytest.raises(TypeError): F.SurfaceQuantity( field=1, - surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + surface=F.SurfaceSubdomain1D(id=1, x=0), ) @@ -116,7 +101,7 @@ def test_writer(tmp_path): my_export = F.SurfaceQuantity( filename=os.path.join(tmp_path, "my_export.csv"), field=F.Species("test"), - surface_subdomain=F.SurfaceSubdomain1D(id=1, x=0), + surface=F.SurfaceSubdomain1D(id=1, x=0), ) my_export.value = 2.0 From 9852021726c6e94e48576c367b52d96e22abe33d Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 17:41:04 -0400 Subject: [PATCH 48/62] updated doc strings --- festim/hydrogen_transport_problem.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 2610d43a8..da28d81e3 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -53,6 +53,10 @@ class HydrogenTransportProblem: that is used to update the temperature_fenics temperature_time_dependent (bool): True if the temperature is time dependent + V_DG_0 (dolfinx.fem.FunctionSpace): A DG function space of degree 0 + over domain + V_DG_1 (dolfinx.fem.FunctionSpace): A DG function space of degree 1 + over domain Usage: From 8ddc1615c5abede19e38167ba569244057e277de Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 18:23:19 -0400 Subject: [PATCH 49/62] use temperature_fenics --- 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 da28d81e3..5681a9e7e 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -298,8 +298,9 @@ def define_D_global(self, species): # create global D function D = fem.Function(self.V_DG_1) + expr = D_0 * ufl.exp( - -E_D / F.as_fenics_constant(F.k_B, self.mesh.mesh) / self.temperature + -E_D / F.as_fenics_constant(F.k_B, self.mesh.mesh) / self.temperature_fenics ) D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) D.interpolate(D_expr) From e1abfabff2f75f23c3ae701a139c281c73ad75a6 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 18:23:46 -0400 Subject: [PATCH 50/62] test defining D_global --- test/test_h_transport_problem.py | 78 ++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 55a5c07e2..1b51ba1af 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -226,3 +226,81 @@ def test_initialise_exports_find_species_with_one_field(): # TEST with pytest.raises(ValueError, match="Species J not found in list of species"): my_model.initialise_exports() + + +def test_define_D_global_different_temperatures(): + D_0, E_D = 1.5, 0.1 + my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") + surf = F.SurfaceSubdomain1D(id=1, x=0) + H = F.Species("H") + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(np.linspace(0, 4, num=101)), + subdomains=[ + F.VolumeSubdomain1D(id=1, borders=[0, 2], material=my_mat), + F.VolumeSubdomain1D(id=2, borders=[2, 4], material=my_mat), + ], + species=[H], + temperature=lambda x: 100.0 * x[0] + 50, + exports=[ + F.SurfaceFlux( + field=H, + surface=surf, + ), + ], + ) + + my_model.define_function_spaces() + my_model.define_markers_and_measures() + my_model.define_temperature() + + D_computed, D_expr = my_model.define_D_global(H) + + computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] + + D_analytical_left = D_0 * np.exp(-E_D / (F.k_B * 50)) + D_analytical_right = D_0 * np.exp(-E_D / (F.k_B * 450)) + + expected_values = [D_analytical_left, D_analytical_right] + + assert np.isclose(computed_values, expected_values).all() + + +def test_define_D_global_different_materials(): + D_0_left, E_D_left = 1.0, 0.1 + D_0_right, E_D_right = 2.0, 0.2 + my_mat_L = F.Material(D_0=D_0_left, E_D=E_D_left, name="my_mat_L") + my_mat_R = F.Material(D_0=D_0_right, E_D=E_D_right, name="my_mat_R") + surf = F.SurfaceSubdomain1D(id=1, x=0) + H = F.Species("H") + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(np.linspace(0, 4, num=101)), + subdomains=[ + F.VolumeSubdomain1D(id=1, borders=[0, 2], material=my_mat_L), + F.VolumeSubdomain1D(id=2, borders=[2, 4], material=my_mat_R), + ], + species=[F.Species("H"), F.Species("D")], + temperature=500, + exports=[ + F.SurfaceFlux( + field=H, + surface=surf, + ), + ], + ) + + my_model.define_function_spaces() + my_model.define_markers_and_measures() + my_model.define_temperature() + + D_computed, D_expr = my_model.define_D_global(H) + + computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] + + D_analytical_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature)) + D_analytical_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature)) + + expected_values = [D_analytical_left, D_analytical_right] + + assert np.isclose(computed_values, expected_values).all() From 8248d7888ae54dbd00160e451ac28b433704e546 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 19:13:41 -0400 Subject: [PATCH 51/62] test surface setter --- test/test_surface_quantity.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index 25b8ad782..0ff67f19e 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -113,3 +113,16 @@ def test_writer(tmp_path): expected_value = 10 assert computed_value == expected_value + + +def test_surface_setter_raises_TypeError(): + """Test that a TypeError is raised when the surface is not a + F.SurfaceSubdomain1D""" + + with pytest.raises( + TypeError, match="surface should be an int or F.SurfaceSubdomain1D" + ): + F.SurfaceQuantity( + field=F.Species("H"), + surface="1", + ) From b003bde119a1ea552f63a4186f3481263c18f201 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 2 Nov 2023 19:30:13 -0400 Subject: [PATCH 52/62] test global D with same species --- test/test_h_transport_problem.py | 39 ++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 1b51ba1af..36a3df1ce 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -304,3 +304,42 @@ def test_define_D_global_different_materials(): expected_values = [D_analytical_left, D_analytical_right] assert np.isclose(computed_values, expected_values).all() + + +def test_initialise_exports_multiple_exports_same_species(): + D_0, E_D = 1.5, 0.1 + my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") + surf_1 = F.SurfaceSubdomain1D(id=1, x=0) + surf_2 = F.SurfaceSubdomain1D(id=1, x=4) + H = F.Species("H") + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(np.linspace(0, 4, num=101)), + subdomains=[ + F.VolumeSubdomain1D(id=1, borders=[0, 2], material=my_mat), + F.VolumeSubdomain1D(id=2, borders=[2, 4], material=my_mat), + ], + species=[H], + temperature=500, + exports=[ + F.SurfaceFlux( + field=H, + surface=surf_1, + ), + F.SurfaceFlux( + field=H, + surface=surf_2, + ), + ], + ) + + my_model.define_function_spaces() + my_model.define_markers_and_measures() + my_model.define_temperature() + my_model.initialise_exports() + + Ds = [] + for export in my_model.exports: + Ds.append(export.D) + + assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) From 387832428a19e25fcc8a43ebd0697b591280e11c Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 4 Nov 2023 10:43:05 -0400 Subject: [PATCH 53/62] test post processing --- test/test_h_transport_problem.py | 47 ++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 36a3df1ce..9b22598f8 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -343,3 +343,50 @@ def test_initialise_exports_multiple_exports_same_species(): Ds.append(export.D) assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) + + +def test_post_processing_update_D_global(): + """Test that the D_global attribute is updated at each time + step when temperture is time dependent""" + my_mesh = F.Mesh1D(np.linspace(0, 1, num=11)) + my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat") + surf = F.SurfaceSubdomain1D(id=1, x=1) + + # create species and interpolate solution + H = F.Species("H") + V = fem.FunctionSpace(my_mesh.mesh, ("CG", 1)) + u = fem.Function(V) + u.interpolate(lambda x: 2 * x[0] ** 2 + 1) + H.solution = u + + my_export = F.SurfaceFlux( + field=H, + surface=surf, + ) + + # Build the model + my_model = F.HydrogenTransportProblem( + mesh=my_mesh, + subdomains=[F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat), surf], + species=[H], + temperature=lambda t: 500 * t, + exports=[my_export], + ) + + my_model.define_function_spaces() + my_model.define_markers_and_measures() + my_model.t = fem.Constant(my_model.mesh.mesh, 1.0) + my_model.define_temperature() + my_model.initialise_exports() + + # RUN + my_model.post_processing() + value_t_1 = my_export.D.x.array[-1] + + my_model.t = fem.Constant(my_model.mesh.mesh, 2.0) + my_model.update_time_dependent_values() + my_model.post_processing() + value_t_2 = my_export.D.x.array[-1] + + # TEST + assert value_t_1 != value_t_2 From bda7ad845b077ad6dcd31edbe82b320e47e5e548 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Sat, 4 Nov 2023 10:44:05 -0400 Subject: [PATCH 54/62] Apply suggestions from code review 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_h_transport_problem.py | 4 +--- test/test_surface_quantity.py | 13 +++++++------ 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 9b22598f8..045c9b529 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -338,9 +338,7 @@ def test_initialise_exports_multiple_exports_same_species(): my_model.define_temperature() my_model.initialise_exports() - Ds = [] - for export in my_model.exports: - Ds.append(export.D) + Ds = [export.D for export in my_model_exports] assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index 0ff67f19e..3d0420342 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -30,11 +30,11 @@ def surface_flux_export_compute(): # give function to species V = fem.FunctionSpace(my_mesh.mesh, ("CG", 1)) - u = fem.Function(V) - u.interpolate(lambda x: 2 * x[0] ** 2 + 1) + c = fem.Function(V) + c.interpolate(lambda x: 2 * x[0] ** 2 + 1) my_species = F.Species("H") - my_species.solution = u + my_species.solution = c my_export = F.SurfaceFlux( filename="my_surface_flux.csv", @@ -47,6 +47,7 @@ def surface_flux_export_compute(): my_export.compute(n=my_mesh.n, ds=ds) # TEST + # flux = -D grad(c)_ \cdot n = -D dc/dx = -D * 4 * x expected_value = -D * 4 * dummy_surface.x computed_value = my_export.value @@ -108,11 +109,11 @@ def test_writer(tmp_path): for i in range(10): my_export.write(i) - computed_value = len(np.genfromtxt(my_export.filename, delimiter=",")) + file_length = len(np.genfromtxt(my_export.filename, delimiter=",")) - expected_value = 10 + expected_length = i + 1 - assert computed_value == expected_value + assert file_length == expected_length def test_surface_setter_raises_TypeError(): From 52da3d007267bd31477da0dafd82b6edb721cbae Mon Sep 17 00:00:00 2001 From: J Dark Date: Sat, 4 Nov 2023 10:59:22 -0400 Subject: [PATCH 55/62] fixes --- test/test_h_transport_problem.py | 2 +- test/test_surface_quantity.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 045c9b529..985e2d713 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -338,7 +338,7 @@ def test_initialise_exports_multiple_exports_same_species(): my_model.define_temperature() my_model.initialise_exports() - Ds = [export.D for export in my_model_exports] + Ds = [export.D for export in my_model.exports] assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index 3d0420342..4bcc192b6 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -99,19 +99,19 @@ def test_field_setter_raises_TypeError(): def test_writer(tmp_path): """Test that the writes values at each timestep""" - my_export = F.SurfaceQuantity( + my_export = F.SurfaceFlux( filename=os.path.join(tmp_path, "my_export.csv"), field=F.Species("test"), surface=F.SurfaceSubdomain1D(id=1, x=0), ) my_export.value = 2.0 + my_export.initialise_export() for i in range(10): my_export.write(i) - file_length = len(np.genfromtxt(my_export.filename, delimiter=",")) - expected_length = i + 1 + expected_length = i + 2 assert file_length == expected_length From 96d47e03d52d193a287cf72deed8f2b77dd2eafc Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 11:28:19 -0500 Subject: [PATCH 56/62] new test and doc strings --- test/test_h_transport_problem.py | 56 ++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 985e2d713..13eb718bd 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -229,6 +229,8 @@ def test_initialise_exports_find_species_with_one_field(): def test_define_D_global_different_temperatures(): + """Test that the D_global attribute is correctly defined when the temperature + is different in the volume subdomains""" D_0, E_D = 1.5, 0.1 my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") surf = F.SurfaceSubdomain1D(id=1, x=0) @@ -267,6 +269,8 @@ def test_define_D_global_different_temperatures(): def test_define_D_global_different_materials(): + """Test that the D_global attribute is correctly defined when the material + is different in the volume subdomains""" D_0_left, E_D_left = 1.0, 0.1 D_0_right, E_D_right = 2.0, 0.2 my_mat_L = F.Material(D_0=D_0_left, E_D=E_D_left, name="my_mat_L") @@ -307,6 +311,8 @@ def test_define_D_global_different_materials(): def test_initialise_exports_multiple_exports_same_species(): + """Test that the D attribute is the same for multiple exports of the same species, + and that D_global is only created once per species""" D_0, E_D = 1.5, 0.1 my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") surf_1 = F.SurfaceSubdomain1D(id=1, x=0) @@ -343,6 +349,56 @@ def test_initialise_exports_multiple_exports_same_species(): assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) +def test_define_D_global_multispecies(): + """Test that the D_global attribute is correctly defined when there are multiple + species in one subdomain""" + A = F.Species("A") + B = F.Species("B") + + D_0_A, D_0_B = 1.0, 2.0 + 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}, name="my_mat" + ) + surf = F.SurfaceSubdomain1D(id=1, x=1) + + my_model = F.HydrogenTransportProblem( + mesh=F.Mesh1D(np.linspace(0, 1, num=101)), + subdomains=[ + F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat), + ], + species=[F.Species("A"), F.Species("B")], + temperature=500, + exports=[ + F.SurfaceFlux( + field=A, + surface=surf, + ), + F.SurfaceFlux( + field=A, + surface=surf, + ), + ], + ) + + my_model.define_function_spaces() + my_model.define_markers_and_measures() + my_model.define_temperature() + + D_A_computed, D_A_expr = my_model.define_D_global(A) + D_B_computed, D_B_expr = my_model.define_D_global(B) + + computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]] + + D_analytical_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature)) + D_analytical_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature)) + + expected_values = [D_analytical_A, D_analytical_B] + + assert np.isclose(computed_values, expected_values).all() + + def test_post_processing_update_D_global(): """Test that the D_global attribute is updated at each time step when temperture is time dependent""" From cf41ff464ee6def4657b608a162f28869e0c8955 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Sun, 5 Nov 2023 15:38:17 -0500 Subject: [PATCH 57/62] Apply suggestions from code review 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/exports/surface_flux.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index 1a3fdad88..3585344de 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -18,12 +18,12 @@ class SurfaceFlux(F.SurfaceQuantity): filename (str): name of the file to which the surface flux is exported t (list): list of time values data (list): list of surface flux values - title (str): title of the exported data + title (str): title of the exported data in the csv file """ def __init__( self, - field, + field: F.Species, surface: F.SurfaceSubdomain1D, filename: str = None, ) -> None: From e0d568a7ed63be28679bccb6a3d1558170162acf Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 16:24:57 -0500 Subject: [PATCH 58/62] changes from reveiw --- festim/exports/surface_flux.py | 16 ++++----- festim/exports/surface_quantity.py | 7 ++-- test/test_h_transport_problem.py | 53 +++++++++--------------------- test/test_permeation_problem.py | 1 + test/test_surface_quantity.py | 23 ++++++++----- 5 files changed, 44 insertions(+), 56 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index 1a3fdad88..c0710948d 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -16,9 +16,6 @@ class SurfaceFlux(F.SurfaceQuantity): field (festim.Species): species for which the surface flux is computed surface (festim.SurfaceSubdomain1D): surface subdomain filename (str): name of the file to which the surface flux is exported - t (list): list of time values - data (list): list of surface flux values - title (str): title of the exported data """ def __init__( @@ -29,10 +26,13 @@ def __init__( ) -> None: super().__init__(field, surface, filename) - self.t = [] - self.data = [] - def compute(self, n, ds): + """Computes the value of the surface flux at the surface + + Args: + n (ufl.geometry.FacetNormal): normal vector to the surface + ds (ufl.Measure): surface measure of the model + """ self.value = fem.assemble_scalar( fem.form( -self.D @@ -43,9 +43,9 @@ def compute(self, n, ds): self.data.append(self.value) def initialise_export(self): - self.title = "Flux surface {}: {}".format(self.surface.id, self.field.name) + title = "Flux surface {}: {}".format(self.surface.id, self.field.name) if self.filename is not None: with open(self.filename, mode="w", newline="") as file: writer = csv.writer(file) - writer.writerow(["t(s)", f"{self.title}"]) + writer.writerow(["t(s)", f"{title}"]) diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index 3045d6bab..c03e2bfc2 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -21,6 +21,9 @@ def __init__(self, field, surface, filename: str = None) -> None: self.surface = surface self.filename = filename + self.t = [] + self.data = [] + @property def filename(self): return self._filename @@ -31,8 +34,8 @@ def filename(self, value): self._filename = None elif not isinstance(value, str): raise TypeError("filename must be of type str") - elif not value.endswith(".csv"): - raise ValueError("filename must end with .csv") + elif not value.endswith(".csv") and not value.endswith(".txt"): + raise ValueError("filename must end with .csv or .txt") self._filename = value @property diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 13eb718bd..e3fc55073 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -229,7 +229,7 @@ def test_initialise_exports_find_species_with_one_field(): def test_define_D_global_different_temperatures(): - """Test that the D_global attribute is correctly defined when the temperature + """Test that the D_global object is correctly defined when the temperature is different in the volume subdomains""" D_0, E_D = 1.5, 0.1 my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") @@ -244,12 +244,6 @@ def test_define_D_global_different_temperatures(): ], species=[H], temperature=lambda x: 100.0 * x[0] + 50, - exports=[ - F.SurfaceFlux( - field=H, - surface=surf, - ), - ], ) my_model.define_function_spaces() @@ -269,13 +263,12 @@ def test_define_D_global_different_temperatures(): def test_define_D_global_different_materials(): - """Test that the D_global attribute is correctly defined when the material + """Test that the D_global object is correctly defined when the material is different in the volume subdomains""" D_0_left, E_D_left = 1.0, 0.1 D_0_right, E_D_right = 2.0, 0.2 my_mat_L = F.Material(D_0=D_0_left, E_D=E_D_left, name="my_mat_L") my_mat_R = F.Material(D_0=D_0_right, E_D=E_D_right, name="my_mat_R") - surf = F.SurfaceSubdomain1D(id=1, x=0) H = F.Species("H") my_model = F.HydrogenTransportProblem( @@ -284,14 +277,8 @@ def test_define_D_global_different_materials(): F.VolumeSubdomain1D(id=1, borders=[0, 2], material=my_mat_L), F.VolumeSubdomain1D(id=2, borders=[2, 4], material=my_mat_R), ], - species=[F.Species("H"), F.Species("D")], + species=[H], temperature=500, - exports=[ - F.SurfaceFlux( - field=H, - surface=surf, - ), - ], ) my_model.define_function_spaces() @@ -302,17 +289,19 @@ def test_define_D_global_different_materials(): computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] - D_analytical_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature)) - D_analytical_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature)) + D_expected_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature)) + D_expected_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature)) - expected_values = [D_analytical_left, D_analytical_right] + expected_values = [D_expected_left, D_expected_right] assert np.isclose(computed_values, expected_values).all() def test_initialise_exports_multiple_exports_same_species(): - """Test that the D attribute is the same for multiple exports of the same species, - and that D_global is only created once per species""" + """Test that the diffusion coefficient within the D_global object function is the same + for multiple exports of the same species, and that D_global object is only + created once per species""" + D_0, E_D = 1.5, 0.1 my_mat = F.Material(D_0=D_0, E_D=E_D, name="my_mat") surf_1 = F.SurfaceSubdomain1D(id=1, x=0) @@ -346,11 +335,11 @@ def test_initialise_exports_multiple_exports_same_species(): Ds = [export.D for export in my_model.exports] - assert np.isclose(Ds[0].x.array[0], Ds[1].x.array[0]) + assert Ds[0].x.array[0] == Ds[1].x.array[0] def test_define_D_global_multispecies(): - """Test that the D_global attribute is correctly defined when there are multiple + """Test that the D_global object is correctly defined when there are multiple species in one subdomain""" A = F.Species("A") B = F.Species("B") @@ -370,16 +359,6 @@ def test_define_D_global_multispecies(): ], species=[F.Species("A"), F.Species("B")], temperature=500, - exports=[ - F.SurfaceFlux( - field=A, - surface=surf, - ), - F.SurfaceFlux( - field=A, - surface=surf, - ), - ], ) my_model.define_function_spaces() @@ -391,16 +370,16 @@ def test_define_D_global_multispecies(): computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]] - D_analytical_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature)) - D_analytical_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature)) + D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature)) + D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature)) - expected_values = [D_analytical_A, D_analytical_B] + expected_values = [D_expected_A, D_expected_B] assert np.isclose(computed_values, expected_values).all() def test_post_processing_update_D_global(): - """Test that the D_global attribute is updated at each time + """Test that the D_global object is updated at each time step when temperture is time dependent""" my_mesh = F.Mesh1D(np.linspace(0, 1, num=11)) my_mat = F.Material(D_0=1.5, E_D=0.1, name="my_mat") diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index eb90c22c0..624ce635e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -62,6 +62,7 @@ def test_permeation_problem(mesh_size=1001): ), ] outgassing_flux = F.SurfaceFlux( + filename="outgassing_flux.txt", field=mobile_H, surface=right_surface, ) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index 4bcc192b6..dee0ea5f1 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -37,7 +37,6 @@ def surface_flux_export_compute(): my_species.solution = c my_export = F.SurfaceFlux( - filename="my_surface_flux.csv", field=my_species, surface=dummy_surface, ) @@ -54,15 +53,20 @@ def surface_flux_export_compute(): assert np.isclose(computed_value, expected_value, rtol=1e-2) -def test_title_generation(tmp_path): - """Test that the title is made to be written to the header""" +@pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"]) +def test_title_generation(tmp_path, value): + """Test that the title is made to be written to the header in a csv or txt file""" my_export = F.SurfaceFlux( - filename=os.path.join(tmp_path, "my_export.csv"), + filename=os.path.join(tmp_path, f"{value}"), field=F.Species("TEST"), surface=F.SurfaceSubdomain1D(id=35, x=1), ) my_export.initialise_export() - assert my_export.title == "Flux surface 35: TEST" + title = np.genfromtxt(my_export.filename, delimiter=",", max_rows=1, dtype=str) + + expected_title = "Flux surface 35: TEST" + + assert title[1] == expected_title def test_filename_setter_raises_TypeError(): @@ -77,7 +81,7 @@ def test_filename_setter_raises_TypeError(): def test_filename_setter_raises_ValueError(tmp_path): - """Test that a ValueError is raised when the filename does not end with .csv""" + """Test that a ValueError is raised when the filename does not end with .csv or .txt""" with pytest.raises(ValueError): F.SurfaceQuantity( @@ -97,10 +101,11 @@ def test_field_setter_raises_TypeError(): ) -def test_writer(tmp_path): - """Test that the writes values at each timestep""" +@pytest.mark.parametrize("value", ["my_export.csv", "my_export.txt"]) +def test_writer(tmp_path, value): + """Test that the writes values at each timestep to either a csv or txt file""" my_export = F.SurfaceFlux( - filename=os.path.join(tmp_path, "my_export.csv"), + filename=os.path.join(tmp_path, f"{value}"), field=F.Species("test"), surface=F.SurfaceSubdomain1D(id=1, x=0), ) From 4b20b0542b51308ebe1059ecf41544cd9f7fa707 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 16:25:11 -0500 Subject: [PATCH 59/62] test new materials methods --- festim/material.py | 4 +-- test/test_material.py | 69 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 2 deletions(-) diff --git a/festim/material.py b/festim/material.py index 057760434..401baa01f 100644 --- a/festim/material.py +++ b/festim/material.py @@ -59,7 +59,7 @@ def get_D_0(self, species=None): raise ValueError(f"{species} is not in D_0 keys") else: - raise ValueError("D_0 must be either a float or a dict") + raise TypeError("D_0 must be either a float, int or a dict") def get_E_D(self, species=None): """Returns the activation energy of the diffusion coefficient @@ -85,7 +85,7 @@ def get_E_D(self, species=None): raise ValueError(f"{species} is not in E_D keys") else: - raise ValueError("E_D must be either a float or a dict") + raise TypeError("E_D must be either a float, int or a dict") def get_diffusion_coefficient(self, mesh, temperature, species=None): """Defines the diffusion coefficient diff --git a/test/test_material.py b/test/test_material.py index 265bd9afb..c897a7daa 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -98,6 +98,7 @@ def test_multispecies_dict_different_keys(): def test_D_0_type_raises_error(): """Test that a value error is raised in the get_diffusion_coefficient function""" + # TODO remove this when material class is updated A = F.Species("A") my_mat = F.Material(D_0=[1, 1], E_D=0.1) @@ -128,3 +129,71 @@ def test_error_raised_when_species_not_not_in_D_0_dict(): with pytest.raises(ValueError, match="J is not in D_0 keys"): my_mat.get_diffusion_coefficient(test_mesh.mesh, 500, species=J) + + +def test_D_0_raises_ValueError_if_species_not_provided_in_dict(): + """Test that a value error is raised in the get_diffusion_coefficient + function""" + # TODO remove this when material class is updated + A = F.Species("A") + B = F.Species("B") + my_mat = F.Material(D_0={A: 1, B: 2}, E_D=1) + + with pytest.raises(ValueError, match="species must be provided if D_0 is a dict"): + my_mat.get_D_0() + + +def test_D_0_raises_ValueError_if_species_given_not_in_dict_keys(): + """Test that a value error is raised in the get_diffusion_coefficient + function""" + # TODO remove this when material class is updated + A = F.Species("A") + B = F.Species("B") + J = F.Species("J") + my_mat = F.Material(D_0={A: 1, B: 2}, E_D=1) + + with pytest.raises(ValueError, match="J is not in D_0 keys"): + my_mat.get_D_0(species=J) + + +def test_raises_TypeError_when_D_0_is_not_correct_type(): + """Test that a TypeError is raised when D_0 is not a float or a dict""" + + my_mat = F.Material(D_0=[1, 2], E_D=1) + + with pytest.raises(TypeError, match="D_0 must be either a float, int or a dict"): + my_mat.get_D_0() + + +def test_E_D_raises_ValueError_if_species_not_provided_in_dict(): + """Test that a value error is raised in the get_diffusion_coefficient + function""" + # TODO remove this when material class is updated + A = F.Species("A") + B = F.Species("B") + my_mat = F.Material(D_0=1, E_D={A: 1, B: 2}) + + with pytest.raises(ValueError, match="species must be provided if E_D is a dict"): + my_mat.get_E_D() + + +def test_E_D_raises_ValueError_if_species_given_not_in_dict_keys(): + """Test that a value error is raised in the get_diffusion_coefficient + function""" + # TODO remove this when material class is updated + A = F.Species("A") + B = F.Species("B") + J = F.Species("J") + my_mat = F.Material(D_0=1, E_D={A: 1, B: 2}) + + with pytest.raises(ValueError, match="J is not in E_D keys"): + my_mat.get_E_D(species=J) + + +def test_raises_TypeError_when_E_D_is_not_correct_type(): + """Test that a TypeError is raised when E_D is not a float or a dict""" + + my_mat = F.Material(D_0=1, E_D=[1, 2]) + + with pytest.raises(TypeError, match="E_D must be either a float, int or a dict"): + my_mat.get_E_D() From ee155e91498c4dec2984105a7a5208f887115514 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 19:04:22 -0500 Subject: [PATCH 60/62] export one first write --- festim/exports/surface_flux.py | 8 -------- festim/exports/surface_quantity.py | 9 +++++++++ festim/hydrogen_transport_problem.py | 2 -- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index 96403b8c6..c21db46ba 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -41,11 +41,3 @@ def compute(self, n, ds): ) ) self.data.append(self.value) - - def initialise_export(self): - title = "Flux surface {}: {}".format(self.surface.id, self.field.name) - - if self.filename is not None: - with open(self.filename, mode="w", newline="") as file: - writer = csv.writer(file) - writer.writerow(["t(s)", f"{title}"]) diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index c03e2bfc2..3ff9b6a56 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -1,5 +1,6 @@ import festim as F import csv +import os class SurfaceQuantity: @@ -63,6 +64,14 @@ def field(self, value): self._field = value def write(self, t): + if not os.path.isfile(self.filename): + title = "Flux surface {}: {}".format(self.surface.id, self.field.name) + + if self.filename is not None: + with open(self.filename, mode="w", newline="") as file: + writer = csv.writer(file) + writer.writerow(["t(s)", f"{title}"]) + with open(self.filename, mode="a", newline="") as file: writer = csv.writer(file) writer.writerow([t, self.value]) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 5681a9e7e..7307ab7f1 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -250,8 +250,6 @@ def initialise_exports(self): export.define_writer(MPI.COMM_WORLD) if isinstance(export, F.XDMFExport): export.writer.write_mesh(self.mesh.mesh) - elif isinstance(export, F.SurfaceQuantity): - export.initialise_export() # compute diffusivity function for surface fluxes From 849de0419cf2a707df5322d6f9891e2c0a5d3de6 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 19:04:33 -0500 Subject: [PATCH 61/62] update test --- test/test_surface_quantity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index dee0ea5f1..b0801f862 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -61,7 +61,8 @@ def test_title_generation(tmp_path, value): field=F.Species("TEST"), surface=F.SurfaceSubdomain1D(id=35, x=1), ) - my_export.initialise_export() + my_export.value = 2.0 + my_export.write(0) title = np.genfromtxt(my_export.filename, delimiter=",", max_rows=1, dtype=str) expected_title = "Flux surface 35: TEST" @@ -110,7 +111,6 @@ def test_writer(tmp_path, value): surface=F.SurfaceSubdomain1D(id=1, x=0), ) my_export.value = 2.0 - my_export.initialise_export() for i in range(10): my_export.write(i) From 3f818283ad56962755188ce33cc786896de81340 Mon Sep 17 00:00:00 2001 From: J Dark Date: Sun, 5 Nov 2023 19:26:15 -0500 Subject: [PATCH 62/62] updated docs --- festim/exports/surface_flux.py | 1 - festim/exports/surface_quantity.py | 5 +++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/festim/exports/surface_flux.py b/festim/exports/surface_flux.py index c21db46ba..93831d490 100644 --- a/festim/exports/surface_flux.py +++ b/festim/exports/surface_flux.py @@ -1,7 +1,6 @@ from dolfinx import fem import festim as F import ufl -import csv class SurfaceFlux(F.SurfaceQuantity): diff --git a/festim/exports/surface_quantity.py b/festim/exports/surface_quantity.py index 3ff9b6a56..148cb0ba4 100644 --- a/festim/exports/surface_quantity.py +++ b/festim/exports/surface_quantity.py @@ -15,6 +15,8 @@ class SurfaceQuantity: field (festim.Species): species for which the surface flux is computed surface (festim.SurfaceSubdomain1D): surface subdomain filename (str): name of the file to which the surface flux is exported + t (list): list of time values + data (list): list of values of the surface quantity """ def __init__(self, field, surface, filename: str = None) -> None: @@ -64,6 +66,9 @@ def field(self, value): self._field = value def write(self, t): + """If the filename doesnt exist yet, create it and write the header, + then append the time and value to the file""" + if not os.path.isfile(self.filename): title = "Flux surface {}: {}".format(self.surface.id, self.field.name)