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 
jgs 
107 
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 
artak 
1971 
\begin{figure} [h!] 
37 


\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 
artak 
1971 
\begin{figure}[h] 
132 
gross 
599 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}} 
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 


For more details we refer the reader to the literature, for instance 
157 


\Ref{Zienc,NumHand}. 
158 



159 


\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}. 
160 


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


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


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


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


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


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


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


\finley method \method{Rectangle} 
168 


\begin{python} 
169 
ksteube 
1316 
from esys.finley import Rectangle 
170 


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


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


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


other \Domain generators within the \finley module, 
177 


see \Chap{CHAPTER ON FINLEY}. 
178 



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


\begin{python} 
182 
ksteube 
1316 
from esys.escript.linearPDEs import Poisson 
183 


mypde = Poisson(mydomain) 
184 


mypde.setValue(f=1) 
185 
jgs 
102 
\end{python} 
186 


We have not specified any boundary condition but the 
187 


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


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


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


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


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


which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set 
194 
jgs 
102 
and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, 
195 
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 
196 
ksteube 
1316 
an object \var{x} which contains the coordinates of the nodes in the domain use 
197 
jgs 
102 
\begin{python} 
198 
ksteube 
1316 
x=mydomain.getX() 
199 
jgs 
102 
\end{python} 
200 
gross 
660 
The method \method{getX} of the \Domain \var{mydomain} 
201 
jgs 
107 
gives access to locations 
202 
ksteube 
1316 
in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be 
203 


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

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


\begin{python} 
208 
ksteube 
1316 
print "rank ",x.getRank(),", shape ",x.getShape() 
209 
gross 
565 
\end{python} 
210 
ksteube 
1316 
This will print something like 
211 
gross 
565 
\begin{python} 
212 
ksteube 
1316 
rank 1, shape (2,) 
213 
gross 
565 
\end{python} 
214 


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


\FunctionSpace of the object. For instance 
216 


\begin{python} 
217 
ksteube 
1316 
print x.getFunctionSpace() 
218 
gross 
565 
\end{python} 
219 


will print 
220 


\begin{python} 
221 
ksteube 
1316 
Function space type: Finley_Nodes on FinleyMesh 
222 
gross 
565 
\end{python} 
223 
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. 
224 
gross 
565 
To get the $x\hackscore{0}$ coordinates of the locations we use the 
225 


statement 
226 


\begin{python} 
227 
ksteube 
1316 
x0=x[0] 
228 
gross 
565 
\end{python} 
229 


Object \var{x0} 
230 


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


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


\begin{python} 
233 
ksteube 
1316 
print x0.getRank(),x0.getShape(),x0.getFunctionSpace() 
234 
gross 
565 
\end{python} 
235 


will print 
236 


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


of the domain with 
241 
gross 
565 
\begin{python} 
242 
ksteube 
1316 
from esys.escript import whereZero 
243 


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

246 


\code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero 
247 
jgs 
107 
and $0$ elsewhere. 
248 
gross 
565 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
249 
jgs 
107 
equal to zero and $0$ elsewhere. 
250 
gross 
565 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
251 
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. 
252 
gross 
565 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
253 


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


one gets 
257 


\begin{python} 
258 
ksteube 
1316 
0 () Function space type: Finley_Nodes on FinleyMesh 
259 
gross 
565 
\end{python} 
260 
ksteube 
1316 
An additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the 
261 
jgs 
102 
characteristic function \index{characteristic function} of the locations 
262 


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


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


\begin{python} 
265 
ksteube 
1316 
from esys.linearPDEs import Poisson 
266 


x = mydomain.getX() 
267 


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


mypde = Poisson(domain=mydomain) 
269 


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



275 
ksteube 
1316 
Now we can write the script to solve our Poisson problem 
276 
jgs 
102 
\begin{python} 
277 
ksteube 
1316 
from esys.escript import * 
278 


from esys.escript.linearPDEs import Poisson 
279 


from esys.finley import Rectangle 
280 


# generate domain: 
281 


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


# define characteristic function of Gamma^D 
283 


x = mydomain.getX() 
284 


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


# define PDE and get its solution u 
286 


mypde = Poisson(domain=mydomain) 
287 


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


u = mypde.getSolution() 
289 


# write u to an external file 
290 


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

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


\VTK file format. 
296 


Now you may run the script and visualize the solution using \mayavi: 
297 


\begin{verbatim} 
298 


python poisson.py 
299 


mayavi d u.xml m SurfaceMap 
300 


\end{verbatim} 
301 


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



303 
jgs 
102 
\begin{figure} 
304 
gross 
599 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}} 
305 
lkettle 
573 
\caption{Visualization of the Poisson Equation Solution for $f=1$} 
306 
jgs 
102 
\label{fig:FirstSteps.3} 
307 


\end{figure} 
308 


