From f3b304bd5a40f582667dc784956c6247bee86607 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 27 Oct 2025 12:13:11 -0400 Subject: [PATCH 1/5] allow D to be func for multi vol subdomains --- src/festim/hydrogen_transport_problem.py | 62 +++++++++++++++++++++--- test/test_h_transport_problem.py | 15 ++---- 2 files changed, 57 insertions(+), 20 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 133b9fd78..ab2c2fc09 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -569,6 +569,34 @@ def define_D_global(self, species): D.interpolate(D_expr) return D, D_expr + def define_D_mult_subdomains( + self, species, volume, mesh=None, temperature=None, field=None + ): + """Defines the global diffusion coefficient for a given species + + Args: + species (F.Species): the species + volume (F.VolumeSubdomain): the volume + + 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, _species.Species) + + # if diffusion coeffient has been given as a function, use that + global_D = {} + if volume.material.D: + global_D[volume.id] = volume.material.D + else: + global_D[volume.id] = ( + volume.material.volume.material.get_diffusion_coefficient( + self.mesh, self.temperature, self.field + ) + ) + return global_D + def define_function_spaces(self, element_degree=1): """Creates the function space of the model, creates a mixed element if model is multispecies. Creates the main solution and previous solution @@ -1597,13 +1625,17 @@ def mixed_term(u, v, n): interface.id ) - 0.5 * mixed_term( v_b, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface(interface.id) + ) * dInterface( + interface.id + ) F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface( interface.id ) - 0.5 * mixed_term( v_t, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface(interface.id) + ) * dInterface( + interface.id + ) F_0 += ( 2 * gamma @@ -1706,12 +1738,26 @@ def initialise_exports(self): # compute diffusivity function for surface fluxes # for the discontinuous case, we don't use D_global as in # HydrogenTransportProblem + + vol_to_D_global = {} + for export in self.exports: if isinstance(export, exports.SurfaceQuantity): volume = self.surface_to_volume[export.surface] - D = volume.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics, export.field - ) + if volume.id in vol_to_D_global: + D = vol_to_D_global[volume.id] + else: + D = self.define_D_mult_subdomains( + species=export.field, + volume=volume, + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + ) + vol_to_D_global[volume.id] = D + + # D = volume.material.get_diffusion_coefficient( + # self.mesh.mesh, self.temperature_fenics, export.field + # ) # NOTE: maybe we need to make sure there are no functionspace clashes? export.D = D @@ -1819,9 +1865,9 @@ def post_processing(self): export.write(t=float(self.t)) elif isinstance(export, exports.Profile1DExport): - assert export.subdomain, ( - "Profile1DExport requires a subdomain to be set" - ) + assert ( + export.subdomain + ), "Profile1DExport requires a subdomain to be set" u = export.subdomain.u if export._dofs is None: index = self.subdomain_to_species[export.subdomain].index( diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 6a92ef348..567ac8a84 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1253,8 +1253,8 @@ def test_create_flux_values_fenics_multispecies(): def test_not_implemented_error_raised_with_D_as_function(): - """Test that NotImplementedError is raised when D is a function and more than one - volume subdomain has been defined""" + """Test that if a function is given as diffusion coeff of a material for multiple subdomain case, + that it is passed to D_global_mult""" # BUILD test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101)) @@ -1275,16 +1275,7 @@ def test_not_implemented_error_raised_with_D_as_function(): my_model.define_function_spaces() my_model.assign_functions_to_species() - - # TEST - with pytest.raises( - NotImplementedError, - match=( - "Giving the diffusion coefficient as a function is currently only " - "supported for a single volume subdomain case" - ), - ): - my_model.define_D_global(H) + D = my_model.define_D_mult_subdomains(H, vol_1) def test_define_D_global_with_D_as_function(): From 2f5da14d79df87a8c4f3e4ed0525e3fea1b7a0eb Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 27 Oct 2025 12:18:25 -0400 Subject: [PATCH 2/5] test correctly --- src/festim/hydrogen_transport_problem.py | 21 +++++++++------------ test/test_h_transport_problem.py | 21 ++++++++++++++------- 2 files changed, 23 insertions(+), 19 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ab2c2fc09..044f9b2b6 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -546,11 +546,11 @@ def define_D_global(self, species): # if diffusion coeffient has been given as a function, use that if self.volume_subdomains[0].material.D: - if len(self.volume_subdomains) > 1: - raise NotImplementedError( - "Giving the diffusion coefficient as a function is currently " - "only supported for a single volume subdomain case" - ) + # if len(self.volume_subdomains) > 1: + # raise NotImplementedError( + # "Giving the diffusion coefficient as a function is currently " + # "only supported for a single volume subdomain case" + # ) return self.volume_subdomains[0].material.D, None D_0 = fem.Function(self.V_DG_0) @@ -586,16 +586,13 @@ def define_D_mult_subdomains( assert isinstance(species, _species.Species) # if diffusion coeffient has been given as a function, use that - global_D = {} if volume.material.D: - global_D[volume.id] = volume.material.D + D = volume.material.D else: - global_D[volume.id] = ( - volume.material.volume.material.get_diffusion_coefficient( - self.mesh, self.temperature, self.field - ) + D = volume.material.volume.material.get_diffusion_coefficient( + self.mesh, self.temperature, self.field ) - return global_D + return D def define_function_spaces(self, element_degree=1): """Creates the function space of the model, creates a mixed element if diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 567ac8a84..05d961fc1 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1259,12 +1259,15 @@ def test_not_implemented_error_raised_with_D_as_function(): # BUILD test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101)) V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) - D_func = fem.Function(V) - D_func.x.array[:] = 7.0 - - my_mat = F.Material(D=D_func) - vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=my_mat) - vol_2 = F.VolumeSubdomain1D(id=1, borders=[0.5, 1], material=my_mat) + D_func1 = fem.Function(V) + D_func1.x.array[:] = 7.0 + D_func2 = fem.Function(V) + D_func2.x.array[:] = 1.0 + + my_mat1 = F.Material(D=D_func1) + my_mat2 = F.Material(D=D_func2) + vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=my_mat1) + vol_2 = F.VolumeSubdomain1D(id=1, borders=[0.5, 1], material=my_mat2) H = F.Species("H") my_model = F.HydrogenTransportProblem( mesh=test_mesh, @@ -1275,7 +1278,11 @@ def test_not_implemented_error_raised_with_D_as_function(): my_model.define_function_spaces() my_model.assign_functions_to_species() - D = my_model.define_D_mult_subdomains(H, vol_1) + D1 = my_model.define_D_mult_subdomains(H, vol_1) + D2 = my_model.define_D_mult_subdomains(H, vol_2) + + assert D1 == D_func1 + assert D2 == D_func2 def test_define_D_global_with_D_as_function(): From 68a6a903141c47b64bbc509740156eaf451eb00a Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 27 Oct 2025 12:19:12 -0400 Subject: [PATCH 3/5] ruff --- src/festim/hydrogen_transport_problem.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 044f9b2b6..056ae3c72 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1622,17 +1622,13 @@ def mixed_term(u, v, n): interface.id ) - 0.5 * mixed_term( v_b, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface( - interface.id - ) + ) * dInterface(interface.id) F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface( interface.id ) - 0.5 * mixed_term( v_t, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface( - interface.id - ) + ) * dInterface(interface.id) F_0 += ( 2 * gamma @@ -1862,9 +1858,9 @@ def post_processing(self): export.write(t=float(self.t)) elif isinstance(export, exports.Profile1DExport): - assert ( - export.subdomain - ), "Profile1DExport requires a subdomain to be set" + assert export.subdomain, ( + "Profile1DExport requires a subdomain to be set" + ) u = export.subdomain.u if export._dofs is None: index = self.subdomain_to_species[export.subdomain].index( From ae4db79489246790cec29d005e8d0c760aaede40 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 27 Oct 2025 12:21:26 -0400 Subject: [PATCH 4/5] fix func args --- src/festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 056ae3c72..df006ea0b 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -590,7 +590,7 @@ def define_D_mult_subdomains( D = volume.material.D else: D = volume.material.volume.material.get_diffusion_coefficient( - self.mesh, self.temperature, self.field + mesh, temperature, field ) return D From 9a95372c99caf54f63ff0e8a9786bd66052051a5 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 27 Oct 2025 12:37:47 -0400 Subject: [PATCH 5/5] fix & test D_mult for non func cases --- src/festim/hydrogen_transport_problem.py | 24 +++++++---- test/test_h_transport_problem.py | 52 ++++++++++++++++++++++-- 2 files changed, 64 insertions(+), 12 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index df006ea0b..f4c21c413 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -570,7 +570,11 @@ def define_D_global(self, species): return D, D_expr def define_D_mult_subdomains( - self, species, volume, mesh=None, temperature=None, field=None + self, + species, + volume, + mesh=None, + temperature=None, ): """Defines the global diffusion coefficient for a given species @@ -589,9 +593,7 @@ def define_D_mult_subdomains( if volume.material.D: D = volume.material.D else: - D = volume.material.volume.material.get_diffusion_coefficient( - mesh, temperature, field - ) + D = volume.material.get_diffusion_coefficient(mesh, temperature, species) return D def define_function_spaces(self, element_degree=1): @@ -1622,13 +1624,17 @@ def mixed_term(u, v, n): interface.id ) - 0.5 * mixed_term( v_b, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface(interface.id) + ) * dInterface( + interface.id + ) F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface( interface.id ) - 0.5 * mixed_term( v_t, (u_b / K_b - u_t / K_t), n_0 - ) * dInterface(interface.id) + ) * dInterface( + interface.id + ) F_0 += ( 2 * gamma @@ -1858,9 +1864,9 @@ def post_processing(self): export.write(t=float(self.t)) elif isinstance(export, exports.Profile1DExport): - assert export.subdomain, ( - "Profile1DExport requires a subdomain to be set" - ) + assert ( + export.subdomain + ), "Profile1DExport requires a subdomain to be set" u = export.subdomain.u if export._dofs is None: index = self.subdomain_to_species[export.subdomain].index( diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 05d961fc1..3f1678cdb 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -10,6 +10,8 @@ import festim as F +from festim import k_B + test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) dummy_mat = F.Material(D_0=1, E_D=1, name="dummy_mat") @@ -1252,7 +1254,7 @@ def test_create_flux_values_fenics_multispecies(): assert np.isclose(float(my_model.boundary_conditions[1].value_fenics), 11) -def test_not_implemented_error_raised_with_D_as_function(): +def test_D_as_function_multiple_subdomains(): """Test that if a function is given as diffusion coeff of a material for multiple subdomain case, that it is passed to D_global_mult""" @@ -1278,13 +1280,57 @@ def test_not_implemented_error_raised_with_D_as_function(): my_model.define_function_spaces() my_model.assign_functions_to_species() - D1 = my_model.define_D_mult_subdomains(H, vol_1) - D2 = my_model.define_D_mult_subdomains(H, vol_2) + D1 = my_model.define_D_mult_subdomains( + species=H, + volume=vol_1, + ) + D2 = my_model.define_D_mult_subdomains( + species=H, + volume=vol_2, + ) assert D1 == D_func1 assert D2 == D_func2 +def test_D_not_function_multiple_subdomains(): + """Test that if a float or int is given as diffusion coeff of a material for multiple subdomain case, + that it is passed to D_global_mult and defined correctly""" + + # BUILD + test_mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=101)) + V = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) + D_01 = 5 + D_02 = 7 + + my_mat1 = F.Material(D_0=D_01, E_D=1) + my_mat2 = F.Material(D_0=D_02, E_D=1) + vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=my_mat1) + vol_2 = F.VolumeSubdomain1D(id=1, borders=[0.5, 1], material=my_mat2) + H = F.Species("H") + my_model = F.HydrogenTransportProblem( + mesh=test_mesh, + temperature=10, + subdomains=[vol_1, vol_2], + species=[F.Species("H")], + ) + + my_model.define_function_spaces() + my_model.assign_functions_to_species() + D1 = my_model.define_D_mult_subdomains( + species=H, volume=vol_1, mesh=test_mesh.mesh, temperature=my_model.temperature + ) + D2 = my_model.define_D_mult_subdomains( + species=H, volume=vol_2, mesh=test_mesh.mesh, temperature=my_model.temperature + ) + + D_analytical_1 = D_01 * np.exp(-1 / F.k_B / my_model.temperature) + D_analytical_2 = D_02 * np.exp(-1 / F.k_B / my_model.temperature) + + assert np.isclose(float(D1), D_analytical_1) + assert np.isclose(float(D2), D_analytical_2) + + def test_define_D_global_with_D_as_function(): """Test that if a function is given as diffusion coeff of a material, that it is passed to D_global"""