From fe8051d43980c73255294f34e1fcefa33ac09ffa Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 15:36:37 -0400 Subject: [PATCH 01/41] new material class --- festim/material.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 festim/material.py diff --git a/festim/material.py b/festim/material.py new file mode 100644 index 000000000..9aaa3f56e --- /dev/null +++ b/festim/material.py @@ -0,0 +1,35 @@ +import ufl +from dolfinx import fem +import festim as F + + +class Material: + """ + Material class + """ + + def __init__(self, D_0=None, E_D=None, name=None) -> None: + """Inits Material + Args: + D_0 (float or fem.Constant): the diffusion coefficient at 0 K + E_D (float or fem.Constant): the activation energy for diffusion + name (str): the name of the material + """ + self.D_0 = D_0 + self.E_D = E_D + self.name = name + + def define_diffusion_coefficient(self, mesh, temperature): + """Defines the diffusion coefficient + Args: + temperature (dolfinx.fem.Constant): the temperature + Returns: + float: the diffusion coefficient + """ + + if isinstance(self.D_0, float): + self.D_0 = fem.Constant(mesh, self.D_0) + if isinstance(self.E_D, float): + self.E_D = fem.Constant(mesh, self.E_D) + + return self.D_0 * ufl.exp(-self.E_D / F.k_B / temperature) From 29446dd1f0ab54a3c89a5b86b4804afb8455831a Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 15:37:10 -0400 Subject: [PATCH 02/41] pull material properties from volum subdomain --- festim/hydrogen_transport_problem.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index e5dfd56b2..036cba7c4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -97,6 +97,17 @@ def initialise(self): self.define_function_space() self.define_markers_and_measures() self.assign_functions_to_species() + + self.volume_subdomains = [] + for sub_dom in self.subdomains: + if isinstance(sub_dom, F.VolumeSubdomain1D): + self.volume_subdomains.append(sub_dom) + self.D = sub_dom.material.define_diffusion_coefficient( + self.mesh.mesh, self.temperature + ) + if len(self.volume_subdomains) > 1: + raise NotImplementedError("Multiple volume subdomains not implemented yet") + self.create_formulation() def define_function_space(self): @@ -159,16 +170,9 @@ def create_formulation(self): if len(self.species) > 1: raise NotImplementedError("Multiple species not implemented yet") - # TODO expose D_0 and E_D as parameters of a Material class - D_0 = fem.Constant(self.mesh.mesh, 1.9e-7) - E_D = fem.Constant(self.mesh.mesh, 0.2) - - D = D_0 * exp(-E_D / F.k_B / self.temperature) - # TODO expose dt as parameter of the model dt = fem.Constant(self.mesh.mesh, 1 / 20) - self.D = D # TODO remove this self.dt = dt # TODO remove this for spe in self.species: @@ -176,7 +180,7 @@ def create_formulation(self): u_n = spe.prev_solution v = spe.test_function - formulation = dot(D * grad(u), grad(v)) * self.dx + formulation = dot(self.D * grad(u), grad(v)) * self.dx formulation += ((u - u_n) / dt) * v * self.dx # add sources From ff7e7aabf4d7132d2b872ef55b9e7d6a38694b28 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 15:37:18 -0400 Subject: [PATCH 03/41] update test --- test/test_permeation_problem.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 87492dec4..44ae480d2 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -8,8 +8,7 @@ form, assemble_scalar, ) -from ufl import dot, grad, exp, FacetNormal, ds, Measure -from dolfinx import log +from ufl import dot, grad, exp, FacetNormal import numpy as np import tqdm.autonotebook @@ -26,7 +25,8 @@ def test_permeation_problem(): my_model = F.HydrogenTransportProblem() my_model.mesh = my_mesh - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=None) + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) left_surface = F.SurfaceSubdomain1D(id=1, x=0) right_surface = F.SurfaceSubdomain1D(id=2, x=L) my_model.subdomains = [my_subdomain, left_surface, right_surface] From 4dd4c50bf846d13d8b9ae0517b34cfd210e7b3af Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 15:37:24 -0400 Subject: [PATCH 04/41] update init --- festim/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/festim/__init__.py b/festim/__init__.py index de4a2652c..c871ab014 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -16,12 +16,14 @@ R = 8.314462618 # Gas constant J.mol-1.K-1 k_B = 8.6173303e-5 # Boltzmann constant eV.K-1 +from .material import Material + from .mesh.mesh import Mesh from .mesh.mesh_1d import Mesh1D +from .hydrogen_transport_problem import HydrogenTransportProblem + from .species import Species, Trap from .subdomain.surface_subdomain import SurfaceSubdomain1D from .subdomain.volume_subdomain import VolumeSubdomain1D - -from .hydrogen_transport_problem import HydrogenTransportProblem From 17aa043e94be3cc019749a879aa24f4b2eb26b20 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 15:39:38 -0400 Subject: [PATCH 05/41] process temperature --- festim/hydrogen_transport_problem.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 036cba7c4..976d64022 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -18,7 +18,7 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float or dolfinx.Function): the temperature of the model + temperature (float or fem.Constant): the temperature of the model sources (list of festim.Source): the hydrogen sources of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model @@ -28,7 +28,7 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float or dolfinx.Function): the temperature of the model + temperature (float or fem.Constant): the temperature of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model exports (list of festim.Export): the exports of the model @@ -98,6 +98,17 @@ def initialise(self): self.define_markers_and_measures() self.assign_functions_to_species() + if isinstance(self.temperature, float): + self.temperature = fem.Constant(self.mesh.mesh, self.temperature) + elif isinstance(self.temperature, fem.Constant): + pass + elif self.temperature is None: + raise ValueError("Temperature not defined") + else: + raise TypeError( + f"Temperature must be float or dolfinx.Constant, not {type(self.temperature)}" + ) + self.volume_subdomains = [] for sub_dom in self.subdomains: if isinstance(sub_dom, F.VolumeSubdomain1D): From 4c5e208af8a391d379685b8ea434efa3c2d256ed Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:09:53 -0400 Subject: [PATCH 06/41] tests for temperature --- test/test_temperature.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 test/test_temperature.py diff --git a/test/test_temperature.py b/test/test_temperature.py new file mode 100644 index 000000000..b445d4958 --- /dev/null +++ b/test/test_temperature.py @@ -0,0 +1,38 @@ +import festim as F +from dolfinx import fem +import numpy as np + + +def test_temperature_value(): + # Test that the temperature value is correctly set + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.array([0, 1, 3, 4])) + my_model.species = [F.Species("H")] + my_mat = F.Material(1, 1, "1") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) + my_model.subdomains = [my_subdomain] + + my_model.temperature = fem.Constant(my_model.mesh.mesh, 23.0) + my_model.initialise() + + assert float(my_model.temperature) == 23.0 + + +def test_temperature_type(): + # Test that the temperature type is correctly set + my_mesh = F.Mesh1D(vertices=np.array([0, 1, 3, 4])) + my_mesh.generate_mesh() + values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] + + for value in values: + my_model = F.HydrogenTransportProblem() + my_model.mesh = my_mesh + my_model.species = [F.Species("H")] + my_mat = F.Material(1, 1, "1") + my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) + my_model.subdomains = [my_subdomain] + + my_model.temperature = value + my_model.initialise() + + assert isinstance(my_model.temperature, fem.Constant) From 70be3a9d44b4084f2a3f920b0de46efa010e5417 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:24:57 -0400 Subject: [PATCH 07/41] test temperature types --- test/test_temperature.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index b445d4958..7a9b8df2b 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -1,6 +1,8 @@ import festim as F from dolfinx import fem import numpy as np +import sympy as sp +from pytest import raises def test_temperature_value(): @@ -24,15 +26,22 @@ def test_temperature_type(): my_mesh.generate_mesh() values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] - for value in values: + def model(value): my_model = F.HydrogenTransportProblem() my_model.mesh = my_mesh my_model.species = [F.Species("H")] my_mat = F.Material(1, 1, "1") my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) my_model.subdomains = [my_subdomain] - my_model.temperature = value my_model.initialise() + return my_model + + for value in values: + my_model = model(value) assert isinstance(my_model.temperature, fem.Constant) + + with raises(TypeError): + x = sp.Symbol("x") + model(sp.sin(sp.pi * x)) From d9b41e693c88a9cb57b3a06013c8b8c8756e0317 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:31:01 -0400 Subject: [PATCH 08/41] use ufl not sympy --- test/test_temperature.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index 7a9b8df2b..c78cea49f 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -1,8 +1,8 @@ import festim as F from dolfinx import fem import numpy as np -import sympy as sp from pytest import raises +import ufl def test_temperature_value(): @@ -43,5 +43,5 @@ def model(value): assert isinstance(my_model.temperature, fem.Constant) with raises(TypeError): - x = sp.Symbol("x") - model(sp.sin(sp.pi * x)) + x = ufl.SpatialCoordinate(my_mesh.mesh) + model(2 * x[0]) From d38963999799557a1a1fe96e330d2f38bac6a75f Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:38:27 -0400 Subject: [PATCH 09/41] test materials class --- test/test_material.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 test/test_material.py diff --git a/test/test_material.py b/test/test_material.py new file mode 100644 index 000000000..74edc6efa --- /dev/null +++ b/test/test_material.py @@ -0,0 +1,17 @@ +import festim as F +import numpy as np + + +def test_define_diffusion_coefficient(): + """Test that the diffusion coefficient is correctly defined""" + T = 10 + D_0 = 1.2 + E_D = 0.5 + + my_mesh = F.Mesh1D(vertices=np.array([0, 1, 2])) + my_mat = F.Material(D_0=D_0, E_D=E_D) + D = my_mat.define_diffusion_coefficient(my_mesh.mesh, T) + + D_analytical = D_0 * np.exp(-E_D / F.k_B / T) + + assert float(D) == D_analytical From 0539c2e090bb9c48407a0e23d8864ac5faec3f76 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:38:43 -0400 Subject: [PATCH 10/41] allow integers --- festim/hydrogen_transport_problem.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 976d64022..8f9f9d342 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -98,8 +98,8 @@ def initialise(self): self.define_markers_and_measures() self.assign_functions_to_species() - if isinstance(self.temperature, float): - self.temperature = fem.Constant(self.mesh.mesh, self.temperature) + if isinstance(self.temperature, (float, int)): + self.temperature = fem.Constant(self.mesh.mesh, float(self.temperature)) elif isinstance(self.temperature, fem.Constant): pass elif self.temperature is None: From f5032934f043a65202b500242d79938e422451f4 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:38:59 -0400 Subject: [PATCH 11/41] diffusion values required --- festim/material.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/material.py b/festim/material.py index 9aaa3f56e..af3b2dcb9 100644 --- a/festim/material.py +++ b/festim/material.py @@ -8,7 +8,7 @@ class Material: Material class """ - def __init__(self, D_0=None, E_D=None, name=None) -> None: + def __init__(self, D_0, E_D, name=None) -> None: """Inits Material Args: D_0 (float or fem.Constant): the diffusion coefficient at 0 K From 27156916443ebc7b310b590aca71e4cef7527058 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:44:50 -0400 Subject: [PATCH 12/41] doc strings and floats --- test/test_temperature.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index c78cea49f..20a635ea4 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -6,9 +6,9 @@ def test_temperature_value(): - # Test that the temperature value is correctly set + """Test that the temperature value is correctly set""" my_model = F.HydrogenTransportProblem() - my_model.mesh = F.Mesh1D(vertices=np.array([0, 1, 3, 4])) + my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) my_model.species = [F.Species("H")] my_mat = F.Material(1, 1, "1") my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) @@ -21,8 +21,8 @@ def test_temperature_value(): def test_temperature_type(): - # Test that the temperature type is correctly set - my_mesh = F.Mesh1D(vertices=np.array([0, 1, 3, 4])) + """Test that the temperature type is correctly set""" + my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) my_mesh.generate_mesh() values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] From 7a9600d0c5c7a53c70f4de46c241f34e96e7c73a Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:44:58 -0400 Subject: [PATCH 13/41] use floats --- test/test_material.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_material.py b/test/test_material.py index 74edc6efa..48167937e 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -8,7 +8,7 @@ def test_define_diffusion_coefficient(): D_0 = 1.2 E_D = 0.5 - my_mesh = F.Mesh1D(vertices=np.array([0, 1, 2])) + my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0])) my_mat = F.Material(D_0=D_0, E_D=E_D) D = my_mat.define_diffusion_coefficient(my_mesh.mesh, T) From d14ca16e76f512d1028e6933d994132610a295a8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:45:14 -0400 Subject: [PATCH 14/41] mini test for type of input --- festim/material.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/festim/material.py b/festim/material.py index af3b2dcb9..6d2c06bda 100644 --- a/festim/material.py +++ b/festim/material.py @@ -27,9 +27,17 @@ def define_diffusion_coefficient(self, mesh, temperature): float: the diffusion coefficient """ - if isinstance(self.D_0, float): - self.D_0 = fem.Constant(mesh, self.D_0) - if isinstance(self.E_D, float): - self.E_D = fem.Constant(mesh, self.E_D) + if isinstance(self.D_0, (float, int)): + self.D_0 = fem.Constant(mesh, float(self.D_0)) + else: + raise TypeError( + f"D_0 must be float, int or dolfinx.fem.Constant, not {type(self.D_0)}" + ) + if isinstance(self.E_D, (float, int)): + self.E_D = fem.Constant(mesh, float(self.E_D)) + else: + raise TypeError( + f"E_D must be float, int or dolfinx.fem.Constant, not {type(self.E_D)}" + ) return self.D_0 * ufl.exp(-self.E_D / F.k_B / temperature) From 7e8ac31009e8138709494579be89e68f390f0854 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 11 Oct 2023 16:53:11 -0400 Subject: [PATCH 15/41] test material type --- festim/material.py | 4 ++++ test/test_material.py | 24 ++++++++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/festim/material.py b/festim/material.py index 6d2c06bda..8a674edd5 100644 --- a/festim/material.py +++ b/festim/material.py @@ -29,12 +29,16 @@ def define_diffusion_coefficient(self, mesh, temperature): if isinstance(self.D_0, (float, int)): self.D_0 = fem.Constant(mesh, float(self.D_0)) + elif isinstance(self.D_0, fem.Constant): + pass else: raise TypeError( f"D_0 must be float, int or dolfinx.fem.Constant, not {type(self.D_0)}" ) if isinstance(self.E_D, (float, int)): self.E_D = fem.Constant(mesh, float(self.E_D)) + elif isinstance(self.E_D, fem.Constant): + pass else: raise TypeError( f"E_D must be float, int or dolfinx.fem.Constant, not {type(self.E_D)}" diff --git a/test/test_material.py b/test/test_material.py index 48167937e..78fdfaaf7 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -1,5 +1,8 @@ import festim as F import numpy as np +from dolfinx import fem +from pytest import raises +import ufl def test_define_diffusion_coefficient(): @@ -15,3 +18,24 @@ def test_define_diffusion_coefficient(): D_analytical = D_0 * np.exp(-E_D / F.k_B / T) assert float(D) == D_analytical + + +def test_diffusion_values_type(): + """Test that the diffusion coefficient values types are correctly + processed""" + my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) + my_mesh.generate_mesh() + values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] + + def model(value): + my_mat = F.Material(value, value, "1") + my_mat.define_diffusion_coefficient( + my_mesh.mesh, fem.Constant(my_mesh.mesh, 1.0) + ) + + for value in values: + model(value) + + with raises(TypeError): + x = ufl.SpatialCoordinate(my_mesh.mesh) + model(2 * x[0]) From 6ffb0def685591d25c2f2c4429ecabaad19da47f Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 12:15:04 -0400 Subject: [PATCH 16/41] use setters for temperature --- festim/hydrogen_transport_problem.py | 28 ++++++++------ test/test_temperature.py | 57 ++++++++++------------------ 2 files changed, 37 insertions(+), 48 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 8f9f9d342..6c1075e3d 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -89,6 +89,23 @@ def __init__( self.volume_meshtags = None self.formulation = None + @property + def temperature(self): + return self._temperature + + @temperature.setter + def temperature(self, value): + if value is None: + self._temperature = value + elif isinstance(value, (float, int)): + self._temperature = fem.Constant(self.mesh.mesh, float(value)) + elif isinstance(value, fem.Constant): + self._temperature = value + else: + raise TypeError( + f"Temperature must be float or dolfinx.Constant, not {type(value)}" + ) + def initialise(self): """Initialise the model. Creates suitable function spaces, facet and volume tags... @@ -98,17 +115,6 @@ def initialise(self): self.define_markers_and_measures() self.assign_functions_to_species() - if isinstance(self.temperature, (float, int)): - self.temperature = fem.Constant(self.mesh.mesh, float(self.temperature)) - elif isinstance(self.temperature, fem.Constant): - pass - elif self.temperature is None: - raise ValueError("Temperature not defined") - else: - raise TypeError( - f"Temperature must be float or dolfinx.Constant, not {type(self.temperature)}" - ) - self.volume_subdomains = [] for sub_dom in self.subdomains: if isinstance(sub_dom, F.VolumeSubdomain1D): diff --git a/test/test_temperature.py b/test/test_temperature.py index 20a635ea4..caa76f012 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -1,47 +1,30 @@ import festim as F from dolfinx import fem import numpy as np -from pytest import raises -import ufl +import pytest - -def test_temperature_value(): - """Test that the temperature value is correctly set""" - my_model = F.HydrogenTransportProblem() - my_model.mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) - my_model.species = [F.Species("H")] - my_mat = F.Material(1, 1, "1") - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) - my_model.subdomains = [my_subdomain] - - my_model.temperature = fem.Constant(my_model.mesh.mesh, 23.0) - my_model.initialise() - - assert float(my_model.temperature) == 23.0 +test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) -def test_temperature_type(): +@pytest.mark.parametrize( + "test_type", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou"] +) +def test_temperature_type(test_type): """Test that the temperature type is correctly set""" - my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) - my_mesh.generate_mesh() - values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] - - def model(value): - my_model = F.HydrogenTransportProblem() - my_model.mesh = my_mesh - my_model.species = [F.Species("H")] - my_mat = F.Material(1, 1, "1") - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=my_mat) - my_model.subdomains = [my_subdomain] - my_model.temperature = value - my_model.initialise() + my_model = F.HydrogenTransportProblem() + if not isinstance(test_type, (fem.Constant, int, float)): + with pytest.raises(TypeError): + my_model.temperature = test_type - return my_model - for value in values: - my_model = model(value) - assert isinstance(my_model.temperature, fem.Constant) +@pytest.mark.parametrize( + "test_value", + [1, fem.Constant(test_mesh.mesh, 1.0), 1.0], +) +def test_temperature_value_processing(test_value): + """Test that the temperature type is correctly processed""" + my_model = F.HydrogenTransportProblem() + my_model.mesh = test_mesh + my_model.temperature = test_value - with raises(TypeError): - x = ufl.SpatialCoordinate(my_mesh.mesh) - model(2 * x[0]) + assert isinstance(my_model.temperature, fem.Constant) From 408460115eec8dbc52671b6ca0447cf26b9b7dd4 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 12:39:38 -0400 Subject: [PATCH 17/41] use parametrised tests --- test/test_material.py | 39 ++++++++++++++++++--------------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index 78fdfaaf7..67f5ed515 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -1,41 +1,38 @@ import festim as F import numpy as np from dolfinx import fem -from pytest import raises +import pytest import ufl +test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) +x = ufl.SpatialCoordinate(test_mesh.mesh) + def test_define_diffusion_coefficient(): """Test that the diffusion coefficient is correctly defined""" - T = 10 - D_0 = 1.2 - E_D = 0.5 + T, D_0, E_D = 10, 1.2, 0.5 - my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0])) my_mat = F.Material(D_0=D_0, E_D=E_D) - D = my_mat.define_diffusion_coefficient(my_mesh.mesh, T) + D = my_mat.define_diffusion_coefficient(test_mesh.mesh, T) D_analytical = D_0 * np.exp(-E_D / F.k_B / T) assert float(D) == D_analytical -def test_diffusion_values_type(): +@pytest.mark.parametrize( + "test_type", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]] +) +def test_diffusion_values_type(test_type): """Test that the diffusion coefficient values types are correctly processed""" - my_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) - my_mesh.generate_mesh() - values = [int(1), fem.Constant(my_mesh.mesh, 1.0), float(1.0)] - - def model(value): - my_mat = F.Material(value, value, "1") + my_mat = F.Material(test_type, test_type, "1") + if isinstance(test_type, (fem.Constant, int, float)): my_mat.define_diffusion_coefficient( - my_mesh.mesh, fem.Constant(my_mesh.mesh, 1.0) + test_mesh.mesh, fem.Constant(test_mesh.mesh, 1.0) ) - - for value in values: - model(value) - - with raises(TypeError): - x = ufl.SpatialCoordinate(my_mesh.mesh) - model(2 * x[0]) + else: + with pytest.raises(TypeError): + my_mat.define_diffusion_coefficient( + test_mesh.mesh, fem.Constant(test_mesh.mesh, 1.0) + ) From 4ac713212ef1498c266eefc9744ceb79791cbb89 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 12:47:48 -0400 Subject: [PATCH 18/41] added comments --- festim/material.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/festim/material.py b/festim/material.py index 8a674edd5..aee057b6e 100644 --- a/festim/material.py +++ b/festim/material.py @@ -22,11 +22,13 @@ def __init__(self, D_0, E_D, name=None) -> None: def define_diffusion_coefficient(self, mesh, temperature): """Defines the diffusion coefficient Args: + mesh (dolfinx.mesh.Mesh): the domain mesh temperature (dolfinx.fem.Constant): the temperature Returns: - float: the diffusion coefficient + ufl.algebra.Product: the diffusion coefficient """ + # check type of values and convert to fem.Constant if needed if isinstance(self.D_0, (float, int)): self.D_0 = fem.Constant(mesh, float(self.D_0)) elif isinstance(self.D_0, fem.Constant): From 600c582c38b9b258bc417257b45523d4ce0d62d3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:23:55 -0400 Subject: [PATCH 19/41] use helper, new func name --- festim/material.py | 25 +++++-------------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/festim/material.py b/festim/material.py index aee057b6e..0bd03a3b1 100644 --- a/festim/material.py +++ b/festim/material.py @@ -1,5 +1,4 @@ import ufl -from dolfinx import fem import festim as F @@ -19,7 +18,7 @@ def __init__(self, D_0, E_D, name=None) -> None: self.E_D = E_D self.name = name - def define_diffusion_coefficient(self, mesh, temperature): + def get_diffusion_coefficient(self, mesh, temperature): """Defines the diffusion coefficient Args: mesh (dolfinx.mesh.Mesh): the domain mesh @@ -28,22 +27,8 @@ def define_diffusion_coefficient(self, mesh, temperature): ufl.algebra.Product: the diffusion coefficient """ - # check type of values and convert to fem.Constant if needed - if isinstance(self.D_0, (float, int)): - self.D_0 = fem.Constant(mesh, float(self.D_0)) - elif isinstance(self.D_0, fem.Constant): - pass - else: - raise TypeError( - f"D_0 must be float, int or dolfinx.fem.Constant, not {type(self.D_0)}" - ) - if isinstance(self.E_D, (float, int)): - self.E_D = fem.Constant(mesh, float(self.E_D)) - elif isinstance(self.E_D, fem.Constant): - pass - else: - raise TypeError( - f"E_D must be float, int or dolfinx.fem.Constant, not {type(self.E_D)}" - ) + # check type of values + D_0 = F.as_fenics_constant(self.D_0, mesh) + E_D = F.as_fenics_constant(self.E_D, mesh) - return self.D_0 * ufl.exp(-self.E_D / F.k_B / temperature) + return D_0 * ufl.exp(-E_D / F.k_B / temperature) From 395d8fd233c92ccfa9ece371c80e3df41a34edaa Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:24:09 -0400 Subject: [PATCH 20/41] removed test, new func name --- test/test_material.py | 22 +--------------------- 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index 67f5ed515..6de2865c7 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -1,7 +1,5 @@ import festim as F import numpy as np -from dolfinx import fem -import pytest import ufl test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) @@ -13,26 +11,8 @@ def test_define_diffusion_coefficient(): T, D_0, E_D = 10, 1.2, 0.5 my_mat = F.Material(D_0=D_0, E_D=E_D) - D = my_mat.define_diffusion_coefficient(test_mesh.mesh, T) + D = my_mat.get_diffusion_coefficient(test_mesh.mesh, T) D_analytical = D_0 * np.exp(-E_D / F.k_B / T) assert float(D) == D_analytical - - -@pytest.mark.parametrize( - "test_type", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]] -) -def test_diffusion_values_type(test_type): - """Test that the diffusion coefficient values types are correctly - processed""" - my_mat = F.Material(test_type, test_type, "1") - if isinstance(test_type, (fem.Constant, int, float)): - my_mat.define_diffusion_coefficient( - test_mesh.mesh, fem.Constant(test_mesh.mesh, 1.0) - ) - else: - with pytest.raises(TypeError): - my_mat.define_diffusion_coefficient( - test_mesh.mesh, fem.Constant(test_mesh.mesh, 1.0) - ) From 9577c8f819f01231715bc52231861ea1aa118ed1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:27:23 -0400 Subject: [PATCH 21/41] use helpers, D not attribute of sim --- festim/hydrogen_transport_problem.py | 55 ++++++++++++++-------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 6c1075e3d..3851c75c4 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -28,7 +28,7 @@ class HydrogenTransportProblem: mesh (festim.Mesh): the mesh of the model subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model - temperature (float or fem.Constant): the temperature of the model + temperature (fem.Constant): the temperature of the model boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model solver_parameters (dict): the solver parameters of the model exports (list of festim.Export): the exports of the model @@ -97,10 +97,8 @@ def temperature(self): def temperature(self, value): if value is None: self._temperature = value - elif isinstance(value, (float, int)): - self._temperature = fem.Constant(self.mesh.mesh, float(value)) - elif isinstance(value, fem.Constant): - self._temperature = value + elif isinstance(value, (float, int, fem.Constant)): + self._temperature = F.as_fenics_constant(value, self.mesh.mesh) else: raise TypeError( f"Temperature must be float or dolfinx.Constant, not {type(value)}" @@ -119,11 +117,6 @@ def initialise(self): for sub_dom in self.subdomains: if isinstance(sub_dom, F.VolumeSubdomain1D): self.volume_subdomains.append(sub_dom) - self.D = sub_dom.material.define_diffusion_coefficient( - self.mesh.mesh, self.temperature - ) - if len(self.volume_subdomains) > 1: - raise NotImplementedError("Multiple volume subdomains not implemented yet") self.create_formulation() @@ -193,24 +186,30 @@ def create_formulation(self): self.dt = dt # TODO remove this for spe in self.species: - u = spe.solution - u_n = spe.prev_solution - v = spe.test_function - - formulation = dot(self.D * grad(u), grad(v)) * self.dx - formulation += ((u - u_n) / dt) * v * self.dx - - # add sources - for source in self.sources: - # f = Constant(my_mesh.mesh, (PETSc.ScalarType(0))) - if source.species == spe: - formulation += source * v * self.dx - # add fluxes - # TODO implement this - # for bc in self.boundary_conditions: - # pass - # if bc.species == spe and bc.type != "dirichlet": - # formulation += bc * v * self.ds + for vol in self.volume_subdomains: + u = spe.solution + u_n = spe.prev_solution + v = spe.test_function + + D = vol.material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature + ) + + formulation = dot(D * grad(u), grad(v)) * self.dx + formulation += ((u - u_n) / dt) * v * self.dx + + # add sources + # TODO implement this + # for source in self.sources: + # # f = Constant(my_mesh.mesh, (PETSc.ScalarType(0))) + # if source.species == spe: + # formulation += source * v * self.dx + # add fluxes + # TODO implement this + # for bc in self.boundary_conditions: + # pass + # if bc.species == spe and bc.type != "dirichlet": + # formulation += bc * v * self.ds self.formulation = formulation From a0dd726ab1437c76a8e61e0804725b133de8712e Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:27:34 -0400 Subject: [PATCH 22/41] helpers --- festim/helpers.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 festim/helpers.py diff --git a/festim/helpers.py b/festim/helpers.py new file mode 100644 index 000000000..cb622d4f7 --- /dev/null +++ b/festim/helpers.py @@ -0,0 +1,13 @@ +import festim as F +from dolfinx import fem + + +def as_fenics_constant(value, mesh): + if isinstance(value, (float, int)): + return fem.Constant(mesh, float(value)) + elif isinstance(value, fem.Constant): + return value + else: + raise TypeError( + f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" + ) From 7b60d4724634d5b641a841cfd38f16e9e112584a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:27:43 -0400 Subject: [PATCH 23/41] merge tests --- test/test_temperature.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/test/test_temperature.py b/test/test_temperature.py index caa76f012..827a6960c 100644 --- a/test/test_temperature.py +++ b/test/test_temperature.py @@ -2,29 +2,23 @@ from dolfinx import fem import numpy as np import pytest +import ufl test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) +x = ufl.SpatialCoordinate(test_mesh.mesh) @pytest.mark.parametrize( - "test_type", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou"] + "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]] ) -def test_temperature_type(test_type): +def test_temperature_type_and_processing(value): """Test that the temperature type is correctly set""" my_model = F.HydrogenTransportProblem() - if not isinstance(test_type, (fem.Constant, int, float)): - with pytest.raises(TypeError): - my_model.temperature = test_type - - -@pytest.mark.parametrize( - "test_value", - [1, fem.Constant(test_mesh.mesh, 1.0), 1.0], -) -def test_temperature_value_processing(test_value): - """Test that the temperature type is correctly processed""" - my_model = F.HydrogenTransportProblem() my_model.mesh = test_mesh - my_model.temperature = test_value - assert isinstance(my_model.temperature, fem.Constant) + if not isinstance(value, (fem.Constant, int, float)): + with pytest.raises(TypeError): + my_model.temperature = value + else: + my_model.temperature = value + assert isinstance(my_model.temperature, fem.Constant) From 7f3be9a1c724923a73a7eb203aa0f3ac50789aa8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:27:53 -0400 Subject: [PATCH 24/41] D not attribute of model --- test/test_permeation_problem.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 44ae480d2..92269b6b0 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -39,6 +39,8 @@ def test_permeation_problem(): my_model.initialise() + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) + V = my_model.function_space u = mobile_H.solution @@ -98,7 +100,7 @@ def siverts_law(T, S_0, E_S, pressure): mobile_xdmf.write_function(u, t) - surface_flux = form(my_model.D * dot(grad(u), n) * my_model.ds(2)) + surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) @@ -109,13 +111,12 @@ def siverts_law(T, S_0, E_S, pressure): # analytical solution S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability = float(my_model.D) * S + permeability = float(D) * S times = np.array(times) n_array = np.arange(1, 10000)[:, np.newaxis] summation = np.sum( - (-1) ** n_array - * np.exp(-((np.pi * n_array) ** 2) * float(my_model.D) / L**2 * times), + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), axis=0, ) analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) From 56c432198a02d9ff69e38f4e1c9ab096dba977ea Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:28:00 -0400 Subject: [PATCH 25/41] test helpers --- test/test_helpers.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 test/test_helpers.py diff --git a/test/test_helpers.py b/test/test_helpers.py new file mode 100644 index 000000000..14b269bda --- /dev/null +++ b/test/test_helpers.py @@ -0,0 +1,21 @@ +import festim as F +from dolfinx import fem +import numpy as np +import pytest +import ufl + +test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) +x = ufl.SpatialCoordinate(test_mesh.mesh) + + +@pytest.mark.parametrize( + "value", [1, fem.Constant(test_mesh.mesh, 1.0), 1.0, "coucou", 2 * x[0]] +) +def test_temperature_type_and_processing(value): + """Test that the temperature type is correctly set""" + + if not isinstance(value, (fem.Constant, int, float)): + with pytest.raises(TypeError): + F.as_fenics_constant(value, test_mesh.mesh) + else: + assert isinstance(F.as_fenics_constant(value, test_mesh.mesh), fem.Constant) From 6a31bd9f605f5a7466bb2cebaea29f98f95f21fc Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 14:28:07 -0400 Subject: [PATCH 26/41] import helpers --- festim/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/festim/__init__.py b/festim/__init__.py index c871ab014..0edcfb9ea 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -16,6 +16,8 @@ R = 8.314462618 # Gas constant J.mol-1.K-1 k_B = 8.6173303e-5 # Boltzmann constant eV.K-1 +from .helpers import as_fenics_constant + from .material import Material from .mesh.mesh import Mesh From 5a9f771e6a203bc152fd393a4b0f347250e25ee9 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 15:24:03 -0400 Subject: [PATCH 27/41] support and new test for multi-vol --- festim/hydrogen_transport_problem.py | 16 ++-- test/test_permeation_problem.py | 133 ++++++++++++++++++++++++++- 2 files changed, 138 insertions(+), 11 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3851c75c4..7de02dc2c 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -88,6 +88,7 @@ def __init__( self.facet_meshtags = None self.volume_meshtags = None self.formulation = None + self.volume_subdomains = [] @property def temperature(self): @@ -112,12 +113,6 @@ def initialise(self): self.define_function_space() self.define_markers_and_measures() self.assign_functions_to_species() - - self.volume_subdomains = [] - for sub_dom in self.subdomains: - if isinstance(sub_dom, F.VolumeSubdomain1D): - self.volume_subdomains.append(sub_dom) - self.create_formulation() def define_function_space(self): @@ -143,6 +138,7 @@ def define_markers_and_measures(self): tags_facets.append(sub_dom.id) if isinstance(sub_dom, F.VolumeSubdomain1D): # find all cells in subdomain and mark them as sub_dom.id + self.volume_subdomains.append(sub_dom) entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim) tags_volumes[entities] = sub_dom.id @@ -185,6 +181,8 @@ def create_formulation(self): self.dt = dt # TODO remove this + self.formulation = 0 + for spe in self.species: for vol in self.volume_subdomains: u = spe.solution @@ -195,8 +193,8 @@ def create_formulation(self): self.mesh.mesh, self.temperature ) - formulation = dot(D * grad(u), grad(v)) * self.dx - formulation += ((u - u_n) / dt) * v * self.dx + self.formulation += dot(D * grad(u), grad(v)) * self.dx(vol.id) + self.formulation += ((u - u_n) / dt) * v * self.dx(vol.id) # add sources # TODO implement this @@ -211,8 +209,6 @@ def create_formulation(self): # if bc.species == spe and bc.type != "dirichlet": # formulation += bc * v * self.ds - self.formulation = formulation - def create_solver(self): """Creates the solver of the model""" problem = fem.petsc.NonlinearProblem( diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 92269b6b0..64995a394 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -134,5 +134,136 @@ def siverts_law(T, S_0, E_S, pressure): assert error < 0.01 +def test_permeation_problem_multi_volume(): + L = 3e-04 + vertices = np.linspace(0, L, num=1001) + + my_mesh = F.Mesh1D(vertices) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = my_mesh + + my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") + my_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, L / 4], material=my_mat) + my_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[L / 4, L / 2], material=my_mat) + my_subdomain_3 = F.VolumeSubdomain1D( + id=3, borders=[L / 2, 3 * L / 4], material=my_mat + ) + my_subdomain_4 = F.VolumeSubdomain1D(id=4, borders=[3 * L / 4, L], material=my_mat) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=L) + my_model.subdomains = [ + my_subdomain_1, + my_subdomain_2, + my_subdomain_3, + my_subdomain_4, + left_surface, + right_surface, + ] + + mobile_H = F.Species("H") + my_model.species = [mobile_H] + + temperature = Constant(my_mesh.mesh, 500.0) + my_model.temperature = temperature + + my_model.initialise() + + D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) + + V = my_model.function_space + u = mobile_H.solution + + # TODO this should be a property of Mesh + n = FacetNormal(my_mesh.mesh) + + def siverts_law(T, S_0, E_S, pressure): + S = S_0 * exp(-E_S / F.k_B / T) + return S * pressure**0.5 + + fdim = my_mesh.mesh.topology.dim - 1 + left_facets = my_model.facet_meshtags.find(1) + left_dofs = locate_dofs_topological(V, fdim, left_facets) + right_facets = my_model.facet_meshtags.find(2) + right_dofs = locate_dofs_topological(V, fdim, right_facets) + + S_0 = 4.02e21 + E_S = 1.04 + P_up = 100 + surface_conc = siverts_law(T=temperature, S_0=S_0, E_S=E_S, pressure=P_up) + bc_sieverts = dirichletbc( + Constant(my_mesh.mesh, PETSc.ScalarType(surface_conc)), left_dofs, V + ) + bc_outgas = dirichletbc(Constant(my_mesh.mesh, PETSc.ScalarType(0)), right_dofs, V) + my_model.boundary_conditions = [bc_sieverts, bc_outgas] + my_model.create_solver() + + my_model.solver.convergence_criterion = "incremental" + my_model.solver.rtol = 1e-10 + my_model.solver.atol = 1e10 + + my_model.solver.report = True + ksp = my_model.solver.krylov_solver + opts = PETSc.Options() + option_prefix = ksp.getOptionsPrefix() + opts[f"{option_prefix}ksp_type"] = "cg" + opts[f"{option_prefix}pc_type"] = "gamg" + opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" + ksp.setFromOptions() + + mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") + mobile_xdmf.write_mesh(my_model.mesh.mesh) + + final_time = 50 + + flux_values = [] + times = [] + t = 0 + progress = tqdm.autonotebook.tqdm( + desc="Solving H transport problem", total=final_time + ) + while t < final_time: + progress.update(float(my_model.dt)) + t += float(my_model.dt) + + my_model.solver.solve(u) + + mobile_xdmf.write_function(u, t) + + surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(t) + + mobile_H.prev_solution.x.array[:] = u.x.array[:] + + mobile_xdmf.close() + + # analytical solution + S = S_0 * exp(-E_S / F.k_B / float(temperature)) + permeability = float(D) * S + times = np.array(times) + + n_array = np.arange(1, 10000)[:, np.newaxis] + summation = np.sum( + (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), + axis=0, + ) + analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) + + analytical_flux = np.abs(analytical_flux) + flux_values = np.array(np.abs(flux_values)) + + relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) + + relative_error = relative_error[ + np.where(analytical_flux > 0.01 * np.max(analytical_flux)) + ] + error = relative_error.mean() + + assert error < 0.01 + + if __name__ == "__main__": - test_permeation_problem() + # test_permeation_problem() + test_permeation_problem_multi_volume() From 9c2fcf32a2ca7f7d8b1fb8aa0530d9a2b60f6af8 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:26:15 -0400 Subject: [PATCH 28/41] added doc strings --- festim/helpers.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/festim/helpers.py b/festim/helpers.py index cb622d4f7..41d1bac8b 100644 --- a/festim/helpers.py +++ b/festim/helpers.py @@ -3,6 +3,18 @@ def as_fenics_constant(value, mesh): + """Converts a value to a dolfinx.Constant + + Args: + value (float, int or dolfinx.Constant): the value to convert + mesh (dolfinx.mesh.Mesh): the mesh of the domiain + + Returns: + dolfinx.Constant: the converted value + + Raises: + TypeError: if the value is not a float, an int or a dolfinx.Constant + """ if isinstance(value, (float, int)): return fem.Constant(mesh, float(value)) elif isinstance(value, fem.Constant): From bd6fc9286687e58531b8b33db2382929615c66ec Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:26:48 -0400 Subject: [PATCH 29/41] moved doc strings --- festim/material.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/festim/material.py b/festim/material.py index 0bd03a3b1..574bdcd0d 100644 --- a/festim/material.py +++ b/festim/material.py @@ -5,15 +5,19 @@ class Material: """ Material class + + Args: + D_0 (float or fem.Constant): the diffusion coefficient at 0 K + E_D (float or fem.Constant): the activation energy for diffusion + name (str): the name of the material + + Attributes: + D_0 (float or fem.Constant): the diffusion coefficient at 0 K + E_D (float or fem.Constant): the activation energy for diffusion + name (str): the name of the material """ def __init__(self, D_0, E_D, name=None) -> None: - """Inits Material - Args: - D_0 (float or fem.Constant): the diffusion coefficient at 0 K - E_D (float or fem.Constant): the activation energy for diffusion - name (str): the name of the material - """ self.D_0 = D_0 self.E_D = E_D self.name = name From ed612d3c24ec7de566ca90ccb7f7c96f2cb61f52 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:27:02 -0400 Subject: [PATCH 30/41] check borders within domain --- festim/subdomain/volume_subdomain.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py index a3398993c..c97d848b0 100644 --- a/festim/subdomain/volume_subdomain.py +++ b/festim/subdomain/volume_subdomain.py @@ -37,9 +37,20 @@ def locate_subdomain_entities(self, mesh, vdim): Returns: entities (np.array): the entities of the subdomain """ + self.check_borders_within_domain(mesh) + entities = locate_entities( - mesh, + mesh.mesh, vdim, lambda x: np.logical_and(x[0] >= self.borders[0], x[0] <= self.borders[1]), ) return entities + + def check_borders_within_domain(self, mesh): + """Checks that the borders of the subdomain are within the domain + + Returns: + bool: True if borders are within domain, False otherwise + """ + if self.borders[0] < mesh.vertices[0] or self.borders[1] > mesh.vertices[-1]: + raise ValueError("borders of subdomain are outside of domain") From b0b717eef5f1c24c066a5b9e5996b9c9b94ec60a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:27:21 -0400 Subject: [PATCH 31/41] remove lines not used, use np.close --- test/test_material.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/test_material.py b/test/test_material.py index 6de2865c7..e706cfe88 100644 --- a/test/test_material.py +++ b/test/test_material.py @@ -1,9 +1,7 @@ import festim as F import numpy as np -import ufl test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) -x = ufl.SpatialCoordinate(test_mesh.mesh) def test_define_diffusion_coefficient(): @@ -15,4 +13,4 @@ def test_define_diffusion_coefficient(): D_analytical = D_0 * np.exp(-E_D / F.k_B / T) - assert float(D) == D_analytical + assert np.isclose(float(D), D_analytical) From 7143a00ccb583e202b29258a51095acbeb5a2de6 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:27:36 -0400 Subject: [PATCH 32/41] remove unecessary lines --- test/test_permeation_problem.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 64995a394..b8534107e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -262,8 +262,3 @@ def siverts_law(T, S_0, E_S, pressure): error = relative_error.mean() assert error < 0.01 - - -if __name__ == "__main__": - # test_permeation_problem() - test_permeation_problem_multi_volume() From 0084de6d82666bd1b1126cf53b896a7a204607e1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 16:27:52 -0400 Subject: [PATCH 33/41] refactoring --- festim/hydrogen_transport_problem.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 7de02dc2c..304261454 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -98,12 +98,8 @@ def temperature(self): def temperature(self, value): if value is None: self._temperature = value - elif isinstance(value, (float, int, fem.Constant)): - self._temperature = F.as_fenics_constant(value, self.mesh.mesh) else: - raise TypeError( - f"Temperature must be float or dolfinx.Constant, not {type(value)}" - ) + self._temperature = F.as_fenics_constant(value, self.mesh.mesh) def initialise(self): """Initialise the model. Creates suitable function @@ -139,7 +135,7 @@ def define_markers_and_measures(self): if isinstance(sub_dom, F.VolumeSubdomain1D): # find all cells in subdomain and mark them as sub_dom.id self.volume_subdomains.append(sub_dom) - entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim) + entities = sub_dom.locate_subdomain_entities(self.mesh, vdim) tags_volumes[entities] = sub_dom.id # dofs and tags need to be in np.in32 format for meshtags @@ -184,11 +180,11 @@ def create_formulation(self): self.formulation = 0 for spe in self.species: - for vol in self.volume_subdomains: - u = spe.solution - u_n = spe.prev_solution - v = spe.test_function + u = spe.solution + u_n = spe.prev_solution + v = spe.test_function + for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( self.mesh.mesh, self.temperature ) From 7e3692f461cb77081ffd72fb6f1e79f11cda2e26 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:02:09 -0400 Subject: [PATCH 34/41] additional tests for volume subdomain borders --- test/test_subdomains.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/test/test_subdomains.py b/test/test_subdomains.py index c446804fc..89bb37113 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -1,5 +1,6 @@ import numpy as np import festim as F +import pytest def test_different_surface_ids(): @@ -40,9 +41,37 @@ def test_different_volume_ids(): vol_subdomain_3 = F.VolumeSubdomain1D( id=vol_subdomains_ids[2], borders=[sub_dom_3[0], sub_dom_3[-1]], material=None ) - my_model.subdomains = [vol_subdomain_1, vol_subdomain_2, vol_subdomain_3] + my_model.subdomains = [vol_subdomain_1, vol_subdomain_3, vol_subdomain_2] my_model.define_markers_and_measures() for vol_id in vol_subdomains_ids: assert vol_id in np.array(my_model.volume_meshtags.values) + + +def test_non_matching_volume_borders(): + my_model = F.HydrogenTransportProblem() + + sub_dom_1 = np.linspace(0, 1e-04, num=5) + sub_dom_2 = np.linspace(8e-05, 2e-04, num=5) + my_model.mesh = F.Mesh1D(np.unique(np.concatenate([sub_dom_1, sub_dom_2]))) + + vol_subdomains_ids = [1, 2] + vol_subdomain_1 = F.VolumeSubdomain1D( + id=vol_subdomains_ids[0], borders=[sub_dom_1[0], sub_dom_1[-1]], material=None + ) + vol_subdomain_2 = F.VolumeSubdomain1D( + id=vol_subdomains_ids[1], borders=[sub_dom_2[0], sub_dom_2[-1]], material=None + ) + my_model.subdomains = [vol_subdomain_1, vol_subdomain_2] + + with pytest.raises(ValueError): + my_model.define_markers_and_measures() + + +def test_borders_out_of_domain(): + test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=5)) + sub_dom = F.VolumeSubdomain1D(id=1, borders=[1, 1.5], material=None) + + with pytest.raises(ValueError): + sub_dom.locate_subdomain_entities(test_mesh, 1) From 28e4131f811bcddb083446c24764d1c3fe9476b0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:04:00 -0400 Subject: [PATCH 35/41] test borders of subdomains --- festim/hydrogen_transport_problem.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 304261454..5dc849bf5 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -3,8 +3,8 @@ import ufl from mpi4py import MPI from dolfinx.fem import Function -from dolfinx.mesh import meshtags, locate_entities -from ufl import TestFunction, dot, grad, exp, Measure +from dolfinx.mesh import meshtags +from ufl import TestFunction, dot, grad, Measure import numpy as np import festim as F @@ -138,6 +138,17 @@ def define_markers_and_measures(self): entities = sub_dom.locate_subdomain_entities(self.mesh, vdim) tags_volumes[entities] = sub_dom.id + # check that subdomains are connected + if len(self.volume_subdomains) > 1: + all_borders = [] + for vol in self.volume_subdomains: + for border in vol.borders: + all_borders.append(border) + sorted_borders = np.sort(all_borders).reshape(int(len(all_borders) / 2), 2) + for i in range(0, len(sorted_borders) - 1): + if sorted_borders[i][1] != sorted_borders[i + 1][0]: + raise ValueError("Subdomain borders don't match to each other ") + # dofs and tags need to be in np.in32 format for meshtags dofs_facets = np.array(dofs_facets, dtype=np.int32) tags_facets = np.array(tags_facets, dtype=np.int32) From 653af2488945a8556307c5a2a5d97f75f2d13b06 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:11:39 -0400 Subject: [PATCH 36/41] typo --- festim/subdomain/volume_subdomain.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py index c97d848b0..b073f3a0e 100644 --- a/festim/subdomain/volume_subdomain.py +++ b/festim/subdomain/volume_subdomain.py @@ -49,8 +49,8 @@ def locate_subdomain_entities(self, mesh, vdim): def check_borders_within_domain(self, mesh): """Checks that the borders of the subdomain are within the domain - Returns: - bool: True if borders are within domain, False otherwise + Raises: + Value error: if borders outside the domain """ if self.borders[0] < mesh.vertices[0] or self.borders[1] > mesh.vertices[-1]: raise ValueError("borders of subdomain are outside of domain") From b835bee275e8696a8ea763f33f7d30ae4af93b1a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:18:50 -0400 Subject: [PATCH 37/41] update doc strings --- festim/hydrogen_transport_problem.py | 21 ++++++++++++--------- festim/mesh/mesh.py | 7 +++---- festim/mesh/mesh_1d.py | 9 +++------ festim/species.py | 5 ----- 4 files changed, 18 insertions(+), 24 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 5dc849bf5..3384e422f 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -20,7 +20,8 @@ class HydrogenTransportProblem: species (list of festim.Species): the species of the model temperature (float or fem.Constant): the temperature of the model sources (list of festim.Source): the hydrogen sources of the model - boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model + boundary_conditions (list of festim.BoundaryCondition): the boundary + conditions of the model solver_parameters (dict): the solver parameters of the model exports (list of festim.Export): the exports of the model @@ -29,14 +30,17 @@ class HydrogenTransportProblem: subdomains (list of festim.Subdomain): the subdomains of the model species (list of festim.Species): the species of the model temperature (fem.Constant): the temperature of the model - boundary_conditions (list of festim.BoundaryCondition): the boundary conditions of the model + boundary_conditions (list of festim.BoundaryCondition): the boundary + conditions of the model solver_parameters (dict): the solver parameters of the model exports (list of festim.Export): the exports of the model dx (dolfinx.fem.dx): the volume measure of the model ds (dolfinx.fem.ds): the surface measure of the model - function_space (dolfinx.fem.FunctionSpace): the function space of the model + function_space (dolfinx.fem.FunctionSpace): the function space of the + model facet_meshtags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model - volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model + volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the + model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model @@ -102,10 +106,6 @@ def temperature(self, value): self._temperature = F.as_fenics_constant(value, self.mesh.mesh) def initialise(self): - """Initialise the model. Creates suitable function - spaces, facet and volume tags... - """ - self.define_function_space() self.define_markers_and_measures() self.assign_functions_to_species() @@ -116,6 +116,9 @@ def define_function_space(self): self.function_space = fem.FunctionSpace(self.mesh.mesh, elements) def define_markers_and_measures(self): + """Defines the markers and measures of the model, and check borders of + VolumeSubdomains match to each other""" + dofs_facets, tags_facets = [], [] # TODO this should be a property of mesh @@ -147,7 +150,7 @@ def define_markers_and_measures(self): sorted_borders = np.sort(all_borders).reshape(int(len(all_borders) / 2), 2) for i in range(0, len(sorted_borders) - 1): if sorted_borders[i][1] != sorted_borders[i + 1][0]: - raise ValueError("Subdomain borders don't match to each other ") + raise ValueError("Subdomain borders don't match to each other") # dofs and tags need to be in np.in32 format for meshtags dofs_facets = np.array(dofs_facets, dtype=np.int32) diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py index 90f7746bc..5f6be3dac 100644 --- a/festim/mesh/mesh.py +++ b/festim/mesh/mesh.py @@ -5,15 +5,14 @@ class Mesh: """ Mesh class + Args: + mesh (dolfinx.mesh.Mesh, optional): the mesh. Defaults to None. + Attributes: mesh (dolfinx.mesh.Mesh): the mesh """ def __init__(self, mesh=None): - """Inits Mesh - Args: - mesh (dolfinx.mesh.Mesh, optional): the mesh. Defaults to None. - """ self.mesh = mesh if self.mesh is not None: diff --git a/festim/mesh/mesh_1d.py b/festim/mesh/mesh_1d.py index 6a29e4c54..f7de96bee 100644 --- a/festim/mesh/mesh_1d.py +++ b/festim/mesh/mesh_1d.py @@ -9,17 +9,14 @@ class Mesh1D(Mesh): """ 1D Mesh + Args: + vertices (list): the mesh x-coordinates (m) + Attributes: vertices (list): the mesh x-coordinates (m) """ def __init__(self, vertices, **kwargs) -> None: - """Inits Mesh1D - - Args: - vertices (list): the mesh x-coordinates (m) - """ - self.vertices = vertices mesh = self.generate_mesh() diff --git a/festim/species.py b/festim/species.py index 18ed7683d..eda76f2e1 100644 --- a/festim/species.py +++ b/festim/species.py @@ -22,11 +22,6 @@ class Species: """ def __init__(self, name: str = None) -> None: - """_summary_ - - Args: - name (str, optional): a name given to the species. Defaults to None. - """ self.name = name self.solution = None From 3a545bee4bf33c476d4a4c3c3303541bec29fbfa Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:18:57 -0400 Subject: [PATCH 38/41] line not needed --- festim/helpers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/festim/helpers.py b/festim/helpers.py index 41d1bac8b..9b17fb92b 100644 --- a/festim/helpers.py +++ b/festim/helpers.py @@ -1,4 +1,3 @@ -import festim as F from dolfinx import fem From d9266f5191600ed773860ee97ac3386aba944109 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:22:05 -0400 Subject: [PATCH 39/41] updated doc strings --- festim/material.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/festim/material.py b/festim/material.py index 574bdcd0d..494181e99 100644 --- a/festim/material.py +++ b/festim/material.py @@ -7,13 +7,17 @@ class Material: Material class Args: - D_0 (float or fem.Constant): the diffusion coefficient at 0 K - E_D (float or fem.Constant): the activation energy for diffusion + D_0 (float or fem.Constant): the pre-exponential factor of the + diffusion coefficient (m2/s) + E_D (float or fem.Constant): the activation energy of the diffusion + coeficient (eV) name (str): the name of the material Attributes: - D_0 (float or fem.Constant): the diffusion coefficient at 0 K - E_D (float or fem.Constant): the activation energy for diffusion + D_0 (float or fem.Constant): the pre-exponential factor of the + diffusion coefficient (m2/s) + E_D (float or fem.Constant): the activation energy of the diffusion + coeficient (eV) name (str): the name of the material """ From 484264d3b63c0160f5036ebdeb1256405a7fb39b Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 12 Oct 2023 17:44:37 -0400 Subject: [PATCH 40/41] remove border checks for another PR --- festim/hydrogen_transport_problem.py | 14 +------------ festim/subdomain/volume_subdomain.py | 13 +----------- test/test_subdomains.py | 31 +--------------------------- 3 files changed, 3 insertions(+), 55 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3384e422f..df2abea22 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -116,8 +116,7 @@ def define_function_space(self): self.function_space = fem.FunctionSpace(self.mesh.mesh, elements) def define_markers_and_measures(self): - """Defines the markers and measures of the model, and check borders of - VolumeSubdomains match to each other""" + """Defines the markers and measures of the model""" dofs_facets, tags_facets = [], [] @@ -141,17 +140,6 @@ def define_markers_and_measures(self): entities = sub_dom.locate_subdomain_entities(self.mesh, vdim) tags_volumes[entities] = sub_dom.id - # check that subdomains are connected - if len(self.volume_subdomains) > 1: - all_borders = [] - for vol in self.volume_subdomains: - for border in vol.borders: - all_borders.append(border) - sorted_borders = np.sort(all_borders).reshape(int(len(all_borders) / 2), 2) - for i in range(0, len(sorted_borders) - 1): - if sorted_borders[i][1] != sorted_borders[i + 1][0]: - raise ValueError("Subdomain borders don't match to each other") - # dofs and tags need to be in np.in32 format for meshtags dofs_facets = np.array(dofs_facets, dtype=np.int32) tags_facets = np.array(tags_facets, dtype=np.int32) diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py index b073f3a0e..6fa4b1695 100644 --- a/festim/subdomain/volume_subdomain.py +++ b/festim/subdomain/volume_subdomain.py @@ -30,27 +30,16 @@ def locate_subdomain_entities(self, mesh, vdim): """Locates all cells in subdomain borders within domain Args: - mesh (dolfinx.cpp.mesh.Mesh): the mesh of the model + mesh (festim.Mesh): the mesh of the model vdim (int): the dimension of the volumes of the mesh, for 1D this is always 1 Returns: entities (np.array): the entities of the subdomain """ - self.check_borders_within_domain(mesh) - entities = locate_entities( mesh.mesh, vdim, lambda x: np.logical_and(x[0] >= self.borders[0], x[0] <= self.borders[1]), ) return entities - - def check_borders_within_domain(self, mesh): - """Checks that the borders of the subdomain are within the domain - - Raises: - Value error: if borders outside the domain - """ - if self.borders[0] < mesh.vertices[0] or self.borders[1] > mesh.vertices[-1]: - raise ValueError("borders of subdomain are outside of domain") diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 89bb37113..c446804fc 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -1,6 +1,5 @@ import numpy as np import festim as F -import pytest def test_different_surface_ids(): @@ -41,37 +40,9 @@ def test_different_volume_ids(): vol_subdomain_3 = F.VolumeSubdomain1D( id=vol_subdomains_ids[2], borders=[sub_dom_3[0], sub_dom_3[-1]], material=None ) - my_model.subdomains = [vol_subdomain_1, vol_subdomain_3, vol_subdomain_2] + my_model.subdomains = [vol_subdomain_1, vol_subdomain_2, vol_subdomain_3] my_model.define_markers_and_measures() for vol_id in vol_subdomains_ids: assert vol_id in np.array(my_model.volume_meshtags.values) - - -def test_non_matching_volume_borders(): - my_model = F.HydrogenTransportProblem() - - sub_dom_1 = np.linspace(0, 1e-04, num=5) - sub_dom_2 = np.linspace(8e-05, 2e-04, num=5) - my_model.mesh = F.Mesh1D(np.unique(np.concatenate([sub_dom_1, sub_dom_2]))) - - vol_subdomains_ids = [1, 2] - vol_subdomain_1 = F.VolumeSubdomain1D( - id=vol_subdomains_ids[0], borders=[sub_dom_1[0], sub_dom_1[-1]], material=None - ) - vol_subdomain_2 = F.VolumeSubdomain1D( - id=vol_subdomains_ids[1], borders=[sub_dom_2[0], sub_dom_2[-1]], material=None - ) - my_model.subdomains = [vol_subdomain_1, vol_subdomain_2] - - with pytest.raises(ValueError): - my_model.define_markers_and_measures() - - -def test_borders_out_of_domain(): - test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=5)) - sub_dom = F.VolumeSubdomain1D(id=1, borders=[1, 1.5], material=None) - - with pytest.raises(ValueError): - sub_dom.locate_subdomain_entities(test_mesh, 1) From ce48583664aed98029a0ad51a60512cd280cfae0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 13 Oct 2023 08:17:27 -0400 Subject: [PATCH 41/41] pass dolfin mesh to volume subdomain --- festim/hydrogen_transport_problem.py | 2 +- festim/subdomain/volume_subdomain.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index df2abea22..9898f4ba2 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -137,7 +137,7 @@ def define_markers_and_measures(self): if isinstance(sub_dom, F.VolumeSubdomain1D): # find all cells in subdomain and mark them as sub_dom.id self.volume_subdomains.append(sub_dom) - entities = sub_dom.locate_subdomain_entities(self.mesh, vdim) + entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim) tags_volumes[entities] = sub_dom.id # dofs and tags need to be in np.in32 format for meshtags diff --git a/festim/subdomain/volume_subdomain.py b/festim/subdomain/volume_subdomain.py index 6fa4b1695..dd3d1bf24 100644 --- a/festim/subdomain/volume_subdomain.py +++ b/festim/subdomain/volume_subdomain.py @@ -30,7 +30,7 @@ def locate_subdomain_entities(self, mesh, vdim): """Locates all cells in subdomain borders within domain Args: - mesh (festim.Mesh): the mesh of the model + mesh (dolfinx.mesh.Mesh): the mesh of the model vdim (int): the dimension of the volumes of the mesh, for 1D this is always 1 @@ -38,7 +38,7 @@ def locate_subdomain_entities(self, mesh, vdim): entities (np.array): the entities of the subdomain """ entities = locate_entities( - mesh.mesh, + mesh, vdim, lambda x: np.logical_and(x[0] >= self.borders[0], x[0] <= self.borders[1]), )