16. Poisson equation with periodic boundary conditions

This demo is implemented in a single Python file, demo_periodic.py, which contains both the variational form and the solver.

This demo illustrates how to:

  • Solve a linear partial differential equation
  • Read mesh and subdomains from file
  • Create and apply Dirichlet and periodic boundary conditions

The solution for u in this demo will look as follows:

../../../../_images/plot_u1.png

16.1. Equation and problem definition

The Poisson equation is the canonical elliptic partial differential equation. For a domain \(\Omega \subset \mathbb{R}^n\) with boundary \(\partial \Omega = \Gamma_{D} \cup \Gamma_{P}\), the Poisson equation with particular boundary conditions reads:

\[\begin{split}- \nabla \cdot (\nabla u) &= f \quad {\rm in} \ \Omega, \\ u &= 0 \quad {\rm on} \ \Gamma_{D}, \\ u(0, y) &= u(1, y) \quad {\rm on} \ \Gamma_{P}. \\\end{split}\]

Here, \(f\) is a given source function. The most standard variational form of Poisson equation reads: find \(u \in V\) such that

\[a(u, v) = L(v) \quad \forall \ v \in V,\]

where \(V\) is a suitable function space and

\[\begin{split}a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x, \\ L(v) &= \int_{\Omega} f v \, {\rm d} x.\end{split}\]

The expression \(a(u, v)\) is the bilinear form and \(L(v)\) is the linear form. It is assumed that all functions in \(V\) satisfy the Dirichlet boundary conditions (\(u = 0 \ {\rm on} \ \Gamma_{D}\)).

In this demo, we shall consider the following definitions of the input functions, the domain, and the boundaries:

  • \(\Omega = [0,1] \times [0,1]\) (a unit square)
  • \(\Gamma_{D} = \{(x, 0) \cup (x, 1) \subset \partial \Omega\}\) (Dirichlet boundary)
  • \(\Gamma_{P} = \{(0, y) \cup (1, y) \subset \partial \Omega\}\) (Periodic boundary)
  • \(f = x \sin(5.0 \pi y) + \exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)\) (source term)

16.2. Implementation

This demo is implemented in a single file, demo_periodic.py.

First, the dolfin module is imported

from dolfin import *

A subclass of Expression, Source, is created for the source term f. The function eval() returns values for a function at the given point x.

# Source term
class Source(Expression):
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02)

To define the boundaries, we create subclasses of the class SubDomain. A simple Python function, returning a boolean, can be used to define the subdomain for the Dirichlet boundary condition (\(\Gamma_D\)). The function should return True for those points inside the subdomain and False for the points outside. In our case, we want to say that the points \((x, y)\) such that \(y = 0\) or \(y = 1\) are inside of \(\Gamma_D\). (Note that because of round-off errors, it is often wise to instead specify \(y < \epsilon\) or \(y > 1 - \epsilon\) where \(\epsilon\) is a small number (such as machine precision).)

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool((x[1] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS)) and on_boundary)

The periodic boundary is defined by PeriodicBoundary and we define what is inside the boundary in the same way as in DirichletBoundary. The function map maps a coordinate x in domain H to a coordinate y in the domain G, it is used for periodic boundary conditions, so that the right boundary of the domain is mapped to the left boundary. When the class is defined, we create the boundary by making an instance of the class.

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]

# Create periodic boundary condition
pbc = PeriodicBoundary()

A 2D mesh is created using the built-in class UnitSquareMesh, and we define a finite element function space relative to this space. Notice the fourth argument of FunctionSpace. It specifies that all functions in V have periodic boundaries defined by pbc.

# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

Now, we create the Dirichlet boundary condition using the class DirichletBC. A DirichletBC takes three arguments: the function space the boundary condition applies to, the value of the boundary condition, and the part of the boundary on which the condition applies. In our example, the function space is V, the value of the boundary condition (0.0) can be represented using a Constant and the Dirichlet boundary is defined by the class DirichletBoundary. The definition of the Dirichlet boundary condition then looks as follows:

# Create Dirichlet boundary condition
u0 = Constant(0.0)
dbc = DirichletBoundary()
bc0 = DirichletBC(V, u0, dbc)

When all boundary conditions are defined and created we can collect them in a list:

# Collect boundary conditions
bcs = [bc0]

Here only the Dirichlet boundary condition is put into the list because the periodic boundary condition is already applied in the definition of the function space. Next, we want to express the variational problem. First, we need to specify the trial function u and the test function v, both living in the function space V. We do this by defining a TrialFunction and a TestFunction on the previously defined FunctionSpace V. The source function f is created by making an instance of Source. With these ingredients, we can write down the bilinear form a and the linear form L (using UFL operators). In summary, this reads

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Source()
a = dot(grad(u), grad(v))*dx
L = f*v*dx

Now, we have specified the variational form and can consider the solution of the variational problem. First, we need to define a Function u to represent the solution. (Upon initialization, it is simply set to the zero function.) A Function represents a function living in a finite element function space. Next, we can call the solve function with the arguments a == L, u and bcs as follows:

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

The function u will be modified during the call to solve. The default settings for solving a variational problem have been used. However, the solution process can be controlled in much more detail if desired.

A Function can be manipulated in various ways, in particular, it can be plotted and saved to file. Here, we output the solution to a VTK file (using the suffix .pvd) for later visualization and also plot it using the plot command:

# Save solution to file
file = File("periodic.pvd")
file << u

# Plot solution
plot(u, interactive=True)

16.3. Complete code


from dolfin import *

# Source term
class Source(Expression):
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) \
                    + 1.0*exp(-(dx*dx + dy*dy)/0.02)

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool((x[1] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS)) \
                    and on_boundary)

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]

# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())

# Create Dirichlet boundary condition
u0 = Constant(0.0)
dbc = DirichletBoundary()
bc0 = DirichletBC(V, u0, dbc)

# Collect boundary conditions
bcs = [bc0]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Source()
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Save solution to file
file = File("periodic.pvd")
file << u

# Plot solution
plot(u, interactive=True)