1 


2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
% 
4 
% $Id$ 
% Copyright (c) 20032009 by University of Queensland 
5 
% 
% Earth Systems Science Computational Center (ESSCC) 
6 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% http://www.uq.edu.au/esscc 

% 


% Copyright 20032007 by ACceSS MNRF 


% Copyright 2007 by University of Queensland 


% 


% http://esscc.uq.edu.au 


% Primary Business: Queensland, Australia 


% Licensed under the Open Software License version 3.0 


% http://www.opensource.org/licenses/osl3.0.php 

7 
% 
% 
8 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Primary Business: Queensland, Australia 
9 

% Licensed under the Open Software License version 3.0 
10 

% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
% 
12 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 


14 


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


18 
\begin{figure} 


\centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}} 


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


\label{fig:FirstSteps.1} 


\end{figure} 

19 


20 
In this chapter we give an introduction how to use \escript to solve 
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} 
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} 
33 
\label{eq:FirstSteps.1b} 
\label{eq:FirstSteps.1b} 
34 
\end{equation} 
\end{equation} 
35 
The domain is shown in \fig{fig:FirstSteps.1}. 
The domain is shown in \fig{fig:FirstSteps.1}. 
36 

\begin{figure} [ht] 
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 


42 
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by 
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by 
43 
\begin{equation} 
\begin{equation} 
59 
\end{equation*} 
\end{equation*} 
60 
} 
} 
61 
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect 
62 
to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$ 
to the index. To get a more compact form we will write $u\hackscore{,ij}=(u\hackscore {,i})\hackscore{,j}$ 
63 
which leads to 
which leads to 
64 
\begin{equation} 
\begin{equation} 
65 
\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii} 
\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii} 
98 
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} + 
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} + 
99 
n\hackscore{1} u\hackscore{,1}$. 
n\hackscore{1} u\hackscore{,1}$. 
100 
\footnote{Some readers may familiar with the notation 
\footnote{Some readers may familiar with the notation 
101 
\begin{equation*} 
$ 
102 
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i} 
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i} 
103 
\end{equation*} 
$ 
104 
for the normal derivative.} 
for the normal derivative.} 
105 
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the 
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: 
set $\Gamma^N$ which is the top and right edge of the domain: 
125 
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so 
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so 
126 
called boundary value 
called boundary value 
127 
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for 
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for 
128 
the unknown 
the unknown function~$u$. 

function $u$. 

129 


130 


131 
\begin{figure} 
\begin{figure}[ht] 
132 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}} 
\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}} 
133 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here 
134 
each element is a quadrilateral and described by four nodes, namely 
each element is a quadrilateral and described by four nodes, namely 
135 
the corner points. The solution is interpolated by a bilinear 
the corner points. The solution is interpolated by a bilinear 
153 
method!mesh}. \fig{fig:FirstSteps.2} shows an 
method!mesh}. \fig{fig:FirstSteps.2} shows an 
154 
example of a FEM mesh with four elements in the $x_0$ and four elements 
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. 
in the $x_1$ direction over the unit square. 
156 
For more details we refer the reader to the literature, for instance 
For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}. 
157 
\Ref{Zienc,NumHand}. 

158 

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 favorite 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 familiar with python. 
180 


181 
\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}. 
\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} 
(We will discuss a more general form of a PDE \index{partial differential equation!PDE} 
270 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is 
271 
equal to zero and $0$ elsewhere. 
equal to zero and $0$ elsewhere. 
272 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])} 
273 
gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero. 
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 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from 
275 
\begin{python} 
\begin{python} 
276 
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace() 
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace() 
315 


316 
The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in 
The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in 
317 
\VTK file format. 
\VTK file format. 
318 
Now you may run the script and visualize the solution using \mayavi: 
Now you may run the script using the \escript environment 
319 

and visualize the solution using \mayavi: 
320 
\begin{verbatim} 
\begin{verbatim} 
321 
python poisson.py 
escript poisson.py 
322 
mayavi d u.xml m SurfaceMap 
mayavi d u.xml m SurfaceMap 
323 
\end{verbatim} 
\end{verbatim} 
324 
See \fig{fig:FirstSteps.3}. 
See \fig{fig:FirstSteps.3}. 
325 


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