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

trunk/esys2/doc/user/firststep.tex revision 104 by jgs, Fri Dec 17 07:43:12 2004 UTC trunk/doc/user/firststep.tex revision 565 by gross, Mon Feb 27 00:04:17 2006 UTC
# Line 1  Line 1
1  % $Id$  % $Id$
2
3  \chapter{The First Steps}  \section{The First Steps}
4  \label{FirstSteps}  \label{FirstSteps}
5
6  \begin{figure}  \begin{figure}
# Line 9  Line 9
9  \label{fig:FirstSteps.1}  \label{fig:FirstSteps.1}
10  \end{figure}  \end{figure}
11
12  In this chapter we will show the basics of how to use \escript to solve  In this chapter we will give an introduction how to use \escript to solve
13  a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). The reader should be familiar with  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}
14  the basics of Python. The knowledge presented the Python tutorial at \url{http://docs.python.org/tut/tut.html}  is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.
is sufficient. It is helpful if the reader has some basic knowledge on PDEs \index{PDE}.
15
16  The \index{PDE} we want 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}
17
18  -\Delta u =f  -\Delta u =f
19  \label{eq:FirstSteps.1}  \label{eq:FirstSteps.1}
20
21  for the solution $u$. The domain of interest which we denote by $\Omega$  for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$
22  is the unit square  is the unit square
23
24  \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 \}
25  \label{eq:FirstSteps.1b}  \label{eq:FirstSteps.1b}
26
27  The domain is shown in Figure~\fig{fig:FirstSteps.1}.  The domain is shown in \fig{fig:FirstSteps.1}.
28
29  $\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by  $\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by
30
31  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}  \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
32  \label{eq:FirstSteps.1.1}  \label{eq:FirstSteps.1.1}
33
34  where for any function $w$ and any direction $i$ $u\hackscore{,i}$  where, for any function $w$ and any direction $i$, $u\hackscore{,i}$
35  denotes the partial derivative \index{partial derivative} of $w$ with respect to $i$.    denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.
37  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
38  as $\nabla^2$, and written in the form  as $\nabla^2$, and written in the form
39  \begin{equation*}  \begin{equation*}
40  \nabla^2 u = \frac{\partial^2 u}{\partial x\hackscore 0^2}  \nabla^2 u = \nabla^t \cdot \nabla u =  \frac{\partial^2 u}{\partial x\hackscore 0^2}
41  + \frac{\partial^2 u}{\partial  x\hackscore 1^2}  + \frac{\partial^2 u}{\partial  x\hackscore 1^2}
42  \end{equation*}  \end{equation*}
43  and \eqn{eq:FirstSteps.1} as  and \eqn{eq:FirstSteps.1} as
# Line 46  and \eqn{eq:FirstSteps.1} as Line 45  and \eqn{eq:FirstSteps.1} as
45  -\nabla^2 u = f  -\nabla^2 u = f
46  \end{equation*}  \end{equation*}
47  }  }
48  Basically in the subindex of a function any index left to 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
49  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 $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
51
54
55  In some cases, and we will see examples for this in the next chapter,  In some cases, and we will see examples for this in the next chapter,
56  the usage of the nested $\sum$ symbols blows up the formulas and therefore  the usage of the nested $\sum$ symbols blows up the formulas and therefore
57  it is convenient to use Einstein summation convention \index{summation convention} which  it is convenient to use the Einstein summation convention \index{summation convention}. This
58  says that $\sum$ sign is dropped and a summation over a repeated index is performed  drops the $\sum$ sign and assumes that a summation over a repeated index is performed
59  ("repeated index means summation"). For instance we write  ("repeated index means summation"). For instance we write
60  \begin{eqnarray}  \begin{eqnarray}
61  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}   \\
62  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}   \\
63  u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\  u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\
64    x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j}   \\
65  \label{eq:FirstSteps.1.1c}  \label{eq:FirstSteps.1.1c}
66  \end{eqnarray}  \end{eqnarray}
67  With the summation convention we can write the Poisson equation \index{Poisson equation} as  With the summation convention we can write the Poisson equation \index{Poisson equation} as
# Line 80  n\hackscore{i} u\hackscore{,i}= 0 \;. Line 80  n\hackscore{i} u\hackscore{,i}= 0 \;.
80  $n=(n\hackscore{i})$ denotes the outer normal field  $n=(n\hackscore{i})$ denotes the outer normal field
81  of the domain, see \fig{fig:FirstSteps.1}. Remember that we  of the domain, see \fig{fig:FirstSteps.1}. Remember that we
82  are applying the Einstein summation convention \index{summation convention}, i.e  are applying the Einstein summation convention \index{summation convention}, i.e
83  $n\hackscore{i} u\hackscore{,i}= n\hackscore1 u\hackscore{,1} +$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
84  n\hackscore2 u\hackscore{,2}$. n\hackscore{1} u\hackscore{,1}$.
85  \footnote{Some readers may familiar with the notation  \footnote{Some readers may familiar with the notation
86  \begin{equation*}  \begin{equation*}
87  \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}  \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
# Line 93  set $\Gamma^N$ which is the top and righ Line 93  set $\Gamma^N$ which is the top and righ
93  \Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1  \}  \Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1  \}
94  \label{eq:FirstSteps.2b}  \label{eq:FirstSteps.2b}
95
96  On the bottom and the left edge of the domain which defined  On the bottom and the left edge of the domain which is defined
97  as  as
98
99  \Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0  \}  \Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0  \}
# Line 104  the solution shall be identically zero: Line 104  the solution shall be identically zero:
104  u=0 \; .  u=0 \; .
105  \label{eq:FirstSteps.2d}  \label{eq:FirstSteps.2d}
106
107  The kind of boundary condition is called a homogeneous Dirichlet boundary condition  This kind of boundary condition is called a homogeneous Dirichlet boundary condition
108  \index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together  \index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
109  with the Neumann boundary condition \eqn{eq:FirstSteps.2} and  with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
110  Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so  Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
111  called boundary value  called boundary value
112  problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for unknown  problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
113    the unknown
114  function $u$.  function $u$.
115
116
# Line 127  methods have to be used construct an app Line 128  methods have to be used construct an app
128  $u$. Here we will use the finite element method\index{finite element  $u$. Here we will use the finite element method\index{finite element
129  method} (FEM\index{finite element  method} (FEM\index{finite element
130  method!FEM}). The basic idea is to fill the domain with a  method!FEM}). The basic idea is to fill the domain with a
131  set of points, so called nodes. The solution is approximated by its  set of points called nodes. The solution is approximated by its
132  values on the nodes\index{finite element  values on the nodes\index{finite element
133  method!nodes}. Moreover, the domain is subdivide into small  method!nodes}. Moreover, the domain is subdivided into small,
134  sub-domain, so-called elements \index{finite element  sub-domain called elements \index{finite element
135  method!element}. On each element the solution is  method!element}. On each element the solution is
136  represented by a polynomial of a certain degree through its values at  represented by a polynomial of a certain degree through its values at
137  the nodes located in the element. The nodes and its connection through  the nodes located in the element. The nodes and its connection through
138  elements is called a mesh\index{finite element  elements is called a mesh\index{finite element
139  method!mesh}. Figure~\fig{fig:FirstSteps.2} shows an  method!mesh}. \fig{fig:FirstSteps.2} shows an
140  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
141  in the $x_1$ direction over the unit square.    in the $x_1$ direction over the unit square.
142  For more details we refer the reader to the literature, for instance  For more details we refer the reader to the literature, for instance
# Line 151  to solve the PDE. Here we are using the Line 152  to solve the PDE. Here we are using the
152  method}. The following statements create the \Domain object \var{mydomain} from the  method}. The following statements create the \Domain object \var{mydomain} from the
153  \finley method \method{Rectangle}  \finley method \method{Rectangle}
154  \begin{python}  \begin{python}
155  import finley  from esys.finley import Rectangle
156  mydomain = finley.Rectangle(l0=1.,l1=1.,n0=40, n1=20)  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
157  \end{python}  \end{python}
158  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
159  the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.  the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.
160  The arguments \var{l0} and \var{l1} define the number of elements in $x\hackscore{0}$ and  The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and
161  $x\hackscore{1}$-direction respectively. For more details on \method{Rectangle} and  $x\hackscore{1}$-direction respectively. For more details on \method{Rectangle} and
162  other \Domain generators within the \finley module,  other \Domain generators within the \finley module,
163  see \Chap{CHAPTER ON FINLEY}.  see \Chap{CHAPTER ON FINLEY}.
164
165  The following statements define the \Poisson object \var{mypde} with domain var{mydomain} and  The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
166  the right hand side $f$ of the PDE to constant $1$:  the right hand side $f$ of the PDE to constant $1$:
167  \begin{python}  \begin{python}
168  import escript  from esys.escript.linearPDEs import Poisson
169  mypde = escript.Poisson(domain=mydomain,f=1)  mypde = Poisson(mydomain)
170    mypde.setValue(f=1)
171  \end{python}  \end{python}
172  We have not specified any boundary condition but the  We have not specified any boundary condition but the
173  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
# Line 173  boundary condition!homogeneous} defined Line 175  boundary condition!homogeneous} defined
175  condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$  condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$
176  and any constant $C$ the function $u+C$ becomes a solution as well. We have to add  and any constant $C$ the function $u+C$ becomes a solution as well. We have to add
177  a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done  a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done
178  by defines a characteristic function \index{characteristic function}  by defining a characteristic function \index{characteristic function}
179  which has a positive values at locations $(x_0,x_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
180  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},
181  we need a function which is positive for the cases $x_0=0$ or $x_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
182    an object \var{x} which represents locations in the domain one uses
183  \begin{python}  \begin{python}
184  x=mydomain.getX()  x=mydomain.getX() \;.
185  \end{python}  \end{python}
186  In first statement returns, the method \method{getX} of the \Domain \var{mydomain} access to the locations  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
187  in the domain defined by \var{mydomain}. The object \var{x} is actually an \Data object  that \var{x} has a
which we will learn more about later. \code{x[0]} returns the $x_0$ coordinates of the locations and
\code{x[0].whereZero()} creates function which equals $1$ where \code{x[0]} is (nearly) equal to zero
and $0$ elsewhere. 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_0$ or $x_1$ is equal to zero.
188
189  The additional parameter \var{q} of the \Poisson object creater defines the  In the first statement, the method \method{getX} of the \Domain \var{mydomain}
191    in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which is
192    discussed in Chpater\ref{X} in more details. What we need to know here is that
193    \var{x} has \Rank (=number of dimensions) and a \Shape (=tuple of dimensions) which can be checked by
194    calling the \method{getRank} and \method{getShape} methods:
195    \begin{python}
196    print "rank ",x.getRank(),", shape ",x.getShape()
197    \end{python}
198    will print something like
199    \begin{python}
200    rank 1, shape (2,)
201    \end{python}
202    The \Data object also maintains type information which is represented by the
203    \FunctionSpace of the object. For instance
204    \begin{python}
205    print x.getFunctionSpace()
206    \end{python}
207    will print
208    \begin{python}
209    Function space type: Finley_Nodes on FinleyMesh
210    \end{python}
211    which tells us that the coordinates are stored on the nodes of a \finley mesh.
212    To get the  $x\hackscore{0}$ coordinates of the locations we use the
213    statement
214    \begin{python}
215    x0=x[0]
216    \end{python}
217    Object \var{x0}
218    is again a \Data object now with \Rank $0$ and
219    \Shape $()$. It inherits the \FunctionSpace from \var{x}:
220    \begin{python}
221    print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
222    \end{python}
223    will print
224    \begin{python}
225    0 () Function space type: Finley_Nodes on FinleyMesh
226    \end{python}
227    We can now construct the function \var{gammaD} by
228    \begin{python}
230    \end{python}
231    where
232    \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (allmost) equal to zero
233    and $0$ elsewhere.
234    Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is
235    equal to zero and $0$ elsewhere.
236    The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}
237    gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
238    Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
239    \begin{python}
241    \end{python}
242    one gets
243    \begin{python}
244    0 () Function space type: Finley_Nodes on FinleyMesh
245    \end{python}
246    The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the
247  characteristic function \index{characteristic function} of the locations  characteristic function \index{characteristic function} of the locations
248  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}
249  are set. The complete definition of our example is now:  are set. The complete definition of our example is now:
250  \begin{python}  \begin{python}
251  from linearPDEs import Poisson  from esys.linearPDEs import Poisson
252  x = mydomain.getX()  x = mydomain.getX()
254  mypde = Poisson(domain=mydomain,f=1,q=gammaD)  mypde = Poisson(domain=mydomain)
256  \end{python}  \end{python}
257  The first statement imports the \Poisson class definition form the \linearPDEsPack module which is part of the \escript module.  The first statement imports the \Poisson class definition form the \linearPDEs module \escript package.
258  To get the solution of the Poisson equation defines 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
259  \method{getSolution}.  \method{getSolution}.
260
261  Now we can write the script to solve our test problem (Remember that  Now we can write the script to solve our test problem (Remember that
262  lines starting with '\#' are commend lines in Python) (available as \file{mypoisson.py}  lines starting with '\#' are comment lines in Python) (available as \file{mypoisson.py}
263  in the \ExampleDirectory):  in the \ExampleDirectory):
264  \begin{python}  \begin{python}
265  import esys.finley  from esys.escript import *
266  from esys.linearPDEs import Poisson  from linearPDEs import Poisson
267    from esys.finley import Rectangle
268  # generate domain:  # generate domain:
269  mydomain = esys.finley.Rectangle(l0=1.,l1=1.,n0=40, n1=20)  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
270  # define characteristic function of Gamma^D  # define characteristic function of Gamma^D
271  x = mydomain.getX()  x = mydomain.getX()
273  # define PDE and get its solution u  # define PDE and get its solution u
274  mypde = Poisson(domain=mydomain,f=1,q=gammaD)  mypde = Poisson(domain=mydomain)
276  u = mypde.getSolution()  u = mypde.getSolution()
277  # write u to an external file  # write u to an external file
278  u.saveDX("u.dx")  saveVTK("u.xml",sol=u)
279  \end{python}  \end{python}
280  The last statement writes the solution to the external file \file{u.dx} in  The last statement writes the solution tagged with the name "sol" to the external file \file{u.xml} in
281  \OpenDX file format. \OpenDX is a software package  \VTK file format. \VTK is a software library
282  for the visualization of scientific, engineering and analytical data and is freely available  for the visualization of scientific, engineering and analytical data and is freely available
283  from \url{http://www.opendx.org}.  from \url{http://www.vtk.org}. There are a variaty of graphical user interfaces
284    for \VTK available, for instance \mayavi which can be downloaded from \url{http://mayavi.sourceforge.net/} but is also available on most
285    \LINUX distributions.
286
287  \begin{figure}  \begin{figure}
288  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}
289  \caption{\OpenDX visualization of the Possion equation soluition for $f=1$}  \caption{Visualization of the Possion Equation Solution for $f=1$}
290  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
291  \end{figure}  \end{figure}
292
293  You can edit this script using your favourite text editor (or the Integrated DeveLopment Environment IDLE  You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE
294  for Python). If the script file has the name \file{mypoisson.py} \index{scripts!\file{mypoisson.py}} you can run the  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
295  script from any shell using the command:  script from any shell using the command:
296  \begin{verbatim}  \begin{python}
297  python mypoisson.py  python mypoisson.py
298  \end{verbatim}  \end{python}
299  After the script has (hopefully successfully) been completed you will find the file \file{u.dx} in the current  After the script has (hopefully successfully) been completed you will find the file \file{u.dx} in the current
300  directory. An easy way to visualize the results is the command  directory. An easy way to visualize the results is the command
301  \begin{verbatim}  \begin{python}
302  dx -prompter  mayavi -d u.xml -m SurfaceMap &
303  \end{verbatim}  \end{python}
304  to start the generic data visualization interface of \OpenDX. \fig{fig:FirstSteps.3} shows the result.  to show the results, see \fig{fig:FirstSteps.3}.

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