Your first code using SymPDE & PsyDAC

Your first code using SymPDE & PsyDAC#

Author: Ahmed Ratnani

We start by writing our first example using SymPDE. Let \(\Omega := (0,1)^2\). We consider the Poisson problem with homogeneous Dirichlet boundary conditions.

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

Variational Formulation#

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

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

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

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

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']]

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

# linear form
f = 2*pi**2*sin(pi*x)*sin(pi*y)
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)

This very simple Python code reflects well the abstract mathematical framework needed for variational formulations. The structure of the code is as follows,

  1. Create a domain.

  2. Create a space of scalar functions over the domain.

  3. Create elements from this function space. These elements will denote the test and trial functions.

  4. Create the Bilinear and Linear forms, \(a\) and \(l\) respectively.

  5. Create Essential Boundary Conditions.

  6. Create the variational problem.

Most of the time, you will need to follow the same steps, with some minor variants depending on the problem you’re considering.

Discretization#

We shall need the discretize function from PsyDAC.

from psydac.api.discretization import discretize

degree = [3,3]
ncells = [16,16]

# 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#

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#

ue = sin(pi*x)*sin(pi*y)

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)

9.49756716724664e-07

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)

9.768702410721585e-05

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)

9.769164089397408e-05
Contents