First mixed formulation of the Stokes problem

First mixed formulation of the Stokes problem#

For the variational formulation, we take the dot product of the first equation with \(v\) and integrate over the whole domain

\[ \int_{\Omega} (-\Delta \mathbf{u} + \nabla p)\cdot \mathbf{v} ~\mathrm{d} \mathbf{x} = \int_{\Omega}\nabla \mathbf{u}:\nabla \mathbf{v} ~\mathrm{d} \mathbf{x} + \int_{\Omega} \nabla p \cdot \mathbf{v} ~\mathrm{d} \mathbf{x} = \int_{\Omega} \mathbf{f}\cdot \mathbf{v} ~\mathrm{d} \mathbf{x} \]

The integration by parts is performed component by component. We impose the essential boundary condition \(\mathbf{v}=0\) on \(\partial\Omega\), and we denote by

\[ \int_{\Omega}\nabla \mathbf{u}:\nabla \mathbf{v} ~\mathrm{d} \mathbf{x} =\sum_{i=1}^3 \int_{\Omega}\nabla u_i\cdot\nabla v_i ~\mathrm{d} \mathbf{x} =\sum_{i,j=1}^3 \int_{\Omega}\partial_j u_i\partial_j v_i ~\mathrm{d} \mathbf{x} . \]

We now need to deal with the constraint \(\nabla\cdot \mathbf{u}=0\). The theoretical framework for saddle point problems requires that the corresponding bilinear form is the same as the second one appearing in the first part of the variational formulation. To this aim we multiply \(\nabla\cdot \mathbf{u}=0\) by a scalar test function (which will be associated to \(p\)) and integrate on the full domain, with an integration by parts in order to get the same bilinear form as in the first equation

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

using that \(q=0\) on the boundary as an essential boundary condition. We finally obtain the mixed variational formulation:

Find \(( \mathbf{u},p)\in H^1_0(\Omega)^3\times H^1_0(\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} \nabla p \cdot \mathbf{v} ~\mathrm{d} \mathbf{x} &= \int_{\Omega} \mathbf{f}\cdot \mathbf{v} ~\mathrm{d} \mathbf{x}, & \forall \mathbf{v}\in H_0^1(\Omega)^3 \\ \int_{\Omega} \mathbf{u}\cdot \nabla q~\mathrm{d} \mathbf{x} & &=0, & \forall p\in H^1_0(\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='H1')
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)) + dot(grad(p), v) + dot(u, grad(q)) ) )

# 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 =  5.852908166372369e-05
>>> norm-l2 ph =  0.01411215698700372
Contents