First mixed formulation of the Poisson problem

First mixed formulation of the Poisson problem#

Instead of having one unknown, we now have two, along with the above two equations. In order to get a mixed variational formulation, we first take the dot product of the first one by \( \mathbf{v}\) and integrate by parts

\[ \begin{align} \int_{\Omega} \mathbf{u}\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} \mathbf{u}\cdot \mathbf{v}~\mathrm{d} \mathbf{x} -\int_{\Omega} p ~ \nabla\cdot \mathbf{v}~\mathrm{d} \mathbf{x}=0, \end{align} \]

using \(p=0\) as a natural boundary condition. Then multiplying the second equation by \(q\) and integrating yields

\[ \begin{align} \int_{\Omega} \nabla\cdot\mathbf{u} ~ q ~\mathrm{d} \mathbf{x} = \int_{\Omega} f q ~\mathrm{d} \mathbf{x}. \end{align} \]

No integration by parts is necessary here. And we thus get the following mixed variational formulation:

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

\[\begin{split} \begin{align} \left\{ \begin{array}{llll} \int_{\Omega} \mathbf{u}\cdot \mathbf{v}~\mathrm{d} \mathbf{x} &- \int_{\Omega} p ~ \nabla\cdot \mathbf{v}~\mathrm{d} \mathbf{x} &=0, & \forall \mathbf{v}\in H(\mathrm{div},\Omega) \\ \int_{\Omega} \nabla\cdot\mathbf{u} ~ q ~\mathrm{d} \mathbf{x} & &= \int_{\Omega} f q ~\mathrm{d} \mathbf{x}, & \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,
                             element_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

from psydac.api.discretization import discretize

domain = Square()

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

x,y = domain.coordinates

# rhs
f = 2*pi**2*sin(pi*x)*sin(pi*y)
# exact solution
pe = sin(pi*x)*sin(pi*y)
#ue = -grad(pe) # not working when computing the norm
ue = Tuple(-pi*cos(pi*x)*sin(pi*y), -pi*sin(pi*x)*cos(pi*y))

u,v = [element_of(V1, name=i) for i in ['u', 'v']]
p,q = [element_of(V2, name=i) for i in ['p', 'q']]

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

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

# Variational problem
equation = find([u,p], forall=[v,q], lhs=a((u,p),(v,q)), rhs=l(v,q))

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 =  0.012985071480053038
>>> norm-l2 ph =  0.000641861830583991
Contents