1 
% 
2 
% $Id$ 
3 
% 
4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
5 
% 
6 
% Copyright 20032007 by ACceSS MNRF 
7 
% Copyright 2007 by University of Queensland 
8 
% 
9 
% http://esscc.uq.edu.au 
10 
% Primary Business: Queensland, Australia 
11 
% Licensed under the Open Software License version 3.0 
12 
% http://www.opensource.org/licenses/osl3.0.php 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 
% 
16 

17 
\section{The First Steps} 
18 
\label{FirstSteps} 
19 

20 
\begin{figure} 
21 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}} 
22 
\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.} 
23 
\label{fig:FirstSteps.1} 
24 
\end{figure} 
25 

26 
In this chapter we give an introduction how to use \escript to solve 
27 
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} 
28 
is more than sufficient. 
29 

30 
The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation} 
31 
\begin{equation} 
32 
\Delta u =f 
33 
\label{eq:FirstSteps.1} 
34 
\end{equation} 
35 
for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$, 
36 
is the unit square 
37 
\begin{equation} 
38 
\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 \} 
39 
\label{eq:FirstSteps.1b} 
40 
\end{equation} 
41 
The domain is shown in \fig{fig:FirstSteps.1}. 
42 

43 
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by 
44 
\begin{equation} 
45 
\Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1} 
46 
\label{eq:FirstSteps.1.1} 
47 
\end{equation} 
48 
where, for any function $u$ and any direction $i$, $u\hackscore{,i}$ 
49 
denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$. 
50 
\footnote{You 
51 
may be more familiar with the Laplace operator\index{Laplace operator} being written 
52 
as $\nabla^2$, and written in the form 
53 
\begin{equation*} 
54 
\nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x\hackscore 0^2} 
55 
+ \frac{\partial^2 u}{\partial x\hackscore 1^2} 
56 
\end{equation*} 
57 
and \eqn{eq:FirstSteps.1} as 
58 
\begin{equation*} 
59 
\nabla^2 u = f 
60 
\end{equation*} 
61 
} 
62 
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
63 
to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$ 
64 
which leads to 
65 
\begin{equation} 
66 
\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii} 
67 
\label{eq:FirstSteps.1.1b} 
68 
\end{equation} 
69 
We often find that use 
70 
of nested $\sum$ symbols makes formulas cumbersome, and we use the more 
71 
convenient Einstein summation convention \index{summation convention}. This 
72 
drops the $\sum$ sign and assumes that a summation is performed over any repeated index. 
73 
For instance we write 
74 
\begin{eqnarray} 
75 
x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i} \\ 
76 
x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i} \\ 
77 
u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\ 
78 
x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j} \\ 
79 
\label{eq:FirstSteps.1.1c} 
80 
\end{eqnarray} 
81 
With the summation convention we can write the Poisson equation \index{Poisson equation} as 
82 
\begin{equation} 
83 
 u\hackscore{,ii} =1 
84 
\label{eq:FirstSteps.1.sum} 
85 
\end{equation} 
86 
where $f=1$ in this example. 
87 

88 
On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$ 
89 
of the solution $u$ shall be zero, ie. $u$ shall fulfill 
90 
the homogeneous Neumann boundary condition\index{Neumann 
91 
boundary condition!homogeneous} 
92 
\begin{equation} 
93 
n\hackscore{i} u\hackscore{,i}= 0 \;. 
94 
\label{eq:FirstSteps.2} 
95 
\end{equation} 
96 
$n=(n\hackscore{i})$ denotes the outer normal field 
97 
of the domain, see \fig{fig:FirstSteps.1}. Remember that we 
98 
are applying the Einstein summation convention \index{summation convention}, i.e 
99 
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} + 
100 
n\hackscore{1} u\hackscore{,1}$. 
101 
\footnote{Some readers may familiar with the notation 
102 
\begin{equation*} 
103 
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i} 
104 
\end{equation*} 
105 
for the normal derivative.} 
106 
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the 
107 
set $\Gamma^N$ which is the top and right edge of the domain: 
108 
\begin{equation} 
109 
\Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega  x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1 \} 
110 
\label{eq:FirstSteps.2b} 
111 
\end{equation} 
112 
On the bottom and the left edge of the domain which is defined 
113 
as 
114 
\begin{equation} 
115 
\Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega  x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0 \} 
116 
\label{eq:FirstSteps.2c} 
117 
\end{equation} 
118 
the solution shall be identically zero: 
119 
\begin{equation} 
120 
u=0 \; . 
121 
\label{eq:FirstSteps.2d} 
122 
\end{equation} 
123 
This kind of boundary condition is called a homogeneous Dirichlet boundary condition 
124 
\index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together 
125 
with the Neumann boundary condition \eqn{eq:FirstSteps.2} and 
126 
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so 
127 
called boundary value 
128 
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for 
129 
the unknown 
130 
function $u$. 
131 

132 

133 
\begin{figure} 
134 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}} 
135 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
136 
each element is a quadrilateral and described by four nodes, namely 
137 
the corner points. The solution is interpolated by a bilinear 
138 
polynomial.} 
139 
\label{fig:FirstSteps.2} 
140 
\end{figure} 
141 

142 
In general the BVP\index{boundary value problem!BVP} cannot be solved analytically and numerical 
143 
methods have to be used construct an approximation of the solution 
144 
$u$. Here we will use the finite element method\index{finite element 
145 
method} (FEM\index{finite element 
146 
method!FEM}). The basic idea is to fill the domain with a 
147 
set of points called nodes. The solution is approximated by its 
148 
values on the nodes\index{finite element 
149 
method!nodes}. Moreover, the domain is subdivided into smaller 
150 
subdomains called elements \index{finite element 
151 
method!element}. On each element the solution is 
152 
represented by a polynomial of a certain degree through its values at 
153 
the nodes located in the element. The nodes and its connection through 
154 
elements is called a mesh\index{finite element 
155 
method!mesh}. \fig{fig:FirstSteps.2} shows an 
156 
example of a FEM mesh with four elements in the $x_0$ and four elements 
157 
in the $x_1$ direction over the unit square. 
158 
For more details we refer the reader to the literature, for instance 
159 
\Ref{Zienc,NumHand}. 
160 

161 
\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}. 
162 
(We will discuss a more general form of a PDE \index{partial differential equation!PDE} 
163 
that can be defined through the \LinearPDE class later). The instantiation of 
164 
a \Poisson class object requires the specification of the domain $\Omega$. In \escript 
165 
the \Domain class objects are used to describe the geometry of a domain but it also 
166 
contains information about the discretization methods and the actual solver which is used 
167 
to solve the PDE. Here we are using the FEM library \finley \index{finite element 
168 
method}. The following statements create the \Domain object \var{mydomain} from the 
169 
\finley method \method{Rectangle} 
170 
\begin{python} 
171 
from esys.finley import Rectangle 
172 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
173 
\end{python} 
174 
In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and 
175 
the right, upper corner at $(\var{l0},\var{l1})=(1,1)$. 
176 
The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and 
177 
$x\hackscore{1}$direction respectively. For more details on \method{Rectangle} and 
178 
other \Domain generators within the \finley module, 
179 
see \Chap{CHAPTER ON FINLEY}. 
180 

181 
The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and 
182 
the right hand side $f$ of the PDE to constant $1$: 
183 
\begin{python} 
184 
from esys.escript.linearPDEs import Poisson 
185 
mypde = Poisson(mydomain) 
186 
mypde.setValue(f=1) 
187 
\end{python} 
188 
We have not specified any boundary condition but the 
189 
\Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann 
190 
boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary 
191 
condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$ 
192 
and any constant $C$ the function $u+C$ becomes a solution as well. We have to add 
193 
a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done 
194 
by defining a characteristic function \index{characteristic function} 
195 
which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set 
196 
and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, 
197 
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 
198 
an object \var{x} which contains the coordinates of the nodes in the domain use 
199 
\begin{python} 
200 
x=mydomain.getX() 
201 
\end{python} 
202 
The method \method{getX} of the \Domain \var{mydomain} 
203 
gives access to locations 
204 
in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be 
205 
discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that 
206 

207 
\var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by 
208 
calling the \method{getRank} and \method{getShape} methods: 
209 
\begin{python} 
210 
print "rank ",x.getRank(),", shape ",x.getShape() 
211 
\end{python} 
212 
This will print something like 
213 
\begin{python} 
214 
rank 1, shape (2,) 
215 
\end{python} 
216 
The \Data object also maintains type information which is represented by the 
217 
\FunctionSpace of the object. For instance 
218 
\begin{python} 
219 
print x.getFunctionSpace() 
220 
\end{python} 
221 
will print 
222 
\begin{python} 
223 
Function space type: Finley_Nodes on FinleyMesh 
224 
\end{python} 
225 
which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a \finley mesh. 
226 
To get the $x\hackscore{0}$ coordinates of the locations we use the 
227 
statement 
228 
\begin{python} 
229 
x0=x[0] 
230 
\end{python} 
231 
Object \var{x0} 
232 
is again a \Data object now with \Rank $0$ and 
233 
\Shape $()$. It inherits the \FunctionSpace from \var{x}: 
234 
\begin{python} 
235 
print x0.getRank(),x0.getShape(),x0.getFunctionSpace() 
236 
\end{python} 
237 
will print 
238 
\begin{python} 
239 
0 () Function space type: Finley_Nodes on FinleyMesh 
240 
\end{python} 
241 
We can now construct a function \var{gammaD} which is only nonzero on the bottom and left edges 
242 
of the domain with 
243 
\begin{python} 
244 
from esys.escript import whereZero 
245 
gammaD=whereZero(x[0])+whereZero(x[1]) 
246 
\end{python} 
247 

248 
\code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero 
249 
and $0$ elsewhere. 
250 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
251 
equal to zero and $0$ elsewhere. 
252 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
253 
gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero. 
254 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
255 
\begin{python} 
256 
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace() 
257 
\end{python} 
258 
one gets 
259 
\begin{python} 
260 
0 () Function space type: Finley_Nodes on FinleyMesh 
261 
\end{python} 
262 
An additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the 
263 
characteristic function \index{characteristic function} of the locations 
264 
of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous} 
265 
are set. The complete definition of our example is now: 
266 
\begin{python} 
267 
from esys.linearPDEs import Poisson 
268 
x = mydomain.getX() 
269 
gammaD = whereZero(x[0])+whereZero(x[1]) 
270 
mypde = Poisson(domain=mydomain) 
271 
mypde.setValue(f=1,q=gammaD) 
272 
\end{python} 
273 
The first statement imports the \Poisson class definition from the \linearPDEs module \escript package. 
274 
To get the solution of the Poisson equation defined by \var{mypde} we just have to call its 
275 
\method{getSolution}. 
276 

277 
Now we can write the script to solve our Poisson problem 
278 
\begin{python} 
279 
from esys.escript import * 
280 
from esys.escript.linearPDEs import Poisson 
281 
from esys.finley import Rectangle 
282 
# generate domain: 
283 
mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20) 
284 
# define characteristic function of Gamma^D 
285 
x = mydomain.getX() 
286 
gammaD = whereZero(x[0])+whereZero(x[1]) 
287 
# define PDE and get its solution u 
288 
mypde = Poisson(domain=mydomain) 
289 
mypde.setValue(f=1,q=gammaD) 
290 
u = mypde.getSolution() 
291 
# write u to an external file 
292 
saveVTK("u.xml",sol=u) 
293 
\end{python} 
294 
The entire code is available as \file{poisson.py} in the \ExampleDirectory 
295 

296 
The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in 
297 
\VTK file format. 
298 
Now you may run the script and visualize the solution using \mayavi: 
299 
\begin{verbatim} 
300 
python poisson.py 
301 
mayavi d u.xml m SurfaceMap 
302 
\end{verbatim} 
303 
See \fig{fig:FirstSteps.3}. 
304 

305 
\begin{figure} 
306 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}} 
307 
\caption{Visualization of the Poisson Equation Solution for $f=1$} 
308 
\label{fig:FirstSteps.3} 
309 
\end{figure} 
310 
