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

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

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

revision 106 by jgs, Fri Dec 17 07:43:12 2004 UTC revision 107 by jgs, Thu Jan 27 06:21:48 2005 UTC
# 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  \begin{equation}  \begin{equation}
18  -\Delta u =f  -\Delta u =f
19  \label{eq:FirstSteps.1}  \label{eq:FirstSteps.1}
20  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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$.  
36  \footnote{Some readers  \footnote{Some readers
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 \cot \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}$
50  which leads to  which leads to
51  \begin{equation}  \begin{equation}
# Line 55  which leads to Line 54  which leads to
54  \end{equation}  \end{equation}
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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
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  \end{equation}  \end{equation}
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 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 a function which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$:  
182  \begin{python}  \begin{python}
183  x=mydomain.getX()  x=mydomain.getX()
184  gammaD=x[0].whereZero()+x[1].whereZero()  gammaD=x[0].whereZero()+x[1].whereZero()
185  \end{python}  \end{python}
186  In first statement returns, the method \method{getX} of the \Domain \var{mydomain} access to the locations  In the first statement, the method \method{getX} of the \Domain \var{mydomain}
187  in the domain defined by \var{mydomain}. The object \var{x} is actually an \Data object  gives access to locations
188  which we will learn more about later. \code{x[0]} returns the $x_0$ coordinates of the locations and  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object
189    which we will learn more about later. \code{x[0]} returns the $x\hackscore{0}$ coordinates of the locations and
190  \code{x[0].whereZero()} creates function which equals $1$ where \code{x[0]} is (nearly) equal to zero  \code{x[0].whereZero()} creates function which equals $1$ where \code{x[0]} is (nearly) equal to zero
191  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.  and $0$ elsewhere.
192    Similarly, \code{x[1].whereZero()} creates function which equals $1$ where \code{x[1]} is
193    equal to zero and $0$ elsewhere.
194    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.
195    
196  The additional parameter \var{q} of the \Poisson object creater defines the  The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the
197  characteristic function \index{characteristic function} of the locations  characteristic function \index{characteristic function} of the locations
198  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}
199  are set. The complete definition of our example is now:  are set. The complete definition of our example is now:
200  \begin{python}  \begin{python}
201  from linearPDEs import Poisson  from esys.linearPDEs import Poisson
202  x = mydomain.getX()  x = mydomain.getX()
203  gammaD = x[0].whereZero()+x[1].whereZero()  gammaD = x[0].whereZero()+x[1].whereZero()
204  mypde = Poisson(domain=mydomain,f=1,q=gammaD)  mypde = Poisson(domain=mydomain)
205    mypde = setValue(f=1,q=gammaD)
206  \end{python}  \end{python}
207  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 \linearPDEsPack module which is part of the \ESyS package.
208  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
209  \method{getSolution}.  \method{getSolution}.
210    
211  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
212  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}
213  in the \ExampleDirectory):  in the \ExampleDirectory):
214  \begin{python}  \begin{python}
215  import esys.finley  from esys.finley import Rectangle
216  from esys.linearPDEs import Poisson  from esys.linearPDEs import Poisson
217  # generate domain:  # generate domain:
218  mydomain = esys.finley.Rectangle(l0=1.,l1=1.,n0=40, n1=20)  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
219  # define characteristic function of Gamma^D  # define characteristic function of Gamma^D
220  x = mydomain.getX()  x = mydomain.getX()
221  gammaD = x[0].whereZero()+x[1].whereZero()  gammaD = x[0].whereZero()+x[1].whereZero()
# Line 225  from \url{http://www.opendx.org}. Line 232  from \url{http://www.opendx.org}.
232    
233  \begin{figure}  \begin{figure}
234  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}  \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}
235  \caption{\OpenDX visualization of the Possion equation soluition for $f=1$}  \caption{\OpenDX Visualization of the Possion Equation Solution for $f=1$}
236  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
237  \end{figure}  \end{figure}
238    
239  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
240  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
241  script from any shell using the command:  script from any shell using the command:
242  \begin{verbatim}  \begin{verbatim}
243  python mypoisson.py  python mypoisson.py
# Line 238  python mypoisson.py Line 245  python mypoisson.py
245  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
246  directory. An easy way to visualize the results is the command  directory. An easy way to visualize the results is the command
247  \begin{verbatim}  \begin{verbatim}
248  dx -prompter  dx -prompter &
249  \end{verbatim}  \end{verbatim}
250  to start the generic data visualization interface of \OpenDX. \fig{fig:FirstSteps.3} shows the result.  to start the generic data visualization interface of \OpenDX. \fig{fig:FirstSteps.3} shows the result.
   
   
   
   
Legend:
Removed from v.106  
changed lines
  Added in v.107
  ViewVC Help
Powered by ViewVC 1.1.26