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

revision 565 by gross, Mon Feb 27 00:04:17 2006 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1
1  % $Id$
2    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3    %
4    % Copyright (c) 2003-2008 by University of Queensland
5    % Earth Systems Science Computational Center (ESSCC)
6    % http://www.uq.edu.au/esscc
7    %
8    % Primary Business: Queensland, Australia
11    %
12    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14
15  \section{The First Steps}  \section{The First Steps}
16  \label{FirstSteps}  \label{FirstSteps}
17
18  \begin{figure}  \begin{figure}
19  \centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
20  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
21  \label{fig:FirstSteps.1}  \label{fig:FirstSteps.1}
22  \end{figure}  \end{figure}
23
24  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
25  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}
26  is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.  is more than sufficient.
27
28  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}
29  \begin{equation}  \begin{equation}
30  -\Delta u =f  -\Delta u =f
31  \label{eq:FirstSteps.1}  \label{eq:FirstSteps.1}
32  \end{equation}  \end{equation}
33  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$,
34  is the unit square  is the unit square
35  \begin{equation}  \begin{equation}
36  \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 38  is the unit square
38  \end{equation}  \end{equation}
39  The domain is shown in \fig{fig:FirstSteps.1}.  The domain is shown in \fig{fig:FirstSteps.1}.
40
41  $\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by  $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
42  \begin{equation}  \begin{equation}
43  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
44  \label{eq:FirstSteps.1.1}  \label{eq:FirstSteps.1.1}
45  \end{equation}  \end{equation}
46  where, for any function $w$ and any direction $i$, $u\hackscore{,i}$  where, for any function $u$ and any direction $i$, $u\hackscore{,i}$
47  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$.
49  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
50  as $\nabla^2$, and written in the form  as $\nabla^2$, and written in the form
51  \begin{equation*}  \begin{equation*}
64  \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}
65  \label{eq:FirstSteps.1.1b}  \label{eq:FirstSteps.1.1b}
66  \end{equation}  \end{equation}
67  In some cases, and we will see examples for this in the next chapter,  We often find that use
68  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
69  it is convenient to use the Einstein summation convention \index{summation convention}. This  convenient Einstein summation convention \index{summation convention}. This
70  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.
71  ("repeated index means summation"). For instance we write  For instance we write
72  \begin{eqnarray}  \begin{eqnarray}
73  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}   \\
74  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 81  With the summation convention we can wri
81  - u\hackscore{,ii} =1  - u\hackscore{,ii} =1
82  \label{eq:FirstSteps.1.sum}  \label{eq:FirstSteps.1.sum}
83  \end{equation}  \end{equation}
84    where $f=1$ in this example.
85
86  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}$
87  of the solution $u$ shall be zero, ie. $u$ shall fulfill  of the solution $u$ shall be zero, ie. $u$ shall fulfill
88  the homogeneous Neumann boundary condition\index{Neumann  the homogeneous Neumann boundary condition\index{Neumann
# Line 115  function $u$. Line 129  function $u$.
129
130
131  \begin{figure}  \begin{figure}
132  \centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}}
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 bi-linear  the corner points. The solution is interpolated by a bi-linear
# Line 130  method} (FEM\index{finite element Line 144  method} (FEM\index{finite element
144  method!FEM}). The basic idea is to fill the domain with a  method!FEM}). The basic idea is to fill the domain with a
145  set of points called nodes. The solution is approximated by its  set of points called nodes. The solution is approximated by its
146  values on the nodes\index{finite element  values on the nodes\index{finite element
147  method!nodes}. Moreover, the domain is subdivided into small,  method!nodes}. Moreover, the domain is subdivided into smaller
148  sub-domain called elements \index{finite element  sub-domains called elements \index{finite element
149  method!element}. On each element the solution is  method!element}. On each element the solution is
150  represented by a polynomial of a certain degree through its values at  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  the nodes located in the element. The nodes and its connection through
# Line 152  to solve the PDE. Here we are using the Line 166  to solve the PDE. Here we are using the
166  method}. The following statements create the \Domain object \var{mydomain} from the  method}. The following statements create the \Domain object \var{mydomain} from the
167  \finley method \method{Rectangle}  \finley method \method{Rectangle}
168  \begin{python}  \begin{python}
169  from esys.finley import Rectangle    from esys.finley import Rectangle
170  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
171  \end{python}  \end{python}
172  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
173  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 179  see \Chap{CHAPTER ON FINLEY}.
179  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
180  the right hand side $f$ of the PDE to constant $1$:  the right hand side $f$ of the PDE to constant $1$:
181  \begin{python}  \begin{python}
182  from esys.escript.linearPDEs import Poisson    from esys.escript.linearPDEs import Poisson
183  mypde = Poisson(mydomain)    mypde = Poisson(mydomain)
184  mypde.setValue(f=1)    mypde.setValue(f=1)
185  \end{python}  \end{python}
186  We have not specified any boundary condition but the  We have not specified any boundary condition but the
187  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
# Line 179  by defining a characteristic function \i Line 193  by defining a characteristic function \i
193  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
194  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},
195  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
196  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
197  \begin{python}  \begin{python}
198  x=mydomain.getX() \;.    x=mydomain.getX()
199  \end{python}  \end{python}
200  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}
202  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
203  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
204  \var{x} has \Rank (=number of dimensions) and a \Shape (=tuple of dimensions) which can be checked by
205    \var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by
206  calling the \method{getRank} and \method{getShape} methods:  calling the \method{getRank} and \method{getShape} methods:
207  \begin{python}  \begin{python}
208  print "rank ",x.getRank(),", shape ",x.getShape()    print "rank ",x.getRank(),", shape ",x.getShape()
209  \end{python}  \end{python}
210  will print something like  This will print something like
211  \begin{python}  \begin{python}
212  rank 1, shape (2,)    rank 1, shape (2,)
213  \end{python}  \end{python}
214  The \Data object also maintains type information which is represented by the  The \Data object also maintains type information which is represented by the
215  \FunctionSpace of the object. For instance  \FunctionSpace of the object. For instance
216  \begin{python}  \begin{python}
217  print x.getFunctionSpace()    print x.getFunctionSpace()
218  \end{python}  \end{python}
219  will print  will print
220  \begin{python}  \begin{python}
221  Function space type: Finley_Nodes on FinleyMesh    Function space type: Finley_Nodes on FinleyMesh
222  \end{python}  \end{python}
223  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.
224  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
225  statement  statement
226  \begin{python}  \begin{python}
227  x0=x    x0=x
228  \end{python}  \end{python}
229  Object \var{x0}  Object \var{x0}
230  is again a \Data object now with \Rank $0$ and  is again a \Data object now with \Rank $0$ and
231  \Shape $()$. It inherits the \FunctionSpace from \var{x}:  \Shape $()$. It inherits the \FunctionSpace from \var{x}:
232  \begin{python}  \begin{python}
233  print x0.getRank(),x0.getShape(),x0.getFunctionSpace()    print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
234  \end{python}  \end{python}
235  will print  will print
236  \begin{python}  \begin{python}
237  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
238  \end{python}  \end{python}
239  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
240    of the domain with
241  \begin{python}  \begin{python}
242  gammaD=whereZero(x)+whereZero(x)    from esys.escript import whereZero
244  \end{python}  \end{python}
245  where
246  \code{whereZero(x)} creates function which equals $1$ where \code{x} is (allmost) equal to zero  \code{whereZero(x)} creates function which equals $1$ where \code{x} is (almost) equal to zero
247  and $0$ elsewhere.  and $0$ elsewhere.
248  Similarly, \code{whereZero(x)} creates function which equals $1$ where \code{x} is  Similarly, \code{whereZero(x)} creates function which equals $1$ where \code{x} is
249  equal to zero and $0$ elsewhere.  equal to zero and $0$ elsewhere.
# Line 237  The sum of the results of \code{whereZer Line 251  The sum of the results of \code{whereZer
251  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.
252  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
253  \begin{python}  \begin{python}
255  \end{python}  \end{python}
256  one gets  one gets
257  \begin{python}  \begin{python}
258  0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
259  \end{python}  \end{python}
260  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
261  characteristic function \index{characteristic function} of the locations  characteristic function \index{characteristic function} of the locations
262  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}
263  are set. The complete definition of our example is now:  are set. The complete definition of our example is now:
264  \begin{python}  \begin{python}
265  from esys.linearPDEs import Poisson    from esys.linearPDEs import Poisson
266  x = mydomain.getX()    x = mydomain.getX()
268  mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
270  \end{python}  \end{python}
271  The first statement imports the \Poisson class definition form the \linearPDEs module \escript package.  The first statement imports the \Poisson class definition from the \linearPDEs module \escript package.
272  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
273  \method{getSolution}.  \method{getSolution}.
274
275  Now we can write the script to solve our test problem (Remember that  Now we can write the script to solve our Poisson problem
276  lines starting with '\#' are comment lines in Python) (available as \file{mypoisson.py}  \begin{python}
277  in the \ExampleDirectory):    from esys.escript import *
278  \begin{python}    from esys.escript.linearPDEs import Poisson
279  from esys.escript import *    from esys.finley import Rectangle
280  from linearPDEs import Poisson    # generate domain:
281  from esys.finley import Rectangle    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
282  # generate domain:    # define characteristic function of Gamma^D
283  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    x = mydomain.getX()
284  # define characteristic function of Gamma^D    gammaD = whereZero(x)+whereZero(x)
285  x = mydomain.getX()    # define PDE and get its solution u
286  gammaD = whereZero(x)+whereZero(x)    mypde = Poisson(domain=mydomain)
287  # define PDE and get its solution u    mypde.setValue(f=1,q=gammaD)
288  mypde = Poisson(domain=mydomain)    u = mypde.getSolution()
289  mypde.setValue(f=1,q=gammaD)    # write u to an external file
290  u = mypde.getSolution()    saveVTK("u.xml",sol=u)
291  # write u to an external file  \end{python}
292  saveVTK("u.xml",sol=u)  The entire code is available as \file{poisson.py} in the \ExampleDirectory
293  \end{python}
294  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
295  \VTK file format. \VTK is a software library  \VTK file format.
296  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:
297  from \url{http://www.vtk.org}. There are a variaty of graphical user interfaces  \begin{verbatim}
298  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
299  \LINUX distributions.    mayavi -d u.xml -m SurfaceMap
300    \end{verbatim}
301    See \fig{fig:FirstSteps.3}.
302
303  \begin{figure}  \begin{figure}
304  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}
305  \caption{Visualization of the Possion Equation Solution for $f=1$}  \caption{Visualization of the Poisson Equation Solution for $f=1$}
306  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
307  \end{figure}  \end{figure}
308
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{python}
python mypoisson.py
\end{python}
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{python}
mayavi -d u.xml -m SurfaceMap &
\end{python}
to show the results, see \fig{fig:FirstSteps.3}.

Legend:
 Removed from v.565 changed lines Added in v.1811