fix Nitsche symmetry term on cylindrical and spherical meshes

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
 

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):

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

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

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

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",
+        )