This page was generated from unit-2.5-mixed/mixed.ipynb.

2.5 Mixed formulation for second order equations

Motivation:

• exact flux conservation

• useful for a posteriori error estimates

• model problem for \(4^{th}\) order problems, Stokes, …

We consider the diffusion equation

\[\begin{split}\DeclareMathOperator{\Div}{div} \begin{array}{rcll} -\Div \lambda \nabla u & = & f & \text{ in } \Omega \\ u & = & u_D & \text{ on } \Gamma_D \\ \lambda \frac{\partial u}{\partial n} & = & g & \text{ on } \Gamma_N \end{array}\end{split}\]

Primal variational formulation

Find \(u \in H^1, u = u_D\) on \(\Gamma_D\) such that

\[\int_\Omega \lambda \nabla u \nabla v = \int_\Omega f v + \int_{\Gamma_N} g v \quad \forall v \in H^1, v = 0 \text{ on } \Gamma_D\]

First order system

Find the scalar \(u\) and the flux \(\sigma\) such that

\[\lambda^{-1} \sigma = \nabla u, \quad \Div \sigma = -f\]

with boundary conditions

\[\sigma \cdot n = g \text{ on } \Gamma_N, \quad \text{ and } \quad u = u_D \text{ on } \Gamma_D\]

Mixed variational formulation

Find \(\sigma \in H(\Div)\) and \(u \in L_2\) such that \(\sigma_n = g\) on \(\Gamma_N\) and

\begin{eqnarray*} \int_\Omega \lambda^{-1} \sigma \tau + \int_\Omega \Div \tau u & = & \int_{\Gamma_D} u_D \tau_n \\ \int_\Omega \Div \sigma v & = & -\int_\Omega f v \end{eqnarray*}

for all test functions \(\tau \in H(\Div)\) and \(v \in L_2\) with \(\tau_n = 0\) on \(\Gamma_N\).

Here \(\sigma_n\) is the normal trace \(\sigma \cdot n\) on the boundary, which is meaningful in \(H(\Div)\). The roles of the boundary conditions have changed: The Neumann condition is now an essential boundary condition for the flux, while the Dirichlet condition has become a natural one, and enters the right hand side.

A compact notation is the single-liner

Find \((\sigma, u) \in H(\Div) \times L_2\) such that \(\sigma_n = g\) on \(\Gamma_N\) and

\[\int_\Omega \lambda^{-1} \sigma \tau + \Div \sigma v + \Div \tau u = -\int_\Omega f v + \int_{\Gamma_D} u_D \tau_n\]

for all test functions \((\tau, v) \in H(\Div) \times L_2\) with \(\tau_n = 0\) on \(\Gamma_N\).

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

source = sin(pi*x)
ud = CF(5)
g = mesh.BoundaryCF( {"left" : y*(1-y)}, default=0)
lam = 10

Set up and solve the primal problem:

fesp = H1(mesh, order=4)
up, vp = fesp.TnT()

gfup = Solve(lam*grad(up)*grad(vp)*dx == source*vp*dx + g*vp*ds,
             up[BND("bottom")] == ud)
Draw (gfup)
Draw (lam * grad(gfup), mesh, "flux-primal");

Solving the mixed problem

Define the spaces for the mixed problem. Multiplying spaces generates a product space. The flag lowest_order_wb is explained below.

order_flux=1
V = HDiv(mesh, order=order_flux)
Q = L2(mesh, order=order_flux-1, lowest_order_wb=True)
fesm = V*Q

The product space provides a tuple of trial functions, and a tuple of test functions:

(sigma, u), (tau, v) = fesm.TnT()

The normal vector is provided as a special coefficient function (which may be used only at the boundary). The orientation depends on the definition of the geometry. In 2D, it is the tangential vector rotated to the right, which is the outer normal vector in our case. Since every CoefficientFunction must know its vector dimension, we have to provide the space dimension:

normal = specialcf.normal(mesh.dim)
print (normal)

coef normal vector, real, dim=2

Define the variational equation on the product space. The proxy functions know to which component of the product space they belong. For the boundary term we have to use the Trace operator, which restricts the \(H(\Div)\) function to its normal component on the boundary.

lhs = (1/lam*sigma*tau + div(sigma)*v + div(tau)*u)*dx
rhs = -source*v*dx + ud*(tau.Trace()*normal)*ds("bottom")

The essential boundary conditions are now set for the flux, on the opposite part of the boundary.

Solve eliminates element-local unknowns by static condensation. All unknowns of the scalar \(u\) are element-local, but the piecewise constants cannot be eliminated on the element level, since they are coupled only to the normal fluxes on the element boundary. Thus we have declared the lowest order functions of the \(L_2\) space as wirebasket dofs by the flag lowest_order_wb, and they stay in the global system. The remaining saddle point problem is solved by a sparse direct solver.

gfm = Solve(lhs==rhs, sigma[BND("right|top|left")] == g*normal)

The solution lives in the product space. We access its components via the trial functions (or as gfm.components). They don’t have their own coefficient vectors, but refer to ranges of the big coefficient vector.

gfsigma, gfu = gfm[sigma], gfm[u]
Draw (gfsigma, mesh, "flux-mixed")
Draw (gfu, mesh, "u-mixed");

Calculate the difference to the primal solution:

print ("err-u:   ", sqrt(Integrate( (gfup-gfu)**2, mesh)))
errflux = lam * grad(gfup) - gfsigma
print ("err-flux:", sqrt(Integrate(errflux*errflux, mesh)))

err-u:    0.0010281098221311617
err-flux: 0.00165350563604929

Post-processing for the scalar variable

The scalar variable is approximated one order lower than the flux, which is (up to the factor \(\lambda\)) its gradient. Knowing the gradient of a function more accurately, and knowing its mean value, one can recover the function itself. For this post-processing trick we refer to [Arnold+Brezzi 85].

Find \(\widehat u \in P^{k+1, dc}\) and \(\widehat \lambda \in P^{0, dc}\) such that

\[\begin{split}\begin{array}{ccccl} \int \lambda \nabla \widehat u \nabla \widehat v & + & \int \widehat \lambda \widehat v & = & \int \sigma \nabla \widehat v & \forall \, \widehat v\\ \int \widehat u \widehat \mu & & & = & \int u \widehat \mu & \forall \, \widehat \mu \end{array}\end{split}\]

All unknowns are element-local, and the local problems are regular. Here Solve computes everything element by element.

fespost_u = L2(mesh, order=order_flux+1)
fespost_lam = L2(mesh, order=0)
fes_post = fespost_u*fespost_lam

(uh,lah),(vh,muh) = fes_post.TnT()

gfpost = Solve((lam*grad(uh)*grad(vh)+lah*vh+muh*uh)*dx == (gfsigma*grad(vh)+gfu*muh)*dx)

Draw (gfpost[uh], mesh, "upost")
print ("err-upost:   ", sqrt(Integrate( (gfup-gfpost[uh])**2, mesh)))

err-upost:    nan