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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: application/x-tex
File size: 13704 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26