The Module esys.escript

Figure 3.1: Dependency of Function Spaces in Finley. An arrow indicates that a function in the function space at the starting point can be interpolated to the function space of the arrow target. All functionspaces on the left side can be interpolated to any of the functionspaces on the right.
\includegraphics[width=\textwidth]{figures/EscriptDiagram1}

esys.escript is a Python module that allows you to represent the values of a function at points in a Domain in such a way that the function will be useful for the Finite Element Method (FEM) simulation. It also provides what we call a function space that describes how the data is used in the simulation. Stored along with the data is information about the elements and nodes which will be used by esys.finley.

In order to understand what we mean by the term 'function space', consider that the solution of a partial differential equation (PDE) is a function on a domain $ \Omega$. When solving a PDE using FEM, the solution is piecewise-differentiable but, in general, its gradient is discontinuous. To reflect these different degrees of smoothness, different function spaces are used. For instance, in FEM, the displacement field is represented by its values at the nodes of the mesh, and so is continuous. The strain, which is the symmetric part of the gradient of the displacement field, is stored on the element centers, and so is considered to be discontinuous.

A function space is described by a FunctionSpace object. The following statement generates the object solution_space which is a FunctionSpace object and provides access to the function space of PDE solutions on the Domain mydomain:


\begin{python}
solution_space=Solution(mydomain)
\end{python}
The following generators for function spaces on a Domain mydomain are available:

The reduced smoothness for PDE solution is often used to fulfill the Ladyzhenskaya–-Babuska–-Brezzi condition [11] when solving saddle point problems , eg. the Stokes equation. A discontinuity is a region within the domain across which functions may be discontinuous. The location of discontinuity is defined in the Domain object. Figure 3.1 shows the dependency between the types of function spaces in Finley (other libraries may have different relationships).

The solution of a PDE is a continuous function. Any continuous function can be seen as a general function on the domain and can be restricted to the boundary as well as to one side of discontinuity (the result will be different depending on which side is chosen). Functions on any side of the discontinuity can be seen as a function on the corresponding other side.

A function on the boundary or on one side of the discontinuity cannot be seen as a general function on the domain as there are no values defined for the interior. For most PDE solver libraries the space of the solution and continuous functions is identical, however in some cases, eg. when periodic boundary conditions are used in esys.finley, a solution fulfills periodic boundary conditions while a continuous function does not have to be periodic.

The concept of function spaces describes the properties of functions and allows abstraction from the actual representation of the function in the context of a particular application. For instance, in the FEM context a function of the general FunctionSpace type (written as Function() in Figure 3.1) is usually represented by its values at the element center, but in a finite difference scheme the edge midpoint of cells is preferred. By changing its function space you can use the same function in a Finite Difference scheme instead of Finite Element scheme. Changing the function space of a particular function will typically lead to a change of its representation. So, when seen as a general function, a continuous function which is typically represented by its values on the node of the FEM mesh or finite difference grid must be interpolated to the element centers or the cell edges, respectively. Interpolation happens automatically in esys.escript whenever it is required.

In esys.escript the class that stores these functions is called Data. The function is represented through its values on data sample points where the data sample points are chosen according to the function space of the function. Data class objects are used to define the coefficients of the PDEs to be solved by a PDE solver library and also to store the solutions of the PDE.

The values of the function have a rank which gives the number of indices, and a shape defining the range of each index. The rank in esys.escript is limited to the range 0 through $ 4$ and it is assumed that the rank and shape is the same for all data sample points. The shape of a Data object is a tuple (list) s of integers. The length of s is the rank of the Data object and the i-th index ranges between 0 and $ \var{s[i]}-1$. For instance, a stress field has rank $ 2$ and shape $ (d,d)$ where $ d$ is the spatial dimension. The following statement creates the Data object mydat representing a continuous function with values of shape $ (2,3)$ and rank $ 2$:
\begin{python}
mydat=Data(value=1,what=ContinuousFunction(myDomain),shape=(2,3))
\end{python}
The initial value is the constant $ 1$ for all data sample points and all components.

Data objects can also be created from any numpy array or any object, such as a list of floating point numbers, that can be converted into a numpy.ndarray [6]. The following two statements create objects which are equivalent to mydat:
\begin{python}
mydat1=Data(value=numpy.ones((2,3)),what=ContinuousFunction(myDo...
...2=Data(value=[[1,1],[1,1],[1,1]],what=ContinuousFunction(myDomain))
\end{python}
In the first case the initial value is numpy.ones((2,3)) which generates a $ 2 \times 3$ matrix as a numpy.ndarray filled with ones. The shape of the created Data object it taken from the shape of the array. In the second case, the creator converts the initial value, which is a list of lists, and converts it into a numpy.ndarray before creating the actual Data object.

For convenience esys.escript provides creators for the most common types of Data objects in the following forms (d defines the spatial dimension):

Here the initial value is 0 but any object that can be converted into a numpy.ndarray and whose shape is consistent with shape of the Data object to be created can be used as the initial value.

Data objects can be manipulated by applying unary operations (eg. cos, sin, log) point and can be combined point-wise by applying arithmetic operations (eg. +, - ,* , /). It is to be emphasized that esys.escript itself does not handle any spatial dependencies as it does not know how values are interpreted by the processing PDE solver library. However esys.escript invokes interpolation if this is needed during data manipulations. Typically, this occurs in binary operation when both arguments belong to different function spaces or when data are handed over to a PDE solver library which requires functions to be represented in a particular way.

The following example shows the usage of Data objects: Assume we have a displacement field $ u$ and we want to calculate the corresponding stress field $ \sigma$ using the linear-elastic isotropic material model

$\displaystyle \sigma\hackscore {ij}=\lambda u\hackscore {k,k} \delta\hackscore {ij} + \mu ( u\hackscore {i,j} + u\hackscore {j,i})$     (64)

where $ \delta\hackscore{ij}$ is the Kronecker symbol and $ \lambda$ and $ \mu$ are the Lame coefficients. The following function takes the displacement u and the Lame coefficients lam and mu as arguments and returns the corresponding stress:
\begin{python}
from esys.escript import *
def getStress(u,lam,mu):
d=u.getDom...
...tress=lam*trace(g)*kronecker(d)+mu*(g+transpose(g))
return stress
\end{python}
The variable d gives the spatial dimension of the domain on which the displacements are defined. kronecker returns the Kronecker symbol with indexes $ i$ and $ j$ running from 0 to d-1. The call grad(u) requires the displacement field u to be in the Solution or continuous FunctionSpace function space. The result g as well as the returned stress will be in the general FunctionSpace function space. If, for example, u is the solution of a PDE then getStress might be called in the following way:
\begin{python}
s=getStress(u,1.,2.)
\end{python}
However getStress can also be called with Data objects as values for lam and mu which, for instance in the case of a temperature dependency, are calculated by an expression. The following call is equivalent to the previous example:
\begin{python}
lam=Scalar(1.,ContinuousFunction(mydomain))
mu=Scalar(2.,Function(mydomain))
s=getStress(u,lam,mu)
\end{python}

The function lam belongs to the continuous FunctionSpace function space but with g the function trace(g) is in the general FunctionSpace function space. In the evaluation of the product lam*trace(g) we have different function spaces (on the nodes versus in the centers) and at first glance we have incompatible data. esys.escript converts the arguments in an appropriate function space according to Table 3.1. In this example that means esys.escript sees lam as a function of the general FunctionSpace function space. In the context of FEM this means the nodal values of lam are interpolated to the element centers. The interpolation is automatic and requires no special handling.

Figure 3.2: Element Tagging. A rectangular mesh over a region with two rock types white and gray. The number in each cell refers to the major rock type present in the cell ($ 1$ for white and $ 2$ for gray).
\includegraphics[width=\textwidth]{figures/EscriptDiagram2}

Material parameters such as the Lame coefficients are typically dependent on rock types present in the area of interest. A common technique to handle these kinds of material parameters is "tagging", which uses storage efficiently. Figure 3.2 shows an example. In this case two rock types white and gray can be found in the domain. The domain is subdivided into triangular shaped cells. Each cell has a tag indicating the rock type predominately found in this cell. Here $ 1$ is used to indicate rock type white and $ 2$ for rock type gray. The tags are assigned at the time when the cells are generated and stored in the Domain class object. To allow easier usage of tags, names can be used instead of numbers. These names are typically defined at the time when the geometry is generated.

The following statements show how, for the example of Figure 3.2, the stress calculation discussed above and tagged values are used for lam:
\begin{python}
lam=Scalar(value=2.,what=Function(mydomain))
insertTaggedValue(lam,white=30.,gray=5000.)
s=getStress(u,lam,2.)
\end{python}
In this example lam is set to $ 30$ for those cells with tag white (=$ 1$) and to $ 5000.$ for those cells with tag gray (=$ 2$_. The initial value $ 2$ of lam is used as a default value for the case when a tag is encountered which has not been linked with a value. The getStress method does not need to be changed now that we are using tags. esys.escript resolves the tags when lam*trace(g) is calculated.

This brings us to a very important point about esys.escript. You can develop a simulation with constant Lame coefficients, and then later switch to tagged Lame coefficients without otherwise changing your python script. In short, you can use the same script to model with different domains and different types of input data.

There are three main ways in which Data objects are represented internally: constant, tagged, and expanded. In the constant case, the same value is used at each sample point and only a single value is stored to save memory. In the expanded case, each sample point has an individual value (such as for the solution of a PDE). This is where your largest data sets will be created because the values are stored as a complete array. The tagged case has already been discussed above.

Expanded data is created when you create a Data object with expanded=True. Tagged data sets are created when you use the insertTaggedValue() method as shown above.

Values are accessed through a sample reference number. Operations on expanded Data objects have to be performed for each sample point individually. When tagged values are used, the values are held in a dictionary. Operations on tagged data require processing the set of tagged values only, rather than processing the value for each individual sample point. esys.escript allows any mixture of constant, tagged and expanded data in a single expression.

Data objects can be written to disk files and read with dump and load, both of which use netCDF [16]. Use these to save data for visualization, checkpoint/restart or simply to save and reuse data that was expensive to compute.

For instance to save the coordinates of the data points of the continuous FunctionSpace to the file x.nc use
\begin{python}
x=ContinuousFunction(mydomain).getX()
x.dump(''x.nc'')
mydomain.dump(\lq dom.nc\lq )
\end{python}
To recover the object x and mydomain was a esys.finley mesh use
\begin{python}
from esys.finley import LoadMesh
mydomain=LoadMesh('dom.nc')
x=load(''x.nc'', mydomain)
\end{python}
It possible to rerun the mechanism that was originally used to generates mydomain to recreate mydomain. However in most cases using dump and load is faster in particular if optimization has been applied. In case that esys.escript is running on more than one MPI processor the dump will create an individual file for each processor containing the local data. In order to avoid conflicts the file name is extended by the MPI processor rank.

The function space of the Data is stored in x.nc, though. If the Data object is expanded, the number of data points in the file and of the Domain for the particular FunctionSpace must match. Moreover, the ordering of the values is checked using the reference identifiers provided by FunctionSpace on the Domain. In some cases, data points will be re-ordered. Take care to be sure you get what you want!



Subsections
esys@esscc.uq.edu.au