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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 3 months ago) by trankine
File MIME type: application/x-tex
File size: 13727 byte(s)
And get the *(&(*&(* name right
1 ksteube 1316 %
2 jgs 102 % $Id$
3 gross 625 %
4 ksteube 1316 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 gross 625 %
6 ksteube 1316 % 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 jgs 102
17 jgs 121 \section{The First Steps}
18 jgs 102 \label{FirstSteps}
19    
20     \begin{figure}
21 gross 599 \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
22 jgs 102 \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
23     \label{fig:FirstSteps.1}
24     \end{figure}
25    
26 ksteube 1316 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}). 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 more than sufficient.
29 jgs 102
30 jgs 107 The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}
31 jgs 102 \begin{equation}
32     -\Delta u =f
33     \label{eq:FirstSteps.1}
34     \end{equation}
35 ksteube 1316 for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$,
36 jgs 102 is the unit square
37     \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 \}
39     \label{eq:FirstSteps.1b}
40     \end{equation}
41 jgs 107 The domain is shown in \fig{fig:FirstSteps.1}.
42 jgs 102
43 ksteube 1316 $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
44 jgs 102 \begin{equation}
45     \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
46     \label{eq:FirstSteps.1.1}
47     \end{equation}
48 ksteube 1316 where, for any function $u$ and any direction $i$, $u\hackscore{,i}$
49 jgs 107 denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.
50 ksteube 1316 \footnote{You
51 jgs 102 may be more familiar with the Laplace operator\index{Laplace operator} being written
52     as $\nabla^2$, and written in the form
53     \begin{equation*}
54 jgs 110 \nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x\hackscore 0^2}
55 jgs 102 + \frac{\partial^2 u}{\partial x\hackscore 1^2}
56     \end{equation*}
57     and \eqn{eq:FirstSteps.1} as
58     \begin{equation*}
59     -\nabla^2 u = f
60     \end{equation*}
61     }
62 jgs 107 Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect
63 jgs 102 to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
64     which leads to
65     \begin{equation}
66     \Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
67     \label{eq:FirstSteps.1.1b}
68     \end{equation}
69 ksteube 1316 We often find that use
70     of nested $\sum$ symbols makes formulas cumbersome, and we use the more
71     convenient Einstein summation convention \index{summation convention}. This
72     drops the $\sum$ sign and assumes that a summation is performed over any repeated index.
73     For instance we write
74 jgs 102 \begin{eqnarray}
75     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} \\
77     u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\
78 jgs 107 x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j} \\
79 jgs 102 \label{eq:FirstSteps.1.1c}
80     \end{eqnarray}
81     With the summation convention we can write the Poisson equation \index{Poisson equation} as
82     \begin{equation}
83     - u\hackscore{,ii} =1
84     \label{eq:FirstSteps.1.sum}
85     \end{equation}
86 lkettle 575 where $f=1$ in this example.
87    
88 jgs 102 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
90     the homogeneous Neumann boundary condition\index{Neumann
91     boundary condition!homogeneous}
92     \begin{equation}
93     n\hackscore{i} u\hackscore{,i}= 0 \;.
94     \label{eq:FirstSteps.2}
95     \end{equation}
96     $n=(n\hackscore{i})$ denotes the outer normal field
97     of the domain, see \fig{fig:FirstSteps.1}. Remember that we
98     are applying the Einstein summation convention \index{summation convention}, i.e
99 jgs 107 $n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
100     n\hackscore{1} u\hackscore{,1}$.
101 jgs 102 \footnote{Some readers may familiar with the notation
102     \begin{equation*}
103     \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
104     \end{equation*}
105     for the normal derivative.}
106     The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
107     set $\Gamma^N$ which is the top and right edge of the domain:
108     \begin{equation}
109     \Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1 \}
110     \label{eq:FirstSteps.2b}
111     \end{equation}
112 jgs 107 On the bottom and the left edge of the domain which is defined
113 jgs 102 as
114     \begin{equation}
115     \Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0 \}
116     \label{eq:FirstSteps.2c}
117     \end{equation}
118     the solution shall be identically zero:
119     \begin{equation}
120     u=0 \; .
121     \label{eq:FirstSteps.2d}
122     \end{equation}
123 jgs 107 This kind of boundary condition is called a homogeneous Dirichlet boundary condition
124 jgs 102 \index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
125     with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
126     Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
127     called boundary value
128 jgs 107 problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
129     the unknown
130 jgs 102 function $u$.
131    
132    
133     \begin{figure}
134 gross 599 \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}}
135 jgs 102 \caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here
136     each element is a quadrilateral and described by four nodes, namely
137     the corner points. The solution is interpolated by a bi-linear
138     polynomial.}
139     \label{fig:FirstSteps.2}
140     \end{figure}
141    
142     In general the BVP\index{boundary value problem!BVP} cannot be solved analytically and numerical
143     methods have to be used construct an approximation of the solution
144     $u$. Here we will use the finite element method\index{finite element
145     method} (FEM\index{finite element
146     method!FEM}). The basic idea is to fill the domain with a
147 jgs 107 set of points called nodes. The solution is approximated by its
148 jgs 102 values on the nodes\index{finite element
149 lkettle 573 method!nodes}. Moreover, the domain is subdivided into smaller
150     sub-domains called elements \index{finite element
151 jgs 102 method!element}. On each element the solution is
152     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
154     elements is called a mesh\index{finite element
155 jgs 107 method!mesh}. \fig{fig:FirstSteps.2} shows an
156 jgs 102 example of a FEM mesh with four elements in the $x_0$ and four elements
157     in the $x_1$ direction over the unit square.
158     For more details we refer the reader to the literature, for instance
159     \Ref{Zienc,NumHand}.
160    
161     \escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.
162     (We will discuss a more general form of a PDE \index{partial differential equation!PDE}
163     that can be defined through the \LinearPDE class later). The instantiation of
164     a \Poisson class object requires the specification of the domain $\Omega$. In \escript
165     the \Domain class objects are used to describe the geometry of a domain but it also
166     contains information about the discretization methods and the actual solver which is used
167     to solve the PDE. Here we are using the FEM library \finley \index{finite element
168     method}. The following statements create the \Domain object \var{mydomain} from the
169     \finley method \method{Rectangle}
170     \begin{python}
171 ksteube 1316 from esys.finley import Rectangle
172     mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
173 jgs 102 \end{python}
174     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)$.
176 jgs 107 The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and
177 jgs 102 $x\hackscore{1}$-direction respectively. For more details on \method{Rectangle} and
178     other \Domain generators within the \finley module,
179     see \Chap{CHAPTER ON FINLEY}.
180    
181 jgs 107 The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
182 jgs 102 the right hand side $f$ of the PDE to constant $1$:
183     \begin{python}
184 ksteube 1316 from esys.escript.linearPDEs import Poisson
185     mypde = Poisson(mydomain)
186     mypde.setValue(f=1)
187 jgs 102 \end{python}
188     We have not specified any boundary condition but the
189     \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
190     boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary
191     condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$
192     and any constant $C$ the function $u+C$ becomes a solution as well. We have to add
193     a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done
194 jgs 107 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
196 jgs 102 and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c},
197 gross 565 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 ksteube 1316 an object \var{x} which contains the coordinates of the nodes in the domain use
199 jgs 102 \begin{python}
200 ksteube 1316 x=mydomain.getX()
201 jgs 102 \end{python}
202 gross 660 The method \method{getX} of the \Domain \var{mydomain}
203 jgs 107 gives access to locations
204 ksteube 1316 in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be
205     discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that
206 gross 660
207 ksteube 1316 \var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by
208 gross 565 calling the \method{getRank} and \method{getShape} methods:
209     \begin{python}
210 ksteube 1316 print "rank ",x.getRank(),", shape ",x.getShape()
211 gross 565 \end{python}
212 ksteube 1316 This will print something like
213 gross 565 \begin{python}
214 ksteube 1316 rank 1, shape (2,)
215 gross 565 \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 ksteube 1316 print x.getFunctionSpace()
220 gross 565 \end{python}
221     will print
222     \begin{python}
223 ksteube 1316 Function space type: Finley_Nodes on FinleyMesh
224 gross 565 \end{python}
225 ksteube 1316 which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a \finley mesh.
226 gross 565 To get the $x\hackscore{0}$ coordinates of the locations we use the
227     statement
228     \begin{python}
229 ksteube 1316 x0=x[0]
230 gross 565 \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 ksteube 1316 print x0.getRank(),x0.getShape(),x0.getFunctionSpace()
236 gross 565 \end{python}
237     will print
238     \begin{python}
239 ksteube 1316 0 () Function space type: Finley_Nodes on FinleyMesh
240 gross 565 \end{python}
241 ksteube 1316 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 gross 565 \begin{python}
244 ksteube 1316 from esys.escript import whereZero
245     gammaD=whereZero(x[0])+whereZero(x[1])
246 gross 565 \end{python}
247 ksteube 1316
248     \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero
249 jgs 107 and $0$ elsewhere.
250 gross 565 Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is
251 jgs 107 equal to zero and $0$ elsewhere.
252 gross 565 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     Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
255     \begin{python}
256 ksteube 1316 print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
257 gross 565 \end{python}
258     one gets
259     \begin{python}
260 ksteube 1316 0 () Function space type: Finley_Nodes on FinleyMesh
261 gross 565 \end{python}
262 ksteube 1316 An additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the
263 jgs 102 characteristic function \index{characteristic function} of the locations
264     of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
265     are set. The complete definition of our example is now:
266     \begin{python}
267 ksteube 1316 from esys.linearPDEs import Poisson
268     x = mydomain.getX()
269     gammaD = whereZero(x[0])+whereZero(x[1])
270     mypde = Poisson(domain=mydomain)
271     mypde.setValue(f=1,q=gammaD)
272 jgs 102 \end{python}
273 lkettle 573 The first statement imports the \Poisson class definition from the \linearPDEs module \escript package.
274 jgs 107 To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
275 jgs 102 \method{getSolution}.
276    
277 ksteube 1316 Now we can write the script to solve our Poisson problem
278 jgs 102 \begin{python}
279 ksteube 1316 from esys.escript import *
280     from esys.escript.linearPDEs import Poisson
281     from esys.finley import Rectangle
282     # generate domain:
283     mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
284     # define characteristic function of Gamma^D
285     x = mydomain.getX()
286     gammaD = whereZero(x[0])+whereZero(x[1])
287     # define PDE and get its solution u
288     mypde = Poisson(domain=mydomain)
289     mypde.setValue(f=1,q=gammaD)
290     u = mypde.getSolution()
291     # write u to an external file
292     saveVTK("u.xml",sol=u)
293 jgs 102 \end{python}
294 ksteube 1316 The entire code is available as \file{poisson.py} in the \ExampleDirectory
295 jgs 102
296 ksteube 1316 The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in
297     \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 jgs 102 \begin{figure}
306 gross 599 \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}
307 lkettle 573 \caption{Visualization of the Poisson Equation Solution for $f=1$}
308 jgs 102 \label{fig:FirstSteps.3}
309     \end{figure}
310    

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26