Taylor-Hood elements for Darcy's equation

Question

I would like to know if Taylor-Hood elements $P_2$-$P_1$ form a stable pair for the mixed approximation of Darcy's equation ( or Poisson's equation) with Dirichlet B.C. In the literature I only find them only in the context of the Stokes problem.

Answer 1

I found some evidence that this method is not optimal, see the introductions in the following papers:

• He, Mingyan; Sun, Pengtao, Error analysis of mixed finite element method for Poisson-Nernst-Planck system, Numer. Methods Partial Differ. Equations 33, No. 6, 1924-1948 (2017). ZBL1383.78029.

• Brezzi, Franco; Fortin, Michel; Marini, L. Donatella, Mixed finite element methods with continuous stresses, Math. Models Methods Appl. Sci. 3, No. 2, 275-287 (1993). ZBL0774.73066.

I tried doing some experimental numerics to check the convergence rates. I'm working with the mixed formulation $$ \begin{aligned} (\sigma, \tau) + (u, \mathrm{div}\, \tau) & = 0 \\ (\mathrm{div}\, \sigma, v) &=-(f,v) \end{aligned} $$ and $f=2 \pi^2 \sin(\pi x) \sin(\pi y)$. Then I'm evaluating $\| \sigma - \sigma_h \|_{0}$ for different elements as a function of the mesh parameter $h$.

Results for $RT_0 - P_0$:

Results for $BDM_1 - P_0$:

Results for $[P_2]^2 - P_1$:

While you would expect $\| \sigma - \sigma_h \|_0 = O(h^2)$ for such basis functions, I'm instead observing $\| \sigma - \sigma_h \|_0 = O(h^{1.5})$ which might be caused by the missing (uniform?) compatibility condition. Despite this, I'm observing a quadratic convergence rate for $\|u - u_h\|_0$. However, I believe the energy norm is $\| \sigma - \sigma_h \|_{div} + \| u - u_h\|_0$ and for the first term I'm observing $\| \sigma - \sigma_h\|_{div} = O(h^{0.5})$ for $[P_2]^2-P_1$ while for the other elements ($RT_0$ and $BDM_1$) I'm observing $O(h)$.

So I guess it depends on what you mean by "stable". There is inf-sup condition and there is ellipticity on the kernel. To get optimal rates you need these to hold uniformly with constants independent of the mesh. I think this is not the case here since we observe suboptimal rates.

Edit: Here is the source code after pip install scikit-fem[all]==6.0.0:

import numpy as np
import matplotlib.pyplot as plt
from skfem import *
from skfem.helpers import dot, div

hs, errors = [], []

for refs in [3, 4, 5]:
  m = MeshTri.init_sqsymmetric().refined(refs)
  e = ElementVector(ElementTriP2()) * ElementTriP1()
  #e = ElementTriRT0() * ElementTriP0()
  #e = ElementTriBDM1() * ElementTriP0()
  basis = Basis(m, e)

  @BilinearForm
  def bilinf(sigma, u, tau, v, w):
      return dot(sigma, tau) + div(sigma) * v + div(tau) * u

  @LinearForm
  def linf(tau, v, w):
      return - 2* np.pi ** 2 * np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1]) *v

  A = asm(bilinf, basis)
  b = asm(linf, basis)
  x = solve(A, b)
  (sigma, sigmabasis), (u, ubasis) = basis.split(x)

  @Functional
  def error(w):
    return (w['u'] - np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1])) ** 2

  @Functional
  def derror(w):
    return (div(w['sigma']) + 2 * np.pi ** 2 * np.sin(np.pi * w.x[0]) * np.sin(np.pi * w.x[1])) ** 2

  hs.append(m.param())
  #errors.append(error.assemble(ubasis, u=u))
  errors.append(derror.assemble(sigmabasis, sigma=sigma))

errors = np.sqrt(np.array(errors))
plt.title(str(e.elems[0].__class__))
plt.loglog(hs, errors, 'bo-')
plt.legend(['rate = {:.2f}'.format(np.polyfit(np.log(hs), np.log(errors), 1)[0])])

Answer 2

Yes, they are. I don't remember a reference to the literature, but the proof goes along the same line as for the Stokes case.