Mixed Methods for the Diffusion equations

5. Mixed Methods for the Diffusion equations#

We consider the mixed system

\[\begin{eqnarray*} a \sigma + \nabla u & = & 0 \\ -\operatorname{div} \sigma & = & f \end{eqnarray*}\]

with boundary conditions

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

5.1. Primal and Dual mixed methods#

5.1.1. Primal mixed method#

Find \(\sigma \in \Sigma := L_2(\Omega)^d\) and \(u \in V := H^1\) s.t. \(u = u_D\) on \(\Gamma_D\) and

\[\begin{split} \begin{array}{ccccll} (a \sigma, \tau)_{L_2} & + & (\nabla u, \tau)_{L_2} & = & 0 & \forall \, \tau \in \Sigma \\ (\nabla v, \sigma)_{L_2} & & & = & \int_\Omega f v + \int_{\Gamma_N} g v & \forall \, v \in V, v = 0 \text{ on } \Gamma_D \end{array} \end{split}\]

5.1.2. Dual mixed method#

Find \(\sigma \in \Sigma := H(\operatorname{div})\) and \(u \in L_2\) s.t. \(\sigma\cdot n = g\) on \(\Gamma_N\) and

\[\begin{split} \begin{array}{ccccll} (a \sigma, \tau)_{L_2} & - & (u, \operatorname{div} \tau)_{L_2} & = & \int_{\Gamma_D} u_D \tau_n & \forall \, \tau \in \Sigma \\ -(v, \operatorname{div} \sigma)_{L_2} & & & = & \int_\Omega f v & \forall \, v \in H_{0,D}^1 \end{array} \end{split}\]

5.1.3. Balkenwaage:#

Two different formulations for the same equation? We shift regularity between primal and dual spaces

Different choices of spaces for the same biform

5.2. Proving Brezzi conditions#

5.2.1. Primal mixed method:#

\(V = H^1, \Sigma = [L_2]^d\)

\(a(.,.)\) is the inner product in \(\Sigma\), thus (kernel) coercivity holds with constant 1.
LBB condition:

\[ \inf_{u \in H^1} \sup_{\sigma \in [L_2]^d} \frac{(\sigma, \nabla u)_{L_2}}{ \| \sigma \|_{L_2} \| \nabla u \|_{L_2} } \; \underset{\sigma := \nabla u} {=} \; \inf_{u \in H^1} \frac{ \| \nabla u \|} { \| \nabla u \| } = 1 \]

The flux is chosen as the gradient of \(u\).

5.2.2. Dual mixed method:#

\(V = L_2, \Sigma = H(\operatorname{div})\), with \(\| \sigma \|_{H(\operatorname{div})}^2 = \| \sigma \|_{L_2}^2 + \| \operatorname{div} \sigma \|_{L_2}^2\).

On \(V_{0h} = \{ \sigma \in \Sigma : \operatorname{div} \sigma = 0 \}\), we have

\[a(\sigma, \sigma) = (\sigma, \sigma)_{L_2} = (\sigma, \sigma)_{L_2} + (\operatorname{div} \sigma, \operatorname{div}\sigma)_{L_2} = \| \sigma \|_{H(\operatorname{div})}^2 \]
LBB condition:

Now we have to find a right-inverse to the \(\operatorname{div}\). We can reduce it to solving a Poisson equation:

\[ \sigma := \nabla w \quad \text{with} \quad \Delta w = u \]

Then \(\| w \|_{H^1} \preceq \| u \|_{H^-1} \preceq \| u \|_{L_2}\), and thus \(\| \sigma \|_{L_2} = \| \nabla w \|_{L_2} \preceq \| u \|_{L_2}\)

And, \(\operatorname{div} \sigma = \Delta w = u\), and thus also \(\| \operatorname{div} \sigma \|_{L_2} \preceq \| u \|_{L_2}\).

We have found a \(\sigma\) such that

\[ \inf_{u \in V} \sup_{\tau \in \Sigma} \frac{ (\operatorname{div} \tau, u) }{ \| \tau \|_{\Sigma} \| u \|_V } \geq \inf_{u \in V} \frac{ (\operatorname{div} \sigma, u) }{ \| \sigma \|_{\Sigma} \| u \|_V } \succeq \frac { \| u \|_{L_2}^2 } { \| u \|_{L_2}^2 } = 1 \]

5.3. Finite elements for \(H^1\) and \(H(\operatorname{div})\)#

For \(H^1\), we choose piecewise polynomials, which are continuous across elements. Thus, these functions have a weak gradient, which is exactly the element-wise gradient.
For \(H(\operatorname{div})\), we choose vectorial piecewise polynomials, whose normal components are continuous across elements. Thus, these functions have a weak divergence, which is exactly the element-wise divergence (Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements).

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.25))

fes = H1(mesh, order=1)
gfu = GridFunction(fes, multidim=fes.ndof)
for i in range(fes.ndof):
    gfu.vecs[i][i] = 1
Draw (gfu, animate=True, deformation=True, scale=0.3, euler_angles=(-65,0,-170));

fes = HDiv(mesh, order=0)  # use RT=True for higher order
gfu = GridFunction(fes, multidim=fes.ndof)
for i in range(fes.ndof):
    gfu.vecs[i][i] = 1
Draw (gfu, animate=True, vectors={"grid_size" : 40 }, scale=0.3);

5.4. Proving discrete LBB conditions#

5.4.1. Primal mixed#

Same as continuous one, choose \(\sigma_h = \nabla u_h\).

5.4.2. Dual mixed#

Use a \(H(\operatorname{div})\)-continuous, commuting interpolation operator

\[ \sigma_h = \Pi_h \sigma, \]

where \(\sigma\) is a candidate for the infinite-dimensional LBB condition. Then

\[ \frac{(\operatorname{div} \sigma_h, u_h)}{ \| \sigma_h \|_{\Sigma} \| u_h \|_{V_h}} \succeq \frac{(\operatorname{div} \sigma, u_h)}{ \| \sigma \|_{\Sigma} \| u_h \|_{V_h}} \geq \beta \]

5.4.3. Dual mixed method with primal mixed analysis#

Discrete norms:

\[\begin{eqnarray*} \| \sigma_h \|_{\Sigma_h} & = & \| \Sigma \|_{L_2} \\[0.5em] \| u_h \|_{V_h}^2 & = & \sum_{\text{cells } T} \| \nabla u_h \|_{L_2(T)}^2 + \sum_{\text{facets } F} \tfrac{1}{h} \| [u_h] \|_{L_2(F)}^2 \end{eqnarray*}\]

Candidate for LBB: Given \(u_h \in V_h\), choose \(\sigma_h \in \Sigma_h\) such that

\[\begin{eqnarray*} \sigma_h \cdot n & = & \tfrac{1}{h} [u_h] \qquad \forall F \\ \int_T \sigma_h \cdot q_h & = & \int_T \nabla u_h \cdot q_h \qquad \forall \, T \forall q_h \in [P^{k'}]^d \end{eqnarray*}\]

proves the discrete \(\inf-\sup\) condition.

5.5. Super-convergence#

For the dual mixed discretizaton we have

\[ \operatorname{div} \Sigma_h = V_h \]

and thus the kernel-inclusion property.

From the discrete \(\inf-\sup\) condition we get:

\[ \| I_h \sigma - \sigma_h \|_{\Sigma_h} + \| P_h u - u_h \|_{V_h} \preceq \sup_{\tau_h, v_h} \frac{a(\sigma - I_h \sigma, \tau_h) + b(\sigma - I_h \sigma, q_h) + b(\tau_h, u - P_h u)}{\| \tau_h \|_{\Sigma} + \| v_h \|_{V_h}} \]

Using the commuting diagram property we get:

\[ \int (\sigma - I_h \sigma) \tau_h + \int \underbrace{(I - P_h) \operatorname{div} \sigma}_{\in V_h^\bot} \;\underbrace{ q_h}_{\in V_h} + \int \underbrace{ \operatorname{div} \tau_h}_{\in V_h} \, \underbrace{ (u-P_h u) }_{\in V_h^\bot} \]

and thus

\[ \| I_h \sigma - \sigma_h \|_{\Sigma_h} + \| I_h u - u_h \|_{V_h} \preceq \frac{a(\sigma - I_h \sigma, \tau_h)}{\| \tau_h \|_{\Sigma_h}} \leq \| \sigma - I_h \sigma \|_{L_2} \]

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

order=1
Sigma = HDiv(mesh,order=order)
V = L2(mesh, order=order-1)
X = Sigma * V

sigma,u = X.TrialFunction()
tau,v = X.TestFunction()

a = BilinearForm(X)
a += (sigma*tau + div(sigma)*v + div(tau)*u)*dx
a.Assemble()

f = LinearForm(10*v*dx).Assemble()

gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse() * f.vec

sol_sigma = gfu.components[0]
sol_u = gfu.components[1]
print ("The flux")
Draw (sol_sigma, mesh, "sigma");
print ("The scalar")
Draw (sol_u, mesh, "u");

The flux

The scalar

5.6. Local post-processing#

Since \(\lambda \nabla u = \sigma\), we can reconstruct a better approximation \(\widetilde u\) by small, element-wise problems:

\[ \widetilde u = \operatorname{arg}\min_{v_h \in P^{k+1} \atop \int_T v_h = \int_T u_h} \| \lambda \nabla v_h - \sigma \|_{L_2, \lambda^{-1}}^2 \]

This optimization problems can be written as a mixed variational problem:

Find: \(\widetilde u \in P^{k+1,dc}\) and \(p_h \in P^0\):

\[\begin{split} \begin{array}{ccccll} \sum_T \int_T \lambda \nabla \widetilde u \nabla \widetilde v & + & \int_{\Omega} \widetilde v_h p_h & = & \sum_T \int_T \sigma_h \nabla \widetilde v_h & \forall \, \widetilde v_h \\ \int_{\Omega} \widetilde u_h q_h & & & = & \int_{\Omega} u_h q_h & \forall q_h \end{array} \end{split}\]

V2 = L2(mesh, order=order+1)
Q2 = L2(mesh, order=0)
X2 = V2*Q2

u2,p2 = X2.TrialFunction()
v2,q2 = X2.TestFunction()
a2 = BilinearForm(X2)
a2 += (grad(u2)*grad(v2) + u2*q2 + v2*p2)*dx
f2 = LinearForm(X2)
f2 += (sol_sigma*grad(v2) + sol_u * q2)*dx

a2.Assemble()
f2.Assemble()
gfu2 = GridFunction(X2)
gfu2.vec.data = a2.mat.Inverse() * f2.vec

Draw(gfu2.components[0], mesh, "upost");

5.7. Hybridization Techniques#

Discretizing an elliptic equation by a primal method leads to a positive definite matrix. Using a mixed method leads to a saddle point problem. This disadvantage can be overcome by the so called hybridization technique:

The idea is as follows: Break the normal-continuity of the RT - space, and reinforce it by another Lagrange parameter. This new Lagrange parameter is a polynomial on every mesh edge (or face, in 3D).

The variational formulation is as follows:

Find: \(\sigma_h \in RT_k^{dc}, u_h \in P^k, \widehat u_h \in P^k(\cup E)\) such that

\[\begin{split} \begin{array}{ccccccl} \int_\Omega \sigma_h \tau_h & + & \sum_T \int_T \operatorname{div} \tau_h \, u_h & - & \sum_E \int_E [\![\tau_{h,n}]\!] \widehat u_h & = & 0 & \forall \, \tau_h \\ \sum_T \int_T \operatorname{div} \sigma_h \, v_h & & & & & = & \int f v_h & \forall \, v_h \\ -\sum_E \int_E [\![\sigma_{h,n}]\!] \widehat v_h & & & & & = & 0 & \forall \, \widehat v_h \end{array} \end{split}\]

We used the definition of the jump term \([\![ \tau_n ]\!] = \tau_1 n_1 + \tau_2 n_2\)

This formulation gives the same solution as the mixed formulation, so we don’t need an extra error analysis

The physical meaning of the Lagrange paramter \(\widehat u_h\) is the primal variable, what can be seen by integration by parts of the first equation.

Now, Dirichlet boundary conditions are set by constraining the \(\widehat u\) on the Dirichlet boundary, and Neumann boundary conditions are formulated by \(\int_{\Gamma_N} g \widehat v\). The hybridized formulation is thus similar to a primal method.

The discretization system has the form

\[\begin{split} \left( \begin{array}{ccc} A & B_1^T & B_2^T \\ B_1 & & \\ B_2 & & \end{array} \right) \left( \begin{array}{c} \sigma \\ u_1 \\ u_2 \end{array} \right) = \left( \begin{array}{c} 0 \\ f \\ 0 \end{array} \right) \end{split}\]

The submatrix

\[\begin{split} \left( \begin{array}{cc} A & B_1^T \\ B_1 & 0 \end{array} \right) \end{split}\]

is regular and block-diagonal, each block corresponds to an element. Thus it can be cheaply inverted.

If \(u_2\) were known, we could compute the first two variables from \(u_2\):

\[\begin{split} \left( \begin{array}{c} \sigma \\ u_1 \end{array} \right) = \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left[ \left( \begin{array}{c} 0 \\ f \end{array} \right) - \left( \begin{array}{c} B_2^T \\ 0 \end{array} \right) u_2 \right] \end{split}\]

Plugging this term into the third equation \(B_2 \sigma = 0\), we obtain the system

\[\begin{split} \left( B_2 \; \; 0 \right) \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left( \begin{array}{c} B_2^T \\ 0 \end{array} \right) \; \; u_2 = \left( B_2 \; \; 0 \right) \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left( \begin{array}{c} 0 \\ f \end{array} \right) \end{split}\]

The matrix on the left hand side is symmetric positive definite. It behaves like a standard system matrix (e.g. for condition number), and standard iterative methods and preconditioners can be used for solution.

The lowest order hybrid method has the same degrees of freedom as the non-conforming \(P^1\) element. Here, the extension to higher order is straight forward.

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

order=2
Sigma = Discontinuous (HDiv(mesh, order=order, RT=True))
V = L2(mesh, order=order)
Vhat = FacetFESpace(mesh, order=order, dirichlet="left|bottom")

X = Sigma*V*Vhat

sigma,u,uhat = X.TrialFunction()
tau,v,vhat = X.TestFunction()

n = specialcf.normal(mesh.dim)

a = BilinearForm(X, eliminate_internal = True)
a += (sigma*tau + div(sigma)*v + div(tau)*u)*dx
a += (-sigma*n*vhat-tau*n*uhat) * dx(element_boundary=True)

c = Preconditioner(a, "bddc")
a.Assemble()

f = LinearForm(X)
f += -v*dx
f.Assemble()

gfu = GridFunction(X)

f.vec.data += a.harmonic_extension_trans * f.vec

# gfu.vec.data = a.mat.Inverse(X.FreeDofs(True)) * f.vec
from ngsolve.solvers import CG
CG (mat=a.mat, pre=c.mat, rhs=f.vec, sol=gfu.vec, 
    printrates=True, maxsteps=200)

gfu.vec.data += a.harmonic_extension * gfu.vec
gfu.vec.data += a.inner_solve * f.vec

Draw (gfu.components[0], mesh, "sigma")
Draw (gfu.components[1], mesh, "u");

CG iteration 1, residual = 0.37591294565277006     
CG iteration 2, residual = 0.0757921849074072     
CG iteration 3, residual = 0.013587579030952717     
CG iteration 4, residual = 0.0044402342217845295     
CG iteration 5, residual = 0.0016170166562525501     
CG iteration 6, residual = 0.0006636169222860922     
CG iteration 7, residual = 0.0001874144178565463     
CG iteration 8, residual = 6.0290886204633915e-05     
CG iteration 9, residual = 2.0777211605223626e-05     
CG iteration 10, residual = 6.153358671770975e-06     
CG iteration 11, residual = 2.0768303974751285e-06     
CG iteration 12, residual = 5.592765774424452e-07     
CG iteration 13, residual = 2.1404260555987917e-07     
CG iteration 14, residual = 5.581193089556126e-08     
CG iteration 15, residual = 2.1521755208779676e-08     
CG iteration 16, residual = 6.3948425814477085e-09     
CG iteration 17, residual = 2.257408527935346e-09     
CG iteration 18, residual = 7.646470690282719e-10     
CG iteration 19, residual = 2.296191954510052e-10     
CG iteration 20, residual = 8.967090603271594e-11     
CG iteration 21, residual = 2.1711429021155548e-11     
CG iteration 22, residual = 7.154168881453398e-12     
CG iteration 23, residual = 1.8130000866380565e-12     
CG iteration 24, residual = 6.554662617866787e-13     
CG iteration 25, residual = 1.9067207689610788e-13     

5.8. Does it make a difference?#

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

r1 = MoveTo(0,0).RectangleC(2,2).Face()
r1.edges.name = "outer"
r2 = MoveTo(0,0).RectangleC(1.6,0.1).Face()
outer = r1-r2
outer.faces.name="outer"
r2.faces.name="slit"
r2.edges.name="interface"
shape = Glue([outer,r2])
mesh = shape.GenerateMesh(maxh=0.1, dim=2)
Draw (mesh);

Exercise: Replace the thin domain by a curve, and make a lower dimensional model. What happens for small or large conductivities in the slit?

5.8.1. Domains with holes:#

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
shape = Circle((0,0,),1).Face()-Circle((0,0,),0.6).Face()
mesh = shape.GenerateMesh(maxh=0.1, dim=2)
Draw (mesh);

fes = HDiv(mesh, order=0, dirichlet=".*") * L2(mesh, order=0) 
(sigma,u),(tau,v) = fes.TnT()
a = BilinearForm(sigma*tau*dx + div(sigma)*v*dx + div(tau)*u*dx - u*v*dx).Assemble()
f = LinearForm( CF((y,-x))*tau*dx).Assemble()
gf = GridFunction(fes)
gf.vec[:] = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
gfsigma,gfu = gf.components
Draw (gfsigma, vectors=True);

fes = Discontinuous(HDiv(mesh, order=0)) * L2(mesh, order=0) * FacetFESpace(mesh, order=0)
(sigma,u,uhat),(tau,v,vhat) = fes.TnT()
n = specialcf.normal(2)
a = BilinearForm(sigma*tau*dx - sigma*grad(v)*dx-tau*grad(u)*dx + (sigma*n*(v-vhat)+tau*n*(u-uhat))*dx(element_boundary=True)-u*v*dx).Assemble()
f = LinearForm( CF((y,-x))*tau*dx).Assemble()
gf = GridFunction(fes)
gf.vec[:] = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
gfsigma,gfu,gfuhat = gf.components
Draw (gfsigma, vectors=True);

fes = H1(mesh, order=1)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx+1e-6*u*v*dx).Assemble()
f = LinearForm(CF((y,-x))*grad(v)*dx)
gfu = GridFunction(fes)
gfu.vec[:] = a.mat.Inverse()*f.vec
Draw (CF((y,-x))-grad(gfu), mesh);