From 96c523d6bf85971c5e9c167545956d3648abd136 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 12 Feb 2026 14:53:09 -0500 Subject: [PATCH 1/4] interpolate into prev sol if function --- src/festim/hydrogen_transport_problem.py | 52 +++++++++++++++--------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 0cf83440f..f9c9a3281 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -37,6 +37,7 @@ as_fenics_constant, get_interpolation_points, is_it_time_to_export, + nmm_interpolate, ) from festim.material import SolubilityLaw @@ -1221,20 +1222,27 @@ def create_initial_conditions(self): the previous solution of the condition's species""" for condition in self.initial_conditions: - V = condition.species.subdomain_to_function_space[condition.volume] + idx = self.species.index(condition.species) - condition.create_expr_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, - function_space=V, - ) + # if the value given is a function, then directly interpolate it on the + # previous solution of the species + if isinstance(condition.value, fem.Function): + nmm_interpolate(condition.volume.u_n.sub(idx), condition.value) - # assign to previous solution of species - entities = self.volume_meshtags.find(condition.volume.id) - idx = self.species.index(condition.species) - condition.volume.u_n.sub(idx).interpolate( - condition.expr_fenics, cells1=entities - ) + else: + V = condition.species.subdomain_to_function_space[condition.volume] + + condition.create_expr_fenics( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + function_space=V, + ) + + # assign to previous solution of species + entities = self.volume_meshtags.find(condition.volume.id) + condition.volume.u_n.sub(idx).interpolate( + condition.expr_fenics, cells1=entities + ) def define_function_spaces( self, subdomain: _subdomain.VolumeSubdomain, element_degree=1 @@ -1652,9 +1660,9 @@ def initialise_exports(self): engine="BP5", ) else: - raise NotImplementedError( - f"Export type {type(export)} not implemented for " - f"mixed-domain approach" + adios4dolfinx.write_mesh( + export.filename, + mesh=functions[0].function_space.mesh, ) # compute diffusivity function for surface fluxes @@ -1704,11 +1712,17 @@ def post_processing(self): ) if isinstance(export, exports.VTXSpeciesExport): if export._checkpoint: - raise NotImplementedError( - f"Export type {type(export)} not implemented " - f"for mixed-domain approach" + post_processing_solution = export.field[ + 0 + ].subdomain_to_post_processing_solution[export._subdomain] + adios4dolfinx.write_function( + export.filename, + post_processing_solution, + time=float(self.t), + name=export.field[0].name, ) - export.writer.write(float(self.t)) + else: + export.writer.write(float(self.t)) # handle derived quantities if isinstance(export, exports.SurfaceQuantity): From 81fa63353b4a503d727df7e5a89cace97942e69a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 12 Feb 2026 16:13:42 -0500 Subject: [PATCH 2/4] allows several fields per vtx checkpoint --- src/festim/hydrogen_transport_problem.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index f9c9a3281..39b6f1350 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1712,15 +1712,18 @@ def post_processing(self): ) if isinstance(export, exports.VTXSpeciesExport): if export._checkpoint: - post_processing_solution = export.field[ - 0 - ].subdomain_to_post_processing_solution[export._subdomain] - adios4dolfinx.write_function( - export.filename, - post_processing_solution, - time=float(self.t), - name=export.field[0].name, - ) + for species in export.field: + post_processing_solution = ( + species.subdomain_to_post_processing_solution[ + export._subdomain + ] + ) + adios4dolfinx.write_function( + export.filename, + post_processing_solution, + time=float(self.t), + name=species.name, + ) else: export.writer.write(float(self.t)) From bfa39db225f9948495a34bdeea968fb228bf1813 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 17 Feb 2026 09:48:12 -0500 Subject: [PATCH 3/4] additional tests --- test/test_initial_condition.py | 101 +++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 116bf0359..c8154f12a 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -434,3 +434,104 @@ def test_initial_condition_continuous_multimaterial(): # assert value of spe1 is not all the same across the domain assert not np.allclose(my_model.u_n.x.array[:], intial_cond_value) + + +def test_initial_condition_mixed_domain(): + """Test the initial condition in a multi-material discontinous case""" + + my_model = F.HydrogenTransportProblemDiscontinuous() + + vertices_left = np.linspace(0, 0.5, 500) + vertices_right = np.linspace(0.5, 1, 500) + vertices = np.concatenate((vertices_left, vertices_right)) + my_model.mesh = F.Mesh1D(vertices) + + # assumes the same diffusivity for all species + material_left = F.Material(D_0=1e-01, E_D=0, K_S_0=1, E_K_S=0) + material_right = F.Material(D_0=1e-01, E_D=0, K_S_0=2, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=material_left) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=material_right) + my_model.subdomains = [vol1, vol2] + + spe1 = F.Species("1", mobile=True, subdomains=[vol1, vol2]) + my_model.species = [spe1] + + my_model.temperature = 300 + + intial_cond_value = 100 + V = fem.functionspace(my_model.mesh.mesh, ("CG", 1)) + u = fem.Function(V) + u.x.array[:] = intial_cond_value + + my_model.initial_conditions = [ + F.InitialConcentration(value=u, species=spe1, volume=vol1), + ] + + dt = F.Stepsize(0.1) + my_model.settings = F.Settings( + atol=1e-10, rtol=1e-10, final_time=5, transient=True, stepsize=dt + ) + + my_model.initialise() + + # assert value of spe1 in vol1 + left_spe1_un = vol1.u.sub(0) + + assert not np.allclose(left_spe1_un.x.array[:], intial_cond_value) + + +def test_initial_condition_mixed_domain_multispecies(): + """Test the initial condition in a multispecies multi-material discontinous case""" + + my_model = F.HydrogenTransportProblemDiscontinuous() + + vertices_left = np.linspace(0, 0.5, 500) + vertices_right = np.linspace(0.5, 1, 500) + vertices = np.concatenate((vertices_left, vertices_right)) + my_model.mesh = F.Mesh1D(vertices) + + # assumes the same diffusivity for all species + material_left = F.Material(D_0=1e-01, E_D=0, K_S_0=1, E_K_S=0) + material_right = F.Material(D_0=1e-01, E_D=0, K_S_0=2, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=material_left) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=material_right) + my_model.subdomains = [vol1, vol2] + + spe1 = F.Species("1", mobile=True, subdomains=[vol1, vol2]) + spe2 = F.Species("2", mobile=True, subdomains=[vol1, vol2]) + my_model.species = [spe1, spe2] + + my_model.temperature = 300 + + intial_cond_value_1 = 100 + intial_cond_value_2 = 200 + V = fem.functionspace(my_model.mesh.mesh, ("CG", 1)) + u = fem.Function(V) + u.x.array[:] = intial_cond_value_1 + + V2 = fem.functionspace(my_model.mesh.mesh, ("CG", 1)) + u2 = fem.Function(V2) + u2.x.array[:] = intial_cond_value_2 + + my_model.initial_conditions = [ + F.InitialConcentration(value=u, species=spe1, volume=vol1), + F.InitialConcentration(value=u2, species=spe2, volume=vol1), + ] + + dt = F.Stepsize(0.1) + my_model.settings = F.Settings( + atol=1e-10, rtol=1e-10, final_time=5, transient=True, stepsize=dt + ) + + my_model.initialise() + + spe1_left, spe1_to_vol1 = vol1.u_n.function_space.sub(0).collapse() + spe2_left, spe2_to_vol1 = vol1.u_n.function_space.sub(1).collapse() + + prev_solution_spe1_left = vol1.u_n.x.array[spe1_to_vol1] + prev_solution_spe2_left = vol1.u_n.x.array[spe2_to_vol1] + + assert np.allclose(prev_solution_spe1_left, intial_cond_value_1) + assert np.allclose(prev_solution_spe2_left, intial_cond_value_2) From df45ecf7fe0dd09af52c5f97262f42e417955465 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 17 Feb 2026 09:55:05 -0500 Subject: [PATCH 4/4] format ruff --- test/system_tests/test_multi_mat_penalty.py | 10 ++++------ test/system_tests/test_multi_material.py | 20 ++++++++----------- .../test_multi_material_change_of_var.py | 20 ++++++++----------- 3 files changed, 20 insertions(+), 30 deletions(-) diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py index efd505081..116da435e 100644 --- a/test/system_tests/test_multi_mat_penalty.py +++ b/test/system_tests/test_multi_mat_penalty.py @@ -71,13 +71,11 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = ( - lambda x: 3 - + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) - + 0.1 * ufl.cos(2 * ufl.pi * x[0]) + c_exact_top_ufl = lambda x: ( + 3 + ufl.sin(ufl.pi * (2 * x[1] + 0.5)) + 0.1 * ufl.cos(2 * ufl.pi * x[0]) ) - c_exact_top_np = ( - lambda x: 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) + c_exact_top_np = lambda x: ( + 3 + np.sin(np.pi * (2 * x[1] + 0.5)) + 0.1 * np.cos(2 * np.pi * x[0]) ) def c_exact_bot_ufl(x): diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index 167efbf11..86526d7b4 100644 --- a/test/system_tests/test_multi_material.py +++ b/test/system_tests/test_multi_material.py @@ -71,15 +71,15 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = ( - lambda x: 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) + c_exact_top_ufl = lambda x: ( + 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) ) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) - c_exact_top_np = ( - lambda x: 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) + c_exact_top_np = lambda x: ( + 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) ) def c_exact_bot_np(x): @@ -126,15 +126,11 @@ def c_exact_bot_np(x): F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H), ] - source_top_val = ( - lambda x: 8 - * ufl.pi**2 - * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + source_top_val = lambda x: ( + 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) - source_bottom_val = ( - lambda x: 40 - * ufl.pi**2 - * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + source_bottom_val = lambda x: ( + 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val), diff --git a/test/system_tests/test_multi_material_change_of_var.py b/test/system_tests/test_multi_material_change_of_var.py index fa0674461..8091a2663 100644 --- a/test/system_tests/test_multi_material_change_of_var.py +++ b/test/system_tests/test_multi_material_change_of_var.py @@ -156,15 +156,15 @@ def test_2_materials_2d_mms(tmpdir): K_S_bot = 6.0 D_top = 2.0 D_bot = 5.0 - c_exact_top_ufl = ( - lambda x: 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) + c_exact_top_ufl = lambda x: ( + 1 + ufl.sin(ufl.pi * (2 * x[0] + 0.5)) + ufl.cos(2 * ufl.pi * x[1]) ) def c_exact_bot_ufl(x): return K_S_bot / K_S_top * c_exact_top_ufl(x) - c_exact_top_np = ( - lambda x: 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) + c_exact_top_np = lambda x: ( + 1 + np.sin(np.pi * (2 * x[0] + 0.5)) + np.cos(2 * np.pi * x[1]) ) def c_exact_bot_np(x): @@ -201,15 +201,11 @@ def c_exact_bot_np(x): F.FixedConcentrationBC(bottom_surface, value=c_exact_bot_ufl, species=H), ] - source_top_val = ( - lambda x: 8 - * ufl.pi**2 - * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + source_top_val = lambda x: ( + 8 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) - source_bottom_val = ( - lambda x: 40 - * ufl.pi**2 - * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) + source_bottom_val = lambda x: ( + 40 * ufl.pi**2 * (ufl.cos(2 * ufl.pi * x[0]) + ufl.cos(2 * ufl.pi * x[1])) ) my_model.sources = [ F.ParticleSource(volume=top_domain, species=H, value=source_top_val),