/[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

trunk/esys2/doc/user/firststep.tex revision 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/doc/user/firststep.tex revision 1316 by ksteube, Tue Sep 25 03:18:30 2007 UTC
# Line 1  Line 1 
1    %
2  % $Id$  % $Id$
3    %
4    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5    %
6    %           Copyright 2003-2007 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/osl-3.0.php
13    %
14    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15    %
16    
17  \section{The First Steps}  \section{The First Steps}
18  \label{FirstSteps}  \label{FirstSteps}
19    
20  \begin{figure}  \begin{figure}
21  \centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
22  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
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 26  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 52  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 69  With the summation convention we can wri Line 83  With the summation convention we can wri
83  - u\hackscore{,ii} =1  - u\hackscore{,ii} =1
84  \label{eq:FirstSteps.1.sum}  \label{eq:FirstSteps.1.sum}
85  \end{equation}  \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}$  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  of the solution $u$ shall be zero, ie. $u$ shall fulfill
90  the homogeneous Neumann boundary condition\index{Neumann  the homogeneous Neumann boundary condition\index{Neumann
# Line 115  function $u$. Line 131  function $u$.
131    
132    
133  \begin{figure}  \begin{figure}
134  \centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}}
135  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here
136  each element is a quadrilateral and described by four nodes, namely  each element is a quadrilateral and described by four nodes, namely
137  the corner points. The solution is interpolated by a bi-linear  the corner points. The solution is interpolated by a bi-linear
# Line 130  method} (FEM\index{finite element Line 146  method} (FEM\index{finite element
146  method!FEM}). The basic idea is to fill the domain with a  method!FEM}). The basic idea is to fill the domain with a
147  set of points called nodes. The solution is approximated by its  set of points called nodes. The solution is approximated by its
148  values on the nodes\index{finite element  values on the nodes\index{finite element
149  method!nodes}. Moreover, the domain is subdivided into small,  method!nodes}. Moreover, the domain is subdivided into smaller
150  sub-domain called elements \index{finite element  sub-domains called elements \index{finite element
151  method!element}. On each element the solution is  method!element}. On each element the solution is
152  represented by a polynomial of a certain degree through its values at  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  the nodes located in the element. The nodes and its connection through
# Line 152  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 165  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 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 178  a Dirichlet boundary condition \index{Di Line 194  a Dirichlet boundary condition \index{Di
194  by defining a characteristic function \index{characteristic function}  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  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 a function which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$:    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}  \begin{python}
200  x=mydomain.getX()    x=mydomain.getX()
 gammaD=x[0].whereZero()+x[1].whereZero()  
201  \end{python}  \end{python}
202  In the first statement, the method \method{getX} of the \Domain \var{mydomain}  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  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be
205  which we will learn more about later. \code{x[0]} returns the $x\hackscore{0}$ coordinates of the locations and  discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that
206  \code{x[0].whereZero()} creates function which equals $1$ where \code{x[0]} is (nearly) equal to zero  
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 non-zero 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.  and $0$ elsewhere.
250  Similarly, \code{x[1].whereZero()} 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.
252  The sum of the results of \code{x[0].whereZero()} and \code{x[1].whereZero()} gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.  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  The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the  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  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 = x[0].whereZero()+x[1].whereZero()    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 form the \linearPDEsPack module which is part of the \ESyS 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{mypoisson.py}  \begin{python}
279  in the \ExampleDirectory):    from esys.escript import *
280  \begin{python}    from esys.escript.linearPDEs import Poisson
281  from esys.finley import Rectangle    from esys.finley import Rectangle
282  from esys.linearPDEs import Poisson    # generate domain:
283  # generate domain:    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
284  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    # define characteristic function of Gamma^D
285  # define characteristic function of Gamma^D    x = mydomain.getX()
286  x = mydomain.getX()    gammaD = whereZero(x[0])+whereZero(x[1])
287  gammaD = x[0].whereZero()+x[1].whereZero()    # define PDE and get its solution u
288  # define PDE and get its solution u    mypde = Poisson(domain=mydomain)
289  mypde = Poisson(domain=mydomain,f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
290  u = mypde.getSolution()    u = mypde.getSolution()
291  # write u to an external file    # write u to an external file
292  u.saveDX("u.dx")    saveVTK("u.xml",sol=u)
293  \end{python}  \end{python}
294  The last statement writes the solution to the external file \file{u.dx} in  The entire code is available as \file{poisson.py} in the \ExampleDirectory
295  \OpenDX file format. \OpenDX is a software package  
296  for the visualization of scientific, engineering and analytical data and is freely available  The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in
297  from \url{http://www.opendx.org}.  \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}  \begin{figure}
306  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}
307  \caption{\OpenDX Visualization of the Possion Equation Solution for $f=1$}  \caption{Visualization of the Poisson Equation Solution for $f=1$}
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{mypoisson.py} \index{scripts!\file{mypoisson.py}} you can run the  
 script from any shell using the command:  
 \begin{verbatim}  
 python mypoisson.py  
 \end{verbatim}  
 After the script has (hopefully successfully) been completed you will find the file \file{u.dx} in the current  
 directory. An easy way to visualize the results is the command  
 \begin{verbatim}  
 dx -prompter &  
 \end{verbatim}  
 to start the generic data visualization interface of \OpenDX. \fig{fig:FirstSteps.3} shows the result.  

Legend:
Removed from v.121  
changed lines
  Added in v.1316

  ViewVC Help
Powered by ViewVC 1.1.26