Mixed formulation for the Poisson equation

This demo illustrates how to solve Poisson equation using a mixed (two-field) formulation. In particular, it illustrates how to

• Use mixed and non-continuous finite element spaces.

• Set essential boundary conditions for subspaces and $H(\mathrm{div})$ spaces.

Python script.
Jupyter notebook.

Equation and problem definition

An alternative formulation of Poisson equation can be formulated by introducing an additional (vector) variable, namely the (negative) flux: $\sigma =\mathrm{\nabla }u$. The partial differential equations then read

$$\begin{array}{rl}\sigma -\mathrm{\nabla }u& =0\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ \mathrm{\nabla }\cdot \sigma & =-f\mathrm{i}\mathrm{n}\mathrm{\Omega },\end{array}$$

with boundary conditions

$$\begin{array}{r}u=u_0\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_D,\\ \sigma \cdot n=g\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_N.\end{array}$$

The same equations arise in connection with flow in porous media, and are also referred to as Darcy flow. Here $n$ denotes the outward pointing normal vector on the boundary. Looking at the variational form, we see that the boundary condition for the flux ($\sigma \cdot n=g$) is now an essential boundary condition (which should be enforced in the function space), while the other boundary condition ($u=u_0$) is a natural boundary condition (which should be applied to the variational form). Inserting the boundary conditions, this variational problem can be phrased in the general form: find $(\sigma ,u)\in {\mathrm{\Sigma }}_g\times V$ such that

$$a((\sigma ,u),(\tau ,v))=L((\tau ,v))\mathrm{\forall }(\tau ,v)\in {\mathrm{\Sigma }}_0\times V,$$

where the variational forms $a$ and $L$ are defined as

$$\begin{array}{rl}a((\sigma ,u),(\tau ,v))& ={\int }_{\mathrm{\Omega }}\sigma \cdot \tau +\mathrm{\nabla }\cdot \tau u+\mathrm{\nabla }\cdot \sigma v\mathrm{d}x,\\ L((\tau ,v))& =-{\int }_{\mathrm{\Omega }}fv\mathrm{d}x+{\int }_{{\mathrm{\Gamma }}_D}{u}_0\tau \cdot n\mathrm{d}s,\end{array}$$

and ${\mathrm{\Sigma }}_g=\{\tau \in H(\mathrm{d}\mathrm{i}\mathrm{v})$ such that $\tau \cdot n{|}_{{\mathrm{\Gamma }}_N}=g\}$ and $V=L^2(\mathrm{\Omega })$.

To discretize the above formulation, two discrete function spaces ${\mathrm{\Sigma }}_h\subset \mathrm{\Sigma }$ and $V_h\subset V$ are needed to form a mixed function space ${\mathrm{\Sigma }}_h\times V_h$. A stable choice of finite element spaces is to let ${\mathrm{\Sigma }}_h$ be the Brezzi-Douglas-Fortin-Marini elements of polynomial order $k$ and let $V_h$ be discontinuous elements of polynomial order $k-1$.

We will use the same definitions of functions and boundaries as in the demo for the Poisson equation. These are:

• $\mathrm{\Omega }=[0,1]\times [0,1]$ (a unit square)

• ${\mathrm{\Gamma }}_D=\{(0,y)\cup (1,y)\in \partial \mathrm{\Omega }\}$

• ${\mathrm{\Gamma }}_N=\{(x,0)\cup (x,1)\in \partial \mathrm{\Omega }\}$

• $u_0=0$

• $g=\sin (5x)$ (flux)

• $f=10\exp (-((x-0.5{)}^2+(y-0.5{)}^2)/0.02)$ (source term)

Implementation

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx


fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)


def f1(x):
    values = np.zeros((2, x.shape[1]))
    values[1, :] = np.sin(5 * x[0])
    return values


f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))


facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)


def f2(x):
    values = np.zeros((2, x.shape[1]))
    values[1, :] = -np.sin(5 * x[0])
    return values


f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))


bcs = [bc_top, bc_bottom]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                      "pc_factor_mat_solver_type": "mumps"})
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, "out_mixed_poisson/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_h)