Add advective terms

docs/api/models.md

@@ -4,7 +4,7 @@
 .. automodule:: sumfem.models
 ```
 
-## Accumulation, diffusion and sources
+## Accumulation, diffusion, advection and sources
 
 ```{eval-rst}
 .. autoclass:: Accumulation
@@ -13,6 +13,15 @@
 .. autoclass:: Diffusion
    :members:
 
+.. autoclass:: Advection
+   :members:
+
+.. autoclass:: ConvectiveAdvection
+   :members:
+
+.. autoclass:: ConservativeAdvection
+   :members:
+
 .. autoclass:: Source
    :members:
 

src/sumfem/models/__init__.py

@@ -1,8 +1,9 @@
 """Models built on :mod:`sumfem`, to be summed and extended.
 
-- Accumulation, diffusion and sources: :class:`Accumulation`,
-  :class:`Diffusion`, :class:`Source` and :class:`Poisson`, and Ohmic
-  :class:`Conduction`.
+- Accumulation, diffusion, advection and sources: :class:`Accumulation`,
+  :class:`Diffusion`, :class:`ConvectiveAdvection` and
+  :class:`ConservativeAdvection`, :class:`Source` and :class:`Poisson`, and
+  Ohmic :class:`Conduction`.
 - Boundaries: :class:`BoundaryFlux`, :class:`Outflow`, and constraints
   with Lagrange multipliers, :class:`BoundaryConstraint` and
   :class:`DirichletMultiplier`.
@@ -25,6 +26,7 @@
 Models of electrochemistry are in :mod:`sumfem.electrochemistry`.
 """
 
+from .advection import Advection, ConservativeAdvection, ConvectiveAdvection
 from .boundary import BoundaryConstraint, BoundaryFlux, DirichletMultiplier, Outflow
 from .charge import Conduction
 from .diffusion import Accumulation, Diffusion, Poisson, Source
@@ -60,6 +62,9 @@
     "RateReaction",
     "Reaction",
     "Accumulation",
+    "Advection",
+    "ConservativeAdvection",
+    "ConvectiveAdvection",
     "BoundaryConstraint",
     "BoundaryFlux",
     "Conduction",

src/sumfem/models/advection.py

@@ -0,0 +1,151 @@
+"""Advection of a field by a given velocity."""
+
+from __future__ import annotations
+
+import numpy.typing as npt
+import ufl
+
+from sumfem.symbolic import Field, Residual, as_expr
+
+__all__ = ["Advection", "ConservativeAdvection", "ConvectiveAdvection"]
+
+
+class Advection(Residual):
+    r"""Advection of a field by a velocity, in a form of a subclass.
+
+    The velocity is any expression, e.g. a field solved by a flow model or
+    data transferred from another mesh. The advected quantity is the field,
+    or an expression of it, e.g. an enthalpy of a temperature. The
+    subclasses write the term in the convective form,
+    :class:`ConvectiveAdvection`, or the conservative one,
+    :class:`ConservativeAdvection`.
+
+    For a continuous field :math:`c` and test function :math:`v`,
+    integration by parts on each cell :math:`K` relates the two forms,
+
+    .. math::
+       :label: advection-forms
+
+       \begin{aligned}
+         -(c \mathbf{u}, \nabla v)
+           &= (\mathbf{u} \cdot \nabla c, v)
+           + \sum_K (c \nabla \cdot \mathbf{u}, v)_K
+           - \sum_F ([\![\mathbf{u} \cdot \mathbf{n}]\!] c, v)_F
+           - (\mathbf{u} \cdot \mathbf{n} c, v)_{\partial \Omega},
+       \end{aligned}
+
+    with the jump :math:`[\![\mathbf{u} \cdot \mathbf{n}]\!]
+    = \mathbf{u}^+ \cdot \mathbf{n}^+ + \mathbf{u}^- \cdot \mathbf{n}^-` of
+    the normal velocity across the interior facets :math:`F`, and the
+    outward normal :math:`\mathbf{n}`.
+
+    The term holds only the integrals over the cells. A field that is not
+    continuous, e.g. with a discontinuous element or one mapped by a Piola
+    map, needs a numerical flux on the interior facets as well, e.g. an
+    upwind flux, added as a model of its own.
+
+    Args:
+        name: Name of the model.
+        field: The field ``c`` whose equation has the term.
+        velocity: The velocity ``u``.
+        markers: Cell markers restricting the term. Empty for the whole
+            domain.
+        quantity: The advected quantity, in place of ``c``. Defaults to
+            the field.
+    """
+
+    _field: Field
+    """The advected field."""
+
+    _velocity: ufl.core.expr.Expr
+    """The velocity."""
+
+    _markers: tuple[str, ...]
+    """Cell markers restricting the term."""
+
+    _quantity: ufl.core.expr.Expr
+    """The advected quantity."""
+
+    def __init__(
+        self,
+        name: str,
+        field: Field,
+        velocity: ufl.core.expr.Expr | npt.ArrayLike,
+        markers: tuple[str, ...] = (),
+        quantity: ufl.core.expr.Expr | None = None,
+    ):
+        super().__init__(name)
+        self._field = field
+        self._velocity = as_expr(velocity)
+        self._markers = tuple(markers)
+        self._quantity = field if quantity is None else as_expr(quantity)
+
+    @property
+    def field(self) -> Field:
+        """The advected field."""
+        return self._field
+
+    @property
+    def velocity(self) -> ufl.core.expr.Expr:
+        """The velocity."""
+        return self._velocity
+
+    @property
+    def quantity(self) -> ufl.core.expr.Expr:
+        """The advected quantity."""
+        return self._quantity
+
+    def flux(self) -> ufl.core.expr.Expr:
+        """The advective flux ``c u``, the outer product for a vector ``c``."""
+        c = self._quantity
+        return ufl.outer(c, self._velocity) if c.ufl_shape else c * self._velocity
+
+
+class ConvectiveAdvection(Advection):
+    """Advection in convective form, ``(grad(c) u, v)``.
+
+    It has no boundary term. It misses the jumps of the normal velocity in
+    :eq:`advection-forms`, so it equals :class:`ConservativeAdvection`, up
+    to the boundary term, only for a velocity that is divergence free in
+    each cell with a continuous normal component, e.g. one in H(div).
+
+    Args:
+        name: Name of the model.
+        field: The field ``c`` whose equation has the term.
+        velocity: The velocity ``u``.
+        markers: Cell markers restricting the term. Empty for the whole
+            domain.
+        quantity: The advected quantity, in place of ``c``. Defaults to
+            the field.
+    """
+
+    def residual(self) -> ufl.Form:
+        c, v = self._quantity, self._field.v
+        dx = self._field.domain.dx(*self._markers)
+        return ufl.inner(ufl.dot(ufl.grad(c), self._velocity), v) * dx
+
+
+class ConservativeAdvection(Advection):
+    """Advection in conservative form, ``-(c u, grad(v))``.
+
+    Its boundary term is the advective flux ``(u . n c, v)``:
+    :class:`~sumfem.models.Outflow` on outlets, and the inflow on inlets.
+    It holds the jumps of the normal velocity in :eq:`advection-forms`, so
+    it needs no facet term, and conserves the advected quantity for any
+    velocity: with ``v = 1`` the interior terms cancel. Use it with a
+    velocity whose normal component jumps across facets, e.g. the gradient
+    of a Lagrange pressure, as in :class:`~sumfem.models.Darcy` flow.
+
+    Args:
+        name: Name of the model.
+        field: The field ``c`` whose equation has the term.
+        velocity: The velocity ``u``.
+        markers: Cell markers restricting the term. Empty for the whole
+            domain.
+        quantity: The advected quantity, in place of ``c``. Defaults to
+            the field.
+    """
+
+    def residual(self) -> ufl.Form:
+        dx = self._field.domain.dx(*self._markers)
+        return -ufl.inner(self.flux(), ufl.grad(self._field.v)) * dx

tests/models/test_advection.py

@@ -0,0 +1,91 @@
+"""Advection-diffusion with manufactured solutions."""
+
+from mpi4py import MPI
+
+import dolfinx.mesh
+import numpy as np
+import pytest
+import ufl
+
+from sumfem.discrete import DomainMaps, MeshData, discretize
+from sumfem.models import (
+    Advection,
+    BoundaryFlux,
+    ConservativeAdvection,
+    ConvectiveAdvection,
+    Diffusion,
+    Outflow,
+    Source,
+)
+from sumfem.post_processing import error
+from sumfem.solvers import Solver
+from sumfem.symbolic import DirichletBC, Domain, Field
+
+
+def square(n: int) -> MeshData:
+    """A unit square whose outlet is the side ``x = 1``, and inlet the rest.
+
+    A boundary facet is marked if all its vertices are, so the inlet's
+    locator takes the corners of the outlet too, for the facets beside them.
+    """
+    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)
+    return MeshData.from_locators(
+        mesh,
+        {
+            1: {
+                "Inlet": lambda x: (
+                    ~np.isclose(x[0], 1.0)
+                    | np.isclose(x[1], 0.0)
+                    | np.isclose(x[1], 1.0)
+                ),
+                "Outlet": lambda x: np.isclose(x[0], 1.0),
+            }
+        },
+        boundary=True,
+    )
+
+
+def exact(x):
+    return 1 + ufl.sin(np.pi * x[0]) * x[1] ** 2
+
+
+def solve(n: int, form: type[Advection], degree: int = 2) -> float:
+    """Advection-diffusion with Dirichlet data on the inlet and a free outlet.
+
+    The equation is ``-div(D grad c) + div(u c) = f`` for the conservative
+    form and ``-div(D grad c) + u . grad(c) = f`` for the convective one. The
+    velocity has a divergence, so the two differ. The outlet, ``x = 1``,
+    gets the exact diffusive flux, and in conservative form the advective
+    one by :class:`Outflow`.
+    """
+    domain = Domain("square", "triangle")
+    c = Field("c", domain, ("Lagrange", degree))
+    x = ufl.SpatialCoordinate(domain)
+    normal = ufl.FacetNormal(domain)
+    u = ufl.as_vector((1 + x[0], 0.5 * x[1]))
+    D = 0.1
+    c_exact = exact(x)
+    if form is ConservativeAdvection:
+        advection = ufl.div(u * c_exact)
+    else:
+        advection = ufl.dot(u, ufl.grad(c_exact))
+    model = (
+        Diffusion("diffusion", c, coefficient=D)
+        + form("advection", c, u)
+        + Source("source", c, source=-ufl.div(D * ufl.grad(c_exact)) + advection)
+        + BoundaryFlux(
+            "diffusive_flux", c, D * ufl.dot(ufl.grad(c_exact), normal), ("Outlet",)
+        )
+    )
+    if form is ConservativeAdvection:
+        model = model + Outflow("outflow", c, u, ("Outlet",))
+    disc = discretize([model], DomainMaps({"square": square(n)}))
+    (solution,) = Solver(disc[model], [DirichletBC(c, "Inlet", c_exact)]).solve()
+    return error(solution, exact)
+
+
+@pytest.mark.parametrize("form", [ConvectiveAdvection, ConservativeAdvection])
+def test_advection_diffusion_converges(form):
+    errors = np.array([solve(n, form) for n in (4, 8)])
+    rate = np.log2(errors[0] / errors[1])
+    assert rate > 2.8, (errors, rate)

tests/models/test_models_from_ufl.py

@@ -20,6 +20,8 @@
     BoundaryFlux,
     Brinkman,
     BrinkmanForchheimer,
+    ConservativeAdvection,
+    ConvectiveAdvection,
     Darcy,
     Diffusion,
     DirichletMultiplier,
@@ -148,6 +150,34 @@ def test_boundary_flux_and_outflow(cell, ufl_functionspace, box, compare):
     compare(plain, [plain.field], model, [c], box(cell), values)
 
 
+@pytest.mark.parametrize("form", [ConvectiveAdvection, ConservativeAdvection])
+@pytest.mark.parametrize("cell", CELLS)
+def test_advection(cell, form, ufl_functionspace, box, compare):
+    """Advection of a quantity of the field by a velocity that varies in
+    space, in either form, restricted to the cells of a marker.
+    """
+    gdim = ufl.Cell(cell).topological_dimension
+    mesh, V = ufl_functionspace(cell, degree=2)
+    ch, v = ufl.Coefficient(V), ufl.TestFunction(V)
+    xh = ufl.SpatialCoordinate(mesh)
+    velocity = ufl.as_vector([1.0 + i + xh[i] for i in range(gdim)])
+    quantity = ch + ch**2
+    if form is ConservativeAdvection:
+        F = -ufl.inner(quantity * velocity, ufl.grad(v)) * ufl.dx(1)
+    else:
+        F = ufl.dot(ufl.grad(quantity), velocity) * v * ufl.dx(1)
+    plain = UFLResidual(
+        "plain", F, ch, names={DOMAIN: mesh}, markers={"dx": {1: "LeftHalf"}}
+    )
+
+    domain = Domain(DOMAIN, cell)
+    c = Field("c", domain, ("Lagrange", 2))
+    x = ufl.SpatialCoordinate(domain)
+    u = ufl.as_vector([1.0 + i + x[i] for i in range(gdim)])
+    model = form("a", c, u, markers=("LeftHalf",), quantity=c + c**2)
+    compare(plain, [plain.field], model, [c], box(cell))
+
+
 @pytest.mark.parametrize("cell", CELLS)
 def test_dirichlet_multiplier(cell, ufl_functionspace, box, compare):
     """A Dirichlet condition with a Lagrange multiplier, the two fields on a