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

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.

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2015 by The University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 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. 21 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} 39 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. 84 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$. 124 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} 133 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}. 150 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. 176 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}. 197 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} 264 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. 291 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}. 319 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. 325 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. 335 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. 368 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} 374 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}. 410 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.}. 416 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} 422 423 \subsection{Visualization using export files} 424 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. 429 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. 444 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. 464 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}. 472