1 
ksteube 
1811 

2 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
ksteube 
1316 
% 
4 
ksteube 
1811 
% Copyright (c) 20032008 by University of Queensland 
5 


% Earth Systems Science Computational Center (ESSCC) 
6 


% http://www.uq.edu.au/esscc 
7 
gross 
625 
% 
8 
ksteube 
1811 
% Primary Business: Queensland, Australia 
9 


% Licensed under the Open Software License version 3.0 
10 


% http://www.opensource.org/licenses/osl3.0.php 
11 
gross 
625 
% 
12 
ksteube 
1811 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 
jgs 
102 

14 
ksteube 
1811 

15 
jgs 
121 
\section{The First Steps} 
16 
jgs 
102 
\label{FirstSteps} 
17 



18 



19 
artak 
1971 

20 
ksteube 
1316 
In this chapter we give an introduction how to use \escript to solve 
21 


a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). We assume you are at least a little familiar with Python. The knowledge presented at the Python tutorial at \url{http://docs.python.org/tut/tut.html} 
22 


is more than sufficient. 
23 
jgs 
102 

24 
gross 
2363 
The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation} 
25 
jgs 
102 
\begin{equation} 
26 


\Delta u =f 
27 


\label{eq:FirstSteps.1} 
28 


\end{equation} 
29 
ksteube 
1316 
for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$, 
30 
jgs 
102 
is the unit square 
31 


\begin{equation} 
32 


\Omega=[0,1]^2=\{ (x\hackscore 0;x\hackscore 1)  0\le x\hackscore{0} \le 1 \mbox{ and } 0\le x\hackscore{1} \le 1 \} 
33 


\label{eq:FirstSteps.1b} 
34 


\end{equation} 
35 
jgs 
107 
The domain is shown in \fig{fig:FirstSteps.1}. 
36 
gross 
2371 
\begin{figure} [ht] 
37 
artak 
1971 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}} 
38 


\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.} 
39 


\label{fig:FirstSteps.1} 
40 


\end{figure} 
41 
jgs 
102 

42 
ksteube 
1316 
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by 
43 
jgs 
102 
\begin{equation} 
44 


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


\label{eq:FirstSteps.1.1} 
46 


\end{equation} 
47 
ksteube 
1316 
where, for any function $u$ and any direction $i$, $u\hackscore{,i}$ 
48 
jgs 
107 
denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$. 
49 
ksteube 
1316 
\footnote{You 
50 
jgs 
102 
may be more familiar with the Laplace operator\index{Laplace operator} being written 
51 


as $\nabla^2$, and written in the form 
52 


\begin{equation*} 
53 
jgs 
110 
\nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x\hackscore 0^2} 
54 
jgs 
102 
+ \frac{\partial^2 u}{\partial x\hackscore 1^2} 
55 


\end{equation*} 
56 


and \eqn{eq:FirstSteps.1} as 
57 


\begin{equation*} 
58 


\nabla^2 u = f 
59 


\end{equation*} 
60 


} 
61 
jgs 
107 
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
62 
artak 
1971 
to the index. To get a more compact form we will write $u\hackscore{,ij}=(u\hackscore {,i})\hackscore{,j}$ 
63 
jgs 
102 
which leads to 
64 


\begin{equation} 
65 


\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii} 
66 


\label{eq:FirstSteps.1.1b} 
67 


\end{equation} 
68 
ksteube 
1316 
We often find that use 
69 


of nested $\sum$ symbols makes formulas cumbersome, and we use the more 
70 


convenient Einstein summation convention \index{summation convention}. This 
71 


drops the $\sum$ sign and assumes that a summation is performed over any repeated index. 
72 


For instance we write 
73 
jgs 
102 
\begin{eqnarray} 
74 


x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i} \\ 
75 


x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i} \\ 
76 


u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\ 
77 
jgs 
107 
x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j} \\ 
78 
jgs 
102 
\label{eq:FirstSteps.1.1c} 
79 


\end{eqnarray} 
80 


With the summation convention we can write the Poisson equation \index{Poisson equation} as 
81 


\begin{equation} 
82 


 u\hackscore{,ii} =1 
83 


\label{eq:FirstSteps.1.sum} 
84 


\end{equation} 
85 
lkettle 
575 
where $f=1$ in this example. 
86 



87 
jgs 
102 
On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$ 
88 


of the solution $u$ shall be zero, ie. $u$ shall fulfill 
89 


the homogeneous Neumann boundary condition\index{Neumann 
90 


boundary condition!homogeneous} 
91 


\begin{equation} 
92 


n\hackscore{i} u\hackscore{,i}= 0 \;. 
93 


\label{eq:FirstSteps.2} 
94 


\end{equation} 
95 


$n=(n\hackscore{i})$ denotes the outer normal field 
96 


of the domain, see \fig{fig:FirstSteps.1}. Remember that we 
97 


are applying the Einstein summation convention \index{summation convention}, i.e 
98 
jgs 
107 
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} + 
99 


n\hackscore{1} u\hackscore{,1}$. 
100 
jgs 
102 
\footnote{Some readers may familiar with the notation 
101 
artak 
1971 
$ 
102 
jgs 
102 
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i} 
103 
artak 
1971 
$ 
104 
jgs 
102 
for the normal derivative.} 
105 


The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the 
106 


set $\Gamma^N$ which is the top and right edge of the domain: 
107 


\begin{equation} 
108 


\Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega  x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1 \} 
109 


\label{eq:FirstSteps.2b} 
110 


\end{equation} 
111 
jgs 
107 
On the bottom and the left edge of the domain which is defined 
112 
jgs 
102 
as 
113 


\begin{equation} 
114 


\Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega  x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0 \} 
115 


\label{eq:FirstSteps.2c} 
116 


\end{equation} 
117 


the solution shall be identically zero: 
118 


\begin{equation} 
119 


u=0 \; . 
120 


\label{eq:FirstSteps.2d} 
121 


\end{equation} 
122 
jgs 
107 
This kind of boundary condition is called a homogeneous Dirichlet boundary condition 
123 
jgs 
102 
\index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together 
124 


with the Neumann boundary condition \eqn{eq:FirstSteps.2} and 
125 


Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so 
126 


called boundary value 
127 
jgs 
107 
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for 
128 
artak 
1971 
the unknown function~$u$. 
129 
jgs 
102 

130 



131 
gross 
2371 
\begin{figure}[ht] 
132 
jfenwick 
2335 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}} 
133 
jgs 
102 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
134 


each element is a quadrilateral and described by four nodes, namely 
135 


the corner points. The solution is interpolated by a bilinear 
136 


polynomial.} 
137 


\label{fig:FirstSteps.2} 
138 


\end{figure} 
139 



140 


In general the BVP\index{boundary value problem!BVP} cannot be solved analytically and numerical 
141 


methods have to be used construct an approximation of the solution 
142 


$u$. Here we will use the finite element method\index{finite element 
143 


method} (FEM\index{finite element 
144 


method!FEM}). The basic idea is to fill the domain with a 
145 
jgs 
107 
set of points called nodes. The solution is approximated by its 
146 
jgs 
102 
values on the nodes\index{finite element 
147 
lkettle 
573 
method!nodes}. Moreover, the domain is subdivided into smaller 
148 


subdomains called elements \index{finite element 
149 
jgs 
102 
method!element}. On each element the solution is 
150 


represented by a polynomial of a certain degree through its values at 
151 


the nodes located in the element. The nodes and its connection through 
152 


elements is called a mesh\index{finite element 
153 
jgs 
107 
method!mesh}. \fig{fig:FirstSteps.2} shows an 
154 
jgs 
102 
example of a FEM mesh with four elements in the $x_0$ and four elements 
155 


in the $x_1$ direction over the unit square. 
156 
gross 
2370 
For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}. 
157 
jgs 
102 

158 
gross 
2370 
The \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 
159 


that is more efficient to use scripts which you can edit with your favourite editor. 
160 


To enter the escript environment you use \program{escript} command\footnote{\program{escript} is not available under Windows yet. If you run under windows you can just use the 
161 


\program{python} command and the \env{OMP_NUM_THREADS} environment variable to control the number 
162 


of threads.}: 
163 


\begin{verbatim} 
164 


escript 
165 


\end{verbatim} 
166 


which will pass you on to the python prompt 
167 


\begin{verbatim} 
168 


Python 2.5.2 (r252:60911, Oct 5 2008, 19:29:17) 
169 


[GCC 4.3.2] on linux2 
170 


Type "help", "copyright", "credits" or "license" for more information. 
171 


>>> 
172 


\end{verbatim} 
173 


Here you can use all available python commands and language features, for instance 
174 


\begin{python} 
175 


>>> x=2+3 
176 


>>> print "2+3=",x 
177 


2+3= 5 
178 


\end{python} 
179 


We refer to the python users guide if you not familar with python. 
180 



181 
jgs 
102 
\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}. 
182 


(We will discuss a more general form of a PDE \index{partial differential equation!PDE} 
183 


that can be defined through the \LinearPDE class later). The instantiation of 
184 


a \Poisson class object requires the specification of the domain $\Omega$. In \escript 
185 


the \Domain class objects are used to describe the geometry of a domain but it also 
186 


contains information about the discretization methods and the actual solver which is used 
187 


to solve the PDE. Here we are using the FEM library \finley \index{finite element 
188 


method}. The following statements create the \Domain object \var{mydomain} from the 
189 


\finley method \method{Rectangle} 
190 


\begin{python} 
191 
ksteube 
1316 
from esys.finley import Rectangle 
192 


mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
193 
jgs 
102 
\end{python} 
194 


In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and 
195 


the right, upper corner at $(\var{l0},\var{l1})=(1,1)$. 
196 
jgs 
107 
The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and 
197 
jgs 
102 
$x\hackscore{1}$direction respectively. For more details on \method{Rectangle} and 
198 


other \Domain generators within the \finley module, 
199 


see \Chap{CHAPTER ON FINLEY}. 
200 



201 
jgs 
107 
The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and 
202 
jgs 
102 
the right hand side $f$ of the PDE to constant $1$: 
203 


\begin{python} 
204 
ksteube 
1316 
from esys.escript.linearPDEs import Poisson 
205 


mypde = Poisson(mydomain) 
206 


mypde.setValue(f=1) 
207 
jgs 
102 
\end{python} 
208 


We have not specified any boundary condition but the 
209 


\Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann 
210 


boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary 
211 


condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$ 
212 


and any constant $C$ the function $u+C$ becomes a solution as well. We have to add 
213 


a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done 
214 
jgs 
107 
by defining a characteristic function \index{characteristic function} 
215 


which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set 
216 
jgs 
102 
and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, 
217 
gross 
565 
we need to construct a function \var{gammaD} which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$. To get 
218 
ksteube 
1316 
an object \var{x} which contains the coordinates of the nodes in the domain use 
219 
jgs 
102 
\begin{python} 
220 
ksteube 
1316 
x=mydomain.getX() 
221 
jgs 
102 
\end{python} 
222 
gross 
660 
The method \method{getX} of the \Domain \var{mydomain} 
223 
jgs 
107 
gives access to locations 
224 
ksteube 
1316 
in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be 
225 


discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that 
226 
gross 
660 

227 
ksteube 
1316 
\var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by 
228 
gross 
565 
calling the \method{getRank} and \method{getShape} methods: 
229 


\begin{python} 
230 
ksteube 
1316 
print "rank ",x.getRank(),", shape ",x.getShape() 
231 
gross 
565 
\end{python} 
232 
ksteube 
1316 
This will print something like 
233 
gross 
565 
\begin{python} 
234 
ksteube 
1316 
rank 1, shape (2,) 
235 
gross 
565 
\end{python} 
236 


The \Data object also maintains type information which is represented by the 
237 


\FunctionSpace of the object. For instance 
238 


\begin{python} 
239 
ksteube 
1316 
print x.getFunctionSpace() 
240 
gross 
565 
\end{python} 
241 


will print 
242 


\begin{python} 
243 
ksteube 
1316 
Function space type: Finley_Nodes on FinleyMesh 
244 
gross 
565 
\end{python} 
245 
ksteube 
1316 
which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a \finley mesh. 
246 
gross 
565 
To get the $x\hackscore{0}$ coordinates of the locations we use the 
247 


statement 
248 


\begin{python} 
249 
ksteube 
1316 
x0=x[0] 
250 
gross 
565 
\end{python} 
251 


Object \var{x0} 
252 


is again a \Data object now with \Rank $0$ and 
253 


\Shape $()$. It inherits the \FunctionSpace from \var{x}: 
254 


\begin{python} 
255 
ksteube 
1316 
print x0.getRank(),x0.getShape(),x0.getFunctionSpace() 
256 
gross 
565 
\end{python} 
257 


will print 
258 


\begin{python} 
259 
ksteube 
1316 
0 () Function space type: Finley_Nodes on FinleyMesh 
260 
gross 
565 
\end{python} 
261 
ksteube 
1316 
We can now construct a function \var{gammaD} which is only nonzero on the bottom and left edges 
262 


of the domain with 
263 
gross 
565 
\begin{python} 
264 
ksteube 
1316 
from esys.escript import whereZero 
265 


gammaD=whereZero(x[0])+whereZero(x[1]) 
266 
gross 
565 
\end{python} 
267 
ksteube 
1316 

268 


\code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero 
269 
jgs 
107 
and $0$ elsewhere. 
270 
gross 
565 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
271 
jgs 
107 
equal to zero and $0$ elsewhere. 
272 
gross 
565 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
273 
artak 
1971 
gives a function on the domain \var{mydomain} which is strictly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero. 
274 
gross 
565 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
275 


\begin{python} 
276 
ksteube 
1316 
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace() 
277 
gross 
565 
\end{python} 
278 


one gets 
279 


\begin{python} 
280 
ksteube 
1316 
0 () Function space type: Finley_Nodes on FinleyMesh 
281 
gross 
565 
\end{python} 
282 
ksteube 
1316 
An additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the 
283 
jgs 
102 
characteristic function \index{characteristic function} of the locations 
284 


of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous} 
285 


are set. The complete definition of our example is now: 
286 


\begin{python} 
287 
ksteube 
1316 
from esys.linearPDEs import Poisson 
288 


x = mydomain.getX() 
289 


gammaD = whereZero(x[0])+whereZero(x[1]) 
290 


mypde = Poisson(domain=mydomain) 
291 


mypde.setValue(f=1,q=gammaD) 
292 
jgs 
102 
\end{python} 
293 
lkettle 
573 
The first statement imports the \Poisson class definition from the \linearPDEs module \escript package. 
294 
jgs 
107 
To get the solution of the Poisson equation defined by \var{mypde} we just have to call its 
295 
jgs 
102 
\method{getSolution}. 
296 



297 
ksteube 
1316 
Now we can write the script to solve our Poisson problem 
298 
jgs 
102 
\begin{python} 
299 
ksteube 
1316 
from esys.escript import * 
300 


from esys.escript.linearPDEs import Poisson 
301 


from esys.finley import Rectangle 
302 


# generate domain: 
303 


mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
304 


# define characteristic function of Gamma^D 
305 


x = mydomain.getX() 
306 


gammaD = whereZero(x[0])+whereZero(x[1]) 
307 


# define PDE and get its solution u 
308 


mypde = Poisson(domain=mydomain) 
309 


mypde.setValue(f=1,q=gammaD) 
310 


u = mypde.getSolution() 
311 


# write u to an external file 
312 


saveVTK("u.xml",sol=u) 
313 
jgs 
102 
\end{python} 
314 
ksteube 
1316 
The entire code is available as \file{poisson.py} in the \ExampleDirectory 
315 
jgs 
102 

316 
ksteube 
1316 
The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in 
317 


\VTK file format. 
318 
gross 
2370 
Now you may run the script using the \escript environment 
319 


and visualize the solution using \mayavi: 
320 
ksteube 
1316 
\begin{verbatim} 
321 
gross 
2370 
escript poisson.py 
322 
ksteube 
1316 
mayavi d u.xml m SurfaceMap 
323 


\end{verbatim} 
324 


See \fig{fig:FirstSteps.3}. 
325 



326 
jgs 
102 
\begin{figure} 
327 
jfenwick 
2335 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}} 
328 
lkettle 
573 
\caption{Visualization of the Poisson Equation Solution for $f=1$} 
329 
jgs 
102 
\label{fig:FirstSteps.3} 
330 


\end{figure} 
331 


