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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (14 months, 2 weeks ago) by jfenwick
File MIME type: application/x-tex
File size: 20791 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ksteube 1811
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 gross 625 %
6 ksteube 1811 % Primary Business: Queensland, Australia
7 jfenwick 6112 % Licensed under the Apache License, version 2.0
8     % http://www.apache.org/licenses/LICENSE-2.0
9 gross 625 %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 % Development 2012-2013 by School of Earth Sciences
12     % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3989 %
14     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 jgs 102
16 caltinay 3274 \section{The First Steps}\label{FirstSteps}
17 ksteube 1316 In this chapter we give an introduction how to use \escript to solve
18 caltinay 3274 a partial differential equation\index{partial differential equation} (PDE\index{partial differential equation!PDE}).
19 caltinay 4891 We assume you are at least a little familiar with \PYTHON.
20 caltinay 5295 The knowledge presented in the \PYTHON tutorial at \url{https://docs.python.org/2/tutorial/} is more than sufficient.
21 jgs 102
22 caltinay 3274 The PDE\index{partial differential equation} we wish to solve is the Poisson equation\index{Poisson equation}
23 jgs 102 \begin{equation}
24 caltinay 3274 -\Delta u=f
25     \label{eq:FirstSteps.1}
26 jgs 102 \end{equation}
27 ksteube 1316 for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$,
28 jgs 102 is the unit square
29     \begin{equation}
30 jfenwick 3295 \Omega=[0,1]^2=\{ (x_0;x_1) | 0\le x_{0} \le 1 \mbox{ and } 0\le x_{1} \le 1 \}
31 jgs 102 \label{eq:FirstSteps.1b}
32     \end{equation}
33 jgs 107 The domain is shown in \fig{fig:FirstSteps.1}.
34 caltinay 3274 \begin{figure}[ht]
35 caltinay 3279 \centerline{\includegraphics{FirstStepDomain}}
36 caltinay 3274 \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
37     \label{fig:FirstSteps.1}
38 artak 1971 \end{figure}
39 jgs 102
40 ksteube 1316 $\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
41 jgs 102 \begin{equation}
42 jfenwick 3295 \Delta u = (u_{,0})_{,0}+(u_{,1})_{,1}
43 jgs 102 \label{eq:FirstSteps.1.1}
44     \end{equation}
45 jfenwick 3295 where, for any function $u$ and any direction $i$, $u_{,i}$
46 caltinay 3274 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 jgs 102 \begin{equation*}
51 jfenwick 3295 \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 jgs 102 \end{equation*}
54     and \eqn{eq:FirstSteps.1} as
55     \begin{equation*}
56 caltinay 3274 -\nabla^2 u = f
57 jgs 102 \end{equation*}
58     }
59 gross 4537 Basically, in the subindex of a function, any index to the right of the comma denotes a spatial derivative with respect
60 jfenwick 3295 to the index. To get a more compact form we will write $u_{,ij}=(u_{,i})_{,j}$
61 jgs 102 which leads to
62     \begin{equation}
63 jfenwick 3295 \Delta u = u_{,00}+u_{,11}=\sum_{i=0}^2 u_{,ii}
64 jgs 102 \label{eq:FirstSteps.1.1b}
65     \end{equation}
66 ksteube 1316 We often find that use
67     of nested $\sum$ symbols makes formulas cumbersome, and we use the more
68 jfenwick 3343 compact Einstein summation convention\index{summation convention}. This
69 ksteube 1316 drops the $\sum$ sign and assumes that a summation is performed over any repeated index.
70     For instance we write
71 jgs 102 \begin{eqnarray}
72 jfenwick 3295 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 jgs 102 \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 jfenwick 3295 - u_{,ii} =1
81 jgs 102 \label{eq:FirstSteps.1.sum}
82     \end{equation}
83 lkettle 575 where $f=1$ in this example.
84    
85 jfenwick 3295 On the boundary of the domain $\Omega$ the normal derivative $n_{i} u_{,i}$
86 caltinay 3274 of the solution $u$ shall be zero, i.e. $u$ shall fulfill
87 jgs 102 the homogeneous Neumann boundary condition\index{Neumann
88     boundary condition!homogeneous}
89     \begin{equation}
90 jfenwick 3295 n_{i} u_{,i}= 0 \;.
91 jgs 102 \label{eq:FirstSteps.2}
92     \end{equation}
93 jfenwick 3295 $n=(n_{i})$ denotes the outer normal field
94 jgs 102 of the domain, see \fig{fig:FirstSteps.1}. Remember that we
95 jfenwick 3295 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 jgs 102 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 jfenwick 3295 \Gamma^N=\{(x_0;x_1) \in \Omega | x_{0}=1 \mbox{ or } x_{1}=1 \}
103 caltinay 3274 \label{eq:FirstSteps.2b}
104 jgs 102 \end{equation}
105 jgs 107 On the bottom and the left edge of the domain which is defined
106 jgs 102 as
107     \begin{equation}
108 jfenwick 3295 \Gamma^D=\{(x_0;x_1) \in \Omega | x_{0}=0 \mbox{ or } x_{1}=0 \}
109 caltinay 3274 \label{eq:FirstSteps.2c}
110 jgs 102 \end{equation}
111 caltinay 3274 the solution shall be identical to zero:
112 jgs 102 \begin{equation}
113 caltinay 3274 u=0 \; .
114     \label{eq:FirstSteps.2d}
115 jgs 102 \end{equation}
116 caltinay 3274 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 jgs 102 with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
120 caltinay 3274 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$.
124 jgs 102
125 gross 2371 \begin{figure}[ht]
126 caltinay 3279 \centerline{\includegraphics{FirstStepMesh}}
127 caltinay 3274 \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 jgs 102 \end{figure}
133    
134 caltinay 3274 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 jgs 102 example of a FEM mesh with four elements in the $x_0$ and four elements
148 caltinay 3274 in the $x_1$ direction over the unit square.
149 gross 2370 For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.
150 jgs 102
151 caltinay 4891 The \escript solver we want to use to solve this problem is embedded into the \PYTHON interpreter language.
152 caltinay 3274 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 caltinay 5295 command\footnote{\program{run-escript} is not available under Windows.
156 caltinay 3274 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 gross 2370 \begin{verbatim}
159 caltinay 3274 run-escript
160 gross 2370 \end{verbatim}
161 caltinay 4891 which will pass you on to the \PYTHON prompt
162 gross 2370 \begin{verbatim}
163 caltinay 4891 Python 2.7.6 (default, Mar 22 2014, 15:40:47)
164     [GCC 4.8.2] on linux2
165 gross 2370 Type "help", "copyright", "credits" or "license" for more information.
166     >>>
167     \end{verbatim}
168 caltinay 4891 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 gross 2370 \begin{python}
171 caltinay 3274 >>> x=2+3
172 jfenwick 4853 >>> print("2+3=",x)
173 caltinay 3274 2+3= 5
174 gross 2370 \end{python}
175 caltinay 4891 We refer to the \PYTHON user's guide if you are not familiar with \PYTHON.
176 gross 2370
177 caltinay 3274 \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 caltinay 5295 Here we are using the FEM\index{finite element method} library \finley.
185 caltinay 3274 The following statements create the \Domain object \var{mydomain} from the
186 caltinay 5295 \finley function \method{Rectangle}:
187 jgs 102 \begin{python}
188 ksteube 1316 from esys.finley import Rectangle
189     mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
190 jgs 102 \end{python}
191 caltinay 5295 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 jfenwick 3295 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 caltinay 5295 other \Domain generators see \Chap{chap:finley}, \Chap{chap:ripley}, and
196     \Chap{chap:speckley}.
197 jgs 102
198 jgs 107 The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
199 jgs 102 the right hand side $f$ of the PDE to constant $1$:
200     \begin{python}
201 ksteube 1316 from esys.escript.linearPDEs import Poisson
202     mypde = Poisson(mydomain)
203     mypde.setValue(f=1)
204 jgs 102 \end{python}
205 caltinay 3278 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 jfenwick 3295 which has positive values at locations $x=(x_{0},x_{1})$
214 caltinay 3278 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 jfenwick 3295 construct a function \var{gammaD} which is positive for the cases $x_{0}=0$ or $x_{1}=0$.
217 caltinay 3278 To get an object \var{x} which contains the coordinates of the nodes in the domain use
218 jgs 102 \begin{python}
219 ksteube 1316 x=mydomain.getX()
220 jgs 102 \end{python}
221 caltinay 3278 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 gross 565 \begin{python}
228 jfenwick 4853 print("rank ",x.getRank(),", shape ",x.getShape())
229 gross 565 \end{python}
230 ksteube 1316 This will print something like
231 gross 565 \begin{python}
232 ksteube 1316 rank 1, shape (2,)
233 gross 565 \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 jfenwick 4853 print(x.getFunctionSpace())
238 gross 565 \end{python}
239     will print
240     \begin{python}
241 caltinay 5295 Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
242 gross 565 \end{python}
243 caltinay 3278 which tells us that the coordinates are stored on the nodes of (rather than on
244 caltinay 3331 points in the interior of) a Finley mesh.
245 jfenwick 3295 To get the $x_{0}$ coordinates of the locations we use the statement
246 gross 565 \begin{python}
247 ksteube 1316 x0=x[0]
248 gross 565 \end{python}
249 caltinay 3278 Object \var{x0} is again a \Data object now with \Rank $0$ and \Shape $()$.
250     It inherits the \FunctionSpace from \var{x}:
251 gross 565 \begin{python}
252 jfenwick 4853 print(x0.getRank(), x0.getShape(), x0.getFunctionSpace())
253 gross 565 \end{python}
254     will print
255     \begin{python}
256 caltinay 5295 0 () Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
257 gross 565 \end{python}
258 caltinay 3278 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 gross 565 \begin{python}
261 ksteube 1316 from esys.escript import whereZero
262     gammaD=whereZero(x[0])+whereZero(x[1])
263 gross 565 \end{python}
264 ksteube 1316
265 caltinay 3278 \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 jfenwick 3295 gives a function on the domain \var{mydomain} which is strictly positive where $x_{0}$ or $x_{1}$ is equal to zero.
269 caltinay 3278 Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it.
270     So from
271 gross 565 \begin{python}
272 jfenwick 4853 print(gammaD.getRank(), gammaD.getShape(), gammaD.getFunctionSpace())
273 gross 565 \end{python}
274     one gets
275     \begin{python}
276 caltinay 5295 0 () Finley_Nodes [ContinuousFunction(domain)] on FinleyMesh
277 gross 565 \end{python}
278 caltinay 3278 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 jgs 102 \begin{python}
283 caltinay 3278 from esys.escript.linearPDEs import Poisson
284 ksteube 1316 x = mydomain.getX()
285     gammaD = whereZero(x[0])+whereZero(x[1])
286     mypde = Poisson(domain=mydomain)
287     mypde.setValue(f=1,q=gammaD)
288 jgs 102 \end{python}
289 caltinay 3278 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.
291 jgs 102
292 ksteube 1316 Now we can write the script to solve our Poisson problem
293 jgs 102 \begin{python}
294 ksteube 1316 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 caltinay 4891 mypde.setValue(f=1, q=gammaD)
305 ksteube 1316 u = mypde.getSolution()
306 jgs 102 \end{python}
307 caltinay 3278 The question is what we do with the calculated solution \var{u}.
308 caltinay 5295 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}.
319 jgs 102
320 gross 2574 \subsection{Plotting Using \MATPLOTLIB}
321 caltinay 5295 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.
325 gross 2574
326 caltinay 3278 First we need to create a rectangular grid which is accomplished by the following statements:
327 gross 2574 \begin{python}
328 caltinay 3278 import numpy
329     x_grid = numpy.linspace(0., 1., 50)
330     y_grid = numpy.linspace(0., 1., 50)
331 gross 2574 \end{python}
332 caltinay 3278 \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.
335 gross 2574
336 caltinay 3278 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 gross 2574 \begin{python}
340 caltinay 3278 x=mydomain.getX()[0].toListOfTuples()
341     y=mydomain.getX()[1].toListOfTuples()
342 gross 2574 \end{python}
343 caltinay 3278 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 gross 2574 \begin{python}
348 caltinay 3278 z=interpolate(u, mydomain.getX().getFunctionSpace())
349 gross 2574 \end{python}
350 caltinay 3278 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 gross 2574 \begin{python}
353 caltinay 3278 import matplotlib
354     z_grid = matplotlib.mlab.griddata(x, y, z, xi=x_grid, yi=y_grid)
355 gross 2574 \end{python}
356 caltinay 3278 Now \var{z_grid} gives the values of the PDE solution \var{u} at the grid which can be plotted using \function{contourf}:
357 gross 2574 \begin{python}
358 caltinay 3278 matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
359     matplotlib.pyplot.savefig("u.png")
360 gross 2574 \end{python}
361 caltinay 3278 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 gross 2574 \begin{python}
364 caltinay 3278 matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
365     matplotlib.pyplot.show()
366 gross 2574 \end{python}
367     which gives an interactive browser window.
368    
369     \begin{figure}
370 caltinay 3279 \centerline{\includegraphics[width=\figwidth]{FirstStepResultMATPLOTLIB}}
371 caltinay 3278 \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB}
372 gross 2574 \label{fig:FirstSteps.3b}
373     \end{figure}
374    
375     Now we can write the script to solve our Poisson problem
376     \begin{python}
377 caltinay 3278 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 gross 3666 z=interpolate(u,mydomain.getX().getFunctionSpace()).toListOfTuples()
398 caltinay 3278 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 gross 2574 \end{python}
403 jfenwick 3295 The entire code is available as \file{poisson_matplotlib.py} in the \ExampleDirectory.
404 gross 2574 You can run the script using the {\it escript} environment
405 ksteube 1316 \begin{verbatim}
406 caltinay 3278 run-escript poisson_matplotlib.py
407 ksteube 1316 \end{verbatim}
408 caltinay 5295 This will create a file called \file{u.png}, see \fig{fig:FirstSteps.3b}.
409 gross 2574 For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
410 ksteube 1316
411 caltinay 3278 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.}.
416 gross 2574
417 gross 2580 \begin{figure}
418 caltinay 3279 \centerline{\includegraphics[width=\figwidth]{FirstStepResult}}
419 gross 2580 \caption{Visualization of the Poisson Equation Solution for $f=1$}
420     \label{fig:FirstSteps.3}
421     \end{figure}
422    
423 caltinay 3278 \subsection{Visualization using export files}
424 gross 2574
425 caltinay 3278 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 caltinay 5662 \mayavi\cite{mayavi} and \VisIt~\cite{VisIt}. This method is \MPI safe and
428 caltinay 3278 works with large 2D and 3D problems.
429 gross 2574
430 caltinay 3278 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 gross 2574 \begin{python}
433 caltinay 3348 from esys.weipa import saveVTK
434 caltinay 3278 saveVTK("u.vtu", sol=u)
435 gross 2574 \end{python}
436 caltinay 3278 This file can then be opened in a \VTK compatible visualization tool where the
437 caltinay 3348 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 caltinay 5295 will write \var{u} to a \SILO file if escript was compiled with support for
443     LLNL's \SILO library.
444 gross 2574
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 caltinay 3348 from esys.weipa import saveVTK
451 gross 2574 # 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 caltinay 3278 saveVTK("u.vtu",sol=u)
462 gross 2574 \end{python}
463 jfenwick 3295 The entire code is available as \file{poisson_vtk.py} in the \ExampleDirectory.
464 gross 2574
465 caltinay 3278 You can run the script using the {\it escript} environment and visualize the
466     solution using \mayavi:
467 gross 2574 \begin{verbatim}
468 caltinay 3279 run-escript poisson_vtk.py
469 caltinay 5296 mayavi2 -d u.vtu -m Surface
470 gross 2574 \end{verbatim}
471 caltinay 3278 The result is shown in \fig{fig:FirstSteps.3}.
472    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26