From 50b265c668631dd9f4d651b7558e8970c2432134 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 29 Jan 2026 10:50:25 -0500 Subject: [PATCH 1/8] move decoraters --- src/festim/hydrogen_transport_problem.py | 26 ++++++++++++------------ 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 0cf83440f..d4afbb271 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1090,6 +1090,19 @@ def __init__( self.surface_to_volume = surface_to_volume or {} self.subdomain_to_species = {} # maps subdomain to species defined in it + @property + def method_interface(self): + return self._method_interface + + @method_interface.setter + def method_interface(self, value): + if isinstance(value, InterfaceMethod): + self._method_interface = value + elif isinstance(value, str): + self._method_interface = InterfaceMethod.from_string(value) + else: + raise TypeError("method_interface must be of type str or InterfaceMethod") + def initialise(self): # check that all species have a list of F.VolumeSubdomain as this is # different from F.HydrogenTransportProblem @@ -1159,19 +1172,6 @@ def initialise(self): self.create_solver() self.initialise_exports() - @property - def method_interface(self): - return self._method_interface - - @method_interface.setter - def method_interface(self, value): - if isinstance(value, InterfaceMethod): - self._method_interface = value - elif isinstance(value, str): - self._method_interface = InterfaceMethod.from_string(value) - else: - raise TypeError("method_interface must be of type str or InterfaceMethod") - def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): """ Creates the ``value_fenics`` attribute for a given From aa2fe1d1e9d38ccd5aa2113c1a0bbf3d2631a33c Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 29 Jan 2026 14:22:27 -0500 Subject: [PATCH 2/8] working fix --- .../boundary_conditions/dirichlet_bc.py | 50 ++++--------------- src/festim/hydrogen_transport_problem.py | 30 ++++++++++- src/festim/subdomain/volume_subdomain.py | 18 ++++++- 3 files changed, 56 insertions(+), 42 deletions(-) diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index 1628fb899..f947ab98c 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -243,47 +243,19 @@ def create_value( if "T" in arguments: kwargs["T"] = temperature - try: - self.value_fenics = fem.Function(function_space) - - # store the expression of the boundary condition - # to update the value_fenics later - self.bc_expr = fem.Expression( - self.value(**kwargs), - helpers.get_interpolation_points(function_space.element), - ) - self.value_fenics.interpolate(self.bc_expr) - except RuntimeError: - # if this fails, it is probably because the temperature is a Function - # from the parent mesh and this is used in a mixed domain problem. - # In this case, we need to interpolate the temperature on the submesh - - submesh = mesh - mesh_cell_map = submesh.topology.index_map(submesh.topology.dim) - num_cells_on_proc = ( - mesh_cell_map.size_local + mesh_cell_map.num_ghosts - ) - cells = np.arange(num_cells_on_proc, dtype=np.int32) - parent_functionspace = temperature.function_space - interpolation_data = fem.create_interpolation_data( - function_space, parent_functionspace, cells - ) - - temperature_sub = fem.Function(function_space) - temperature_sub.interpolate_nonmatching( - temperature, - cells=cells, - interpolation_data=interpolation_data, - ) + self.value_fenics = fem.Function(function_space) - # override the kwargs with the temperature_sub - kwargs["T"] = temperature_sub + # store the expression of the boundary condition + # to update the value_fenics later + assert isinstance(self.value(**kwargs), ufl.core.expr.Expr), ( + f"{type(self.value(**kwargs))}" + ) + self.bc_expr = fem.Expression( + self.value(**kwargs), + helpers.get_interpolation_points(function_space.element), + ) - self.bc_expr = fem.Expression( - self.value(**kwargs), - helpers.get_interpolation_points(function_space.element), - ) - self.value_fenics.interpolate(self.bc_expr) + self.value_fenics.interpolate(self.bc_expr) # if K_S is provided, divide the value by K_S (change of variable method) if K_S is not None: diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index d4afbb271..fafc0ff8e 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1172,6 +1172,26 @@ def initialise(self): self.create_solver() self.initialise_exports() + def define_temperature(self): + super().define_temperature() + + # pass temperature function to each subdomain + if isinstance(self.temperature_fenics, fem.Function): + for subdomain in self.volume_subdomains: + element_CG = basix.ufl.element( + basix.ElementFamily.P, + subdomain.submesh.basix_cell(), + 1, # could expose? + basix.LagrangeVariant.equispaced, + ) + V = dolfinx.fem.functionspace(subdomain.submesh, element_CG) + sub_T = dolfinx.fem.Function(V) + from festim.helpers import nmm_interpolate + + nmm_interpolate(f_out=sub_T, f_in=self.temperature_fenics) + + subdomain.sub_T = sub_T + def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): """ Creates the ``value_fenics`` attribute for a given @@ -1190,8 +1210,16 @@ def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC): sub_V = bc.species.subdomain_to_function_space[volume_subdomain] collapsed_V, _ = sub_V.collapse() + # in the discontinuous case, if the temperature is given as a function + # then we can't use the temperature on the parent mesh + # see issue #1007 + if isinstance(self.temperature_fenics, fem.Function): + temp = volume_subdomain.sub_T + else: + temp = self.temperature_fenics + bc.create_value( - temperature=self.temperature_fenics, + temperature=temp, function_space=collapsed_V, t=self.t, ) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 06146759f..fa4476087 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -3,6 +3,7 @@ import dolfinx import numpy as np +from dolfinx import fem from dolfinx.mesh import Mesh, locate_entities from numpy import typing as npt from scifem.mesh import transfer_meshtags_to_submesh @@ -23,8 +24,20 @@ class VolumeSubdomain: Volume subdomain class Args: - id (int): the id of the volume subdomain - material (festim.Material): the material assigned to the subdomain + id: the id of the volume subdomain + submesh: the submesh of the volume subdomain + cell_map: the cell map of the volume subdomain + parent_mesh: the parent mesh of the volume subdomain + parent_to_submesh: the parent to submesh map of the volume subdomain + v_map: the vertex map of the volume subdomain + n_map: the normal map of the volume subdomain + facet_to_parent: the facet to parent map of the volume subdomain + ft: the facet meshtags of the volume subdomain + padded: whether the subdomain is padded (for 0.9 compatibility) + u: the solution function of the subdomain + u_n: the previous solution function of the subdomain + material: the material assigned to the subdomain + sub_T: the sub temperature field in the subdomain """ id: int @@ -40,6 +53,7 @@ class VolumeSubdomain: u: dolfinx.fem.Function u_n: dolfinx.fem.Function material: Material + sub_T: fem.Function | float def __init__(self, id, material, locator: Callable | None = None): self.id = id From d0a587da3d3538fce3c755b48b2aeab5391a85f2 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 30 Jan 2026 12:06:58 -0500 Subject: [PATCH 3/8] removed typoed attribute --- src/festim/hydrogen_transport_problem.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index fafc0ff8e..097373aaf 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -206,7 +206,6 @@ def __init__( self.temperature_fenics = None self._element_for_traps = "DG" - self.petcs_options = petsc_options self._temperature_as_function = None From 94b27157c58e4d95d85f98997984dfb787367532 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 30 Jan 2026 15:25:51 -0500 Subject: [PATCH 4/8] use dolfinx blocked solver instead of scifem --- src/festim/hydrogen_transport_problem.py | 94 +++++++++++++++++++----- 1 file changed, 75 insertions(+), 19 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 097373aaf..4f9cdaf21 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -3,6 +3,7 @@ from enum import Enum from mpi4py import MPI +from petsc4py import PETSc import adios4dolfinx import basix @@ -12,6 +13,8 @@ import tqdm.autonotebook import ufl from dolfinx import fem +from dolfinx.nls.petsc import NewtonSolver +from packaging.version import Version from scifem import BlockedNewtonSolver import festim.boundary_conditions @@ -1636,25 +1639,78 @@ def mixed_term(u, v, n): ) def create_solver(self): - self.solver = BlockedNewtonSolver( - self.forms, - [subdomain.u for subdomain in self.volume_subdomains], - J=self.J, - bcs=self.bc_forms, - petsc_options=self.petsc_options, - ) - self.solver.max_iterations = self.settings.max_iterations - self.solver.convergence_criterion = self.settings.convergence_criterion - self.solver.atol = ( - self.settings.atol - if not callable(self.settings.atol) - else self.settings.atol(float(self.t)) - ) - self.solver.rtol = ( - self.settings.rtol - if not callable(self.settings.rtol) - else self.settings.rtol(float(self.t)) - ) + if Version(dolfinx.__version__) == Version("0.9.0"): + self.solver = BlockedNewtonSolver( + self.forms, + [subdomain.u for subdomain in self.volume_subdomains], + J=self.J, + bcs=self.bc_forms, + petsc_options=self.petsc_options, + ) + self.solver.max_iterations = self.settings.max_iterations + self.solver.convergence_criterion = self.settings.convergence_criterion + self.solver.atol = ( + self.settings.atol + if not callable(self.settings.atol) + else self.settings.atol(float(self.t)) + ) + self.solver.rtol = ( + self.settings.rtol + if not callable(self.settings.rtol) + else self.settings.rtol(float(self.t)) + ) + + elif Version(dolfinx.__version__) > Version("0.9.0"): + from dolfinx.fem.petsc import NonlinearProblem + + if self.petsc_options is None: + # taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28 + use_superlu = ( + PETSc.IntType == np.int64 + ) # or PETSc.ScalarType == np.complex64 + sys = PETSc.Sys() # type: ignore + if sys.hasExternalPackage("mumps") and not use_superlu: + linear_solver = "mumps" + elif sys.hasExternalPackage("superlu_dist"): + linear_solver = "superlu_dist" + else: + linear_solver = "petsc" + + petsc_options = { + "snes_type": "newtonls", + "snes_linesearch_type": "none", + "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) + * 1e-2, + # TODO : make atol and rtol callable + "snes_atol": self.settings.atol, + "snes_rtol": self.settings.rtol, + "snes_max_it": self.settings.max_iterations, + "snes_divergence_tolerance": "PETSC_UNLIMITED", + "ksp_type": "preonly", + "pc_type": "lu", + "snes_monitor": None, + "ksp_monitor": None, + "pc_factor_mat_solver_type": linear_solver, + } + else: + petsc_options = self.petsc_options + + self.solver = NonlinearProblem( + self.forms, + [subdomain.u for subdomain in self.volume_subdomains], + bcs=self.bc_forms, + J=self.J, + petsc_options=petsc_options, + petsc_options_prefix="festim_solver", + ) + + # Delete PETSc options post setting them, ref: + # https://gitlab.com/petsc/petsc/-/issues/1201 + snes = self.solver.solver + prefix = snes.getOptionsPrefix() + opts = PETSc.Options() + for k in petsc_options.keys(): + del opts[f"{prefix}{k}"] def create_flux_values_fenics(self): """For each particle flux create the ``value_fenics`` attribute""" From 0a0422d21a3e705a1177bdc6809d774034f161de Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 30 Jan 2026 16:07:49 -0500 Subject: [PATCH 5/8] fix iterate to work with 0.10 solver --- src/festim/hydrogen_transport_problem.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 4f9cdaf21..000265f61 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1886,7 +1886,15 @@ def iterate(self): self.update_time_dependent_values() # Solve main problem - nb_its, converged = self.solver.solve() + if Version(dolfinx.__version__) == Version("0.9.0"): + nb_its, converged = self.solver.solve(self.u) + elif Version(dolfinx.__version__) > Version("0.9.0"): + _ = self.solver.solve() + converged_reason = self.solver.solver.getConvergedReason() + assert converged_reason > 0, ( + f"Non-linear solver did not converge. Reason code: {converged_reason}. \n See https://petsc.org/release/manualpages/SNES/SNESConvergedReason/ for more information." + ) + nb_its = self.solver.solver.getIterationNumber() # post processing self.post_processing() From 54b8443cbc1c837286246530348e98081d34cd64 Mon Sep 17 00:00:00 2001 From: jhdark Date: Wed, 11 Feb 2026 13:43:41 -0500 Subject: [PATCH 6/8] fix for 0.9 --- 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 000265f61..4ca121026 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1887,7 +1887,7 @@ def iterate(self): # Solve main problem if Version(dolfinx.__version__) == Version("0.9.0"): - nb_its, converged = self.solver.solve(self.u) + nb_its, converged = self.solver.solve() elif Version(dolfinx.__version__) > Version("0.9.0"): _ = self.solver.solve() converged_reason = self.solver.solver.getConvergedReason() From 0459f32c61a4304dc924b33dd649ef36029821e1 Mon Sep 17 00:00:00 2001 From: jhdark Date: Wed, 11 Feb 2026 13:45:21 -0500 Subject: [PATCH 7/8] 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), From 1d28fb0e4bf653c17705d9b22391732fcdd15051 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 17 Feb 2026 10:03:43 -0500 Subject: [PATCH 8/8] test case --- test/test_h_transport_problem.py | 42 ++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 95de032d7..71dbcda84 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1457,3 +1457,45 @@ def test_surface_reaction_BC_discontinuous(): # TEST assert np.isclose(permeation_flux.data[-1], -inlet_flux.data[-1], rtol=1e-5) + + +def test_temperature_as_function_in_discontinuous(): + """Test that a discontinuous problem can be initialised with a temperature field + given as a function""" + + # BUILD + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat) + vol2 = F.VolumeSubdomain1D(id=2, borders=[1, 2], material=mat) + + vertices_left = np.linspace(0, 1, num=100, endpoint=False) + vertices_right = np.linspace(1, 2, num=100) + vertices = np.concatenate([vertices_left, vertices_right]) + my_model.mesh = F.Mesh1D(vertices) + my_model.mesh = F.Mesh1D(vertices) + + my_model.subdomains = [vol1, vol2] + + my_model.method_interface = F.InterfaceMethod.penalty + my_model.interfaces = [ + F.Interface(id=5, subdomains=[vol1, vol2], penalty_term=100), + ] + + H = F.Species("H", subdomains=my_model.volume_subdomains) + my_model.species = [H] + + V = fem.functionspace(my_model.mesh.mesh, ("CG", 1)) + u = fem.Function(V) + u.interpolate(lambda x: 100 + x[0]) + u.x.array[:] = 100 + my_model.temperature = u + + my_model.settings = F.Settings( + atol=1e-10, + rtol=1e-09, + transient=False, + ) + + my_model.initialise()