Second mixed formulation of the Poisson problem

Second mixed formulation of the Poisson problem#

Here, we get an alternative formulation by not integrating by parts, the mixed term in the first formulation but in the second. The first formulation simply becomes

\[ \begin{align} \int_{\Omega} \mathbf{u}\cdot \mathbf{v}~\mathrm{d} \mathbf{x} +\int_{\Omega} \nabla p \cdot \mathbf{v}~\mathrm{d} \mathbf{x}=0, \end{align} \]

and the second, removing immediately the boundary term due to the essential boundary condition \(q=0\)

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

which leads to the variational formulation

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

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

Note that this formulation actually contains the classical variational formulation for the Poisson equation. Indeed for \(q\in H^1_0(\Omega)\), \(\nabla q \in L^2(\Omega)^3\) can be used as a test function in the first equation. And plugging this into the second we get

\[ \int_{\Omega} \nabla p \cdot \nabla q ~\mathrm{d} \mathbf{x} = \int_{\Omega} f q ~\mathrm{d} \mathbf{x}, \quad \forall q\in H^1_0(\Omega). \]

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='L2')
V2 = ScalarFunctionSpace('V2', domain, kind='H1')
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) + dot(grad(p),v) + dot(u,grad(q))) )

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

# Dirichlet boundary conditions
bc = [EssentialBC(p, 0, domain.boundary)]

# Variational problem
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 =  64.20822310868424
>>> norm-l2 ph =  3254706679445337.5
Contents