The First Steps

Figure 2.1: Domain $\Omega =[0,1]^2$ with outer normal field $n$.

In this chapter we will give an introduction how to use esys.escript to solve a partial differential equation (PDE ). The reader should be familiar with Python. The knowledge presented at the Python tutorial at http://docs.python.org/tut/tut.htmlis sufficient. It is helpful if the reader has some basic knowledge of PDEs .

The PDE we wish to solve is the Poisson equation

$$-\Delta u =f$$ (1)

for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$ is the unit square

$$\Omega=[0,1]^2=\{ (x\hackscore 0;x\hackscore 1) \vert 0\le x\hackscore{0} \le 1 0\le x\hackscore{1} \le 1 \}$$ (2)

The domain is shown in Figure 2.1.

$\Delta$ denotes the Laplace operator which is defined by

$$\Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}$$ (3)

where, for any function $w$ and any direction $i$, $u\hackscore{,i}$ denotes the partial derivative of $u$ with respect to $i$. Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$ which leads to

$$\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}$$ (4)

In some cases, and we will see examples for this in the next chapter, the usage of the nested $\sum$ symbols blows up the formulas and therefore it is convenient to use the Einstein summation convention . This drops the $\sum$ sign and assumes that a summation over a repeated index is performed ("repeated index means summation"). For instance we write

$$x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i}$$ (5)

$$x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i}$$ (6)

$$u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii}$$ (7)

$$x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j}$$ (8)

With the summation convention we can write the Poisson equation as

$$- u\hackscore{,ii} =1$$ (9)

where $f=1$ in this example.

On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$ of the solution $u$ shall be zero, ie. $u$ shall fulfill the homogeneous Neumann boundary condition

$$n\hackscore{i} u\hackscore{,i}= 0 \;.$$ (10)

$n=(n\hackscore{i})$ denotes the outer normal field of the domain, see Figure 2.1. Remember that we are applying the Einstein summation convention , i.e $n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} + n\hackscore{1} u\hackscore{,1}$. The Neumann boundary condition of Equation (2.11) should be fulfilled on the set $\Gamma^N$ which is the top and right edge of the domain:

$$\Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega \vert x\hackscore{0}=1 x\hackscore{1}=1 \}$$ (11)

On the bottom and the left edge of the domain which is defined as

$$\Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega \vert x\hackscore{0}=0 x\hackscore{1}=0 \}$$ (12)

the solution shall be identically zero:

$$u=0 \; .$$ (13)

This kind of boundary condition is called a homogeneous Dirichlet boundary condition . The partial differential equation in Equation (2.10) together with the Neumann boundary condition Equation (2.11) and Dirichlet boundary condition in Equation (2.14) form a so called boundary value problem (BVP) for the unknown function $u$.

Figure 2.2: Mesh of $4 \time 4$ elements on a rectangular domain. Here each element is a quadrilateral and described by four nodes, namely the corner points. The solution is interpolated by a bi-linear polynomial.

In general the BVP cannot be solved analytically and numerical methods have to be used construct an approximation of the solution $u$. Here we will use the finite element method (FEM). The basic idea is to fill the domain with a set of points called nodes. The solution is approximated by its values on the nodes. Moreover, the domain is subdivided into smaller sub-domains called elements . On each element the solution is represented by a polynomial of a certain degree through its values at the nodes located in the element. The nodes and its connection through elements is called a mesh. Figure 2.2 shows an example of a FEM mesh with four elements in the $x_0$ and four elements in the $x_1$ direction over the unit square. For more details we refer the reader to the literature, for instance Reference [6,1].

esys.escript provides the class Poisson to define a Poisson equation . (We will discuss a more general form of a PDE that can be defined through the LinearPDE class later). The instantiation of a Poisson class object requires the specification of the domain $\Omega$. In esys.escript the Domain class objects are used to describe the geometry of a domain but it also contains information about the discretization methods and the actual solver which is used to solve the PDE. Here we are using the FEM library esys.finley . The following statements create the Domain object mydomain from the esys.finley method Rectangle

In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and the right, upper corner at $(\var{l0},\var{l1})=(1,1)$. The arguments n0 and n1 define the number of elements in $x\hackscore{0}$ and $x\hackscore{1}$-direction respectively. For more details on Rectangle and other Domain generators within the esys.finley module, see Chapter 5.

The following statements define the Poisson class object mypde with domain mydomain and the right hand side $f$ of the PDE to constant $1$:

We have not specified any boundary condition but the Poisson class implicitly assumes homogeneous Neuman boundary conditions defined by Equation (2.11). With this boundary condition the BVP we have defined has no unique solution. In fact, with any solution $u$ and any constant $C$ the function $u+C$ becomes a solution as well. We have to add a Dirichlet boundary condition . This is done by defining a characteristic function which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set and 0 elsewhere. In our case of $\Gamma^D$ defined by Equation (2.13), we need to construct a function gammaD which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$. To get an object x which represents locations in the domain one uses

The method getX of the Domain mydomain gives access to locations in the domain defined by mydomain. The object x is actually a Data object which is discussed in Chapter 3 in more details. What we need to know here is that

x has rank (=number of dimensions) and a shape (=tuple of dimensions) which can be checked by calling the getRank and getShape methods:

will print something like

The Data object also maintains type information which is represented by the FunctionSpace of the object. For instance

will print

which tells us that the coordinates are stored on the nodes of a esys.finley mesh. To get the $x\hackscore{0}$ coordinates of the locations we use the statement

Object x0 is again a Data object now with rank 0 and shape $()$. It inherits the FunctionSpace from x:

will print

We can now construct the function gammaD by

where whereZero(x[0]) creates function which equals $1$ where x[0] is (allmost) equal to zero and 0 elsewhere. Similarly, whereZero(x[1]) creates function which equals $1$ where x[1] is equal to zero and 0 elsewhere. The sum of the results of whereZero(x[0]) and whereZero(x[1]) gives a function on the domain mydomain which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero. Note that gammaD has the same rank , shape and FunctionSpace like x0 used to define it. So from

one gets

The additional parameter q of the setValue method of the Poisson class defines the characteristic function of the locations of the domain where homogeneous Dirichlet boundary condition are set. The complete definition of our example is now:

The first statement imports the Poisson class definition from the esys.escript.linearPDEs module esys.escript package. To get the solution of the Poisson equation defined by mypde we just have to call its getSolution.

Now we can write the script to solve our test problem (Remember that lines starting with '#' are comment lines in Python) (available as poisson.py in the example directory):

The last statement writes the solution tagged with the name "sol" to the external file u.xml in vtk [4]file format. vtk [4]is a software library for the visualization of scientific, engineering and analytical data and is freely available from http://www.vtk.org. There are a variety of graphical user interfaces for vtk [4]available, for instance mayavi which can be downloaded from http://mayavi.sourceforge.net/ but is also available on most Linux distributions.

Figure 2.3: Visualization of the Poisson Equation Solution for $f=1$

You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE for Python, see http://idlefork.sourceforge.net). If the script file has the name poisson.py you can run the script from any shell using the command:

After the script has (hopefully successfully) been completed you will find the file u.xml in the current directory. An easy way to visualize the results is the command

to show the results, see Figure 2.3.