Second mixed formulation of the Stokes problem

Second mixed formulation of the Stokes problem#

Another possibility to obtained a well posed variational formulation, is to integrate by parts the \(\int_{\Omega} \nabla p \cdot \mathbf{v} ~\mathrm{d} \mathbf{x}\) term in the first formulation:

\[ \int_{\Omega} \nabla p \cdot \mathbf{v} ~\mathrm{d} \mathbf{x} = - \int_{\Omega} p \nabla \cdot \mathbf{v} ~\mathrm{d} \mathbf{x} + \int_{\partial\Omega} p \mathbf{v} \cdot \mathbf{n}~\mathrm{d} \sigma= -\int_{\Omega} p ~ \nabla \cdot \mathbf{v} ~\mathrm{d} \mathbf{x} , \]

using here \(p=0\) as a natural boundary condition. Note that in the other variational formulation the same boundary condition was essential. In this case, for the second variational formulation, we just multiply \(\nabla\cdot \mathbf{u}=0\) by \(q\) and integrate. No integration by parts is needed in this case.

\[ \int_{\Omega} \nabla \cdot \mathbf{u} ~ q ~\mathrm{d} \mathbf{x} =0. \]

This then leads to the following variational formulation:

Find \(( \mathbf{u},p)\in H^1(\Omega)^3\times L^2(\Omega)\) such that

\[\begin{split} \begin{align} \left\{ \begin{array}{llll} \int_{\Omega}\nabla \mathbf{u}:\nabla \mathbf{v} ~\mathrm{d} \mathbf{x} &- \int_{\Omega} p ~ \nabla\cdot \mathbf{v} ~\mathrm{d} \mathbf{x} &= \int_{\Omega} \mathbf{f}\cdot \mathbf{v} ~\mathrm{d} \mathbf{x}, &\forall \mathbf{v}\in H^1(\Omega)^3 \\ \int_{\Omega} \nabla\cdot\mathbf{u} ~ q~\mathrm{d} \mathbf{x} & &=0, &\forall q\in L^2(\Omega) \end{array} \right. \end{align} \end{split}\]

Formal Model#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import (ScalarFunctionSpace, VectorFunctionSpace, Square,
                             elements_of)
from sympde.calculus import grad, dot, div, inner
from sympde.core     import Constant

from psydac.api.discretization import discretize

from sympy import pi, sin, cos, Tuple

domain = Square()
x, y   = domain.coordinates

V1 = VectorFunctionSpace('V1', domain, kind='H1')
V2 = ScalarFunctionSpace('V2', domain, kind='L2')
X = V1*V2

# rhs
fx = -x**2*(x - 1)**2*(24*y - 12) - 4*y*(x**2 + 4*x*(x - 1) + (x - 1)**2)*(2*y**2 - 3*y + 1) - 2*pi*cos(2*pi*x)
fy = 4*x*(2*x**2 - 3*x + 1)*(y**2 + 4*y*(y - 1) + (y - 1)**2) + y**2*(24*x - 12)*(y - 1)**2 + 2*pi*cos(2*pi*y)
f  = Tuple(fx, fy)

# exact solution
ue1 = x**2*(-x + 1)**2*(4*y**3 - 6*y**2 + 2*y)
ue2 =-y**2*(-y + 1)**2*(4*x**3 - 6*x**2 + 2*x)
ue = Tuple(ue1, ue2)
pe = -sin(2*pi*x) + sin(2*pi*y)

u, v = elements_of(V1, names='u, v')
p, q = elements_of(V2, names='p, q')

# bilinear form
a  = BilinearForm(((u, p), (v, q)), integral(domain, inner(grad(u), grad(v)) - div(u)*q - p*div(v)) )

# linear form
l  = LinearForm((v, q), integral(domain, dot(f, v)))

# Dirichlet boundary conditions
bc = EssentialBC(u, 0, domain.boundary)

equation = find((u, p), forall=(v, q), lhs=a((u, p), (v, q)), rhs=l(v, q), bc=bc)

Discretization#

degree = [2,2]
ncells = [8,8]

# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Create discrete spline space
V1h = discretize(V1, domain_h, degree=degree)
V2h = discretize(V2, domain_h, degree=degree)
Xh  = discretize(X, domain_h, degree=degree)

# Discretize equation
equation_h = discretize(equation, domain_h, [Xh, Xh])

Solving the PDE#

equation_h.set_solver('gmres', info=False, tol=1e-8)

#uh, ph = equation_h.solve()

from psydac.fem.basic          import FemField

phi_h = equation_h.solve()

uh = FemField(V1h)
uh.coeffs[0][:] = phi_h.coeffs[0][:]
uh.coeffs[1][:] = phi_h.coeffs[1][:]

ph = FemField(V2h)
ph.coeffs[:] = phi_h.coeffs[2][:]

Computing the error norm#

Computing the \(L^2\) norm#

# L2 error norm of the velocity field
error_u   = [ue[0]-u[0], ue[1]-u[1]]
l2norm_u  = Norm(error_u, domain, kind='l2')
l2norm_uh = discretize(l2norm_u, domain_h, V1h)

# L2 error norm of the pressure, after removing the average value from the field
l2norm_p  = Norm(pe - p, domain, kind='l2')
l2norm_ph = discretize(l2norm_p, domain_h, V2h)

l2norm = l2norm_uh.assemble(u=uh)
print('>>> norm-l2 uh = ', l2norm)

l2norm = l2norm_ph.assemble(p=ph)
print('>>> norm-l2 ph = ', l2norm)

>>> norm-l2 uh =  1.860724404008467e-05
>>> norm-l2 ph =  0.007373764740051068
Contents