The First Steps

In this chapter we give an introduction how to use esys.escript to solve a partial differential equation (PDE ). We assume you are at least a little familiar with Python. The knowledge presented at the Python tutorial at http://docs.python.org/tut/tut.htmlis more than sufficient.

The PDE we wish to solve is the Poisson equation

$\displaystyle -\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

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

The domain is shown in Figure 1.1.
Figure 1.1: Domain $ \Omega =[0,1]^2$ with outer normal field $ n$.
\includegraphics[width=100mm]{figures/FirstStepDomain}

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

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

where, for any function $ u$ 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 $ u\hackscore{,ij}=(u\hackscore {,i})\hackscore{,j}$ which leads to

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

We often find that use of nested $ \sum$ symbols makes formulas cumbersome, and we use the more convenient Einstein summation convention . This drops the $ \sum$ sign and assumes that a summation is performed over any repeated index. For instance we write
$\displaystyle x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i}$     (5)
$\displaystyle x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i}$     (6)
$\displaystyle u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii}$     (7)
$\displaystyle 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

$\displaystyle - 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

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

$ n=(n\hackscore{i})$ denotes the outer normal field of the domain, see Figure 1.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 (1.11) should be fulfilled on the set $ \Gamma^N$ which is the top and right edge of the domain:

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

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

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

the solution shall be identically zero:

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

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

Figure 1.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.
\includegraphics[width=100mm]{figures/FirstStepMesh}

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 1.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 [28,5].

The esys.escript solver we want to use to solve this problem is embedded into the python interpreter language. So you can solve the problem interactively but you will learn quickly that is more efficient to use scripts which you can edit with your favorite editor. To enter the escript environment you use escript command[*]:

  escript
which will pass you on to the python prompt
Python 2.5.2 (r252:60911, Oct  5 2008, 19:29:17) 
[GCC 4.3.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>
Here you can use all available python commands and language features, for instance
\begin{python}
>>> x=2+3
>>> print ''2+3='',x
2+3= 5
\end{python}
We refer to the python users guide if you not familiar with python.

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
\begin{python}
from esys.finley import Rectangle
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
\end{python}
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 7.

The following statements define the Poisson class object mypde with domain mydomain and the right hand side $ f$ of the PDE to constant $ 1$:
\begin{python}
from esys.escript.linearPDEs import Poisson
mypde = Poisson(mydomain)
mypde.setValue(f=1)
\end{python}
We have not specified any boundary condition but the Poisson class implicitly assumes homogeneous Neuman boundary conditions defined by Equation (1.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 (1.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 contains the coordinates of the nodes in the domain use
\begin{python}
x=mydomain.getX()
\end{python}
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 will be discussed in Chapter 3 in more detail. What we need to know here is that

x has rank (number of dimensions) and a shape (list of dimensions) which can be viewed by calling the getRank and getShape methods:
\begin{python}
print ''rank '',x.getRank(),'', shape '',x.getShape()
\end{python}
This will print something like
\begin{python}
rank 1, shape (2,)
\end{python}
The Data object also maintains type information which is represented by the FunctionSpace of the object. For instance
\begin{python}
print x.getFunctionSpace()
\end{python}
will print
\begin{python}
Function space type: Finley_Nodes on FinleyMesh
\end{python}
which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a esys.finley mesh. To get the $ x\hackscore{0}$ coordinates of the locations we use the statement
\begin{python}
x0=x[0]
\end{python}
Object x0 is again a Data object now with rank 0 and shape $ ()$. It inherits the FunctionSpace from x:
\begin{python}
print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
\end{python}
will print
\begin{python}
0 () Function space type: Finley_Nodes on FinleyMesh
\end{python}
We can now construct a function gammaD which is only non-zero on the bottom and left edges of the domain with
\begin{python}
from esys.escript import whereZero
gammaD=whereZero(x[0])+whereZero(x[1])
\end{python}

whereZero(x[0]) creates function which equals $ 1$ where x[0] is (almost) 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 strictly 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
\begin{python}
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
\end{python}
one gets
\begin{python}
0 () Function space type: Finley_Nodes on FinleyMesh
\end{python}
An 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:
\begin{python}
from esys.linearPDEs import Poisson
x = mydomain.getX()
gammaD...
...1])
mypde = Poisson(domain=mydomain)
mypde.setValue(f=1,q=gammaD)
\end{python}
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 Poisson problem
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import ...
...in=mydomain)
mypde.setValue(f=1,q=gammaD)
u = mypde.getSolution()
\end{python}
The question is what we do with the calculated solution u. Besides postprocessing, eg. calculating the gradient or the average value, which will be discussed later, plotting the solution is one one things you might want to do. esys.escript offers two ways to do this, both base on external modules or packages and so data need to converted to hand over the solution. The first option is using the matplotlib module which allows plotting 2D results relatively quickly, see [12]. However, there are limitations when using this tool, eg. in problem size and when solving 3D problems. Therefore esys.escript provides a second options based on VTK files which is especially designed for large scale and 3D problem and which can be read by a variety of software packages such as mayavi [2], VisIt [24].



Subsections
esys@esscc.uq.edu.au