The module esys.escript

extensionescript Data manipulation

Figure 3.1: Dependency of Function Spaces. An arrow indicates that a function in the function space at the starting point can be interpreted as a function in the function space of the arrow target.

esys.escript is an extension of Python to handle functions represented by their values on data sample points for the geometrical region on which the function is defined. The region as well as the method which is used to interpolate value on the data sample points is defined by Domain class objects. For instance when using the finite element method (FEM) Domain object holds the information about the FEM mesh, eg. a table of nodes and a table of elements. Although Domain contains the discretization method to be used esys.escript does not use this information directly. Domain objects are created from a module which want to make use esys.escript, e.g. esys.finley.

The solution of a PDE is a function of its location in the domain of interest $\Omega$. When solving a partial differential equation (PDE) using FEM the solution is (piecewise) differentiable but, in general, its gradient is discontinuous. To reflect these different degrees of smoothness different representations of the functions are used. For instance; in FEM the displacement field is represented by its values at the nodes of the mesh, while the strain, which is the symmetric part of the gradient of the displacement field, is stored on the element centers. To be able to classify functions with respect to their smoothness, esys.escript has the concept of the "function space". 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:

The following generators for function spaces on a Domain mydomain are available:

• Solution(mydomain): solutions of a PDE.
• ReducedSolution(mydomain): solutions of a PDE with a reduced smoothness requirement.
• ContinuousFunction(mydomain): continuous functions, eg. a temperature distribution.
• Function(mydomain): general functions which are not necessarily continuous, eg. a stress field.
• FunctionOnBoundary(mydomain): functions on the boundary of the domain, eg. a surface pressure.
• FunctionOnContact0(mydomain): functions on side 0 of the discontinuity.
• FunctionOnContact1(mydomain): functions on side $1$ of the discontinuity.

The reduced smoothness for PDE solution is often used to fulfill the Ladyzhenskaya–-Babuska–-Brezzi condition [] 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. 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 any side of the 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 fulfils 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 in the general FunctionSpace function space is typically represented by its values at the element center, but in a finite difference scheme the edge midpoint of cells is preferred. Using the concept of function spaces allows the user to run the same script on different PDE solver libraries by just changing the creator of the Domain object. 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.

Data class objects store functions of the location in a domain. 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 to store the returned solutions.

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 s of integers. The length of s is the rank of the Data object and s[i] is the maximum value for the i-th index. 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$:

The initial value is the constant $1$ for all data sample points and all components.

Data objects can also be created from any numarray array or any object, such as a list of floating point numbers, that can be converted into a numarray.NumArray Reference [2]. The following two statements create objects which are equivalent to mydat:

In the first case the initial value is numarray.ones((2,3)) which generates a $2 \times 3$ matrix as a numarray.NumArray 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 numarray.NumArray 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):

• Scalar(0,Function(mydomain)) is the same as Data(0,Function(myDomain),(,)), e.g a temperature field.
• Vector(0,Function(mydomain))is the same as Data(0,Function(myDomain),(d)), e.g a velocity field.
• Tensor(0,Function(mydomain)) is the same as Data(0,Function(myDomain),(d,d)), eg. a stress field.
• Tensor4(0,Function(mydomain)) is the same as Data(0,Function(myDomain),(d,d,d,d)) eg. a Hook tensor field.

Here the initial value is 0 but any object that can be converted into a numarray.NumArray 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 unitary operations (eg. cos, sin, log) and can be combined by applying binary 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

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

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:

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 u is available, eg. by solving a PDE, getStress might be called in the following way:

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:

The function lam belongs to the continuous FunctionSpace function space but with g the function trace(g) is in the general FunctionSpace function space. Therefore the evaluation of the product lam*trace(g) in the stress calculation produces a problem, as both functions are represented differently, eg. in FEM lam by its values on the node, and in trace(g) by its values at the element centers. In the case of inconsistent function spaces of arguments in a binary operation, esys.escript interprets the arguments in the appropriate function space according to the inclusion defined in 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. Behind the scenes esys.escript calls the appropriate function from the PDE solver library.

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).

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". 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. The following statements show how for the example of Figure 3.2 and the stress calculation discussed before tagged values are used for lam:

In this example lam is set to $30$ for those cells with tag $1$ and to $5000.$ for those cells with tag $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. Note that the getStress method is called without modification. esys.escript resolves the tags when lam*trace(g) is calculated.

The Data class provides a transparent interface to various data representations and the translations between them. As shown in the example of stress calculation, this allows the user to develop and test algorithms for a simple case (for instance with the Lame coefficients as constants) and then without further modifications of the program code to apply the algorithm in a more complex application (for instance a definition of the Lame coefficients using tags). As described here, there are three ways in which Data objects are represented internally, constant, tagged, and expanded (other representations will become available in later versions of esys.escript): In the constant case, if the same value is used at each sample point a single value is stored to save memory and compute time. Any operation on this constant data will only be performed on the single value. In the expanded case, each sample point has an individual value, eg. the solution of a PDE, and the values are stored as a complete array. The tagged case has already been discussed above.

Values are accessed through a sample reference number. Operations on expanded Data objects have to be performed for each sample point individually. If tagged values are used 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 use of constant, tagged and expanded data in a single expression.