From b38a95ec9010994d9eb36d7360c97e4c8c4b8b7c Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Thu, 12 Nov 2026 08:59:56 +0000 Subject: [PATCH] Add advective terms --- docs/api/models.md | 11 +- src/sumfem/models/__init__.py | 11 +- src/sumfem/models/advection.py | 151 +++++++++++++++++++++++++++ tests/models/test_advection.py | 91 ++++++++++++++++ tests/models/test_models_from_ufl.py | 30 ++++++ 5 files changed, 290 insertions(+), 4 deletions(-) create mode 100644 src/sumfem/models/advection.py create mode 100644 tests/models/test_advection.py diff --git a/docs/api/models.md b/docs/api/models.md index 86ae22e..8570a64 100644 --- a/docs/api/models.md +++ b/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: diff --git a/src/sumfem/models/__init__.py b/src/sumfem/models/__init__.py index c5cc111..ccb415a 100644 --- a/src/sumfem/models/__init__.py +++ b/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", diff --git a/src/sumfem/models/advection.py b/src/sumfem/models/advection.py new file mode 100644 index 0000000..0f10cb7 --- /dev/null +++ b/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 diff --git a/tests/models/test_advection.py b/tests/models/test_advection.py new file mode 100644 index 0000000..cf169df --- /dev/null +++ b/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) diff --git a/tests/models/test_models_from_ufl.py b/tests/models/test_models_from_ufl.py index 96afc9d..4a1be3e 100644 --- a/tests/models/test_models_from_ufl.py +++ b/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