Mixed formulation for the Poisson equation

We’re considering the Poisson equation \(\nabla^2 u = -f\) using a mixed formulation on two coupled fields. We start by introducing the negative flux \(\sigma = \nabla u\) as an auxiliary vector-valued variable. This leaves us with the PDE on a unit square \(\Omega = [0,1] \times [0,1]\) with boundary \(\Gamma\)

\[ \begin{align}\begin{aligned}\sigma - \nabla u = 0 \ \textrm{on}\ \Omega\\\nabla \cdot \sigma = -f \ \textrm{on}\ \Omega\\u = u_0 \ \textrm{on}\ \Gamma_D\\\sigma \cdot n = g \ \textrm{on}\ \Gamma_N\end{aligned}\end{align} \]

for some known function \(f\). The solution to this equation will be some functions \(u\in V\) and \(\sigma\in \Sigma\) for some suitable function space \(V\) and \(\Sigma\) that satisfy these equations. We multiply by arbitrary test functions \(\tau \in \Sigma\) and \(\nu \in V\), integrate over the domain and then integrate by parts to obtain a weak formulation of the variational problem: find \(\sigma\in \Sigma\) and \(\nu\in V\) such that:

\[ \begin{align}\begin{aligned}\begin{split}\int_{\Omega} (\sigma \cdot \tau + \nabla \cdot \tau \ u) \ {\rm d} x &= \int_{\Gamma} \tau \cdot n \ u \ {\rm d} s \quad \forall \ \tau \in \Sigma, \\\end{split}\\\int_{\Omega} \nabla \cdot \sigma v \ {\rm d} x &= - \int_{\Omega} f \ v \ {\rm d} x \quad \forall \ v \in V.\end{aligned}\end{align} \]

The flux boundary condition \(\sigma \cdot n = g\) becomes an essential boundary condition to be enforced on the function space, while the boundary condition \(u = u_0\) turn into a natural boundary condition which enters into the variational form, such that the variational problem can be written as: find \((\sigma, u)\in \Sigma_g \times V\) such that

\[a((\sigma, u), (\tau, v)) = L((\tau, v)) \quad \forall \ (\tau, v) \in \Sigma_0 \times V\]

with the variational forms \(a\) and \(L\) defined as

\[\begin{split}a((\sigma, u), (\tau, v)) &= \int_{\Omega} \sigma \cdot \tau + \nabla \cdot \tau \ u + \nabla \cdot \sigma \ v \ {\rm d} x \\ L((\tau, v)) &= - \int_{\Omega} f v \ {\rm d} x + \int_{\Gamma_D} u_0 \tau \cdot n \ {\rm d} s\end{split}\]

The essential boundary condition is reflected in function spaces \(\Sigma_g = \{ \tau \in H({\rm div}) \text{ such that } \tau \cdot n|_{\Gamma_N} = g \}\) and \(V = L^2(\Omega)\).

We need to choose a stable combination of discrete function spaces \(\Sigma_h \subset \Sigma\) and \(V_h \subset V\) to form a mixed function space \(\Sigma_h \times V_h\). One such choice is Brezzi-Douglas-Marini elements of polynomial order \(k\) for \(\Sigma_h\) and discontinuous elements of polynomial order \(k-1\) for \(V_h\).

For the remaining functions and boundaries we choose:

\[ \begin{align}\begin{aligned}\Gamma_{D} = \{(0, y) \cup (1, y) \in \partial \Omega\}, \Gamma_{N} = \{(x, 0) \cup (x, 1) \in \partial \Omega\}\\u_0 = 0, g = \sin(5x)\\f = 10~e^{-\frac{(x - 0.5)^2 + (y - 0.5)^2}{0.02}}\end{aligned}\end{align} \]

To produce a numerical solution to this PDE in Firedrake we procede as follows:

The mesh is chosen as a \(32\times32\) element unit square.

from firedrake import *
mesh = UnitSquareMesh(32, 32)

As argued above, a stable choice of function spaces for our problem is the combination of order \(k\) Brezzi-Douglas-Marini (BDM) elements and order \(k - 1\) discontinuous Galerkin elements (DG). We use \(k = 1\) and combine the BDM and DG spaces into a mixed function space W.

BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG

We obtain test and trial functions on the subspaces of the mixed function spaces as follows:

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

Next we declare our source function f over the DG space and initialise it with our chosen right hand side function value.

x, y = SpatialCoordinate(mesh)
f = Function(DG).interpolate(
    10*exp(-(pow(x - 0.5, 2) + pow(y - 0.5, 2)) / 0.02))

After dropping the vanishing boundary term on the right hand side, the bilinear and linear forms of the variational problem are defined as:

a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx

The strongly enforced boundary conditions on the BDM space on the top and bottom of the domain are declared as:

bc0 = DirichletBC(W.sub(0), as_vector([0.0, -sin(5*x)]), 1)
bc1 = DirichletBC(W.sub(0), as_vector([0.0, sin(5*y)]), 2)

Note that it is necessary to apply these boundary conditions to the first subspace of the mixed function space using W.sub(0). This way the association with the mixed space is preserved. Declaring it on the BDM space directly is not the same and would in fact cause the application of the boundary condition during the later solve to fail.

Now we’re ready to solve the variational problem. We define w to be a function to hold the solution on the mixed space.

w = Function(W)

Then we solve the linear variational problem a == L for w under the given boundary conditions bc0 and bc1. Afterwards we extract the components sigma and u on each of the subspaces with split.

solve(a == L, w, bcs=[bc0, bc1])
sigma, u = w.split()

Lastly we write the component of the solution corresponding to the primal variable on the DG space to a file in VTK format for later inspection with a visualisation tool such as ParaView


We could use the built in plot function of firedrake by calling plot to plot a surface graph. Before that, matplotlib.pyplot should be installed and imported:

  import matplotlib.pyplot as plt
  warning("Matplotlib not imported")

except Exception as e:
  warning("Cannot plot figure. Error msg '%s'" % e.message)

Don’t forget to show the image:

except Exception as e:
  warning("Cannot show figure. Error msg '%s'" % e.message)

This demo is based on the corresponding DOLFIN mixed Poisson demo and can be found as a script in poisson_mixed.py.