The Poisson equation

The Poisson equation#

As a first example, we consider the Poisson equation

\[\begin{split} \begin{align} - \nabla^2 u = f \quad &\text{in $\Omega$}, \\ u = 0 \quad &\text{on $\Gamma_0$}, \\ u = g_i \quad &\text{on $\Gamma_I$}, \\ \partial_n u = g_n \quad &\text{on $\Gamma_N := \partial \Omega \setminus \left( \Gamma_0 \cup \Gamma_I \right)$}. \end{align} \end{split}\]

Variational Formulation#

An \(H^1\)-conforming variational formulation of 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^1(\Omega)\),

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

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

Formal Model#

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 grad, dot, laplace
from sympde.topology import NormalVector, Union

from sympy           import pi, sin

from psydac.api.discretization import discretize

domain = Square()
Gamma_0 = domain.get_boundary(axis=0, ext=-1)
Gamma_i = domain.get_boundary(axis=0, ext=1)
Gamma_n = domain.boundary.complement(Union(Gamma_0, Gamma_i))
nn = NormalVector('nn')

V = ScalarFunctionSpace('V', domain)

x,y = domain.coordinates

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

# bilinear form
a = BilinearForm((u,v), integral(domain , dot(grad(v), grad(u))))

# exact solution
ue = sin(pi*x) * (1+y*sin(pi*y/3))**2
L = lambda w: - laplace(w)
f = L(ue)
gi = ue
gn = ue

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

# Boundary term for the Neumann BC
ln = LinearForm(v, integral(Gamma_n, v * dot(grad(gn), nn)))

# Dirichlet boundary conditions
bc  = [EssentialBC(u,  0, Gamma_0)]
bc += [EssentialBC(u, gi, Gamma_i)]

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(v)+ln(v), 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
Vh = discretize(V, domain_h, degree=degree)

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

Solving the PDE#

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

uh = equation_h.solve()

Computing the error norm#

When the analytical solution is available, you might be interested in computing the \(L^2\) norm or \(H^1_0\) semi-norm. SymPDE allows you to do so, by creating the Norm object. In this example, the analytical solution is given by

\[ u_e = \sin(\pi x) \sin(\pi y) \]

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.00042499840633644024

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)

0.025283770887649267

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)

0.025287342563909892
Contents