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

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


 1 ksteube 1811 2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland 4 jfenwick 3989 5 gross 625 % 6 ksteube 1811 % Primary Business: Queensland, Australia 7 jfenwick 6112 % Licensed under the Apache License, version 2.0 8 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