The Biharmonic problem

The Biharmonic problem#

We consider the (inhomogeneous) biharmonic equation with homogeneous essential boundary conditions,

\[ \begin{align} \nabla^2 \nabla^2 u = f \quad \text{in $\Omega$} , \quad \quad u = 0 \quad \text{on $\partial \Omega$}, \quad \quad \nabla u \cdot \mathbf{n} = 0 \quad \text{on $\partial \Omega$}. \end{align} \]

The Variational Formulation#

An \(H^2\)-conforming variational formulation reads

\[ \begin{align} \text{find $u \in V$ such that} \quad a(u,v) = l(v) \quad \forall v \in V, \end{align} \]

where

  • \(V \subset H^2(\Omega)\),

  • \(a(u,v) := \int_{\Omega} \nabla^2 u ~ \nabla^2 v ~ d\Omega\),

  • \(l(v) := \int_{\Omega} f v ~ d\Omega\).

Formal Model#

In this example, the analytical solution is given by

\[ u_e = \sin(\pi x)^2 \sin(\pi y)^2 \]
from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import ScalarFunctionSpace, Square, element_of
from sympde.calculus import laplace, dot

from sympy import pi, sin

domain = Square()

V = ScalarFunctionSpace('V', domain)

x,y = domain.coordinates

u,v = [element_of(V, name=i) for i in ['u', 'v']]

# the analytical solution and its rhs
ue = sin(pi * x)**2 * sin(pi * y)**2
f  = laplace(laplace(ue))

# bilinear form
a = BilinearForm((u,v), integral(domain , laplace(u)*laplace(v)))

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

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

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)

Discretization#

We shall need the discretize function from PsyDAC.

from psydac.api.discretization import discretize

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

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

# Create discrete spline space
Vh = discretize(V, domain_h, degree=degree)

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

Solving the PDE#

uh = equation_h.solve()

Computing the error norm#

Computing the \(L^2\) norm#

u = element_of(V, name='u')

# create the formal Norm object
l2norm = Norm(u - ue, domain, kind='l2')

# discretize the norm
l2norm_h = discretize(l2norm, domain_h, Vh)

# assemble the norm
l2_error = l2norm_h.assemble(u=uh)

# print the result
print(l2_error)

0.5367871671842632

Computing the \(H^1\) semi-norm#

# create the formal Norm object
h1norm = SemiNorm(u - ue, domain, kind='h1')

# discretize the norm
h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
h1_error = h1norm_h.assemble(u=uh)

# print the result
print(h1_error)

2.4343747827631277

Computing the \(H^1\) norm#

# create the formal Norm object
h1norm = Norm(u - ue, domain, kind='h1')

# discretize the norm
h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
h1_error = h1norm_h.assemble(u=uh)

# print the result
print(h1_error)

2.492853996086961

Computing the H^2 semi-norm#

# create the formal Norm object
h2norm = SemiNorm(u - ue, domain, kind='h2')

# discretize the norm
h2norm_h = discretize(h2norm, domain_h, Vh)

# assemble the norm
h2_error = h2norm_h.assemble(u=uh)

# print the result
print(h2_error)

9.115861322114647
Contents