Codim: fix for the flux conservation bug across internal interfaces

docs/source/userguide/subdomains.rst

@@ -81,6 +81,194 @@ For 1D domains, use the :class:`festim.VolumeSubdomain1D` class, which requires
     my_material = Material(D_0=1, E_D=1)
     my_1D_volume = VolumeSubdomain1D(id=1, material=my_material, borders=[0, 1])
 
+Codimensional (manifold) Subdomains
+====================================
+
+A volume subdomain can be a *manifold* embedded in the mesh: a line in a 2D mesh, or a
+surface in a 3D mesh. Such a subdomain carries its own transport equation, with
+diffusion and advection **along** the manifold, and is coupled to the bulk by a flux.
+This is the way to model transport in a grain boundary, along a crack, or in a thin
+surface layer you do not want to resolve with cells.
+
+Pass ``dim`` one less than the dimension of the mesh:
+
+.. testcode::
+
+    import numpy as np
+    from festim import VolumeSubdomain, Material
+
+    gamma = VolumeSubdomain(
+        id=2,
+        material=Material(D_0=1e-6, E_D=0.0),
+        dim=1,  # a line inside a 2D mesh
+        locator=lambda x: np.isclose(x[0], 0.0),
+    )
+
+A manifold subdomain is tagged in the **facet** meshtags rather than the cell meshtags,
+so its ``id`` must be unique among the surface subdomains as well. It can be used
+directly wherever a surface is expected, for instance as the ``subdomain`` of a
+:class:`festim.ParticleFluxBC` — there is no need to declare a separate
+:class:`festim.SurfaceSubdomain` on the same facets.
+
+A manifold may sit on the outer boundary of the domain, or *between* two volume
+subdomains — a grain boundary, or an interface layer with its own trapping — in which
+case it exchanges with both sides.
+
+Coupling to the bulk
+--------------------
+
+The exchange is written twice: once as a flux leaving the bulk, and once as a source
+entering the manifold. Use ``species_dependent_value`` to let each half see both
+concentrations, even though they live on different meshes::
+
+    k = 0.1
+    J = lambda c_bulk, c_man: k * (c_bulk - c_man)
+
+    # the bulk loses J through gamma
+    flux = F.ParticleFluxBC(
+        subdomain=gamma,
+        value=lambda c_man, c_bulk: -J(c_bulk, c_man),
+        species=H_bulk,
+        species_dependent_value={"c_bulk": H_bulk, "c_man": H_manifold},
+    )
+
+    # ... and the manifold gains it
+    source = F.ParticleSource(
+        volume=gamma,
+        value=lambda c_man, c_bulk: J(c_bulk, c_man),
+        species=H_manifold,
+        species_dependent_value={"c_bulk": H_bulk, "c_man": H_manifold},
+    )
+
+.. warning::
+
+    **Mind the units.** A :class:`festim.ParticleFluxBC` value is a *flux* (H/m²/s for a
+    3D bulk) whereas a :class:`festim.ParticleSource` value is a *volumetric rate*
+    (H/m³/s). Writing the same expression on both sides is therefore not generally
+    dimensionally consistent, and the manifold-side source usually needs a conversion
+    factor.
+
+    FESTIM does not impose a unit convention on a manifold species: depending on what
+    you are modelling it may be a volumetric concentration H/m³ (a layer of thickness
+    :math:`\lambda`, in which case the source is :math:`J/\lambda`), an areal density
+    H/m² (an adsorbed layer, source :math:`J`), or a line density H/m (a grain
+    boundary). Keeping the problem dimensionally consistent is up to you.
+
+Manifolds between two subdomains
+--------------------------------
+
+When a manifold separates two volume subdomains, declare **one exchange per side** —
+one :class:`festim.ParticleFluxBC` and one :class:`festim.ParticleSource` each. Both
+name the same manifold as their subdomain; FESTIM works out which side each belongs to
+from the bulk species it reads, so nothing else has to be specified::
+
+    for bulk_species, k in ((H_left, k_left), (H_right, k_right)):
+        J = lambda c_man, c_bulk: k * (c_bulk - c_man)
+        bcs.append(F.ParticleFluxBC(
+            subdomain=gamma, species=bulk_species,
+            value=lambda c_man, c_bulk: -J(c_man, c_bulk),
+            species_dependent_value={"c_bulk": bulk_species, "c_man": H_manifold}))
+        sources.append(F.ParticleSource(
+            volume=gamma, species=H_manifold, value=J,
+            species_dependent_value={"c_bulk": bulk_species, "c_man": H_manifold}))
+
+A single source may not read the bulk concentrations of *both* sides at once: an
+interior manifold is integrated over interior facets, where each term has to be
+restricted to one side. Split such a source in two, as above.
+
+.. note::
+
+    A pair of volume subdomains may be separated either by a
+    :class:`festim.Interface` -- imposing a jump in concentration across a shared
+    boundary -- or by a codim-1 subdomain carrying its own transport equation, but not
+    both. FESTIM raises if an interface and a manifold cover the same facets.
+
+Advection along a manifold
+--------------------------
+
+An :class:`festim.AdvectionTerm` on a manifold subdomain takes an ordinary ambient
+velocity vector — 2 components in a 2D mesh, 3 in a 3D mesh. There is no need to
+project it onto the manifold: the tangential gradient is orthogonal to the normal, so
+:math:`v \cdot \nabla_\Gamma c` automatically ignores the normal component of
+:math:`v`.
+
+Boundary conditions on a manifold
+---------------------------------
+
+A manifold has a boundary of its own — the endpoints of a line in a 2D mesh, the rim of
+a surface in a 3D mesh — and boundary conditions can be applied there. Declare it as a
+:class:`festim.SurfaceSubdomain` with ``dim`` set to the mesh dimension minus **two**,
+just as a manifold is a :class:`festim.VolumeSubdomain` with ``dim`` set to the mesh
+dimension minus one:
+
+.. code-block:: python
+
+    # a 1D fluid running along a 2D pipe wall
+    fluid = F.VolumeSubdomain(id=2, material=..., dim=1,
+                              locator=lambda x: np.isclose(x[1], H))
+
+    # the inlet: one end of that 1D domain
+    inlet = F.SurfaceSubdomain(id=3, dim=0,
+                               locator=lambda x: np.isclose(x[0], 0.0))
+
+    ...
+    boundary_conditions=[
+        F.FixedConcentrationBC(subdomain=inlet, value=c_in, species=c_fluid),
+    ]
+
+The locator is evaluated on the manifold, not on the parent mesh, and must select a
+point on its boundary — a locator matching only interior points raises rather than
+silently doing nothing.
+
+Such a surface carries no meshtag, so its ``id`` does not have to differ from a manifold
+or interface id. Which manifold it bounds is taken from the ``species`` of the boundary
+condition using it, so that species must live on exactly one manifold; the same surface
+object can be reused on several manifolds, one species each.
+
+Without a condition of this kind, the ends of a manifold carry the natural zero-flux
+condition.
+
+Reactions and trapping on a manifold
+------------------------------------
+
+A :class:`festim.Reaction` runs on a manifold like on any other volume subdomain: give
+it ``volume=gamma`` and species that live on ``gamma``. Trapping is written as a
+reaction against :class:`festim.ImplicitSpecies` empty sites:
+
+.. code-block:: python
+
+    trapped = F.Species("trapped", mobile=False, subdomains=[gamma])
+    empty_sites = F.ImplicitSpecies(n=n_trap, others=[trapped])
+
+    trapping = F.Reaction(
+        reactant=[H_manifold, empty_sites], product=trapped,
+        k_0=k_0, E_k=E_k, p_0=p_0, E_p=E_p, volume=gamma,
+    )
+
+Note that the trapped species is declared with its own ``subdomains``, and that
+:class:`festim.Trap` is not a shortcut for this — it builds a species without one.
+
+The density ``n`` of an implicit species consumed on a manifold is a coefficient of an
+integral over that manifold, so FESTIM builds it there. Two consequences: give ``n`` as
+a float or as a callable of ``x`` and ``t`` rather than as a ready-made
+``dolfinx.fem.Function``, which cannot be moved; and declare one implicit species per
+subdomain rather than sharing one between a reaction on a manifold and a reaction
+elsewhere. Both are raised rather than silently mis-assembled.
+
+Limitations
+-----------
+
+* Only codimension 1 is supported (``dim`` must be the mesh dimension minus one), and a
+  manifold must be adjacent to one volume subdomain (on the boundary of the domain) or
+  two (on an interface). A codimension-2 subdomain carrying its own equation is not
+  supported: a bulk field has no well-defined trace on a line in 3D or a point in 2D,
+  so the exchange with it would not be well posed.
+* Boundary conditions on the boundary of a manifold are limited to
+  :class:`festim.FixedConcentrationBC`.
+* Dedicated exports on a manifold or on its boundary are not available;
+  :class:`festim.VTXSpeciesExport` on a manifold works.
+* Cartesian coordinates only.
+
 ----------
 Materials
 ----------

mwe.py

@@ -0,0 +1,92 @@
+from mpi4py import MPI
+
+import numpy as np
+import ufl
+from dolfinx.mesh import create_unit_square
+
+import festim as F
+
+my_mesh = create_unit_square(MPI.COMM_WORLD, 20, 20)
+
+# Create materials for the two subdomains
+mat_omega = F.Material(D_0=1.5, E_D=0.0)
+mat_gamma = F.Material(D_0=0.7, E_D=0.0)
+
+
+# Subdomains. gamma has dim=1: it is a manifold embedded in the 2D mesh, carrying its
+# own transport equation. It is tagged in the facet meshtags and can be used directly
+# wherever a surface is expected.
+omega = F.VolumeSubdomain(
+    id=1, material=mat_omega, locator=lambda x: np.full_like(x[0], True), name="omega"
+)
+gamma = F.VolumeSubdomain(
+    id=2,
+    material=mat_gamma,
+    dim=1,
+    locator=lambda x: np.isclose(x[0], 0.0),
+    name="gamma",
+)
+right = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 1.0))
+
+H_om = F.Species("H_om", subdomains=[omega])
+H_gam = F.Species("H_gam", subdomains=[gamma])
+species = [H_om, H_gam]
+
+k = 2.0
+beta = 1.0 - 1.5 / k
+
+# Volumetric source in the gamma subdomain
+source_in_gamma = F.ParticleSource(
+    value=lambda x: (0.7 * np.pi**2 * beta - 1.5) * ufl.cos(np.pi * x[1]),
+    species=H_gam,
+    volume=gamma,
+)
+
+# Coupling between the two subdomains: the same flux J = k (c_omega - c_gamma) enters
+# the gamma equation as a source and leaves omega as a flux boundary condition.
+coupling_source_gamma = F.ParticleSource(
+    value=lambda c_g, c_o: k * (c_o - c_g),
+    species=H_gam,
+    volume=gamma,
+    species_dependent_value={"c_o": H_om, "c_g": H_gam},
+)
+coupling_flux_omega = F.ParticleFluxBC(
+    subdomain=gamma,
+    value=lambda c_g, c_o: k * (c_g - c_o),
+    species=H_om,
+    species_dependent_value={"c_o": H_om, "c_g": H_gam},
+)
+
+exact_omega = F.FixedConcentrationBC(
+    subdomain=right,
+    value=lambda x: 1 + x[0] ** 2 + (1 + x[0]) * ufl.cos(np.pi * x[1]),
+    species=H_om,
+)
+
+source_omega = F.ParticleSource(
+    value=lambda x: -1.5 * (2 - np.pi**2 * (1 + x[0]) * ufl.cos(np.pi * x[1])),
+    species=H_om,
+    volume=omega,
+)
+
+my_model = F.HydrogenTransportProblemDiscontinuous(
+    mesh=F.Mesh(my_mesh),
+    species=species,
+    subdomains=[omega, gamma, right],
+    sources=[source_in_gamma, coupling_source_gamma, source_omega],
+    boundary_conditions=[coupling_flux_omega, exact_omega],
+    temperature=500,
+    exports=[
+        F.VTXSpeciesExport(filename="H_om.bp", field=H_om, subdomain=omega),
+        F.VTXSpeciesExport(filename="H_gam.bp", field=H_gam, subdomain=gamma),
+    ],
+)
+my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+my_model.initialise()
+my_model.run()
+
+c_om = H_om.subdomain_to_post_processing_solution[omega]
+c_gam = H_gam.subdomain_to_post_processing_solution[gamma]
+print("c_omega range", c_om.x.array.min(), c_om.x.array.max())
+print("c_gamma range", c_gam.x.array.min(), c_gam.x.array.max())
+print("expected c_gamma range", 1 - beta, 1 + beta)

src/festim/advection.py

@@ -149,11 +149,14 @@ def convert_input_value(
                 )
 
         # create vector function space and function
+        # NOTE: the shape is the *geometric* dimension, not the topological one: on a
+        # codim-1 subdomain the mesh is a manifold (eg. a line in 2D) whose cells are
+        # 1D but whose points, and therefore velocities and gradients, are ambient
         v_cg = basix.ufl.element(
             "Lagrange",
             function_space.mesh.topology.cell_name(),
             1,
-            shape=(function_space.mesh.topology.dim,),
+            shape=(function_space.mesh.geometry.dim,),
         )
         self.vector_function_space = fem.functionspace(function_space.mesh, v_cg)
         self.fenics_object = fem.Function(self.vector_function_space)

src/festim/boundary_conditions/dirichlet_bc.py

@@ -49,7 +49,7 @@ class DirichletBCBase:
         | Callable[[np.ndarray, float], np.ndarray]
         | Callable[[float], float]
     )
-    value_fenics: None | fem.Function | fem.Constant | np.ndarray | float
+    value_fenics: fem.Function | fem.Constant | np.ndarray | float | None
     bc_expr: fem.Expression
 
     def __init__(
@@ -74,7 +74,7 @@ def value_fenics(self):
     @value_fenics.setter
     def value_fenics(
         self,
-        value: None | fem.Function | fem.Constant | np.ndarray | ufl.core.expr.Expr,
+        value: fem.Function | fem.Constant | np.ndarray | ufl.core.expr.Expr | None,
     ):
         if value is None:
             self._value_fenics = value

src/festim/boundary_conditions/flux_bc.py

@@ -3,6 +3,7 @@
 from dolfinx import fem
 
 import festim as F
+from festim.helpers import solution_on
 from festim.subdomain.surface_subdomain import SurfaceSubdomain
 
 
@@ -245,9 +246,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                     if species.concentration is not None:
                         kwargs[name] = species.concentration
                     else:  # discontinuous case: one solution per subdomain
-                        kwargs[name] = species.subdomain_to_solution[
-                            self._volume_subdomain
-                        ]
+                        kwargs[name] = solution_on(species, self._volume_subdomain)
 
                 self.value_fenics = self.value(**kwargs)
 

src/festim/enclosure/enclosure.py

@@ -171,7 +171,8 @@ def convert_input_values_to_fenics_objects(self, function_space, t):
             function_space: a function space on the parent mesh
             t: the time, as a fenics Constant
         """
-        # NOTE could we have so guards to make sure temperature cannot be a function of space or temperature?
+        # NOTE could we have so guards to make sure temperature cannot be a function of
+        # space or temperature?
         self.temperature.convert_input_value(function_space=function_space, t=t)
         for opening in self.openings:
             opening.convert_input_values_to_fenics_objects(