From c699c4b5988529a2b54dde2df36806310de66a38 Mon Sep 17 00:00:00 2001 From: huihua Date: Tue, 3 Nov 2026 12:19:08 -0500 Subject: [PATCH 1/2] fix Nitsche symmetry term on cylindrical and spherical meshes --- .../boundary_conditions/dirichlet_bc.py | 27 +++++- src/festim/hydrogen_transport_problem.py | 8 +- test/system_tests/test_enclosure_dirichlet.py | 90 +++++++++++++++++++ 3 files changed, 121 insertions(+), 4 deletions(-) diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index d91ce0d5e..f0153add9 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -12,6 +12,7 @@ from festim import helpers from festim import subdomain as _subdomain +from festim.mesh import CoordinateSystem from festim.species import Species @@ -217,6 +218,7 @@ def weak_formulation( v: ufl.argument.Argument | ufl.indexed.Indexed, ds: ufl.Measure, D: ufl.core.expr.Expr | fem.Function | fem.Constant, + coordinate_system: CoordinateSystem = CoordinateSystem.CARTESIAN, ) -> ufl.core.expr.Expr: """ Returns the symmetric Nitsche weak formulation for the BC @@ -231,9 +233,14 @@ def weak_formulation( D: the diffusion coefficient of the species at this surface. It multiplies the consistency and symmetry terms, without which the scheme is only consistent when ``D == 1``. + coordinate_system: the coordinate system of the mesh. The symmetry term + must be written like the diffusion term of the bulk. Returns: the weak formulation + + Raises: + NotImplementedError: for an unknown coordinate system """ mesh = ds.ufl_domain() n = ufl.FacetNormal(mesh) @@ -243,8 +250,24 @@ def weak_formulation( # scheme conserves is defined in a single place form = self.numerical_flux(u, D, mesh) * v * ds_bc - # symmetry term, making the bilinear form symmetric (and the L2 error optimal) - form += -D * ufl.inner(n, ufl.grad(v)) * (u - self.value_fenics) * ds_bc + # symmetry term, making the bilinear form symmetric (and the L2 error optimal). + # Written like the diffusion term in HydrogenTransportProblem.create_formulation + # and in festim.drift.drift_form: r * grad(v / r) in cylindrical and + # r**2 * grad(v / r**2) in spherical. + match coordinate_system: + case CoordinateSystem.CARTESIAN: + grad_v = ufl.grad(v) + case CoordinateSystem.CYLINDRICAL: + r = ufl.SpatialCoordinate(mesh)[0] + grad_v = r * ufl.grad(v / r) + case CoordinateSystem.SPHERICAL: + r = ufl.SpatialCoordinate(mesh)[0] + grad_v = r**2 * ufl.grad(v / r**2) + case _: + raise NotImplementedError( + f"Unknown coordinate system {coordinate_system!s}" + ) + form += -D * ufl.inner(n, grad_v) * (u - self.value_fenics) * ds_bc return form diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 5a98b59b3..7dafed3a1 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1087,7 +1087,9 @@ def create_formulation(self): D = vol.material.get_diffusion_coefficient( self.mesh.mesh, self.temperature_fenics, bc.species ) - self.formulation += bc.weak_formulation(u, v, self.ds, D) + self.formulation += bc.weak_formulation( + u, v, self.ds, D, self.mesh.coordinate_system + ) for drift_term in self.drift_terms: vol = drift_term.subdomain @@ -2515,7 +2517,9 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): D = subdomain.material.get_diffusion_coefficient( self.mesh.mesh, self.temperature_fenics, bc.species ) - form_coupling += bc.weak_formulation(u, v, self.ds, D) + form_coupling += bc.weak_formulation( + u, v, self.ds, D, self.mesh.coordinate_system + ) # let drift carry the species out where the user says it flows out if isinstance(bc, boundary_conditions.OutflowBC): diff --git a/test/system_tests/test_enclosure_dirichlet.py b/test/system_tests/test_enclosure_dirichlet.py index 43b8d3b1e..140a932b3 100644 --- a/test/system_tests/test_enclosure_dirichlet.py +++ b/test/system_tests/test_enclosure_dirichlet.py @@ -404,3 +404,93 @@ def test_dirichlet_coupling_agrees_with_surface_reaction(): # have settled: the fast-exchange limit is an approximation, so this is not expected # to be exact. It agrees far better than that in practice (~1e-6). assert H2_d.value == pytest.approx(H2_r.value, rel=5e-3) + + +@pytest.mark.parametrize("law, stoichiometry", [("henry", 1), ("sieverts", 2)]) +def test_dirichlet_coupled_enclosure_conserves_particles_cylindrical( + law, stoichiometry +): + """A closed enclosure exchanging with a hollow cylinder through a weakly enforced + Henry's or Sieverts' BC must conserve hydrogen atoms at every timestep. + + The 1D model is per unit axial length, so the invariant is + ``length * int(c 2 pi r) dr + s*P*V/(k*T)``. + """ + r_in, r_out, axial_length = 0.5, 1.5, 2.0 + V_enc, T, P0 = 1e-6, 500.0, 1e5 + + my_model = F.HydrogenTransportProblemDiscontinuous() + my_model.mesh = F.Mesh1D( + vertices=np.linspace(r_in, r_out, 61), coordinate_system="cylindrical" + ) + material = F.Material(name="mat", D_0=1e-1, E_D=0.0) + vol = F.VolumeSubdomain1D(id=1, borders=[r_in, r_out], material=material) + inner = F.SurfaceSubdomain1D(id=1, x=r_in) + outer = F.SurfaceSubdomain1D(id=2, x=r_out) + my_model.subdomains = [vol, inner, outer] + + H = F.Species("H", subdomains=[vol]) + my_model.species = [H] + my_model.temperature = T + + H2 = F.GasSpecies(name="H2", initial_pressure=P0) + my_model.enclosures = [ + F.Enclosure( + volume=V_enc, species=[H2], temperature=T, surfaces={inner: axial_length} + ) + ] + if law == "henry": + coupling_bc = F.HenrysBC( + subdomain=inner, + H_0=1e15, + E_H=0.0, + pressure=H2, + species=H, + enforce_weakly=True, + penalty=100, + ) + else: + coupling_bc = F.SievertsBC( + subdomain=inner, + S_0=1e15, + E_S=0.0, + pressure=H2, + species=H, + enforce_weakly=True, + penalty=100, + ) + my_model.boundary_conditions = [coupling_bc] + my_model.initial_conditions = [ + F.InitialConcentration(value=0.0, species=H, volume=vol) + ] + # same reasoning as in build_slab_with_enclosure for the tolerances + my_model.settings = F.Settings( + atol=1e8, + rtol=1e-8, + transient=True, + final_time=200.0, + stepsize=F.Stepsize(10.0), + ) + my_model.show_progress_bar = False + my_model.initialise() + + r = ufl.SpatialCoordinate(vol.submesh)[0] + + def total_hydrogen_atoms(): + my_model.post_processing() + c = H.subdomain_to_post_processing_solution[vol] + in_solid = axial_length * MPI.COMM_WORLD.allreduce( + dolfinx.fem.assemble_scalar(dolfinx.fem.form(c * 2 * np.pi * r * ufl.dx)), + op=MPI.SUM, + ) + in_gas = stoichiometry * H2.value * V_enc / (F.k_B_SI * T) + return in_solid + in_gas + + initial_inventory = total_hydrogen_atoms() + + while my_model.t.value < my_model.settings.final_time: + my_model.iterate() + assert total_hydrogen_atoms() == pytest.approx(initial_inventory, rel=1e-10) + + # the gas must actually have been absorbed, otherwise the test is vacuous + assert H2.value < P0 From 2e567fff4a906d461e892911a77cdcfa67453897 Mon Sep 17 00:00:00 2001 From: huihua Date: Tue, 3 Nov 2026 13:32:02 -0500 Subject: [PATCH 2/2] test weakly enforced Dirichlet BCs on cylindrical and spherical meshes --- .../test_cylindrical_coordinates.py | 63 +++++++++++++++++++ .../test_spherical_coordinates.py | 63 +++++++++++++++++++ test/test_dirichlet_bc.py | 18 ++++++ 3 files changed, 144 insertions(+) diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py index 13f82d3c1..0c664d932 100644 --- a/test/system_tests/test_cylindrical_coordinates.py +++ b/test/system_tests/test_cylindrical_coordinates.py @@ -4,6 +4,7 @@ import dolfinx import numpy as np +import pytest import ufl import festim as F @@ -331,3 +332,65 @@ def c_exact(x): net_flux = sum(flux.value for flux in fluxes.values()) assert np.isclose(net_flux, -6 * D * volume, rtol=1e-6) + + +@pytest.mark.parametrize( + "model_class", + [F.HydrogenTransportProblem, F.HydrogenTransportProblemDiscontinuous], +) +def test_run_MMS_cylindrical_weak_dirichlet(model_class): + """Tests that festim produces the correct concentration field in cylindrical + coordinates when the Dirichlet BCs are enforced weakly.""" + + my_model = model_class() + my_model.mesh = F.Mesh1D( + vertices=np.linspace(1, 2, 500), coordinate_system="cylindrical" + ) + + def u_exact(x): + return 1 + x[0] ** 2 + + f = -4 + + my_mat = F.Material(D_0=1.0, E_D=0) + + left = F.SurfaceSubdomain1D(id=1, x=1) + right = F.SurfaceSubdomain1D(id=2, x=2) + my_vol = F.VolumeSubdomain1D(id=3, borders=[1, 2], material=my_mat) + my_model.subdomains = [my_vol, left, right] + + H = F.Species("H") + if model_class is F.HydrogenTransportProblemDiscontinuous: + H.subdomains = [my_vol] + my_model.species = [H] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC( + subdomain=surface, + value=u_exact, + species=H, + enforce_weakly=True, + penalty=100, + ) + for surface in [left, right] + ] + + my_model.temperature = 500 + + my_model.sources = [F.ParticleSource(value=f, volume=my_vol, species=H)] + + my_model.settings = F.Settings( + atol=1e-10, rtol=1e-9, max_iterations=50, transient=False + ) + + my_model.initialise() + my_model.run() + + if model_class is F.HydrogenTransportProblemDiscontinuous: + computed_solution = H.subdomain_to_post_processing_solution[my_vol] + else: + computed_solution = H.post_processing_solution + + L2_error = error_L2(computed_solution, u_exact) + + assert L2_error < 1e-6 diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py index 1e94b5595..0a41c5577 100644 --- a/test/system_tests/test_spherical_coordinates.py +++ b/test/system_tests/test_spherical_coordinates.py @@ -1,6 +1,7 @@ import math import numpy as np +import pytest from dolfinx import fem import festim as F @@ -214,3 +215,65 @@ def f_right(x): assert L2_error_l < 1e-06 assert L2_error_r < 1e-06 + + +@pytest.mark.parametrize( + "model_class", + [F.HydrogenTransportProblem, F.HydrogenTransportProblemDiscontinuous], +) +def test_run_MMS_spherical_weak_dirichlet(model_class): + """Tests that festim produces the correct concentration field in spherical + coordinates when the Dirichlet BCs are enforced weakly.""" + + my_model = model_class() + my_model.mesh = F.Mesh1D( + vertices=np.linspace(1, 2, 1000), coordinate_system="spherical" + ) + + def u_exact(x): + return 3 + 2 * x[0] ** 2 + + f = -12 + + my_mat = F.Material(D_0=1.0, E_D=0) + + left = F.SurfaceSubdomain1D(id=1, x=1) + right = F.SurfaceSubdomain1D(id=2, x=2) + my_vol = F.VolumeSubdomain1D(id=3, borders=[1, 2], material=my_mat) + my_model.subdomains = [my_vol, left, right] + + H = F.Species("H") + if model_class is F.HydrogenTransportProblemDiscontinuous: + H.subdomains = [my_vol] + my_model.species = [H] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC( + subdomain=surface, + value=u_exact, + species=H, + enforce_weakly=True, + penalty=100, + ) + for surface in [left, right] + ] + + my_model.temperature = 500 + + my_model.sources = [F.ParticleSource(value=f, volume=my_vol, species=H)] + + my_model.settings = F.Settings( + atol=1e-10, rtol=1e-9, max_iterations=50, transient=False + ) + + my_model.initialise() + my_model.run() + + if model_class is F.HydrogenTransportProblemDiscontinuous: + computed_solution = H.subdomain_to_post_processing_solution[my_vol] + else: + computed_solution = H.post_processing_solution + + L2_error = error_L2(computed_solution, u_exact) + + assert L2_error < 1e-6 diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index f30f8ec66..79a53deed 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -470,3 +470,21 @@ def test_bc_temperature_dependent_attribute(input, expected_value): my_bc = F.DirichletBC(subdomain, input, species) assert my_bc.temperature_dependent is expected_value + + +def test_weak_formulation_unknown_coordinate_system_raises(): + subdomain = F.SurfaceSubdomain1D(1, x=0) + bc = F.FixedConcentrationBC( + subdomain, value=1.0, species="test", enforce_weakly=True, penalty=10 + ) + bc.value_fenics = fem.Constant(mesh, default_scalar_type(1.0)) + V = fem.functionspace(mesh, ("Lagrange", 1)) + + with pytest.raises(NotImplementedError, match="Unknown coordinate system"): + bc.weak_formulation( + u=fem.Function(V), + v=ufl.TestFunction(V), + ds=ufl.Measure("ds", domain=mesh), + D=fem.Constant(mesh, default_scalar_type(1.0)), + coordinate_system="not a coordinate system", + )