From 36e2ea28739604c222d8fd1bee8647bff0434137 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Fri, 12 Dec 2025 14:31:02 -0500 Subject: [PATCH 1/3] fix flux conservation for surface reaction BC --- src/festim/hydrogen_transport_problem.py | 22 +++--- test/test_h_transport_problem.py | 90 ++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 9 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index b971d39bb..c163830be 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1147,6 +1147,7 @@ def initialise(self): self.define_temperature() self.convert_source_input_values_to_fenics_objects() self.convert_advection_term_to_fenics_objects() + self.define_boundary_conditions() self.create_flux_values_fenics() self.create_initial_conditions() @@ -1154,7 +1155,6 @@ def initialise(self): self.create_subdomain_formulation(subdomain) subdomain.u.name = f"u_{subdomain.id}" - self.define_boundary_conditions() self.create_formulation() self.create_solver() self.initialise_exports() @@ -1408,7 +1408,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): ) # add fluxes - for bc in self.boundary_conditions: + for bc in self._unpacked_bcs: if isinstance(bc, boundary_conditions.ParticleFluxBC): # check that the bc is applied on a surface # belonging to this subdomain @@ -1551,13 +1551,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 @@ -1631,7 +1635,7 @@ def create_solver(self): def create_flux_values_fenics(self): """For each particle flux create the ``value_fenics`` attribute""" - for bc in self.boundary_conditions: + for bc in self._unpacked_bcs: if isinstance(bc, boundary_conditions.ParticleFluxBC): volume_subdomain = self.surface_to_volume[bc.subdomain] bc.create_value_fenics( @@ -1773,9 +1777,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( @@ -1908,7 +1912,7 @@ def create_formulation(self): ) # add fluxes - for bc in self.boundary_conditions: + for bc in self._unpacked_bcs: if isinstance(bc, boundary_conditions.ParticleFluxBC): self.formulation -= ( bc.value_fenics diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 393ccb511..73f67730f 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1373,3 +1373,93 @@ def test_traps_with_CG_elements(): # TEST assert my_model.function_space.sub(0).element.basix_element.discontinuous is False + + +def test_surface_reaction_BC_discontinuous(): + """Test that surface reaction BC conserves flux""" + + # TODO: clean this up + # BUILD + model_probe = F.HydrogenTransportProblemDiscontinuous() + + probe_thick = 1 + breeder_thick = 1 + + probe = F.Material( + D_0=1, + E_D=0, + K_S_0=1, + E_K_S=0, + ) + breeder = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0) + + volume_breeder = F.VolumeSubdomain1D( + id=1, borders=[0, breeder_thick], material=breeder + ) + volume_probe = F.VolumeSubdomain1D( + id=2, + borders=[breeder_thick, breeder_thick + probe_thick], + material=probe, + ) + + vertices_left = np.linspace(0, breeder_thick, num=1000, endpoint=False) + vertices_right = np.linspace(breeder_thick, breeder_thick + probe_thick, num=1000) + vertices = np.concatenate([vertices_left, vertices_right]) + model_probe.mesh = F.Mesh1D(vertices) + + inlet = F.SurfaceSubdomain1D(id=3, x=0) + outlet = F.SurfaceSubdomain1D(id=4, x=breeder_thick + probe_thick) + + model_probe.subdomains = [ + volume_breeder, + volume_probe, + inlet, + outlet, + ] + + model_probe.surface_to_volume = { + inlet: volume_breeder, + outlet: volume_probe, + } + + model_probe.method_interface = F.InterfaceMethod.nitsche + model_probe.interfaces = [ + F.Interface(id=5, subdomains=[volume_breeder, volume_probe], penalty_term=10), + ] + + H = F.Species("H", subdomains=model_probe.volume_subdomains) + model_probe.species = [H] + + model_probe.temperature = 600 + + model_probe.boundary_conditions = [ + F.FixedConcentrationBC( + subdomain=inlet, + value=1, + species=H, + ), + F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=0, + k_r0=1, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=outlet, + ), + ] + permeation_flux = F.SurfaceFlux(field=H, surface=outlet) + inlet_flux = F.SurfaceFlux(field=H, surface=inlet) + model_probe.exports = [permeation_flux, inlet_flux] + + model_probe.settings = F.Settings( + atol=1e-10, + rtol=1e-09, + transient=False, + ) + + model_probe.initialise() + model_probe.run() + + # TEST + assert np.isclose(permeation_flux.data[-1], -inlet_flux.data[-1], rtol=1e-5) From 6612b9c59cdde1c006d0a7ba1ad3265a8ce08368 Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Sat, 13 Dec 2025 15:36:45 -0500 Subject: [PATCH 2/3] 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 c163830be..0cf83440f 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1551,17 +1551,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 @@ -1777,9 +1773,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 63d0311eab5d7684b2bc1b63399c2a00caa87c4b Mon Sep 17 00:00:00 2001 From: kaelyndunnell Date: Mon, 15 Dec 2025 10:06:07 -0500 Subject: [PATCH 3/3] clean up test --- test/test_h_transport_problem.py | 62 +++++++++++++++----------------- 1 file changed, 28 insertions(+), 34 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 73f67730f..95de032d7 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1378,61 +1378,55 @@ def test_traps_with_CG_elements(): def test_surface_reaction_BC_discontinuous(): """Test that surface reaction BC conserves flux""" - # TODO: clean this up # BUILD - model_probe = F.HydrogenTransportProblemDiscontinuous() + my_model = F.HydrogenTransportProblemDiscontinuous() - probe_thick = 1 - breeder_thick = 1 - - probe = F.Material( + right_mat = F.Material( D_0=1, E_D=0, K_S_0=1, E_K_S=0, ) - breeder = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0) + left_mat = F.Material(D_0=1, E_D=0, K_S_0=2, E_K_S=0) - volume_breeder = F.VolumeSubdomain1D( - id=1, borders=[0, breeder_thick], material=breeder - ) - volume_probe = F.VolumeSubdomain1D( + left_vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=left_mat) + right_vol = F.VolumeSubdomain1D( id=2, - borders=[breeder_thick, breeder_thick + probe_thick], - material=probe, + borders=[1, 2], + material=right_mat, ) - vertices_left = np.linspace(0, breeder_thick, num=1000, endpoint=False) - vertices_right = np.linspace(breeder_thick, breeder_thick + probe_thick, num=1000) + vertices_left = np.linspace(0, 1, num=1000, endpoint=False) + vertices_right = np.linspace(1, 2, num=1000) vertices = np.concatenate([vertices_left, vertices_right]) - model_probe.mesh = F.Mesh1D(vertices) + my_model.mesh = F.Mesh1D(vertices) inlet = F.SurfaceSubdomain1D(id=3, x=0) - outlet = F.SurfaceSubdomain1D(id=4, x=breeder_thick + probe_thick) + outlet = F.SurfaceSubdomain1D(id=4, x=2) - model_probe.subdomains = [ - volume_breeder, - volume_probe, + my_model.subdomains = [ + left_vol, + right_vol, inlet, outlet, ] - model_probe.surface_to_volume = { - inlet: volume_breeder, - outlet: volume_probe, + my_model.surface_to_volume = { + inlet: left_vol, + outlet: right_vol, } - model_probe.method_interface = F.InterfaceMethod.nitsche - model_probe.interfaces = [ - F.Interface(id=5, subdomains=[volume_breeder, volume_probe], penalty_term=10), + my_model.method_interface = F.InterfaceMethod.nitsche + my_model.interfaces = [ + F.Interface(id=5, subdomains=[left_vol, right_vol], penalty_term=10), ] - H = F.Species("H", subdomains=model_probe.volume_subdomains) - model_probe.species = [H] + H = F.Species("H", subdomains=my_model.volume_subdomains) + my_model.species = [H] - model_probe.temperature = 600 + my_model.temperature = 600 - model_probe.boundary_conditions = [ + my_model.boundary_conditions = [ F.FixedConcentrationBC( subdomain=inlet, value=1, @@ -1450,16 +1444,16 @@ def test_surface_reaction_BC_discontinuous(): ] permeation_flux = F.SurfaceFlux(field=H, surface=outlet) inlet_flux = F.SurfaceFlux(field=H, surface=inlet) - model_probe.exports = [permeation_flux, inlet_flux] + my_model.exports = [permeation_flux, inlet_flux] - model_probe.settings = F.Settings( + my_model.settings = F.Settings( atol=1e-10, rtol=1e-09, transient=False, ) - model_probe.initialise() - model_probe.run() + my_model.initialise() + my_model.run() # TEST assert np.isclose(permeation_flux.data[-1], -inlet_flux.data[-1], rtol=1e-5)