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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6923 - (show annotations)
Mon Dec 2 05:57:07 2019 UTC (5 months, 3 weeks ago) by aellery
File MIME type: application/x-tex
File size: 89629 byte(s)
- Amended the userguide to include ComplexScalar, ComplexVector etc.
- Fixed a bug in run_comm1 & run_comm4 (failure if escript was compiled without scipy)
- Temporarily removed the ability to interpolate from ReducedFunction to Function in dudley and finley


1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2018 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Apache License, version 2.0
8 % http://www.apache.org/licenses/LICENSE-2.0
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 \chapter{The \escript Module}\label{ESCRIPT CHAP}
17
18 \section{Concepts}
19 \escript is a \PYTHON module that allows you to represent the values of
20 a function at points in a \Domain in such a way that the function will
21 be useful for the Finite Element Method (FEM) simulation. It also
22 provides what we call a function space that describes how the data is
23 used in the simulation. Stored along with the data is information
24 about the elements and nodes which will be used by the domain (e.g. \finley).
25
26 \subsection{Function spaces}
27 In order to understand what we mean by the term 'function space',
28 consider that the solution of a partial differential
29 equation\index{partial differential equation} (PDE) is a function on a domain
30 $\Omega$. When solving a PDE using FEM, the solution is
31 piecewise-differentiable but, in general, its gradient is discontinuous.
32 To reflect these different degrees of smoothness, different function spaces
33 are used.
34 For instance, in FEM, the displacement field is represented by its values at
35 the nodes of the mesh, and so is continuous.
36 The strain, which is the symmetric part of the gradient of the displacement
37 field, is stored on the element centers, and so is considered to be
38 discontinuous.
39
40 A function space is described by a \FunctionSpace object.
41 The following statement generates the object \var{solution_space} which is
42 a \FunctionSpace object and provides access to the function space of
43 PDE solutions on the \Domain \var{mydomain}:
44
45 \begin{python}
46 solution_space=Solution(mydomain)
47 \end{python}
48 The following generators for function spaces on a \Domain \var{mydomain} are commonly used:
49 \begin{itemize}
50 \item \var{Solution(mydomain)}: solutions of a PDE
51 \item \var{ReducedSolution(mydomain)}: solutions of a PDE with a reduced
52 smoothness requirement, e.g. using a lower order approximation on the same
53 element or using macro elements\index{macro elements}
54 \item \var{ContinuousFunction(mydomain)}: continuous functions, e.g. a temperature distribution
55 \item \var{Function(mydomain)}: general functions which are not necessarily continuous, e.g. a stress field
56 \item \var{FunctionOnBoundary(mydomain)}: functions on the boundary of the domain, e.g. a surface pressure
57 \item \var{DiracDeltaFunctions(mydomain)}: functions defined on a set of points
58 \item \var{FunctionOnContact0(mydomain)}: functions on side $0$ of a discontinuity
59 \item \var{FunctionOnContact1(mydomain)}: functions on side $1$ of a discontinuity
60 \end{itemize}
61 In some cases under-integration is used. For these cases the user may use a
62 \FunctionSpace from the following list:
63 \begin{itemize}
64 \item \var{ReducedFunction(mydomain)}
65 \item \var{ReducedFunctionOnBoundary(mydomain)}
66 \item \var{ReducedFunctionOnContact0(mydomain)}
67 \item \var{ReducedFunctionOnContact1(mydomain)}
68 \end{itemize}
69 In comparison to the corresponding full version they use a reduced number of
70 integration nodes (typically one only) to represent values.
71
72 \begin{figure}
73 \centering
74 \scalebox{0.97}{\includegraphics{EscriptDiagram1}}
75 \caption{\label{ESCRIPT DEP}Dependency of function spaces in \finley.
76 An arrow indicates that a function in the \FunctionSpace at the starting point
77 can be interpolated to the \FunctionSpace of the arrow target.
78 All function spaces above the dotted line can be interpolated to any of
79 the function spaces below the line. See also \Sec{SEC Projection}.}
80 \end{figure}
81
82 The reduced smoothness for a PDE solution is often used to fulfill the
83 Ladyzhenskaya-Babuska-Brezzi condition~\cite{LBB} when solving saddle point
84 problems\index{saddle point problems}, e.g. the Stokes equation.
85 A discontinuity\index{discontinuity} is a region within the domain across
86 which functions may be discontinuous.
87 The location of a discontinuity is defined in the \Domain object.
88 \fig{ESCRIPT DEP} shows the dependency between the types of function spaces
89 in \finley (other libraries may have different relationships).
90
91 The solution of a PDE is a continuous function. Any continuous function can
92 be seen as a general function on the domain and can be restricted to the
93 boundary as well as to one side of a discontinuity (the result will be
94 different depending on which side is chosen). Functions on any side of the
95 discontinuity can be seen as a function on the corresponding other side.
96
97 A function on the boundary or on one side of the discontinuity cannot be seen
98 as a general function on the domain as there are no values defined for the
99 interior. For most PDE solver libraries the space of the solution and
100 continuous functions is identical, however in some cases, for example when
101 periodic boundary conditions are used in \finley, a solution fulfills periodic
102 boundary conditions while a continuous function does not have to be periodic.
103
104 The concept of function spaces describes the properties of functions and
105 allows abstraction from the actual representation of the function in the
106 context of a particular application. For instance, in the FEM context a
107 function of the \Function type (written as \emph{Function()} in \fig{ESCRIPT DEP})
108 is usually represented by its values at the element center,
109 but in a finite difference scheme the edge midpoint of cells is preferred.
110 By changing its function space you can use the same function in a Finite
111 Difference scheme instead of Finite Element scheme.
112 Changing the function space of a particular function will typically lead to
113 a change of its representation.
114 So, when seen as a general function, a continuous function which is typically
115 represented by its values on the nodes of the FEM mesh or finite difference
116 grid must be interpolated to the element centers or the cell edges,
117 respectively. Interpolation happens automatically in \escript whenever it is
118 required\index{interpolation}. The user needs to be aware that an
119 interpolation is not always possible, see \fig{ESCRIPT DEP} for \finley.
120 An alternative approach to change the representation (=\FunctionSpace) is
121 projection\index{projection}, see \Sec{SEC Projection}.
122
123 \subsection{\Data Objects}
124 In \escript the class that stores these functions is called \Data.
125 The function is represented through its values on \DataSamplePoints where
126 the \DataSamplePoints are chosen according to the function space of the
127 function.
128 \Data class objects are used to define the coefficients of the PDEs to be
129 solved by a PDE solver library and also to store the solutions of the PDE.
130
131 The values of the function have a rank which gives the number of indices,
132 and a \Shape defining the range of each index.
133 The rank in \escript is limited to the range 0 through 4 and it is assumed
134 that the rank and \Shape is the same for all \DataSamplePoints.
135 The \Shape of a \Data object is a tuple (list) \var{s} of integers.
136 The length of \var{s} is the rank of the \Data object and the \var{i}-th
137 index ranges between 0 and $\var{s[i]}-1$.
138 For instance, a stress field has rank 2 and \Shape $(d,d)$ where $d$ is the
139 number of spatial dimensions.
140 The following statement creates the \Data object \var{mydat} representing a
141 continuous function with values of \Shape $(2,3)$ and rank $2$:
142 \begin{python}
143 mydat=Data(value=1, what=ContinuousFunction(myDomain), shape=(2,3))
144 \end{python}
145 The initial value is the constant 1 for all \DataSamplePoints and all
146 components.
147
148 \Data objects can also be created from any \numpy array or any object, such
149 as a list of floating point numbers, that can be converted into
150 a \numpyNDA\cite{NUMPY}.
151 The following two statements create objects which are equivalent
152 to \var{mydat}:
153 \begin{python}
154 mydat1=Data(value=numpy.ones((2,3)), what=ContinuousFunction(myDomain))
155 mydat2=Data(value=[[1,1], [1,1], [1,1]], what=ContinuousFunction(myDomain))
156 \end{python}
157 In the first case the initial value is \var{numpy.ones((2,3))} which generates
158 a $2 \times 3$ matrix as an instance of \numpyNDA filled with ones.
159 The \Shape of the created \Data object is taken from the \Shape of the array.
160 In the second case, the creator converts the initial value, which is a list of
161 lists, into a \numpyNDA before creating the actual \Data object.
162
163 For convenience \escript provides creators for the most common types
164 of \Data objects in the following forms (\var{d} defines the spatial
165 dimensionality):
166 \begin{itemize}
167 \item \code{Scalar(0, Function(mydomain))} is the same as \code{Data(0, Function(myDomain),(,))}\\
168 (each value is a scalar), e.g. a temperature field
169 \item \code{Vector(0, Function(mydomain))} is the same as \code{Data(0, Function(myDomain),(d,))}\\
170 (each value is a vector), e.g. a velocity field
171 \item \code{Tensor(0, Function(mydomain))} equals \code{Data(0, Function(myDomain), (d,d))},
172 e.g. a stress field
173 \item \code{Tensor4(0,Function(mydomain))} equals \code{Data(0,Function(myDomain), (d,d,d,d))},
174 e.g. a Hook tensor field
175 \item \code{ComplexScalar(0+0j, Function(mydomain))} is the same as \\
176 \code{ComplexData(0+0j, Function(myDomain),(,))}
177 (each value is a complex scalar), e.g. a temperature field
178 \item \code{ComplexVector(0+0j, Function(mydomain))} is the same as \\
179 \code{ComplexData(0+0j, Function(myDomain),(d,))}
180 (each value is a complex vector), e.g. a velocity field
181 \item \code{ComplexTensor(0+0j, Function(mydomain))} is the same as \\
182 \code{ComplexData(0+0j, Function(myDomain), (d,d))}, e.g. a stress field
183 \item \code{ComplexTensor4(0+0j,Function(mydomain))} is the same as \\
184 \code{ComplexData(0+0j,Function(myDomain), (d,d,d,d))}, e.g. a Hook tensor field
185 \end{itemize}
186 Here the initial value is 0 but any object that can be converted into
187 a \numpyNDA and whose \Shape is consistent with \Shape of the \Data object to
188 be created can be used as the initial value.
189
190 \Data objects can be manipulated by applying unary operations (e.g. cos, sin,
191 log), and they can be combined point-wise by applying arithmetic operations
192 (e.g. +, - ,* , /).
193 We emphasize that \escript itself does not handle any spatial dependencies as
194 it does not know how values are interpreted by the processing PDE solver library.
195 However \escript invokes interpolation if this is needed during data manipulations.
196 Typically, this occurs in binary operations when the arguments belong to
197 different function spaces or when data are handed over to a PDE solver library
198 which requires functions to be represented in a particular way.
199
200 The following example shows the usage of \Data objects. Assume we have a
201 displacement field $u$ and we want to calculate the corresponding stress field
202 $\sigma$ using the linear-elastic isotropic material model
203 \begin{eqnarray}\label{eq: linear elastic stress}
204 \sigma_{ij}=\lambda u_{k,k} \delta_{ij} + \mu ( u_{i,j} + u_{j,i})
205 \end{eqnarray}
206 where $\delta_{ij}$ is the Kronecker symbol and
207 $\lambda$ and $\mu$ are the Lam\'e coefficients. The following function
208 takes the displacement \var{u} and the Lam\'e coefficients \var{lam} and \var{mu}
209 as arguments and returns the corresponding stress:
210 \begin{python}
211 from esys.escript import *
212 def getStress(u, lam, mu):
213 d=u.getDomain().getDim()
214 g=grad(u)
215 stress=lam*trace(g)*kronecker(d)+mu*(g+transpose(g))
216 return stress
217 \end{python}
218 The variable \var{d} gives the spatial dimensionality of the domain on which
219 the displacements are defined.
220 The \code{kronecker(d)} call, returns the Kronecker symbol with indices $i$ and $j$ running
221 from 0 to \var{d}-1.
222 The \var{grad(u)} call, requires the displacement field \var{u} to be in
223 the \var{Solution} or \ContinuousFunction.
224 The result \var{g} as well as the returned stress will be in the \Function.
225 If, for example, \var{u} is the solution of a PDE then \code{getStress} might
226 be called in the following way:
227 \begin{python}
228 s=getStress(u, 1., 2.)
229 \end{python}
230 However \code{getStress} can also be called with \Data objects as values for
231 \var{lam} and \var{mu} which, for instance in the case of a temperature
232 dependency, are calculated by an expression.
233 The following call is equivalent to the previous example:
234 \begin{python}
235 lam=Scalar(1., ContinuousFunction(mydomain))
236 mu=Scalar(2., Function(mydomain))
237 s=getStress(u, lam, mu)
238 \end{python}
239 %
240 The function \var{lam} belongs to the \ContinuousFunction but with \var{g} the
241 function \var{trace(g)} is in the \Function.
242 In the evaluation of the product \var{lam*trace(g)} we have different function
243 spaces (on the nodes versus in the centers) and at first glance we have incompatible data.
244 \escript converts the arguments into an appropriate function space according
245 to \fig{ESCRIPT DEP}.
246 In this example that means \escript sees \var{lam} as a function of the \Function.
247 In the context of FEM this means the nodal values of \var{lam} are
248 interpolated to the element centers.
249 The interpolation is automatic and requires no special handling.
250
251 \begin{figure}
252 \centering
253 \includegraphics{EscriptDiagram2}
254 \caption{\label{Figure: tag}Element Tagging. A rectangular mesh over a region
255 with two rock types {\it white} and {\it gray} is shown.
256 The number in each cell refers to the major rock type present in the cell
257 ($1$ for {\it white} and $2$ for {\it gray}).}
258 \end{figure}
259
260 \subsection{Tagged, Expanded and Constant Data}
261 Material parameters such as the Lam\'e coefficients are typically dependent on
262 rock types present in the area of interest.
263 A common technique to handle these kinds of material parameters is
264 \emph{tagging}\index{tagging}, which uses storage efficiently.
265 \fig{Figure: tag} shows an example. In this case two rock types {\it white}
266 and {\it gray} can be found in the domain.
267 The domain is subdivided into triangular shaped cells.
268 Each cell has a tag indicating the rock type predominantly found in this cell.
269 Here $1$ is used to indicate rock type {\it white} and $2$ for rock type {\it gray}.
270 The tags are assigned at the time when the cells are generated and stored in
271 the \Domain class object. To allow easier usage of tags, names can be used
272 instead of numbers. These names are typically defined at the time when the
273 geometry is generated.
274
275 The following statements show how to use tagged values for \var{lam} as shown
276 in \fig{Figure: tag} for the stress calculation discussed above:
277 \begin{python}
278 lam=Scalar(value=2., what=Function(mydomain))
279 insertTaggedValue(lam, white=30., gray=5000.)
280 s=getStress(u, lam, 2.)
281 \end{python}
282 In this example \var{lam} is set to $30$ for those cells with tag {\it white}
283 (=$1$) and to $5000$ for cells with tag {\it gray} (=$2$).
284 The initial value $2$ of \var{lam} is used as a default value for the case
285 when a tag is encountered which has not been linked with a value.
286 The \code{getStress} method does not need to be changed now that we are using tags.
287 \escript resolves the tags when \var{lam*trace(g)} is calculated.
288
289 This brings us to a very important point about \escript.
290 You can develop a simulation with constant Lam\'e coefficients, and then later
291 switch to tagged Lam\'e coefficients without otherwise changing your \PYTHON script.
292 In short, you can use the same script for models with different domains and
293 different types of input data.
294
295 There are three main ways in which \Data objects are represented internally --
296 constant, tagged, and expanded.
297 In the constant case, the same value is used at each sample point while only a
298 single value is stored to save memory.
299 In the expanded case, each sample point has an individual value (such as for the solution of a PDE).
300 This is where your largest data sets will be created because the values are
301 stored as a complete array.
302 The tagged case has already been discussed above.
303 Expanded data is created when specifying \code{expanded=True} in the \Data
304 object constructor, while tagged data requires calling the \member{insertTaggedValue}
305 method as shown above.
306
307 Values are accessed through a sample reference number.
308 Operations on expanded \Data objects have to be performed for each sample
309 point individually.
310 When tagged values are used, the values are held in a dictionary.
311 Operations on tagged data require processing the set of tagged values only,
312 rather than processing the value for each individual sample point.
313 \escript allows any mixture of constant, tagged and expanded data in a single expression.
314
315 \subsection{Saving and Restoring Simulation Data}
316 \Data objects can be written to disk files with the \member{dump} method and
317 read back using the \member{load} method, both of which use the
318 \netCDF\cite{NETCDF} file format.
319 Use these to save data for checkpoint/restart or simply to save and reuse data
320 that was expensive to compute.
321 For instance, to save the coordinates of the data points of a
322 \ContinuousFunction to the file \file{x.nc} use
323 \begin{python}
324 x=ContinuousFunction(mydomain).getX()
325 x.dump("x.nc")
326 mydomain.dump("dom.nc")
327 \end{python}
328 To recover the object \var{x}, and you know that \var{mydomain} was an \finley
329 mesh, use
330 \begin{python}
331 from esys.finley import LoadMesh
332 mydomain=LoadMesh("dom.nc")
333 x=load("x.nc", mydomain)
334 \end{python}
335 Obviously, it is possible to execute the same steps that were originally used
336 to generate \var{mydomain} to recreate it. However, in most cases using
337 \member{dump} and \member{load} is faster, particularly if optimization has
338 been applied.
339 If \escript is running on more than one \MPI process \member{dump} will create
340 an individual file for each process containing the local data.
341 In order to avoid conflicts the \MPI processor
342 rank is appended to the file names.
343 That is instead of one file \file{dom.nc} you would get
344 \file{dom.nc.0000}, \file{dom.nc.0001}, etc.
345 You still call \code{LoadMesh("dom.nc")} to load the domain but you have to
346 make sure that the appropriate file is accessible from the corresponding rank,
347 and loading will only succeed if you run with as many processes as were used
348 when calling \member{dump}.
349
350 The function space of the \Data is stored in \file{x.nc}.
351 If the \Data object is expanded, the number of data points in the file and of
352 the \Domain for the particular \FunctionSpace must match.
353 Moreover, the ordering of the values is checked using the reference
354 identifiers provided by the \FunctionSpace on the \Domain.
355 In some cases, data points will be reordered so be aware and confirm that you
356 get what you wanted.
357
358 A more flexible way of saving and restoring \escript simulation data
359 is through an instance of the \class{DataManager} class.
360 It has the advantage of allowing to save and load not only a \Domain and
361 \Data objects but also other values\footnote{The \PYTHON \emph{pickle} module
362 is used for other types.} you compute in your simulation script.
363 Further, \class{DataManager} objects can simultaneously create files for
364 visualization so no extra calls to \code{saveVTK} etc. are needed.
365
366 The following example shows how the \class{DataManager} class can be used.
367 For an explanation of all member functions and options see the class reference
368 Section \ref{sec:datamanager}.
369 \begin{python}
370 from esys.escript import DataManager, Scalar, Function
371 from esys.finley import Rectangle
372
373 dm = DataManager(formats=[DataManager.RESTART, DataManager.VTK])
374 if dm.hasData():
375 mydomain=dm.getDomain()
376 val=dm.getValue("val")
377 t=dm.getValue("t")
378 t_max=dm.getValue("t_max")
379 else:
380 mydomain=Rectangle()
381 val=Function(mydomain).getX()
382 t=0.
383 t_max=2.5
384
385 while t<t_max:
386 t+=.01
387 val=val+t/2
388 dm.addData(val=val, t=t, t_max=t_max)
389 dm.export()
390 \end{python}
391 In the constructor we specify that we want \code{RESTART} (i.e. dump) files
392 and \code{VTK} files to be saved.
393 By default, the constructor will look for previously saved \code{RESTART}
394 files under the current directory and load them.
395 We can then enquire if such files were found by calling the \member{hasData}
396 method. If it returns \True we retrieve the domain and values into local
397 variables. Otherwise the same variables are initialized with appropriate
398 values to start a new simulation.
399 Note, that \var{t} and \var{t_max} are regular floating point values and not
400 \Data objects. Yet they are treated the same way by the \class{DataManager}.
401
402 After this initialization step the script enters the main simulation loop
403 where calculations are performed.
404 When these are finalized for a time step we call the \member{addData} method
405 to let the manager know which variables to store on disk.
406 This does not actually save the data yet and it is allowed to call
407 \member{addData} more than once to add information incrementally, e.g. from
408 separate functions that have access to the \class{DataManager} instance.
409 Once all variables have been added the \member{export} method has to be called
410 to flush all data to disk and clear the manager.
411 In this example, this call dumps \var{mydomain} and \var{val} to files
412 in a restart directory and also stores \var{t} and \var{t_max} on disk.
413 Additionally, it generates a \VTK file for visualization of the data.
414 If the script would stop running before its completion for some reason (e.g.
415 because its runtime limit was exceeded in a batch job environment), you could
416 simply run it again and it would resume at the point it stopped before.
417
418 \section{\escript Classes}
419
420 \subsection{The \Domain class}
421 \begin{classdesc}{Domain}{}
422 A \Domain object is used to describe a geometric region together with
423 a way of representing functions over this region.
424 The \Domain class provides an abstract interface to the domain of \FunctionSpace and \Data objects.
425 \Domain needs to be subclassed in order to provide a complete implementation.
426 \end{classdesc}
427
428 \vspace{1em}\noindent The following methods are available:
429 \begin{methoddesc}[Domain]{getDim}{}
430 returns the number of spatial dimensions of the \Domain.
431 \end{methoddesc}
432 %
433 \begin{methoddesc}[Domain]{dump}{filename}
434 writes the \Domain to the file \var{filename} using the \netCDF file format.
435 \end{methoddesc}
436 %
437 \begin{methoddesc}[Domain]{getX}{}
438 returns the locations in the \Domain. The \FunctionSpace of the returned
439 \Data object is chosen by the \Domain implementation. Typically it will be
440 in the \ContinuousFunction.
441 \end{methoddesc}
442
443 \begin{methoddesc}[Domain]{getNumpyX}{}
444 returns the locations in the \Domain as a \numpy ndarray. The \FunctionSpace
445 of the returned \Data object is chosen by the \Domain implementation.
446 Typically it will be in the \ContinuousFunction.
447
448 Note that it is necessary to load \numpy first in the escript.
449 \end{methoddesc}
450 %
451 \begin{methoddesc}[Domain]{setX}{newX}
452 assigns new locations to the \Domain. \var{newX} has to have \Shape $(d,)$
453 where $d$ is the spatial dimensionality of the domain. Typically \var{newX}
454 must be in the \ContinuousFunction but the space actually to be used
455 depends on the \Domain implementation. Not all domain families support
456 setting locations.
457 \end{methoddesc}
458 %
459 \begin{methoddesc}[Domain]{getNormal}{}
460 returns the surface normals on the boundary of the \Domain as a \Data object.
461 \end{methoddesc}
462 %
463 \begin{methoddesc}[Domain]{getSize}{}
464 returns the local sample size, i.e. the element diameter, as a \Data object.
465 \end{methoddesc}
466 %
467 \begin{methoddesc}[Domain]{setTagMap}{tag_name, tag}
468 defines a mapping of the tag name \var{tag_name} to the \var{tag}.
469 \end{methoddesc}
470 %
471 \begin{methoddesc}[Domain]{getTag}{tag_name}
472 returns the tag associated with the tag name \var{tag_name}.
473 \end{methoddesc}
474 %
475 \begin{methoddesc}[Domain]{isValidTagName}{tag_name}
476 returns \True if \var{tag_name} is a valid tag name.
477 \end{methoddesc}
478 %
479 \begin{methoddesc}[Domain]{__eq__}{arg}
480 (\PYTHON \var{==} operator) returns \True if the \Domain \var{arg}
481 describes the same domain, \False otherwise.
482 \end{methoddesc}
483 %
484 \begin{methoddesc}[Domain]{__ne__}{arg}
485 (\PYTHON \var{!=} operator) returns \True if the \Domain \var{arg} does
486 not describe the same domain, \False otherwise.
487 \end{methoddesc}
488 %
489 \begin{methoddesc}[Domain]{__str__}{}
490 (\PYTHON \var{str()} function) returns a string representation of the
491 \Domain.
492 \end{methoddesc}
493 %
494 \begin{methoddesc}[Domain]{onMasterProcessor}{}
495 returns \True if the process is the master process within the \MPI
496 process group used by the \Domain. This is the process with rank 0.
497 If \MPI support is not enabled the return value is always \True.
498 \end{methoddesc}
499 %
500 \begin{methoddesc}[Domain]{getMPISize}{}
501 returns the number of \MPI processes used for this \Domain. If \MPI
502 support is not enabled 1 is returned.
503 \end{methoddesc}
504 %
505 \begin{methoddesc}[Domain]{getMPIRank}{}
506 returns the rank of the process executing the statement within the
507 \MPI process group used by the \Domain. If \MPI support is not enabled
508 0 is returned.
509 \end{methoddesc}
510 %
511 \begin{methoddesc}[Domain]{MPIBarrier}{}
512 executes barrier synchronization within the \MPI process group used by
513 the \Domain. If \MPI support is not enabled, this command does nothing.
514 \end{methoddesc}
515
516 \subsection{The \FunctionSpace class}
517 \begin{classdesc}{FunctionSpace}{}
518 \FunctionSpace objects, which are instantiated by generator functions, are
519 used to define properties of \Data objects such as continuity.
520 A \Data object in a particular \FunctionSpace is represented by its values at
521 \DataSamplePoints which are defined by the type and the \Domain of the \FunctionSpace.
522 \end{classdesc}
523
524 \vspace{1em}\noindent The following methods are available:
525 %
526 \begin{methoddesc}[FunctionSpace]{getDim}{}
527 returns the spatial dimensionality of the \Domain of the \FunctionSpace.
528 \end{methoddesc}
529 %
530 \begin{methoddesc}[FunctionSpace]{getX}{}
531 returns the location of the \DataSamplePoints.
532 \end{methoddesc}
533 %
534 \begin{methoddesc}[FunctionSpace]{getNormal}{}
535 If the domain of functions in the \FunctionSpace is a hyper-manifold (e.g.
536 the boundary of a domain) the method returns the outer normal at each of
537 the \DataSamplePoints. Otherwise an exception is raised.
538 \end{methoddesc}
539 %
540 \begin{methoddesc}[FunctionSpace]{getSize}{}
541 returns a \Data object measuring the spacing of the \DataSamplePoints.
542 The size may be zero.
543 \end{methoddesc}
544 %
545 \begin{methoddesc}[FunctionSpace]{getDomain}{}
546 returns the \Domain of the \FunctionSpace.
547 \end{methoddesc}
548 %
549 \begin{methoddesc}[FunctionSpace]{setTags}{new_tag, mask}
550 assigns a new tag \var{new_tag} to all data samples where \var{mask} is
551 positive for a least one data point.
552 \var{mask} must be defined on this \FunctionSpace.
553 Use the \var{setTagMap} to assign a tag name to \var{new_tag}.
554 \end{methoddesc}
555 %
556 \begin{methoddesc}[FunctionSpace]{__eq__}{arg}
557 (\PYTHON \var{==} operator) returns \True if the \FunctionSpace \var{arg}
558 describes the same function space, \False otherwise.
559 \end{methoddesc}
560 %
561 \begin{methoddesc}[FunctionSpace]{__ne__}{arg}
562 (\PYTHON \var{!=} operator) returns \True if the \FunctionSpace \var{arg}
563 does not describe the same function space, \False otherwise.
564 \end{methoddesc}
565
566 \begin{methoddesc}[Domain]{__str__}{}
567 (\PYTHON \var{str()} function) returns a string representation of the
568 \FunctionSpace.
569 \end{methoddesc}
570
571 \noindent The following functions provide generators for \FunctionSpace objects:
572
573 \begin{funcdesc}{Function}{domain}
574 returns the \Function on the \Domain \var{domain}. \Data objects in this
575 type of \Function are defined over the whole geometric region defined by
576 \var{domain}.
577 \end{funcdesc}
578 %
579 \begin{funcdesc}{ContinuousFunction}{domain}
580 returns the \ContinuousFunction on the \Domain domain. \Data objects in
581 this type of \Function are defined over the whole geometric region defined
582 by \var{domain} and assumed to represent a continuous function.
583 \end{funcdesc}
584 %
585 \begin{funcdesc}{FunctionOnBoundary}{domain}
586 returns the \FunctionOnBoundary on the \Domain domain. \Data objects in
587 this type of \Function are defined on the boundary of the geometric region
588 defined by \var{domain}.
589 \end{funcdesc}
590 %
591 \begin{funcdesc}{FunctionOnContactZero}{domain}
592 returns the \FunctionOnContactZero the \Domain domain. \Data objects in
593 this type of \Function are defined on side 0 of a discontinuity within
594 the geometric region defined by \var{domain}.
595 The discontinuity is defined when \var{domain} is instantiated.
596 \end{funcdesc}
597 %
598 \begin{funcdesc}{FunctionOnContactOne}{domain}
599 returns the \FunctionOnContactOne on the \Domain domain. \Data objects in
600 this type of \Function are defined on side 1 of a discontinuity within
601 the geometric region defined by \var{domain}.
602 The discontinuity is defined when \var{domain} is instantiated.
603 \end{funcdesc}
604 %
605 \begin{funcdesc}{Solution}{domain}
606 returns the \SolutionFS on the \Domain domain. \Data objects in this type
607 of \Function are defined on the geometric region defined by \var{domain}
608 and are solutions of partial differential equations\index{partial differential equation}.
609 \end{funcdesc}
610 %
611 \begin{funcdesc}{ReducedSolution}{domain}
612 returns the \ReducedSolutionFS on the \Domain domain. \Data objects in
613 this type of \Function are defined on the geometric region defined by
614 \var{domain} and are solutions of partial differential
615 equations\index{partial differential equation} with a reduced smoothness
616 for the solution approximation.
617 \end{funcdesc}
618
619 \subsection{The \Data Class}
620 \label{SEC ESCRIPT DATA}
621
622 The following table shows arithmetic operations that can be performed
623 point-wise on \Data objects:
624 \begin{center}
625 \begin{tabular}{l|l}
626 \textbf{Expression} & \textbf{Description}\\
627 \hline
628 \code{+arg} & identical to \var{arg}\index{+}\\
629 \code{-arg} & negation of \var{arg}\index{-}\\
630 \code{arg0+arg1} & adds \var{arg0} and \var{arg1}\index{+}\\
631 \code{arg0*arg1} & multiplies \var{arg0} and \var{arg1}\index{*}\\
632 \code{arg0-arg1} & subtracts \var{arg1} from \var{arg0}\index{-}\\
633 \code{arg0/arg1} & divides \var{arg0} by \var{arg1}\index{/}\\
634 \code{arg0**arg1} & raises \var{arg0} to the power of \var{arg1}\index{**}\\
635 \end{tabular}
636 \end{center}
637 At least one of the arguments \var{arg0} or \var{arg1} must be a \Data object.
638 Either of the arguments may be a \Data object, a \PYTHON number or a \numpy
639 object.
640 If \var{arg0} or \var{arg1} are not defined on the same \FunctionSpace, then
641 an attempt is made to convert \var{arg0} to the \FunctionSpace of \var{arg1}
642 or to convert \var{arg1} to \var{arg0}'s \FunctionSpace.
643 Both arguments must have the same \Shape or one of the arguments may be of
644 rank 0 (a constant).
645 The returned \Data object has the same \Shape and is defined on
646 the \DataSamplePoints as \var{arg0} or \var{arg1}.
647
648 The following table shows the update operations that can be applied to
649 \Data objects:
650 \begin{center}
651 \begin{tabular}{l|l}
652 \textbf{Expression} & \textbf{Description}\\
653 \hline
654 \code{arg0+=arg1} & adds \var{arg1} to \var{arg0}\index{+}\\
655 \code{arg0*=arg1} & multiplies \var{arg0} by \var{arg1}\index{*}\\
656 \code{arg0-=arg1} & subtracts \var{arg1} from\var{arg0}\index{-}\\
657 \code{arg0/=arg1} & divides \var{arg0} by \var{arg1}\index{/}\\
658 \code{arg0**=arg1} & raises \var{arg0} to the power of \var{arg1}\index{**}\\
659 \end{tabular}
660 \end{center}
661 \var{arg0} must be a \Data object. \var{arg1} must be a \Data object or an
662 object that can be converted into a \Data object.
663 \var{arg1} must have the same \Shape as \var{arg0} or have rank 0.
664 In the latter case it is assumed that the values of \var{arg1} are constant
665 for all components. \var{arg1} must be defined in the same \FunctionSpace as
666 \var{arg0} or it must be possible to interpolate \var{arg1} onto the
667 \FunctionSpace of \var{arg0}.
668
669 The \Data class supports taking slices as well as assigning new values to a
670 slice of an existing \Data object\index{slicing}.
671 The following expressions for taking and setting slices are valid:
672 \begin{center}
673 \begin{tabular}{l|ll}
674 \textbf{Rank of \var{arg}} & \textbf{Slicing expression} & \textbf{\Shape of returned and assigned object}\\
675 \hline
676 0 & no slicing & N/A\\
677 1 & \var{arg[l0:u0]} & (\var{u0}-\var{l0},)\\
678 2 & \var{arg[l0:u0,l1:u1]} & (\var{u0}-\var{l0},\var{u1}-\var{l1})\\
679 3 & \var{arg[l0:u0,l1:u1,l2:u2]} & (\var{u0}-\var{l0},\var{u1}-\var{l1},\var{u2}-\var{l2})\\
680 4 & \var{arg[l0:u0,l1:u1,l2:u2,l3:u3]} & (\var{u0}-\var{l0},\var{u1}-\var{l1},\var{u2}-\var{l2},\var{u3}-\var{l3})\\
681 \end{tabular}
682 \end{center}
683 Let \var{s} be the \Shape of \var{arg}, then
684 \begin{align*}
685 0 \le \var{l0} \le \var{u0} \le \var{s[0]},\\
686 0 \le \var{l1} \le \var{u1} \le \var{s[1]},\\
687 0 \le \var{l2} \le \var{u2} \le \var{s[2]},\\
688 0 \le \var{l3} \le \var{u3} \le \var{s[3]}.
689 \end{align*}
690 Any of the lower indexes \var{l0}, \var{l1}, \var{l2} and \var{l3} may not be
691 present in which case $0$ is assumed.
692 Any of the upper indexes \var{u0}, \var{u1}, \var{u2} and \var{u3} may be
693 omitted, in which case the upper limit for that dimension is assumed.
694 The lower and upper index may be identical in which case the column and the
695 lower or upper index may be dropped.
696 In the returned or in the object assigned to a slice, the corresponding
697 component is dropped, i.e. the rank is reduced by one in comparison to \var{arg}.
698 The following examples show slicing in action:
699 \begin{python}
700 t=Data(1., (4,4,6,6), Function(mydomain))
701 t[1,1,1,0]=9.
702 s=t[:2,:,2:6,5] # s has rank 3
703 s[:,:,1]=1.
704 t[:2,:2,5,5]=s[2:4,1,:2]
705 \end{python}
706
707
708 \subsection{Generation of \Data objects}
709 \begin{classdesc}{Data}{value=0, shape=(,), what=FunctionSpace(), expanded=\False}
710 creates a \Data object with \Shape \var{shape} in the \FunctionSpace \var{what}.
711 The values at all \DataSamplePoints are set to the double value \var{value}.
712 If \var{expanded} is \True the \Data object is represented in expanded form.
713 \end{classdesc}
714
715 \begin{classdesc}{Data}{value, what=FunctionSpace(), expanded=\False}
716 creates a \Data object in the \FunctionSpace \var{what}.
717 The value for each data sample point is set to \var{value}, which could be a
718 \numpy object, \Data object or a dictionary of \numpy or floating point
719 numbers. In the latter case the keys must be integers and are used as tags.
720 The \Shape of the returned object is equal to the \Shape of \var{value}.
721 If \var{expanded} is \True the \Data object is represented in expanded form.
722 \end{classdesc}
723
724 \begin{classdesc}{Data}{}
725 creates an \EmptyData object. The \EmptyData object is used to indicate that
726 an argument is not present where a \Data object is required.
727 \end{classdesc}
728
729 \begin{funcdesc}{Scalar}{value=0., what=FunctionSpace(), expanded=\False}
730 returns a \Data object of rank 0 (a constant) in the \FunctionSpace \var{what}.
731 Values are initialized with \var{value}, a double precision quantity.
732 If \var{expanded} is \True the \Data object is represented in expanded form.
733 \end{funcdesc}
734
735 \begin{funcdesc}{Vector}{value=0., what=FunctionSpace(), expanded=\False}
736 returns a \Data object of \Shape \var{(d,)} in the \FunctionSpace \var{what},
737 where \var{d} is the spatial dimension of the \Domain of \var{what}.
738 Values are initialized with \var{value}, a double precision quantity.
739 If \var{expanded} is \True the \Data object is represented in expanded form.
740 \end{funcdesc}
741
742 \begin{funcdesc}{Tensor}{value=0., what=FunctionSpace(), expanded=\False}
743 returns a \Data object of \Shape \var{(d,d)} in the \FunctionSpace \var{what},
744 where \var{d} is the spatial dimension of the \Domain of \var{what}.
745 Values are initialized with \var{value}, a double precision quantity.
746 If \var{expanded} is \True the \Data object is represented in expanded form.
747 \end{funcdesc}
748
749 \begin{funcdesc}{Tensor3}{value=0., what=FunctionSpace(), expanded=\False}
750 returns a \Data object of \Shape \var{(d,d,d)} in the \FunctionSpace \var{what},
751 where \var{d} is the spatial dimension of the \Domain of \var{what}.
752 Values are initialized with \var{value}, a double precision quantity.
753 If \var{expanded} is \True the \Data object is represented in expanded form.
754 \end{funcdesc}
755
756 \begin{funcdesc}{Tensor4}{value=0., what=FunctionSpace(), expanded=\False}
757 returns a \Data object of \Shape \var{(d,d,d,d)} in the \FunctionSpace \var{what},
758 where \var{d} is the spatial dimension of the \Domain of \var{what}.
759 Values are initialized with \var{value}, a double precision quantity.
760 If \var{expanded} is \True the \Data object is represented in expanded form.
761 \end{funcdesc}
762
763 \begin{funcdesc}{ComplexData}{value, what=FunctionSpace(), expanded=\False}
764 creates a \Data object in the \FunctionSpace \var{what}.
765 The value for each data sample point is set to the complex value \var{value}, which could be a
766 \numpy object, \Data object or a dictionary of \numpy or floating point
767 numbers. In the latter case the keys must be integers and are used as tags.
768 The \Shape of the returned object is equal to the \Shape of \var{value}.
769 If \var{expanded} is \True the \Data object is represented in expanded form.
770 \end{funcdesc}
771
772 \begin{funcdesc}{ComplexData}{value=0, shape=(,), what=FunctionSpace(), expanded=\False}
773 creates a \Data object with \Shape \var{shape} in the \FunctionSpace \var{what}.
774 The values at all \DataSamplePoints are set to the complex value \var{value}.
775 If \var{expanded} is \True the \Data object is represented in expanded form.
776 \end{funcdesc}
777
778 \begin{funcdesc}{ComplexData}{}
779 creates an \EmptyData object with complex values (i.e. with memory allocated to
780 store a complex number). The \EmptyData object is used to indicate that
781 an argument is not present where a \Data object is required.
782 \end{funcdesc}
783
784 \begin{funcdesc}{ComplexScalar}{value=0.+0.j, what=FunctionSpace(), expanded=\False}
785 returns a \Data object of rank 0 (a constant) in the \FunctionSpace \var{what}.
786 Values are initialized with complex \var{value}, a double precision quantity.
787 If \var{expanded} is \True the \Data object is represented in expanded form.
788 \end{funcdesc}
789
790 \begin{funcdesc}{ComplexData}{value=0.+0.j, what=FunctionSpace(), expanded=\False}
791 returns a \Data object of \Shape \var{(d,)} in the \FunctionSpace \var{what},
792 where \var{d} is the spatial dimension of the \Domain of \var{what}.
793 Values are initialized with complex \var{value}, a double precision quantity.
794 If \var{expanded} is \True the \Data object is represented in expanded form.
795 \end{funcdesc}
796
797 \begin{funcdesc}{ComplexTensor}{value=0.+0.j, what=FunctionSpace(), expanded=\False}
798 returns a \Data object of \Shape \var{(d,d)} in the \FunctionSpace \var{what},
799 where \var{d} is the spatial dimension of the \Domain of \var{what}.
800 Values are initialized with complex \var{value}, a double precision quantity.
801 If \var{expanded} is \True the \Data object is represented in expanded form.
802 \end{funcdesc}
803
804 \begin{funcdesc}{ComplexTensor3}{value=0.+0.j, what=FunctionSpace(), expanded=\False}
805 returns a \Data object of \Shape \var{(d,d,d)} in the \FunctionSpace \var{what},
806 where \var{d} is the spatial dimension of the \Domain of \var{what}.
807 Values are initialized with complex \var{value}, a double precision quantity.
808 If \var{expanded} is \True the \Data object is represented in expanded form.
809 \end{funcdesc}
810
811 \begin{funcdesc}{ComplexTensor4}{value=0.+0.j, what=FunctionSpace(), expanded=\False}
812 returns a \Data object of \Shape \var{(d,d,d,d)} in the \FunctionSpace \var{what},
813 where \var{d} is the spatial dimension of the \Domain of \var{what}.
814 Values are initialized with complex \var{value}, a double precision quantity.
815 If \var{expanded} is \True the \Data object is represented in expanded form.
816 \end{funcdesc}
817
818 \begin{funcdesc}{load}{filename, domain}
819 recovers a \Data object on \Domain \var{domain} from the file \var{filename},
820 which was created by \function{dump}.
821 \end{funcdesc}
822
823 \subsection{Generating random \Data objects}
824 A \Data object filled with random values can be produced using the
825 \function{RandomData} function.
826 By default values are drawn uniformly at random from the interval $[0,1]$ (i.e.
827 including end points).
828 The function takes a shape for the data points and a \FunctionSpace for the new
829 \Data as arguments.
830 For example:
831 \begin{python}
832 from esys.finley import *
833 from esys.escript import *
834
835 domain=Rectangle(11,11)
836 fs=ContinuousFunction(domain)
837 d=RandomData((), fs)
838 \end{python}
839 would result in \var{d} being filled with scalar random data since \texttt{()}
840 is an empty tuple.
841
842 \begin{python}
843 from esys.finley import *
844 from esys.escript import *
845
846 domain=Rectangle(11,11)
847 fs=ContinuousFunction(domain)
848 d=RandomData((2,2), fs)
849 \end{python}
850 would give \var{d} the same number of data points, but each point would be a
851 $2\times 2$ matrix instead of a scalar.
852
853 By default, the seed used to generate the random values will be different each
854 time.
855 If required, you can specify a seed to ensure the same sequence is produced.
856 \begin{python}
857 from esys.dudley import *
858 from esys.escript import *
859
860 seed=-17171717
861 domain=Brick(10,10,10)
862 fs=Function(domain)
863 d=RandomData((2,2), fs, seed)
864 \end{python}
865
866 The \var{seed} can be any integer value\footnote{which can be converted to a
867 C++ long} but 0 is special.
868 A seed of zero will cause \escript to use a different seed each time.
869 Also, note that the mechanism used to produce the random values could be
870 different in different releases.
871
872 \noindent\textbf{Note for MPI users:}
873 \textsl{
874 Even if you specify a seed, you will only get the same results if you are running with the same
875 number of ranks.
876 If you change the number of ranks, you will get different values for the same seed.
877 }
878
879 \subsubsection{Smoothed randoms}
880 The \ripley domains (see Chapter \ref{chap:ripley}) support generating random
881 scalars which are smoothed using Gaussian blur.
882 To use this, you need to supply the radius of the filter kernel (in elements)
883 and the \var{sigma} value used in the filter.
884 For example:
885 \begin{python}
886 from esys.ripley import *
887 from esys.escript import *
888
889 fs=ContinuousFunction(Rectangle(11,11, d1=2,d0=2))
890 d=RandomData((), fs, 0, ('gaussian', 1, 0.5))
891 \end{python}
892 will use a filter that uses the immediate neighbours of each point with a sigma
893 value of $0.5$.
894 The random values will be different each time this code is executed due to the
895 seed of $0$.
896
897 Ripley's Gaussian smoothing has the following requirements:
898 \begin{enumerate}
899 \item If \MPI is in use, then each rank must have at least $5$ elements in
900 it \emph{in each dimension}. This value increases as the radius of
901 the blur increases.
902 \item The data being generated must be scalar. (You can generate random
903 data objects for \ripley domains with whatever shape you require, you
904 just can't smooth them unless that shape is scalar).
905 \end{enumerate}
906 An exception will be raised if either of these requirements is not met.
907
908 The components of the matrix used in the kernal for the 2D case are
909 defined\cite{gaussfilter} by:
910
911 \[ G(x,y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} \]
912
913 \noindent For the 3D case, we use:
914
915 \[ G(x,y) = \frac{1}{(\sqrt{2\pi\sigma^2})^3} e^{-\frac{x^2+y^2+z^2}{2\sigma^2}} \]
916
917 All distances ($x$,$y$,$z$) refer to the number of points from the centre point.
918 That is, the closest neighbours have at least one distance of $1$, the next
919 ``ring'' of neighbours have at least one $2$ and so on.
920 The matrix is normalised before use.
921
922 \subsection{\Data methods}
923 These are the most frequently used methods of the \Data class.
924 A complete list of methods can be found in the reference guide,
925 see \ReferenceGuide.
926
927 \begin{methoddesc}[Data]{getFunctionSpace}{}
928 returns the \FunctionSpace of the object.
929 \end{methoddesc}
930
931 \begin{methoddesc}[Data]{getDomain}{}
932 returns the \Domain of the object.
933 \end{methoddesc}
934
935 \begin{methoddesc}[Data]{getShape}{}
936 returns the \Shape of the object as a \class{tuple} of integers.
937 \end{methoddesc}
938
939 \begin{methoddesc}[Data]{getRank}{}
940 returns the rank of the data on each data point\index{rank}.
941 \end{methoddesc}
942
943 \begin{methoddesc}[Data]{isEmpty}{}
944 returns \True if the \Data object is the \EmptyData object, \False otherwise.
945 Note that this is not the same as asking if the object contains no \DataSamplePoints.
946 \end{methoddesc}
947
948 \begin{methoddesc}[Data]{setTaggedValue}{tag_name, value}
949 assigns the \var{value} to all \DataSamplePoints which have the tag
950 assigned to \var{tag_name}. \var{value} must be an object of class
951 \class{numpy.ndarray} or must be convertible into a \class{numpy.ndarray} object.
952 \var{value} (or the corresponding \class{numpy.ndarray} object) must be of
953 rank $0$ or must have the same rank as the object.
954 If a value has already been defined for tag \var{tag_name} within the object
955 it is overwritten by the new \var{value}. If the object is expanded,
956 the value assigned to \DataSamplePoints with tag \var{tag_name} is replaced by
957 \var{value}. If no value is assigned the tag name \var{tag_name}, no value is set.
958 \end{methoddesc}
959
960 \begin{methoddesc}[Data]{dump}{filename}
961 dumps the \Data object to the file \var{filename}. The file stores the
962 function space but not the \Domain. It is the responsibility of the user to
963 save the \Domain in order to be able to recover the \Data object.
964 \end{methoddesc}
965
966 \begin{methoddesc}[Data]{__str__}{}
967 returns a string representation of the object.
968 \end{methoddesc}
969
970 \subsection{Functions of \Data objects}
971 This section lists the most important functions for \Data class objects.
972 A complete list and a more detailed description of the functionality can be
973 found on \ReferenceGuide.
974
975 \begin{funcdesc}{kronecker}{d}
976 returns a \RankTwo in \FunctionSpace \var{d} such that
977 \begin{equation}
978 \code{kronecker(d)}\left[ i,j\right] = \left\{
979 \begin{array}{l l}
980 1 & \quad \text{if $i=j$}\\
981 0 & \quad \text{otherwise}
982 \end{array}
983 \right.
984 \end{equation}
985 If \var{d} is an integer a $(d,d)$ \numpy array is returned.
986 \end{funcdesc}
987
988 \begin{funcdesc}{identityTensor}{d}
989 is a synonym for \code{kronecker} (see above).
990 \end{funcdesc}
991
992 \begin{funcdesc}{identityTensor4}{d}
993 returns a \RankFour in \FunctionSpace \var{d} such that
994 \begin{equation}
995 \code{identityTensor(d)}\left[ i,j,k,l\right] = \left\{
996 \begin{array}{l l}
997 1 & \quad \text{if $i=k$ and $j=l$}\\
998 0 & \quad \text{otherwise}
999 \end{array}
1000 \right.
1001 \end{equation}
1002 If \var{d} is an integer a $(d,d,d,d)$ \numpy array is returned.
1003 \end{funcdesc}
1004
1005 \begin{funcdesc}{unitVector}{i,d}
1006 returns a \RankOne in \FunctionSpace \var{d} such that
1007 \begin{equation}
1008 \code{identityTensor(d)}\left[ j \right] = \left\{
1009 \begin{array}{l l}
1010 1 & \quad \text{if $j=i$}\\
1011 0 & \quad \text{otherwise}
1012 \end{array}
1013 \right.
1014 \end{equation}
1015 If \var{d} is an integer a $(d,)$ \numpy array is returned.
1016 \end{funcdesc}
1017
1018 \begin{funcdesc}{Lsup}{a}
1019 returns the $L^{sup}$ norm of \var{arg}. This is the maximum of the absolute
1020 values over all components and all \DataSamplePoints of \var{a}.
1021 \end{funcdesc}
1022
1023 \begin{funcdesc}{sup}{a}
1024 returns the maximum value over all components and all \DataSamplePoints of \var{a}.
1025 \end{funcdesc}
1026
1027 \begin{funcdesc}{inf}{a}
1028 returns the minimum value over all components and all \DataSamplePoints of \var{a}
1029 \end{funcdesc}
1030
1031 \begin{funcdesc}{minval}{a}
1032 returns at each data sample point the minimum value over all components.
1033 \end{funcdesc}
1034
1035 \begin{funcdesc}{maxval}{a}
1036 returns at each data sample point the maximum value over all components.
1037 \end{funcdesc}
1038
1039 \begin{funcdesc}{length}{a}
1040 returns the Euclidean norm at each data sample point.
1041 For a \RankFour \var{a} this is
1042 \begin{equation}
1043 \code{length(a)}=\sqrt{\sum_{ijkl} \var{a} \left[i,j,k,l\right]^2}
1044 \end{equation}
1045 \end{funcdesc}
1046
1047 \begin{funcdesc}{trace}{a\optional{, axis_offset=0}}
1048 returns the trace of \var{a}. This is the sum over components \var{axis_offset}
1049 and \var{axis_offset+1} with the same index.
1050 For instance, in the case of a \RankTwo this is
1051 \begin{equation}
1052 \code{trace(a)}=\sum_{i} \var{a} \left[i,i\right]
1053 \end{equation}
1054 and for a \RankFour and \code{axis_offset=1} this is
1055 \begin{equation}
1056 \code{trace(a,1)}\left[i,j\right]=\sum_{k} \var{a} \left[i,k,k,j\right]
1057 \end{equation}
1058 \end{funcdesc}
1059
1060 \begin{funcdesc}{transpose}{a\optional{, axis_offset=None}}
1061 returns the transpose of \var{a}. This swaps the first \var{axis_offset}
1062 components of \var{a} with the rest. If \var{axis_offset} is not
1063 present \code{int(r/2)} is used where \var{r} is the rank of \var{a}.
1064 For instance, in the case of a \RankTwo this is
1065 \begin{equation}
1066 \code{transpose(a)}\left[i,j\right]=\var{a} \left[j,i\right]
1067 \end{equation}
1068 and for a \RankFour and \code{axis_offset=1} this is
1069 \begin{equation}
1070 \code{transpose(a,1)}\left[i,j,k,l\right]=\var{a} \left[j,k,l,i\right]
1071 \end{equation}
1072 \end{funcdesc}
1073
1074 \begin{funcdesc}{swap_axes}{a\optional{, axis0=0 \optional{, axis1=1 }}}
1075 returns \var{a} but with swapped components \var{axis0} and \var{axis1}.
1076 The argument \var{a} must be at least of rank 2. For instance, if \var{a}
1077 is a \RankFour, \code{axis0=1} and \code{axis1=2}, the result is
1078 \begin{equation}
1079 \code{swap_axes(a,1,2)}\left[i,j,k,l\right]=\var{a} \left[i,k,j,l\right]
1080 \end{equation}
1081 \end{funcdesc}
1082
1083 \begin{funcdesc}{symmetric}{a}
1084 returns the symmetric part of \var{a}. This is \code{(a+transpose(a))/2}.
1085 \end{funcdesc}
1086
1087 \begin{funcdesc}{nonsymmetric}{a}
1088 returns the non-symmetric part of \var{a}. This is \code{(a-transpose(a))/2}.
1089 \end{funcdesc}
1090
1091 \begin{funcdesc}{inverse}{a}
1092 return the inverse of \var{a} so that
1093 \begin{equation}
1094 \code{matrix_mult(inverse(a),a)=kronecker(d)}
1095 \end{equation}
1096 if \var{a} has shape \code{(d,d)}. The current implementation is restricted to
1097 arguments of shape \code{(2,2)} and \code{(3,3)}.
1098 \end{funcdesc}
1099
1100 \begin{funcdesc}{eigenvalues}{a}
1101 returns the eigenvalues of \var{a} so that
1102 \begin{equation}
1103 \code{matrix_mult(a,V)=e[i]*V}
1104 \end{equation}
1105 where \code{e=eigenvalues(a)} and \var{V} is a suitable non-zero vector.
1106 The eigenvalues are ordered in increasing size.
1107 The argument \var{a} has to be symmetric, i.e. \code{a=symmetric(a)}.
1108 The current implementation is restricted to arguments of shape \code{(2,2)}
1109 and \code{(3,3)}.
1110 \end{funcdesc}
1111
1112 \begin{funcdesc}{eigenvalues_and_eigenvectors}{a}
1113 returns the eigenvalues and eigenvectors of \var{a}.
1114 \begin{equation}
1115 \code{matrix_mult(a,V[:,i])=e[i]*V[:,i]}
1116 \end{equation}
1117 where \code{e,V=eigenvalues_and_eigenvectors(a)}. The eigenvectors \var{V} are
1118 orthogonal and normalized, i.e.
1119 \begin{equation}
1120 \code{matrix_mult(transpose(V),V)=kronecker(d)}
1121 \end{equation}
1122 if \var{a} has shape \code{(d,d)}. The eigenvalues are ordered in increasing
1123 size. The argument \var{a} has to be the symmetric, i.e. \code{a=symmetric(a)}.
1124 The current implementation is restricted to arguments of shape \code{(2,2)}
1125 and \code{(3,3)}.
1126 \end{funcdesc}
1127
1128 \begin{funcdesc}{maximum}{*a}
1129 returns the maximum value over all arguments at all \DataSamplePoints and for each component.
1130 \begin{equation}
1131 \code{maximum(a0,a1)}\left[i,j\right]=max(\var{a0} \left[i,j\right],\var{a1} \left[i,j\right])
1132 \end{equation}
1133 at all \DataSamplePoints.
1134 \end{funcdesc}
1135
1136 \begin{funcdesc}{minimum}{*a}
1137 returns the minimum value over all arguments at all \DataSamplePoints and for each component.
1138 \begin{equation}
1139 \code{minimum(a0,a1)}\left[i,j\right]=min(\var{a0} \left[i,j\right],\var{a1} \left[i,j\right])
1140 \end{equation}
1141 at all \DataSamplePoints.
1142 \end{funcdesc}
1143
1144 \begin{funcdesc}{clip}{a\optional{, minval=0.}\optional{, maxval=1.}}
1145 cuts back \var{a} into the range between \var{minval} and \var{maxval}.
1146 A value in the returned object equals \var{minval} if the corresponding value
1147 of \var{a} is less than \var{minval}, equals \var{maxval} if the corresponding
1148 value of \var{a} is greater than \var{maxval}, or corresponding value of
1149 \var{a} otherwise.
1150 \end{funcdesc}
1151
1152 \begin{funcdesc}{inner}{a0, a1}
1153 returns the inner product of \var{a0} and \var{a1}. For instance in the
1154 case of a \RankTwo:
1155 \begin{equation}
1156 \code{inner(a)}=\sum_{ij}\var{a0} \left[j,i\right] \cdot \var{a1} \left[j,i\right]
1157 \end{equation}
1158 and for a \RankFour:
1159 \begin{equation}
1160 \code{inner(a)}=\sum_{ijkl}\var{a0} \left[i,j,k,l\right] \cdot \var{a1} \left[j,i,k,l\right]
1161 \end{equation}
1162 \end{funcdesc}
1163
1164 \begin{funcdesc}{matrix_mult}{a0, a1}
1165 returns the matrix product of \var{a0} and \var{a1}.
1166 If \var{a1} is a \RankOne this is
1167 \begin{equation}
1168 \code{matrix_mult(a)}\left[i\right]=\sum_{k}\var{a0} \cdot \left[i,k\right]\var{a1} \left[k\right]
1169 \end{equation}
1170 and if \var{a1} is a \RankTwo this is
1171 \begin{equation}
1172 \code{matrix_mult(a)}\left[i,j\right]=\sum_{k}\var{a0} \cdot \left[i,k\right]\var{a1} \left[k,j\right]
1173 \end{equation}
1174 \end{funcdesc}
1175
1176 \begin{funcdesc}{transposed_matrix_mult}{a0, a1}
1177 returns the matrix product of the transposed of \var{a0} and \var{a1}.
1178 The function is equivalent to \code{matrix_mult(transpose(a0),a1)}.
1179 If \var{a1} is a \RankOne this is
1180 \begin{equation}
1181 \code{transposed_matrix_mult(a)}\left[i\right]=\sum_{k}\var{a0} \cdot \left[k,i\right]\var{a1} \left[k\right]
1182 \end{equation}
1183 and if \var{a1} is a \RankTwo this is
1184 \begin{equation}
1185 \code{transposed_matrix_mult(a)}\left[i,j\right]=\sum_{k}\var{a0} \cdot \left[k,i\right]\var{a1} \left[k,j\right]
1186 \end{equation}
1187 \end{funcdesc}
1188
1189 \begin{funcdesc}{matrix_transposed_mult}{a0, a1}
1190 returns the matrix product of \var{a0} and the transposed of \var{a1}.
1191 The function is equivalent to \code{matrix_mult(a0,transpose(a1))}.
1192 If \var{a1} is a \RankTwo this is
1193 \begin{equation}
1194 \code{matrix_transposed_mult(a)}\left[i,j\right]=\sum_{k}\var{a0} \cdot \left[i,k\right]\var{a1} \left[j,k\right]
1195 \end{equation}
1196 \end{funcdesc}
1197
1198 \begin{funcdesc}{outer}{a0, a1}
1199 returns the outer product of \var{a0} and \var{a1}.
1200 For instance, if both, \var{a0} and \var{a1} is a \RankOne then
1201 \begin{equation}
1202 \code{outer(a)}\left[i,j\right]=\var{a0} \left[i\right] \cdot \var{a1}\left[j\right]
1203 \end{equation}
1204 and if \var{a0} is a \RankOne and \var{a1} is a \RankThree:
1205 \begin{equation}
1206 \code{outer(a)}\left[i,j,k\right]=\var{a0} \left[i\right] \cdot \var{a1}\left[j,k\right]
1207 \end{equation}
1208 \end{funcdesc}
1209
1210 \begin{funcdesc}{tensor_mult}{a0, a1}
1211 returns the tensor product of \var{a0} and \var{a1}.
1212 If \var{a1} is a \RankTwo this is
1213 \begin{equation}
1214 \code{tensor_mult(a)}\left[i,j\right]=\sum_{kl}\var{a0}\left[i,j,k,l\right] \cdot \var{a1} \left[k,l\right]
1215 \end{equation}
1216 and if \var{a1} is a \RankFour this is
1217 \begin{equation}
1218 \code{tensor_mult(a)}\left[i,j,k,l\right]=\sum_{mn}\var{a0} \left[i,j,m,n\right] \cdot \var{a1} \left[m,n,k,l\right]
1219 \end{equation}
1220 \end{funcdesc}
1221
1222 \begin{funcdesc}{transposed_tensor_mult}{a0, a1}
1223 returns the tensor product of the transposed of \var{a0} and \var{a1}.
1224 The function is equivalent to \code{tensor_mult(transpose(a0),a1)}.
1225 If \var{a1} is a \RankTwo this is
1226 \begin{equation}
1227 \code{transposed_tensor_mult(a)}\left[i,j\right]=\sum_{kl}\var{a0}\left[k,l,i,j\right] \cdot \var{a1} \left[k,l\right]
1228 \end{equation}
1229 and if \var{a1} is a \RankFour this is
1230 \begin{equation}
1231 \code{transposed_tensor_mult(a)}\left[i,j,k,l\right]=\sum_{mn}\var{a0} \left[m,n,i,j\right] \cdot \var{a1} \left[m,n,k,l\right]
1232 \end{equation}
1233 \end{funcdesc}
1234
1235 \begin{funcdesc}{tensor_transposed_mult}{a0, a1}
1236 returns the tensor product of \var{a0} and the transposed of \var{a1}.
1237 The function is equivalent to \code{tensor_mult(a0,transpose(a1))}.
1238 If \var{a1} is a \RankTwo this is
1239 \begin{equation}
1240 \code{tensor_transposed_mult(a)}\left[i,j\right]=\sum_{kl}\var{a0}\left[i,j,k,l\right] \cdot \var{a1} \left[l,k\right]
1241 \end{equation}
1242 and if \var{a1} is a \RankFour this is
1243 \begin{equation}
1244 \code{tensor_transposed_mult(a)}\left[i,j,k,l\right]=\sum_{mn}\var{a0} \left[i,j,m,n\right] \cdot \var{a1} \left[k,l,m,n\right]
1245 \end{equation}
1246 \end{funcdesc}
1247
1248 \begin{funcdesc}{grad}{a\optional{, where=None}}
1249 returns the gradient of \var{a}. If \var{where} is present the gradient will
1250 be calculated in the \FunctionSpace \var{where}, otherwise a default
1251 \FunctionSpace is used. In case that \var{a} is a \RankTwo one has
1252 \begin{equation}
1253 \code{grad(a)}\left[i,j,k\right]=\frac{\partial \var{a} \left[i,j\right]}{\partial x_{k}}
1254 \end{equation}
1255 \end{funcdesc}
1256
1257 \begin{funcdesc}{integrate}{a\optional{, where=None}}
1258 returns the integral of \var{a} where the domain of integration is defined by
1259 the \FunctionSpace of \var{a}. If \var{where} is present the argument is
1260 interpolated into \FunctionSpace \var{where} before integration.
1261 For instance in the case of a \RankTwo in \ContinuousFunction it is
1262 \begin{equation}
1263 \code{integrate(a)}\left[i,j\right]=\int_{\Omega}\var{a} \left[i,j\right] \; d\Omega
1264 \end{equation}
1265 where $\Omega$ is the spatial domain and $d\Omega$ volume integration.
1266 To integrate over the boundary of the domain one uses
1267 \begin{equation}
1268 \code{integrate(a,where=FunctionOnBoundary(a.getDomain))}\left[i,j\right]=\int_{\partial \Omega} a\left[i,j\right] \; ds
1269 \end{equation}
1270 where $\partial \Omega$ is the surface of the spatial domain and $ds$ area or
1271 line integration.
1272 \end{funcdesc}
1273
1274 \begin{funcdesc}{interpolate}{a, where}
1275 interpolates argument \var{a} into the \FunctionSpace \var{where}.
1276 \end{funcdesc}
1277
1278 \begin{funcdesc}{div}{a\optional{, where=None}}
1279 returns the divergence of \var{a}:
1280 \begin{equation}
1281 \code{div(a)=trace(grad(a),where)}
1282 \end{equation}
1283 \end{funcdesc}
1284
1285 \begin{funcdesc}{jump}{a\optional{, domain=None}}
1286 returns the jump of \var{a} over the discontinuity in its domain or if
1287 \Domain \var{domain} is present in \var{domain}.
1288 \begin{equation}
1289 \begin{array}{rcl}
1290 \code{jump(a)}& = &\code{interpolate(a,FunctionOnContactOne(domain))} \\
1291 & & \hfill - \code{interpolate(a,FunctionOnContactZero(domain))}
1292 \end{array}
1293 \end{equation}
1294 \end{funcdesc}
1295
1296 \begin{funcdesc}{L2}{a}
1297 returns the $L^2$-norm of \var{a} in its \FunctionSpace. This is
1298 \begin{equation}
1299 \code{L2(a)=integrate(length(a)}^2\code{)} \; .
1300 \end{equation}
1301 \end{funcdesc}
1302
1303 \noindent The following functions operate ``point-wise''.
1304 That is, the operation is applied to each component of each point individually.
1305
1306 \begin{funcdesc}{sin}{a}
1307 applies the sine function to \var{a}.
1308 \end{funcdesc}
1309
1310 \begin{funcdesc}{cos}{a}
1311 applies the cosine function to \var{a}.
1312 \end{funcdesc}
1313
1314 \begin{funcdesc}{tan}{a}
1315 applies the tangent function to \var{a}.
1316 \end{funcdesc}
1317
1318 \begin{funcdesc}{asin}{a}
1319 applies the arc (inverse) sine function to \var{a}.
1320 \end{funcdesc}
1321
1322 \begin{funcdesc}{acos}{a}
1323 applies the arc (inverse) cosine function to \var{a}.
1324 \end{funcdesc}
1325
1326 \begin{funcdesc}{atan}{a}
1327 applies the arc (inverse) tangent function to \var{a}.
1328 \end{funcdesc}
1329
1330 \begin{funcdesc}{sinh}{a}
1331 applies the hyperbolic sine function to \var{a}.
1332 \end{funcdesc}
1333
1334 \begin{funcdesc}{cosh}{a}
1335 applies the hyperbolic cosine function to \var{a}.
1336 \end{funcdesc}
1337
1338 \begin{funcdesc}{tanh}{a}
1339 applies the hyperbolic tangent function to \var{a}.
1340 \end{funcdesc}
1341
1342 \begin{funcdesc}{asinh}{a}
1343 applies the arc (inverse) hyperbolic sine function to \var{a}.
1344 \end{funcdesc}
1345
1346 \begin{funcdesc}{acosh}{a}
1347 applies the arc (inverse) hyperbolic cosine function to \var{a}.
1348 \end{funcdesc}
1349
1350 \begin{funcdesc}{atanh}{a}
1351 applies the arc (inverse) hyperbolic tangent function to \var{a}.
1352 \end{funcdesc}
1353
1354 \begin{funcdesc}{exp}{a}
1355 applies the exponential function to \var{a}.
1356 \end{funcdesc}
1357
1358 \begin{funcdesc}{sqrt}{a}
1359 applies the square root function to \var{a}.
1360 \end{funcdesc}
1361
1362 \begin{funcdesc}{log}{a}
1363 takes the natural logarithm of \var{a}.
1364 \end{funcdesc}
1365
1366 \begin{funcdesc}{log10}{a}
1367 takes the base-$10$ logarithm of \var{a}.
1368 \end{funcdesc}
1369
1370 \begin{funcdesc}{sign}{a}
1371 applies the sign function to \var{a}. The result is $1$ where \var{a} is
1372 positive, $-1$ where \var{a} is negative, and $0$ otherwise.
1373 \end{funcdesc}
1374
1375 \begin{funcdesc}{wherePositive}{a}
1376 returns a function which is $1$ where \var{a} is positive and $0$ otherwise.
1377 \end{funcdesc}
1378
1379 \begin{funcdesc}{whereNegative}{a}
1380 returns a function which is $1$ where \var{a} is negative and $0$ otherwise.
1381 \end{funcdesc}
1382
1383 \begin{funcdesc}{whereNonNegative}{a}
1384 returns a function which is $1$ where \var{a} is non-negative and $0$ otherwise.
1385 \end{funcdesc}
1386
1387 \begin{funcdesc}{whereNonPositive}{a}
1388 returns a function which is $1$ where \var{a} is non-positive and $0$ otherwise.
1389 \end{funcdesc}
1390
1391 \begin{funcdesc}{whereZero}{a\optional{, tol=None\optional{, rtol=1.e-8}}}
1392 returns a function which is $1$ where \var{a} equals zero with tolerance
1393 \var{tol} and $0$ otherwise. If \var{tol} is not present, the absolute maximum
1394 value of \var{a} times \var{rtol} is used.
1395 \end{funcdesc}
1396
1397 \begin{funcdesc}{whereNonZero}{a\optional{, tol=None\optional{, rtol=1.e-8}}}
1398 returns a function which is $1$ where \var{a} is non-zero with tolerance
1399 \var{tol} and $0$ otherwise. If \var{tol} is not present, the absolute maximum
1400 value of \var{a} times \var{rtol} is used.
1401 \end{funcdesc}
1402
1403 \subsection{Interpolating Data}
1404 \index{interpolateTable}
1405 \label{sec:interpolation}
1406 In some cases, it may be useful to produce Data objects which fit some user
1407 defined function.
1408 Manually modifying each value in the Data object is not a good idea since it
1409 depends on knowing the location and order of each data point in the domain.
1410 Instead, \escript can use an interpolation table to produce a \Data object.
1411
1412 The following example is available as \file{int_save.py} in the \ExampleDirectory.
1413 We will produce a \Data object which approximates a sine curve.
1414
1415 \begin{python}
1416 from esys.escript import saveDataCSV, sup, interpolateTable
1417 import numpy
1418 from esys.finley import Rectangle
1419
1420 n=4
1421 r=Rectangle(n,n)
1422 x=r.getX()
1423 toobig=100
1424 \end{python}
1425
1426 \noindent First we produce an interpolation table:
1427 \begin{python}
1428 sine_table=[0, 0.70710678118654746, 1, 0.70710678118654746, 0,
1429 -0.70710678118654746, -1, -0.70710678118654746, 0]
1430 \end{python}
1431 %
1432 We wish to identify $0$ and $1$ with the ends of the curve, that is
1433 with the first and eighth value in the table.
1434
1435 \begin{python}
1436 numslices=len(sine_table)-1
1437 minval=0.
1438 maxval=1.
1439 step=sup(maxval-minval)/numslices
1440 \end{python}
1441 %
1442 So the values $v$ from the input lie in the interval
1443 \var{minval} $\leq v <$ \var{maxval}.
1444 \var{step} represents the gap (in the input range) between entries in the table.
1445 By default, values of $v$ outside the table argument range (minval, maxval)
1446 will be pushed back into the range, i.e. if $v <$ \var{minval} the value
1447 \var{minval} will be used to evaluate the table.
1448 Similarly, for values $v>$ \var{maxval} the value \var{maxval} is used.
1449
1450 Now we produce our new \Data object:
1451
1452 \begin{python}
1453 result=interpolateTable(sine_table, x[0], minval, step, toobig)
1454 \end{python}
1455 Any values which interpolate to larger than \var{toobig} will raise an
1456 exception. You can switch on boundary checking by adding
1457 \code{check_boundaries=True} to the argument list.
1458
1459 Now consider a 2D example. We will interpolate from a plane where $\forall x,y\in[0,9]:(x,y)=x+y\cdot10$.
1460
1461 \begin{python}
1462 from esys.escript import whereZero
1463 table2=[]
1464 for y in range(0,10):
1465 r=[]
1466 for x in range(0,10):
1467 r.append(x+y*10)
1468 table2.append(r)
1469 xstep=(maxval-minval)/(10-1)
1470 ystep=(maxval-minval)/(10-1)
1471
1472 xmin=minval
1473 ymin=minval
1474
1475 result2=interpolateTable(table2, x2, (xmin, ymin), (xstep, ystep), toobig)
1476 \end{python}
1477
1478 We can check the values using \function{whereZero}.
1479 For example, for $x=0$:
1480 \begin{python}
1481 print(result2*whereZero(x[0]))
1482 \end{python}
1483
1484 Finally let us look at a 3D example. Note that the parameter tuples should be
1485 $(x,y,z)$ but that in the interpolation table, $x$ is the innermost dimension.
1486 \begin{python}
1487 b=Brick(n,n,n)
1488 x3=b.getX()
1489 toobig=1000000
1490
1491 table3=[]
1492 for z in range(0,10):
1493 face=[]
1494 for y in range(0,10):
1495 r=[]
1496 for x in range(0,10):
1497 r.append(x+y*10+z*100)
1498 face.append(r)
1499 table3.append(face);
1500
1501 zstep=(maxval-minval)/(10-1)
1502
1503 zmin=minval
1504
1505 result3=interpolateTable(table3, x3, (xmin, ymin, zmin),
1506 (xstep, ystep, zstep), toobig)
1507 \end{python}
1508
1509
1510 \subsubsection{Non-uniform Interpolation}
1511 Non-uniform interpolation is also supported for the one dimensional case.
1512 \begin{python}
1513 Data.nonuniformInterpolate(in, out, check_boundaries)
1514 Data.nonuniformSlope(in, out, check_boundaries)
1515 \end{python}
1516
1517 Will produce a new \Data object by mapping the given \Data object through the user-defined function
1518 specified by \texttt{in} and \texttt{out}.
1519 The \ldots Interpolate version gives the value of the function at the specified point and the
1520 \ldots Slope version gives the slope at those points.
1521 The check_boundaries boolean argument specifies what the function should do if the \Data object contains
1522 values outside the range specified by the \texttt{in} parameter.
1523 If the argument is \texttt{False}, then those datapoints will be interpolated to the value of the edge
1524 they are closest to (or assigned a slope of zero).
1525 If the argument is \texttt{True}, then an exception will be thrown if out of bounds values are detected.
1526 Note that the values given by the \texttt{in} parameter must be monotonically increasing.
1527
1528 \noindent For example:\\
1529 If \texttt{d} contains the values \texttt{\{1,2,3,4,5\}}, then
1530 \begin{python}
1531 d.nonuniformInterpolate([1.5, 2, 2.8, 4.6], [4, 5, -1, 1], False)
1532 \end{python}
1533 would produce a \Data object containing \texttt{\{4, 5, -0.7777, 0.3333, 1\}}.\\
1534 A similar call to \texttt{nonuniformSlope} would produce a \Data object containing \texttt{\{0, 2, 1.1111, 1.1111, 0\}}.
1535 %
1536 %
1537 % We will interpolate a surface such that the bottom
1538 % edge is the sine curve described above.
1539 % The amplitude of the curve decreases as we move towards the top edge.
1540 % Our interpolation table will have three rows:
1541 %
1542 % \begin{python}
1543 % st=numpy.array(sine_table)
1544 % table=[st, 0.5*st, 0*st]
1545 % \end{python}
1546 % %
1547 % The use of \numpy and multiplication here is just to save typing.
1548 %
1549 % % result2=x1.interpolateTable(table, 0, 0.55, x0, minval, step, toobig)
1550 % \begin{python}
1551 % result=interpolateTable(table, x (minval,0), (0.55, step), toobig)
1552 % \end{python}
1553 %
1554 % In the 2D case the start and step parameters are tuples $(x,y)$.
1555 % By default, if a point is specified which is outside the boundary, then
1556 % \var{interpolateTable} will operate as if the point was on the boundary.
1557 % Passing \code{check_boundaries=True} will lead to the rejection of any points
1558 % outside the boundaries by \var{interpolateTable}.
1559 %
1560 % This method can also be called with three dimensional tables and \Data objects.
1561 % Tuples should be ordered $(x,y,z)$.
1562
1563 \subsection{The \var{DataManager} Class}
1564 \label{sec:datamanager}
1565
1566 The \var{DataManager} class can be used to conveniently add checkpoint/restart
1567 functionality to \escript simulations.
1568 Once an instance is created \Data objects and other values can be added and
1569 dumped to disk by a single method call.
1570 If required the object can be set up to also save the data in a format suitable
1571 for visualization.
1572 Internally the \var{DataManager} interfaces with \weipa for this.
1573
1574 \begin{classdesc}{DataManager}{formats=[RESTART], work_dir=".", restart_prefix="restart", do_restart=\True}
1575 initializes a new \var{DataManager} object which can be used to save,
1576 restore and export simulation data in a number of formats.
1577 All files and directories saved or restored by this object are located
1578 under the directory specified by \var{work_dir}.
1579 If \var{RESTART} is specified in \var{formats}, the \var{DataManager} will
1580 look for directories whose name starts with \var{restart_prefix}.
1581 In case \var{do_restart} is \True, the last of these directories is used
1582 to restore simulation data while all others are deleted.
1583 If \var{do_restart} is \False, then all of those directories are deleted.
1584 The \var{restart_prefix} and \var{do_restart} parameters are ignored if
1585 \var{RESTART} is not specified in \var{formats}.
1586 \end{classdesc}
1587
1588 \noindent Valid values for the \var{formats} parameter are:
1589 \begin{memberdesc}[DataManager]{RESTART}
1590 enables writing of checkpoint files to be able to continue simulations
1591 as explained in the class description.
1592 \end{memberdesc}
1593 \begin{memberdesc}[DataManager]{SILO}
1594 exports simulation data in the \SILO file format. \escript must have
1595 been compiled with \SILO support for this to work.
1596 \end{memberdesc}
1597 \begin{memberdesc}[DataManager]{VISIT}
1598 enables the \VisIt simulation interface which allows connecting to and
1599 interacting with the running simulation from a compatible \VisIt client.
1600 \escript must have been compiled with \VisIt (version 2) support and the
1601 version of the client has to match the version used at compile time.
1602 In order to connect to the simulation the client needs to have access and
1603 load the file \file{escriptsim.sim2} located under the work directory.
1604 \end{memberdesc}
1605 \begin{memberdesc}[DataManager]{VTK}
1606 exports simulation data in the \VTK file format.
1607 \end{memberdesc}
1608
1609 \noindent The \var{DataManager} class has the following methods:
1610 \begin{methoddesc}[DataManager]{addData}{**data}
1611 adds \Data objects and other data to the manager. Calling this method does
1612 not save or export the data yet so it is allowed to incrementally add data
1613 at various points in the simulation script if required.
1614 Note, that only a single domain is supported so all \Data objects have to
1615 be defined on the same one or an exception is raised.
1616 \end{methoddesc}
1617
1618 \begin{methoddesc}[DataManager]{setDomain}{domain}
1619 explicitly sets the domain for this manager.
1620 It is generally not required to call this method directly.
1621 Instead, the \var{addData} method will set the domain used by the \Data
1622 objects.
1623 An exception is raised if the domain was set to a different domain before
1624 (explicitly or implicitly).
1625 \end{methoddesc}
1626
1627 \begin{methoddesc}[DataManager]{hasData}{}
1628 returns \True if the manager has loaded simulation data for a restart.
1629 \end{methoddesc}
1630
1631 \begin{methoddesc}[DataManager]{getDomain}{}
1632 returns the domain as recovered from a restart.
1633 \end{methoddesc}
1634
1635 \begin{methoddesc}[DataManager]{getValue}{value_name}
1636 returns a \Data object or other value with the name \var{value_name} that
1637 has been recovered after a restart.
1638 \end{methoddesc}
1639
1640 \begin{methoddesc}[DataManager]{getCycle}{}
1641 returns the export cycle, i.e. the number of times \var{export()} has been
1642 called.
1643 \end{methoddesc}
1644
1645 \begin{methoddesc}[DataManager]{setCheckpointFrequency}{freq}
1646 sets the frequency with which checkpoint files are created. This is only
1647 useful if the \var{DataManager} object was created with at least one other
1648 format next to \var{RESTART}. The frequency is 1 by default which means
1649 that checkpoint files are created every time \var{export()} is called.
1650 Unlike visualization output, a simulation checkpoint is usually not
1651 required at every time step. Thus, the frequency can be decreased by
1652 calling this method with $\var{freq}>1$ which would then create restart
1653 files every \var{freq} times \var{export()} is called.
1654 \end{methoddesc}
1655
1656 \begin{methoddesc}[DataManager]{setTime}{time}
1657 sets the simulation time stamp. This floating point number is stored in
1658 the metadata of exported data but not used by \var{RESTART}.
1659 \end{methoddesc}
1660
1661 \begin{methoddesc}[DataManager]{setMeshLabels}{x, y, z=""}
1662 sets labels for the mesh axes. These are currently only used by the \SILO
1663 exporter.
1664 \end{methoddesc}
1665
1666 \begin{methoddesc}[DataManager]{setMeshUnits}{x, y, z=""}
1667 sets units for the mesh axes. These are currently only used by the \SILO
1668 exporter.
1669 \end{methoddesc}
1670
1671 \begin{methoddesc}[DataManager]{setMetadataSchemaString}{schema, metadata=""}
1672 sets metadata namespaces and the corresponding metadata. These are
1673 currently only used by the \VTK exporter.
1674 \var{schema} is a dictionary that maps prefixes to namespace names, e.g.\\
1675 \code{\{"gml": "http://www.opengis.net/gml"\}} and \var{metadata} is a
1676 string with the actual content which will be enclosed in \var{<MetaData>}
1677 tags.
1678 \end{methoddesc}
1679
1680 \begin{methoddesc}[DataManager]{export}{}
1681 executes the actual data export. Depending on the \var{formats} parameter
1682 used in the constructor all data added by \var{addData()} is written to
1683 disk (\var{RESTART,SILO,VTK}) or made available through the \VisIt
1684 simulation interface (\var{VISIT}).
1685 At least the domain must be set for something to be exported.
1686 \end{methoddesc}
1687
1688 \subsection{Saving Data as CSV}
1689 \label{sec:savedatacsv}
1690 \index{saveDataCSV}\index{CSV}
1691 For simple post-processing, \Data objects can be saved in comma separated
1692 value (\emph{CSV}) format.
1693 If \var{mydata1} and \var{mydata2} are scalar data, the command
1694 \begin{python}
1695 saveDataCSV('output.csv', U=mydata1, V=mydata2)
1696 \end{python}
1697 will record the values in \file{output.csv} in the following format:
1698 \begin{verbatim}
1699 U, V
1700 1.0000000e+0, 2.0000000e-1
1701 5.0000000e-0, 1.0000000e+1
1702 ...
1703 \end{verbatim}
1704
1705 The names of the keyword parameters form the names of columns in the output.
1706 If the data objects are over different function spaces, then \var{saveDataCSV}
1707 will attempt to interpolate to a common function space.
1708 If this is not possible, then an exception is raised.
1709
1710 Output can be restricted using a scalar mask as follows:
1711 \begin{python}
1712 saveDataCSV('outfile.csv', U=mydata1, V=mydata2, mask=myscalar)
1713 \end{python}
1714 This command will only output those rows which correspond to to positive
1715 values of \var{myscalar}.
1716 Some aspects of the output can be tuned using additional parameters:
1717 \begin{python}
1718 saveDataCSV('data.csv', refid=True, append=True, sep=' ', csep='/', mask=mymask, e=mat1)
1719 \end{python}
1720
1721 \begin{itemize}
1722 \item \var{refid} -- specifies that the output should include the reference IDs of the elements or nodes
1723 \item \var{append} -- specifies that the output should be written to the end of an existing file
1724 \item \var{sep} -- defines the separator between fields
1725 \item \var{csep} -- defines the separator between components in the header
1726 line. For example between the components of a matrix.
1727 \end{itemize}
1728 %
1729 The above command would produce output like this:
1730 \begin{verbatim}
1731 refid e/0/0 e/1/0 e/0/1 e/1/1
1732 0 1.0000000000e+00 2.0000000000e+00 3.0000000000e+00 4.0000000000e+00
1733 ...
1734 \end{verbatim}
1735
1736 Note that while the order in which rows are output can vary, all the elements
1737 in a given row always correspond to the same input.
1738
1739 \subsection{Converting \Data to a Numpy Array}
1740 \label{sec:getnumpy}
1741 \index{getNumpy}\index{GN}
1742 \Data objects can be converted into a numpy structured array using the commands \var{getNumpy} and \var{convertNumpy}.
1743 \subsubsection{getNumpy}
1744 If \var{mydata1} and \var{mydata2} are scalar \Data, then the command
1745 \begin{python}
1746 a,b = getNumpy(U=mydata1, V=mydata2)
1747 \end{python}
1748 will return two structured ndarrays with the names '\emph{U}' and '\emph{V}'.
1749 \begin{verbatim}
1750 a['U'] = [1.0000000e+0, 2.0000000e-1, ...
1751 b['V'] = [2.0000000e+0, 3.0000000e-1, ...
1752 \end{verbatim}
1753
1754 Up to five \Data objects can be passed to \var{getNumpy} at the time. These objects can be scalar, vector or tensor \Data objects. The names of the keyword parameters form the names of the returned arrays.
1755 If the data objects are over different function spaces, then \var{getNumpy}
1756 will attempt to interpolate to a common function space.
1757 If this is not possible, then an exception is raised.
1758
1759 Output can be restricted using a scalar mask as follows:
1760 \begin{python}
1761 a,b,c = getNumpy(U=mydata1, V=mydata2, W=mydata3, mask=myscalar)
1762 \end{python}
1763 This command will only output those rows which correspond to to positive
1764 values of \var{myscalar}.
1765
1766 Note that while the order in which output rows are output can vary, all the elements
1767 in a given row always correspond to the same input.
1768
1769 \subsubsection{convertNumpy}
1770 \Data objects can also be converted into a numpy structured array using the command \var{convertNumpy}.
1771 If \var{mydata1} is a \Data object, then the command
1772 \begin{python}
1773 a = convertNumpy(mydata1)
1774 \end{python}
1775 will return a structured ndarray containing all of the data in \var{mydata1}. Unlike \var{getNumpy}, this function
1776 does not support the use of masks and does not use MPI.
1777
1778 \subsection{The \Operator Class}
1779 The \Operator class provides an abstract access to operators built
1780 within the \LinearPDE class. \Operator objects are created
1781 when a PDE is handed over to a PDE solver library and handled
1782 by the \LinearPDE object defining the PDE. The user can gain access
1783 to the \Operator of a \LinearPDE object through the \var{getOperator}
1784 method.
1785
1786 \begin{classdesc}{Operator}{}
1787 creates an empty \Operator object.
1788 \end{classdesc}
1789
1790 \begin{methoddesc}[Operator]{isEmpty}{fileName}
1791 returns \True is the object is empty, \False otherwise.
1792 \end{methoddesc}
1793
1794 \begin{methoddesc}[Operator]{resetValues}{}
1795 resets all entries in the operator.
1796 \end{methoddesc}
1797
1798 \begin{methoddesc}[Operator]{solve}{rhs}
1799 returns the solution \var{u} of: operator * \var{u} = \var{rhs}.
1800 \end{methoddesc}
1801
1802 \begin{methoddesc}[Operator]{of}{u}
1803 applies the operator to the \Data object \var{u}, i.e. performs a matrix-vector
1804 multiplication.
1805 \end{methoddesc}
1806
1807 \begin{methoddesc}[Operator]{saveMM}{fileName}\index{Matrix Market}
1808 saves the object to a Matrix Market format file with name \var{fileName}, see
1809 \url{http://math.nist.gov/MatrixMarket}
1810 \end{methoddesc}
1811
1812 \section{Physical Units}
1813 \escript provides support for physical units in the SI system\index{SI units}
1814 including unit conversion. So the user can define variables in the form
1815 \begin{python}
1816 from esys.escript.unitsSI import *
1817 l=20*m
1818 w=30*kg
1819 w2=40*lb
1820 T=100*Celsius
1821 \end{python}
1822 In the two latter cases a conversion from pounds\index{pounds} and degrees
1823 Celsius\index{Celsius} is performed into the appropriate SI units \emph{kg}
1824 and \emph{Kelvin}.
1825 In addition, composed units can be used, for instance
1826 \begin{python}
1827 from esys.escript.unitsSI import *
1828 rho=40*lb/cm**3
1829 \end{python}
1830 defines the density in the units of pounds per cubic centimeter.
1831 The value $40$ will be converted into SI units, in this case kg per cubic
1832 meter. Moreover unit prefixes are supported:
1833 \begin{python}
1834 from esys.escript.unitsSI import *
1835 p=40*Mega*Pa
1836 \end{python}
1837 The pressure \var{p} is set to 40 Mega Pascal. Units can also be converted
1838 back from the SI system into a desired unit, e.g.
1839 \begin{python}
1840 from esys.escript.unitsSI import *
1841 print(p/atm)
1842 \end{python}
1843 can be used print the pressure in units of atmosphere\index{atmosphere}.
1844
1845 The following is an incomplete list of supported physical units:
1846
1847 \begin{datadesc}{km}
1848 unit of kilometer
1849 \end{datadesc}
1850
1851 \begin{datadesc}{m}
1852 unit of meter
1853 \end{datadesc}
1854
1855 \begin{datadesc}{cm}
1856 unit of centimeter
1857 \end{datadesc}
1858
1859 \begin{datadesc}{mm}
1860 unit of millimeter
1861 \end{datadesc}
1862
1863 \begin{datadesc}{sec}
1864 unit of second
1865 \end{datadesc}
1866
1867 \begin{datadesc}{minute}
1868 unit of minute
1869 \end{datadesc}
1870
1871 \begin{datadesc}{h}
1872 unit of hour
1873 \end{datadesc}
1874
1875 \begin{datadesc}{day}
1876 unit of day
1877 \end{datadesc}
1878
1879 \begin{datadesc}{yr}
1880 unit of year
1881 \end{datadesc}
1882
1883 \begin{datadesc}{gram}
1884 unit of gram
1885 \end{datadesc}
1886
1887 \begin{datadesc}{kg}
1888 unit of kilogram
1889 \end{datadesc}
1890
1891 \begin{datadesc}{lb}
1892 unit of pound
1893 \end{datadesc}
1894
1895 \begin{datadesc}{ton}
1896 metric ton
1897 \end{datadesc}
1898
1899 \begin{datadesc}{A}
1900 unit of Ampere
1901 \end{datadesc}
1902
1903 \begin{datadesc}{Hz}
1904 unit of Hertz
1905 \end{datadesc}
1906
1907 \begin{datadesc}{N}
1908 unit of Newton
1909 \end{datadesc}
1910
1911 \begin{datadesc}{Pa}
1912 unit of Pascal
1913 \end{datadesc}
1914
1915 \begin{datadesc}{atm}
1916 unit of atmosphere
1917 \end{datadesc}
1918
1919 \begin{datadesc}{J}
1920 unit of Joule
1921 \end{datadesc}
1922
1923 \begin{datadesc}{W}
1924 unit of Watt
1925 \end{datadesc}
1926
1927 \begin{datadesc}{C}
1928 unit of Coulomb
1929 \end{datadesc}
1930
1931 \begin{datadesc}{V}
1932 unit of Volt
1933 \end{datadesc}
1934
1935 \begin{datadesc}{F}
1936 unit of Farad
1937 \end{datadesc}
1938
1939 \begin{datadesc}{Ohm}
1940 unit of Ohm
1941 \end{datadesc}
1942
1943 \begin{datadesc}{K}
1944 unit of degrees Kelvin
1945 \end{datadesc}
1946
1947 \begin{datadesc}{Celsius}
1948 unit of degrees Celsius
1949 \end{datadesc}
1950
1951 \begin{datadesc}{Fahrenheit}
1952 unit of degrees Fahrenheit
1953 \end{datadesc}
1954
1955 \noindent Supported unit prefixes:
1956
1957 \begin{datadesc}{Yotta}
1958 prefix yotta = $10^{24}$
1959 \end{datadesc}
1960
1961 \begin{datadesc}{Zetta}
1962 prefix zetta = $10^{21}$
1963 \end{datadesc}
1964
1965 \begin{datadesc}{Exa}
1966 prefix exa = $10^{18}$
1967 \end{datadesc}
1968
1969 \begin{datadesc}{Peta}
1970 prefix peta = $10^{15}$
1971 \end{datadesc}
1972
1973 \begin{datadesc}{Tera}
1974 prefix tera = $10^{12}$
1975 \end{datadesc}
1976
1977 \begin{datadesc}{Giga}
1978 prefix giga = $10^9$
1979 \end{datadesc}
1980
1981 \begin{datadesc}{Mega}
1982 prefix mega = $10^6$
1983 \end{datadesc}
1984
1985 \begin{datadesc}{Kilo}
1986 prefix kilo = $10^3$
1987 \end{datadesc}
1988
1989 \begin{datadesc}{Hecto}
1990 prefix hecto = $10^2$
1991 \end{datadesc}
1992
1993 \begin{datadesc}{Deca}
1994 prefix deca = $10^1$
1995 \end{datadesc}
1996
1997 \begin{datadesc}{Deci}
1998 prefix deci = $10^{-1}$
1999 \end{datadesc}
2000
2001 \begin{datadesc}{Centi}
2002 prefix centi = $10^{-2}$
2003 \end{datadesc}
2004
2005 \begin{datadesc}{Milli}
2006 prefix milli = $10^{-3}$
2007 \end{datadesc}
2008
2009 \begin{datadesc}{Micro}
2010 prefix micro = $10^{-6}$
2011 \end{datadesc}
2012
2013 \begin{datadesc}{Nano}
2014 prefix nano = $10^{-9}$
2015 \end{datadesc}
2016
2017 \begin{datadesc}{Pico}
2018 prefix pico = $10^{-12}$
2019 \end{datadesc}
2020
2021 \begin{datadesc}{Femto}
2022 prefix femto = $10^{-15}$
2023 \end{datadesc}
2024
2025 \begin{datadesc}{Atto}
2026 prefix atto = $10^{-18}$
2027 \end{datadesc}
2028
2029 \begin{datadesc}{Zepto}
2030 prefix zepto = $10^{-21}$
2031 \end{datadesc}
2032
2033 \begin{datadesc}{Yocto}
2034 prefix yocto = $10^{-24}$
2035 \end{datadesc}
2036
2037 \section{Utilities}
2038 The \class{FileWriter} class provides a mechanism to write data to a file.
2039 In essence, this class wraps the standard \PYTHON \class{file} class to write
2040 data that are global in \MPI to a file. In fact, data are written on the
2041 processor with \MPI rank 0 only. It is recommended to use \class{FileWriter}
2042 rather than \class{open} in order to write code that will run with and without
2043 \MPI. It is safe to use \class{open} under \MPI to \emph{read} data which are
2044 global under \MPI.
2045
2046 \begin{classdesc}{FileWriter}{fn\optional{,append=\False, \optional{createLocalFiles=\False}})}
2047 Opens a file with name \var{fn} for writing. If \var{append} is set to \True
2048 data are appended at the end of the file.
2049 If running under \MPI, only the first processor (rank==0) will open the file
2050 and write to it.
2051 If \var{createLocalFiles} is set each individual processor will create a file
2052 where for any processor with rank $> 0$ the file name is extended by its rank.
2053 This option is normally used for debugging purposes only.
2054 \end{classdesc}
2055
2056 \vspace{1em}\noindent The following methods are available:
2057 \begin{methoddesc}[FileWriter]{close}{}
2058 closes the file.
2059 \end{methoddesc}
2060 \begin{methoddesc}[FileWriter]{flush}{}
2061 flushes the internal buffer to disk.
2062 \end{methoddesc}
2063 \begin{methoddesc}[FileWriter]{write}{txt}
2064 writes string \var{txt} to the file. Note that a newline is not added.
2065 \end{methoddesc}
2066 \begin{methoddesc}[FileWriter]{writelines}{txts}
2067 writes the list \var{txts} of strings to the file.
2068 Note that newlines are not added.
2069 This method is equivalent to calling \var{write()} for each string.
2070 \end{methoddesc}
2071 \begin{memberdesc}[FileWriter]{closed}
2072 this member is \True if the file is closed.
2073 \end{memberdesc}
2074 \begin{memberdesc}[FileWriter]{mode}
2075 holds the access mode.
2076 \end{memberdesc}
2077 \begin{memberdesc}[FileWriter]{name}
2078 holds the file name.
2079 \end{memberdesc}
2080 \begin{memberdesc}[FileWriter]{newlines}
2081 holds the line separator.
2082 \end{memberdesc}
2083
2084 \noindent The following additional functions are available in the \escript
2085 module:
2086 \begin{funcdesc}{setEscriptParamInt}{name,value}
2087 assigns the integer value \var{value} to the internal Escript parameter
2088 \var{name}. This should be considered an advanced feature and it is generally
2089 not required to call this function. One parameter worth mentioning is
2090 \var{name}="TOO_MANY_LINES" which affects the conversion of \Data objects to a
2091 string. If more than \var{value} lines would be created, a condensed format is
2092 used instead which reports the minimum and maximum values and general
2093 information about the \Data object rather than all values.
2094 \end{funcdesc}
2095
2096 \begin{funcdesc}{getEscriptParamInt}{name}
2097 returns the current value of internal Escript parameter \var{name}.
2098 \end{funcdesc}
2099
2100 \begin{funcdesc}{listEscriptParams}{a}
2101 returns a list of valid Escript parameters and their description.
2102 \end{funcdesc}
2103
2104 \begin{funcdesc}{getMPISizeWorld}{}
2105 returns the number of \MPI processes in use in the \env{MPI_COMM_WORLD}
2106 process group. If \MPI is not used 1 is returned.
2107 \end{funcdesc}
2108
2109 \begin{funcdesc}{getMPIRankWorld}{}
2110 returns the rank of the current process within the \env{MPI_COMM_WORLD}
2111 process group. If \MPI is not used 0 is returned.
2112 \end{funcdesc}
2113
2114 \begin{funcdesc}{MPIBarrierWorld}{}
2115 performs a barrier synchronization across all processes within the
2116 \env{MPI_COMM_WORLD} process group.
2117 \end{funcdesc}
2118
2119 \begin{funcdesc}{getMPIWorldMax}{a}
2120 returns the maximum value of the integer \var{a} across all processes within
2121 \env{MPI_COMM_WORLD}.
2122 \end{funcdesc}
2123
2124 \section{Lazy Evaluation of Data}
2125 \label{sec:lazy}
2126 Constant and Tagged representations of Data are relatively small but Expanded\footnote{Separate values stored for each point of the FunctionSpace.} are larger and
2127 will not entirely fit in CPU cache.
2128
2129 Escript's lazy evaluation features record operations performed on Data objects but do not actually carry them out until the Data is ``resolved''.
2130
2131 Consider the following code:
2132 \begin{python}
2133 from esys.escript import *
2134 from esys.dudley import Rectangle
2135 x=Rectangle(3,3)
2136 x=Rectangle(3,3).getX()
2137 c=Data((1.5, 1), x.getFunctionSpace())
2138 t=Data(((1,1),(0,1)), x.getFunctionSpace())
2139 t.tag()
2140 \end{python}
2141
2142 The variables \var{c}, \var{t}, \var{x} are stored as \texttt{constant}, \texttt{tagged} and \texttt{expanded} Data respectively.
2143 Printing those variables will show the values stored (or if we were to use a larger Rectangle, a summary).
2144
2145 \begin{python}
2146 v = matrix_mult(t,x) + c
2147 print(v.isExpanded())
2148 print(v)
2149 \end{python}
2150
2151 Will output \texttt{True} followed by all of the values for \var{v}.
2152 Now we'll introduce lazy evaluation:
2153
2154 \begin{python}
2155 xx = x.delay()
2156 print(xx.isExpanded(), xx.isLazy())
2157 print(x.isExpanded(), x.isLazy())
2158 print(xx)
2159 \end{python}
2160
2161 The first print will show that \var{xx} is not considered to be ``expanded'', while the second print shows that \var{x} is unaffected.
2162 The last print will produce something like:
2163 \begin{python}
2164 Lazy Data: [depth=0] E@0x55ed512ad760
2165 \end{python}
2166 The \texttt{E} before the \verb|@| shows that this lazy Data is wrapping ``expanded'' Data.
2167 Calling \texttt{.delay()} on constant or tagged Data results in \verb|C@...| and \verb|T@...| respectively.
2168
2169 If an input to an operation is lazy, then the result will be lazy as well\footnote{Matrix inverse is an exception to this.}:
2170 \begin{python}
2171 res = matrix_mult(t,-xx) + c
2172 print(res)
2173 \end{python}
2174 Will produce:
2175 \begin{python}
2176 Lazy Data: [depth=3] (prod(T@0x..., neg(E@...)) + C@0x...)
2177 \end{python}
2178 Depth indicates the largest number of operators from the top of the expression to the bottom.
2179
2180 To actually find the value of this lazy Data object, we need to resolve it:
2181 \begin{python}
2182 res.resolve()
2183 \end{python}
2184 Note that \texttt{resolve()} doesn't return a new object, but transforms the object it is called on.
2185 Printing, \var{res} now will show the values at each point.
2186
2187 \subsection{Lazyness and non-expanded Data}
2188 While it is possible to call delay on constant or tagged Data, escript will not build expressions consisting solely of such Data.
2189 \begin{python}
2190 cx=c.delay()
2191 res=cx+cx
2192 print(res)
2193 \end{python}
2194 would output:
2195 \begin{python}
2196 Lazy Data: [depth=0] C@0x55ed512cc7c0
2197 # Not
2198 Lazy Data: [depth=1] (C@0x... + C@0x...)
2199 \end{python}
2200
2201
2202 \subsection{When to resolve}
2203
2204 You are never \emph{required} to manually resolve lazy Data in \texttt{escript}.
2205 Any operations which need the actual values of an expression will either
2206 \begin{itemize}
2207 \item compute the values without resolving the whole Data object at once (solvers assembling FEM matrices)
2208 \item resolve the data automatically (everthing else)
2209 \end{itemize}
2210
2211 \noindent Escript will automatically resolve lazy Data:
2212 \begin{enumerate}
2213 \item If a matrix inversion operation is applied to the Data.
2214 \item If the expression tree becomes too deep\footnote{At time of writing, this threshold is somewhat arbitrarily set at \texttt{depth>9}, but this is configurable.}.
2215 \end{enumerate}
2216 Note, the second point is important when writing loops like this:
2217 \begin{python}
2218 # x is initial guess
2219 while err > tol:
2220 construct PDE coefficients involving x
2221 solve PDE
2222 calculate err
2223 update x
2224 \end{python}
2225
2226 After a few iterations of the loop, \var{x} may be something like \texttt{x=F(F(F(F(originalX))))}.
2227 So it will probably be better to \texttt{resolve} \var{x} at the end of each loop iteration.
2228 Alternatively, if \var{x} is included in many expressions in the loop, it may be better to resolve it earlier.
2229
2230 \subsection{Options for using lazy evaluation}
2231
2232 There are two ways to enable lazy evaluation:
2233 \begin{enumerate}
2234 \item Any escript script can make use of lazy evaluation by \texttt{delay()}-ing one of its expanded Data variables.
2235 Any expressions including that delayed variable (directly or indirectly) will be lazy until resolved.
2236 \item Setting the \texttt{AUTOLAZY} parameter for \texttt{escript} to \texttt{1}.
2237 In this case, most escript operation which would normally produce extended Data, will produce lazy Data instead.
2238 In general, this option is not recommended for two reasons:
2239 \begin{itemize}
2240 \item AUTOLAZY uses the \texttt{setEscriptParamInt()} which is not guaranteed to have continued support.
2241 \item Making everything lazy instead of just more complex objects is not likely to give significant efficiency improvements.
2242 \end{itemize}
2243 \end{enumerate}
2244
2245 \subsection{When to use lazy evaluation?}
2246 Exactly when using lazy evaluation will be more efficient is still an open question.
2247 When the objects being manipulated are large (eg 4-Tensors in Drucker-Prager), significant memory and runtime improvements can be achieved.
2248 See~\cite{lazyauspdc}.
2249
2250 Our best advice is to experiment with it.

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26