/[escript]/trunk/doc/user/firststep.tex
ViewVC logotype

Diff of /trunk/doc/user/firststep.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 625 by gross, Thu Mar 23 00:41:25 2006 UTC revision 1316 by ksteube, Tue Sep 25 03:18:30 2007 UTC
# Line 1  Line 1 
1    %
2  % $Id$  % $Id$
3  %  %
4  %           Copyright © 2006 by ACcESS MNRF  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5  %               \url{http://www.access.edu.au  %
6  %         Primary Business: Queensland, Australia.  %           Copyright 2003-2007 by ACceSS MNRF
7  %   Licensed under the Open Software License version 3.0  %       Copyright 2007 by University of Queensland
8  %      http://www.opensource.org/licenses/osl-3.0.php  %
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/osl-3.0.php
13    %
14    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15  %  %
   
16    
17  \section{The First Steps}  \section{The First Steps}
18  \label{FirstSteps}  \label{FirstSteps}
# Line 17  Line 23 
23  \label{fig:FirstSteps.1}  \label{fig:FirstSteps.1}
24  \end{figure}  \end{figure}
25    
26  In this chapter we will give an introduction how to use \escript to solve  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}). The reader should be 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}
28  is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.  is more than sufficient.
29    
30  The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}  The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}
31  \begin{equation}  \begin{equation}
32  -\Delta u =f  -\Delta u =f
33  \label{eq:FirstSteps.1}  \label{eq:FirstSteps.1}
34  \end{equation}  \end{equation}
35  for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$  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  is the unit square
37  \begin{equation}  \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 \}  \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 \}
# Line 34  is the unit square Line 40  is the unit square
40  \end{equation}  \end{equation}
41  The domain is shown in \fig{fig:FirstSteps.1}.  The domain is shown in \fig{fig:FirstSteps.1}.
42    
43  $\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by  $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
44  \begin{equation}  \begin{equation}
45  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
46  \label{eq:FirstSteps.1.1}  \label{eq:FirstSteps.1.1}
47  \end{equation}  \end{equation}
48  where, for any function $w$ and any direction $i$, $u\hackscore{,i}$  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$.    denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.  
50  \footnote{Some readers  \footnote{You
51  may be more familiar with the Laplace operator\index{Laplace operator} being written  may be more familiar with the Laplace operator\index{Laplace operator} being written
52  as $\nabla^2$, and written in the form  as $\nabla^2$, and written in the form
53  \begin{equation*}  \begin{equation*}
# Line 60  which leads to Line 66  which leads to
66  \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}
67  \label{eq:FirstSteps.1.1b}  \label{eq:FirstSteps.1.1b}
68  \end{equation}  \end{equation}
69  In some cases, and we will see examples for this in the next chapter,  We often find that use
70  the usage of the nested $\sum$ symbols blows up the formulas and therefore  of nested $\sum$ symbols makes formulas cumbersome, and we use the more
71  it is convenient to use the Einstein summation convention \index{summation convention}. This  convenient Einstein summation convention \index{summation convention}. This
72  drops the $\sum$ sign and assumes that a summation over a repeated index is performed  drops the $\sum$ sign and assumes that a summation is performed over any repeated index.
73  ("repeated index means summation"). For instance we write  For instance we write
74  \begin{eqnarray}  \begin{eqnarray}
75  x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i}   \\  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}   \\  x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i}   \\
# Line 162  to solve the PDE. Here we are using the Line 168  to solve the PDE. Here we are using the
168  method}. The following statements create the \Domain object \var{mydomain} from the  method}. The following statements create the \Domain object \var{mydomain} from the
169  \finley method \method{Rectangle}  \finley method \method{Rectangle}
170  \begin{python}  \begin{python}
171  from esys.finley import Rectangle    from esys.finley import Rectangle
172  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
173  \end{python}  \end{python}
174  In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and  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)$.  the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.
# Line 175  see \Chap{CHAPTER ON FINLEY}. Line 181  see \Chap{CHAPTER ON FINLEY}.
181  The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and  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$:  the right hand side $f$ of the PDE to constant $1$:
183  \begin{python}  \begin{python}
184  from esys.escript.linearPDEs import Poisson    from esys.escript.linearPDEs import Poisson
185  mypde = Poisson(mydomain)    mypde = Poisson(mydomain)
186  mypde.setValue(f=1)    mypde.setValue(f=1)
187  \end{python}  \end{python}
188  We have not specified any boundary condition but the  We have not specified any boundary condition but the
189  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
# Line 189  by defining a characteristic function \i Line 195  by defining a characteristic function \i
195  which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set  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},  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  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 represents locations in the domain one uses  an object \var{x} which contains the coordinates of the nodes in the domain use
199  \begin{python}  \begin{python}
200  x=mydomain.getX()    x=mydomain.getX()
201  \end{python}  \end{python}
202  In fact \var{x} is a \Data object which we will learn more about in Chapter~\ref{X}. At this stage we only have to know  The method \method{getX} of the \Domain \var{mydomain}
 that \var{x} has a  
   
 In the first statement, the method \method{getX} of the \Domain \var{mydomain}  
203  gives access to locations  gives access to locations
204  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which is  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be
205  discussed in Chpater\ref{X} in more details. What we need to know here is that  discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that
206  \var{x} has \Rank (=number of dimensions) and a \Shape (=tuple of dimensions) which can be checked by  
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:  calling the \method{getRank} and \method{getShape} methods:
209  \begin{python}  \begin{python}
210  print "rank ",x.getRank(),", shape ",x.getShape()    print "rank ",x.getRank(),", shape ",x.getShape()
211  \end{python}  \end{python}
212  will print something like  This will print something like
213  \begin{python}  \begin{python}
214  rank 1, shape (2,)    rank 1, shape (2,)
215  \end{python}  \end{python}
216  The \Data object also maintains type information which is represented by the  The \Data object also maintains type information which is represented by the
217  \FunctionSpace of the object. For instance  \FunctionSpace of the object. For instance
218  \begin{python}  \begin{python}
219  print x.getFunctionSpace()    print x.getFunctionSpace()
220  \end{python}  \end{python}
221  will print  will print
222  \begin{python}  \begin{python}
223  Function space type: Finley_Nodes on FinleyMesh    Function space type: Finley_Nodes on FinleyMesh
224  \end{python}  \end{python}
225  which tells us that the coordinates are stored on the nodes of a \finley mesh.  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  To get the  $x\hackscore{0}$ coordinates of the locations we use the
227  statement  statement
228  \begin{python}  \begin{python}
229  x0=x[0]    x0=x[0]
230  \end{python}  \end{python}
231  Object \var{x0}  Object \var{x0}
232  is again a \Data object now with \Rank $0$ and  is again a \Data object now with \Rank $0$ and
233  \Shape $()$. It inherits the \FunctionSpace from \var{x}:  \Shape $()$. It inherits the \FunctionSpace from \var{x}:
234  \begin{python}  \begin{python}
235  print x0.getRank(),x0.getShape(),x0.getFunctionSpace()    print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
236  \end{python}  \end{python}
237  will print  will print
238  \begin{python}  \begin{python}
239  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
240  \end{python}  \end{python}
241  We can now construct the function \var{gammaD} by  We can now construct a function \var{gammaD} which is only non-zero on the bottom and left edges
242    of the domain with
243  \begin{python}  \begin{python}
244  from esys.escript import whereZero    from esys.escript import whereZero
245  gammaD=whereZero(x[0])+whereZero(x[1])    gammaD=whereZero(x[0])+whereZero(x[1])
246  \end{python}  \end{python}
247  where  
248  \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (allmost) equal to zero  \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero
249  and $0$ elsewhere.  and $0$ elsewhere.
250  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
251  equal to zero and $0$ elsewhere.  equal to zero and $0$ elsewhere.
# Line 248  The sum of the results of \code{whereZer Line 253  The sum of the results of \code{whereZer
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.  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  Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
255  \begin{python}  \begin{python}
256  print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()    print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
257  \end{python}  \end{python}
258  one gets  one gets
259  \begin{python}  \begin{python}
260  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
261  \end{python}  \end{python}
262  The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the  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  characteristic function \index{characteristic function} of the locations
264  of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}  of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
265  are set. The complete definition of our example is now:  are set. The complete definition of our example is now:
266  \begin{python}  \begin{python}
267  from esys.linearPDEs import Poisson    from esys.linearPDEs import Poisson
268  x = mydomain.getX()    x = mydomain.getX()
269  gammaD = whereZero(x[0])+whereZero(x[1])    gammaD = whereZero(x[0])+whereZero(x[1])
270  mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
271  mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
272  \end{python}  \end{python}
273  The first statement imports the \Poisson class definition from the \linearPDEs module \escript package.  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  To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
275  \method{getSolution}.  \method{getSolution}.
276    
277  Now we can write the script to solve our test problem (Remember that  Now we can write the script to solve our Poisson problem
278  lines starting with '\#' are comment lines in Python) (available as \file{poisson.py}  \begin{python}
279  in the \ExampleDirectory):    from esys.escript import *
280  \begin{python}    from esys.escript.linearPDEs import Poisson
281  from esys.escript import *    from esys.finley import Rectangle
282  from esys.escript.linearPDEs import Poisson    # generate domain:
283  from esys.finley import Rectangle    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
284  # generate domain:    # define characteristic function of Gamma^D
285  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    x = mydomain.getX()
286  # define characteristic function of Gamma^D    gammaD = whereZero(x[0])+whereZero(x[1])
287  x = mydomain.getX()    # define PDE and get its solution u
288  gammaD = whereZero(x[0])+whereZero(x[1])    mypde = Poisson(domain=mydomain)
289  # define PDE and get its solution u    mypde.setValue(f=1,q=gammaD)
290  mypde = Poisson(domain=mydomain)    u = mypde.getSolution()
291  mypde.setValue(f=1,q=gammaD)    # write u to an external file
292  u = mypde.getSolution()    saveVTK("u.xml",sol=u)
293  # write u to an external file  \end{python}
294  saveVTK("u.xml",sol=u)  The entire code is available as \file{poisson.py} in the \ExampleDirectory
295  \end{python}  
296  The last statement writes the solution tagged with the name "sol" to the external file \file{u.xml} in  The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in
297  \VTK file format. \VTK is a software library  \VTK file format.
298  for the visualization of scientific, engineering and analytical data and is freely available  Now you may run the script and visualize the solution using \mayavi:
299  from \url{http://www.vtk.org}. There are a variety of graphical user interfaces  \begin{verbatim}
300  for \VTK available, for instance \mayavi which can be downloaded from \url{http://mayavi.sourceforge.net/} but is also available on most    python poisson.py
301  \LINUX distributions.    mayavi -d u.xml -m SurfaceMap
302    \end{verbatim}
303    See \fig{fig:FirstSteps.3}.
304    
305  \begin{figure}  \begin{figure}
306  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}
# Line 301  for \VTK available, for instance \mayavi Line 308  for \VTK available, for instance \mayavi
308  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
309  \end{figure}  \end{figure}
310    
 You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE  
 for Python, see \url{http://idlefork.sourceforge.net}). If the script file has the name \file{poisson.py} \index{scripts!\file{poisson.py}} you can run the  
 script from any shell using the command:  
 \begin{python}  
 python poisson.py  
 \end{python}  
 After the script has (hopefully successfully) been completed you will find the file \file{u.xml} in the current  
 directory. An easy way to visualize the results is the command  
 \begin{python}  
 mayavi -d u.xml -m SurfaceMap &  
 \end{python}  
 to show the results, see \fig{fig:FirstSteps.3}.  
Legend:
Removed from v.625  
changed lines
  Added in v.1316
  ViewVC Help
Powered by ViewVC 1.1.26