ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 5662 - (show annotations)
Thu Jun 18 04:40:34 2015 UTC (3 years, 11 months ago) by caltinay
File MIME type: application/x-tex
File size: 20801 byte(s)
Fixes #259 - replaced all references to mayavi by Mayavi2 in the docs.
The arguments were already correct.

2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2015 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Open Software License version 3.0
8 % http://www.opensource.org/licenses/osl-3.0.php
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development 2012-2013 by School of Earth Sciences
12 % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 %
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16 \section{The First Steps}\label{FirstSteps}
17 In this chapter we give an introduction how to use \escript to solve
18 a partial differential equation\index{partial differential equation} (PDE\index{partial differential equation!PDE}).
19 We assume you are at least a little familiar with \PYTHON.
20 The knowledge presented in the \PYTHON tutorial at \url{https://docs.python.org/2/tutorial/} is more than sufficient.
22 The PDE\index{partial differential equation} we wish to solve is the Poisson equation\index{Poisson equation}
23 \begin{equation}
24 -\Delta u=f
25 \label{eq:FirstSteps.1}
26 \end{equation}
27 for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$,
28 is the unit square
29 \begin{equation}
30 \Omega=[0,1]^2=\{ (x_0;x_1) | 0\le x_{0} \le 1 \mbox{ and } 0\le x_{1} \le 1 \}
31 \label{eq:FirstSteps.1b}
32 \end{equation}
33 The domain is shown in \fig{fig:FirstSteps.1}.
34 \begin{figure}[ht]
35 \centerline{\includegraphics{FirstStepDomain}}
36 \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
37 \label{fig:FirstSteps.1}
38 \end{figure}
40 $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
41 \begin{equation}
42 \Delta u = (u_{,0})_{,0}+(u_{,1})_{,1}
43 \label{eq:FirstSteps.1.1}
44 \end{equation}
45 where, for any function $u$ and any direction $i$, $u_{,i}$
46 denotes the partial derivative \index{partial derivative} of $u$ with respect
47 to $i$.\footnote{You may be more familiar with the Laplace
48 operator\index{Laplace operator} being written as $\nabla^2$, and written in
49 the form
50 \begin{equation*}
51 \nabla^2 u = \nabla^t \cdot \nabla u = \frac{\partial^2 u}{\partial x_0^2}
52 + \frac{\partial^2 u}{\partial x_1^2}
53 \end{equation*}
54 and \eqn{eq:FirstSteps.1} as
55 \begin{equation*}
56 -\nabla^2 u = f
57 \end{equation*}
58 }
59 Basically, in the subindex of a function, any index to the right of the comma denotes a spatial derivative with respect
60 to the index. To get a more compact form we will write $u_{,ij}=(u_{,i})_{,j}$
61 which leads to
62 \begin{equation}
63 \Delta u = u_{,00}+u_{,11}=\sum_{i=0}^2 u_{,ii}
64 \label{eq:FirstSteps.1.1b}
65 \end{equation}
66 We often find that use
67 of nested $\sum$ symbols makes formulas cumbersome, and we use the more
68 compact Einstein summation convention\index{summation convention}. This
69 drops the $\sum$ sign and assumes that a summation is performed over any repeated index.
70 For instance we write
71 \begin{eqnarray}
72 x_{i}y_{i}=\sum_{i=0}^2 x_{i}y_{i} \\
73 x_{i}u_{,i}=\sum_{i=0}^2 x_{i}u_{,i} \\
74 u_{,ii}=\sum_{i=0}^2 u_{,ii} \\
75 x_{ij}u_{i,j}=\sum_{j=0}^2\sum_{i=0}^2 x_{ij}u_{i,j} \\
76 \label{eq:FirstSteps.1.1c}
77 \end{eqnarray}
78 With the summation convention we can write the Poisson equation \index{Poisson equation} as
79 \begin{equation}
80 - u_{,ii} =1
81 \label{eq:FirstSteps.1.sum}
82 \end{equation}
83 where $f=1$ in this example.
85 On the boundary of the domain $\Omega$ the normal derivative $n_{i} u_{,i}$
86 of the solution $u$ shall be zero, i.e. $u$ shall fulfill
87 the homogeneous Neumann boundary condition\index{Neumann
88 boundary condition!homogeneous}
89 \begin{equation}
90 n_{i} u_{,i}= 0 \;.
91 \label{eq:FirstSteps.2}
92 \end{equation}
93 $n=(n_{i})$ denotes the outer normal field
94 of the domain, see \fig{fig:FirstSteps.1}. Remember that we
95 are applying the Einstein summation convention \index{summation convention}, i.e. $n_{i} u_{,i}= n_{0} u_{,0} +%
96 n_{1} u_{,1}$.\footnote{Some readers may familiar with the
97 notation $\frac{\partial u}{\partial n} = n_{i} u_{,i}$
98 for the normal derivative.}
99 The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
100 set $\Gamma^N$ which is the top and right edge of the domain:
101 \begin{equation}
102 \Gamma^N=\{(x_0;x_1) \in \Omega | x_{0}=1 \mbox{ or } x_{1}=1 \}
103 \label{eq:FirstSteps.2b}
104 \end{equation}
105 On the bottom and the left edge of the domain which is defined
106 as
107 \begin{equation}
108 \Gamma^D=\{(x_0;x_1) \in \Omega | x_{0}=0 \mbox{ or } x_{1}=0 \}
109 \label{eq:FirstSteps.2c}
110 \end{equation}
111 the solution shall be identical to zero:
112 \begin{equation}
113 u=0 \; .
114 \label{eq:FirstSteps.2d}
115 \end{equation}
116 This kind of boundary condition is called a homogeneous Dirichlet boundary
117 condition\index{Dirichlet boundary condition!homogeneous}.
118 The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
119 with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
120 Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so-called
121 boundary value
122 problem\index{boundary value problem} (BVP\index{boundary value problem!BVP})
123 for the unknown function~$u$.
125 \begin{figure}[ht]
126 \centerline{\includegraphics{FirstStepMesh}}
127 \caption{Mesh of $4 \times 4$ elements on a rectangular domain. Here
128 each element is a quadrilateral and described by four nodes, namely
129 the corner points. The solution is interpolated by a bi-linear
130 polynomial.}
131 \label{fig:FirstSteps.2}
132 \end{figure}
134 In general the BVP\index{boundary value problem!BVP} cannot be solved
135 analytically and numerical methods have to be used to construct an
136 approximation of the solution $u$.
137 Here we will use the finite element method\index{finite element method}
138 (FEM\index{finite element method!FEM}).
139 The basic idea is to fill the domain with a set of points called nodes.
140 The solution is approximated by its values on the nodes\index{finite element method!nodes}.
141 Moreover, the domain is subdivided into smaller sub-domains called
142 elements\index{finite element method!element}.
143 On each element the solution is represented by a polynomial of a certain
144 degree through its values at the nodes located in the element.
145 The nodes and their connection through elements is called a
146 mesh\index{finite element method!mesh}. \fig{fig:FirstSteps.2} shows an
147 example of a FEM mesh with four elements in the $x_0$ and four elements
148 in the $x_1$ direction over the unit square.
149 For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.
151 The \escript solver we want to use to solve this problem is embedded into the \PYTHON interpreter language.
152 So you can solve the problem interactively but you will learn quickly that it
153 is more efficient to use scripts which you can edit with your favorite editor.
154 To enter the escript environment, use the \program{run-escript}
155 command\footnote{\program{run-escript} is not available under Windows.
156 If you run under Windows you can just use the \program{python} command and the
157 \env{OMP_NUM_THREADS} environment variable to control the number of threads.}:
158 \begin{verbatim}
159 run-escript
160 \end{verbatim}
161 which will pass you on to the \PYTHON prompt
162 \begin{verbatim}
163 Python 2.7.6 (default, Mar 22 2014, 15:40:47)
164 [GCC 4.8.2] on linux2
165 Type "help", "copyright", "credits" or "license" for more information.
166 >>>
167 \end{verbatim}
168 Here you can use all available \PYTHON commands and language features\footnote{Throughout our examples, we use the python 3 form of
169 print. That is, print(1) instead of print 1.}, for instance
170 \begin{python}
171 >>> x=2+3
172 >>> print("2+3=",x)
173 2+3= 5
174 \end{python}
175 We refer to the \PYTHON user's guide if you are not familiar with \PYTHON.
177 \escript provides the class \Poisson to define a Poisson equation\index{Poisson equation}.
178 (We will discuss a more general form of a PDE\index{partial differential equation!PDE}
179 that can be defined through the \LinearPDE class later.)
180 The instantiation of a \Poisson class object requires the specification of the domain $\Omega$.
181 In \escript the \Domain class objects are used to describe the geometry of a
182 domain but it also contains information about the discretization methods and
183 the actual solver which is used to solve the PDE.
184 Here we are using the FEM\index{finite element method} library \finley.
185 The following statements create the \Domain object \var{mydomain} from the
186 \finley function \method{Rectangle}:
187 \begin{python}
188 from esys.finley import Rectangle
189 mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
190 \end{python}
191 In this case the domain is a rectangle with the lower left corner at point $(0,0)$
192 and the upper right corner at $(\var{l0},\var{l1})=(1,1)$.
193 The arguments \var{n0} and \var{n1} define the number of elements in $x_{0}$ and
194 $x_{1}$-direction respectively. For more details on \method{Rectangle} and
195 other \Domain generators see \Chap{chap:finley}, \Chap{chap:ripley}, and
196 \Chap{chap:speckley}.
198 The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
199 the right hand side $f$ of the PDE to constant $1$:
200 \begin{python}
201 from esys.escript.linearPDEs import Poisson
202 mypde = Poisson(mydomain)
203 mypde.setValue(f=1)
204 \end{python}
205 We have not specified any boundary condition but the \Poisson class implicitly
206 assumes homogeneous Neuman boundary conditions\index{Neumann boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}.
207 With this boundary condition the BVP\index{boundary value problem!BVP} we have
208 defined has no unique solution.
209 In fact, with any solution $u$ and any constant $C$ the function $u+C$ becomes
210 a solution as well.
211 We have to add a Dirichlet boundary condition\index{Dirichlet boundary condition}.
212 This is done by defining a characteristic function\index{characteristic function}
213 which has positive values at locations $x=(x_{0},x_{1})$
214 where Dirichlet boundary condition is set and $0$ elsewhere.
215 In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, we need to
216 construct a function \var{gammaD} which is positive for the cases $x_{0}=0$ or $x_{1}=0$.
217 To get an object \var{x} which contains the coordinates of the nodes in the domain use
218 \begin{python}
219 x=mydomain.getX()
220 \end{python}
221 The method \method{getX} of the \Domain \var{mydomain} gives access to locations
222 in the domain defined by \var{mydomain}.
223 The object \var{x} is actually a \Data object which will be discussed in
224 \Chap{ESCRIPT CHAP} in more detail.
225 What we need to know here is that \var{x} has \Rank (number of dimensions) and
226 a \Shape (list of dimensions) which can be viewed by calling the \method{getRank} and \method{getShape} methods:
227 \begin{python}
228 print("rank ",x.getRank(),", shape ",x.getShape())
229 \end{python}
230 This will print something like
231 \begin{python}
232 rank 1, shape (2,)
233 \end{python}
234 The \Data object also maintains type information which is represented by the
235 \FunctionSpace of the object. For instance
236 \begin{python}
237 print(x.getFunctionSpace())
238 \end{python}
239 will print
240 \begin{python}
241 Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
242 \end{python}
243 which tells us that the coordinates are stored on the nodes of (rather than on
244 points in the interior of) a Finley mesh.
245 To get the $x_{0}$ coordinates of the locations we use the statement
246 \begin{python}
247 x0=x[0]
248 \end{python}
249 Object \var{x0} is again a \Data object now with \Rank $0$ and \Shape $()$.
250 It inherits the \FunctionSpace from \var{x}:
251 \begin{python}
252 print(x0.getRank(), x0.getShape(), x0.getFunctionSpace())
253 \end{python}
254 will print
255 \begin{python}
256 0 () Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
257 \end{python}
258 We can now construct a function \var{gammaD} which is only non-zero on the
259 bottom and left edges of the domain with
260 \begin{python}
261 from esys.escript import whereZero
262 gammaD=whereZero(x[0])+whereZero(x[1])
263 \end{python}
265 \code{whereZero(x[0])} creates a function which equals $1$ where \code{x[0]} is (almost) equal to zero and $0$ elsewhere.
266 Similarly, \code{whereZero(x[1])} creates a function which equals $1$ where \code{x[1]} is equal to zero and $0$ elsewhere.
267 The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}
268 gives a function on the domain \var{mydomain} which is strictly positive where $x_{0}$ or $x_{1}$ is equal to zero.
269 Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it.
270 So from
271 \begin{python}
272 print(gammaD.getRank(), gammaD.getShape(), gammaD.getFunctionSpace())
273 \end{python}
274 one gets
275 \begin{python}
276 0 () Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
277 \end{python}
278 An additional parameter \var{q} of the \code{setValue} method of the \Poisson
279 class defines the characteristic function\index{characteristic function} of
280 the locations of the domain where the homogeneous Dirichlet boundary condition\index{Dirichlet boundary condition!homogeneous} is set.
281 The complete definition of our example is now:
282 \begin{python}
283 from esys.escript.linearPDEs import Poisson
284 x = mydomain.getX()
285 gammaD = whereZero(x[0])+whereZero(x[1])
286 mypde = Poisson(domain=mydomain)
287 mypde.setValue(f=1,q=gammaD)
288 \end{python}
289 The first statement imports the \Poisson class definition from the \linearPDEs module.
290 To get the solution of the Poisson equation defined by \var{mypde} we just have to call its \method{getSolution} method.
292 Now we can write the script to solve our Poisson problem
293 \begin{python}
294 from esys.escript import *
295 from esys.escript.linearPDEs import Poisson
296 from esys.finley import Rectangle
297 # generate domain:
298 mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
299 # define characteristic function of Gamma^D
300 x = mydomain.getX()
301 gammaD = whereZero(x[0])+whereZero(x[1])
302 # define PDE and get its solution u
303 mypde = Poisson(domain=mydomain)
304 mypde.setValue(f=1, q=gammaD)
305 u = mypde.getSolution()
306 \end{python}
307 The question is what we do with the calculated solution \var{u}.
308 Besides postprocessing, e.g. calculating the gradient or the average value,
309 which will be discussed later, plotting the solution is one of the things you
310 might want to do.
311 \escript offers two ways to do this, both based on external modules or packages.
312 The first option is using the \MATPLOTLIB module which allows plotting 2D
313 results relatively quickly from within the \PYTHON script, see~\cite{matplotlib}.
314 However, there are limitations when using this tool, especially for large
315 problems and when solving three-dimensional problems.
316 Therefore, \escript provides functionality to export data as files which can
317 subsequently be read by third-party software packages such as
318 \mayavi\cite{mayavi} or \VisIt~\cite{VisIt}.
320 \subsection{Plotting Using \MATPLOTLIB}
321 The \MATPLOTLIB module provides a simple and easy-to-use way to visualize PDE
322 solutions (or other \Data objects).
323 To hand over data from \escript to \MATPLOTLIB the values need to be mapped onto
324 a rectangular grid. We will make use of the \numpy module for this.
326 First we need to create a rectangular grid which is accomplished by the following statements:
327 \begin{python}
328 import numpy
329 x_grid = numpy.linspace(0., 1., 50)
330 y_grid = numpy.linspace(0., 1., 50)
331 \end{python}
332 \var{x_grid} is an array defining the x coordinates of the grid while
333 \var{y_grid} defines the y coordinates of the grid.
334 In this case we use $50$ points over the interval $[0,1]$ in both directions.
336 Now the values created by \escript need to be interpolated to this grid.
337 We will use the \MATPLOTLIB \function{mlab.griddata} function to do this.
338 Spatial coordinates are easily extracted as a \var{list} by
339 \begin{python}
340 x=mydomain.getX()[0].toListOfTuples()
341 y=mydomain.getX()[1].toListOfTuples()
342 \end{python}
343 In principle we can apply the same \member{toListOfTuples} method to extract the values from the PDE solution \var{u}.
344 However, we have to make sure that the \Data object we extract the values from
345 uses the same \FunctionSpace as we have used when extracting \var{x} and \var{y}.
346 We apply the \function{interpolation} to \var{u} before extraction to achieve this:
347 \begin{python}
348 z=interpolate(u, mydomain.getX().getFunctionSpace())
349 \end{python}
350 The values in \var{z} are the values at the points with the coordinates given by \var{x} and \var{y}.
351 These values are interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
352 \begin{python}
353 import matplotlib
354 z_grid = matplotlib.mlab.griddata(x, y, z, xi=x_grid, yi=y_grid)
355 \end{python}
356 Now \var{z_grid} gives the values of the PDE solution \var{u} at the grid which can be plotted using \function{contourf}:
357 \begin{python}
358 matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
359 matplotlib.pyplot.savefig("u.png")
360 \end{python}
361 Here we use 5 contours. The last statement writes the plot to the file \file{u.png} in the PNG format.
362 Alternatively, one can use
363 \begin{python}
364 matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
365 matplotlib.pyplot.show()
366 \end{python}
367 which gives an interactive browser window.
369 \begin{figure}
370 \centerline{\includegraphics[width=\figwidth]{FirstStepResultMATPLOTLIB}}
371 \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB}
372 \label{fig:FirstSteps.3b}
373 \end{figure}
375 Now we can write the script to solve our Poisson problem
376 \begin{python}
377 from esys.escript import *
378 from esys.escript.linearPDEs import Poisson
379 from esys.finley import Rectangle
380 import numpy
381 import matplotlib
382 import pylab
383 # generate domain:
384 mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
385 # define characteristic function of Gamma^D
386 x = mydomain.getX()
387 gammaD = whereZero(x[0])+whereZero(x[1])
388 # define PDE and get its solution u
389 mypde = Poisson(domain=mydomain)
390 mypde.setValue(f=1,q=gammaD)
391 u = mypde.getSolution()
392 # interpolate u to a matplotlib grid:
393 x_grid = numpy.linspace(0.,1.,50)
394 y_grid = numpy.linspace(0.,1.,50)
395 x=mydomain.getX()[0].toListOfTuples()
396 y=mydomain.getX()[1].toListOfTuples()
397 z=interpolate(u,mydomain.getX().getFunctionSpace()).toListOfTuples()
398 z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
399 # interpolate u to a rectangular grid:
400 matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
401 matplotlib.pyplot.savefig("u.png")
402 \end{python}
403 The entire code is available as \file{poisson_matplotlib.py} in the \ExampleDirectory.
404 You can run the script using the {\it escript} environment
405 \begin{verbatim}
406 run-escript poisson_matplotlib.py
407 \end{verbatim}
408 This will create a file called \file{u.png}, see \fig{fig:FirstSteps.3b}.
409 For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
411 As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
412 should be used for small problems only.
413 It can not be used under \MPI as the \member{toListOfTuples} method is not
414 safe under \MPI\footnote{The phrase 'safe under \MPI' means that a program
415 will produce correct results when run on more than one processor under \MPI.}.
417 \begin{figure}
418 \centerline{\includegraphics[width=\figwidth]{FirstStepResult}}
419 \caption{Visualization of the Poisson Equation Solution for $f=1$}
420 \label{fig:FirstSteps.3}
421 \end{figure}
423 \subsection{Visualization using export files}
425 As an alternative to \MATPLOTLIB, {\it escript} supports exporting data to
426 \VTK and \SILO files which can be read by visualization tools such as
427 \mayavi\cite{mayavi} and \VisIt~\cite{VisIt}. This method is \MPI safe and
428 works with large 2D and 3D problems.
430 To write the solution \var{u} of the Poisson problem in the \VTK file format
431 to the file \file{u.vtu} one needs to add:
432 \begin{python}
433 from esys.weipa import saveVTK
434 saveVTK("u.vtu", sol=u)
435 \end{python}
436 This file can then be opened in a \VTK compatible visualization tool where the
437 solution is accessible by the name {\it sol}. Similarly,
438 \begin{python}
439 from esys.weipa import saveSilo
440 saveSilo("u.silo", sol=u)
441 \end{python}
442 will write \var{u} to a \SILO file if escript was compiled with support for
443 LLNL's \SILO library.
445 The Poisson problem script is now
446 \begin{python}
447 from esys.escript import *
448 from esys.escript.linearPDEs import Poisson
449 from esys.finley import Rectangle
450 from esys.weipa import saveVTK
451 # generate domain:
452 mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
453 # define characteristic function of Gamma^D
454 x = mydomain.getX()
455 gammaD = whereZero(x[0])+whereZero(x[1])
456 # define PDE and get its solution u
457 mypde = Poisson(domain=mydomain)
458 mypde.setValue(f=1,q=gammaD)
459 u = mypde.getSolution()
460 # write u to an external file
461 saveVTK("u.vtu",sol=u)
462 \end{python}
463 The entire code is available as \file{poisson_vtk.py} in the \ExampleDirectory.
465 You can run the script using the {\it escript} environment and visualize the
466 solution using \mayavi:
467 \begin{verbatim}
468 run-escript poisson_vtk.py
469 mayavi2 -d u.vtu -m Surface
470 \end{verbatim}
471 The result is shown in \fig{fig:FirstSteps.3}.


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

  ViewVC Help
Powered by ViewVC 1.1.26