esys.escript: Scripting Environment for Finite Element Simulations in Earth Sciences

 

Further Information:

• Homepage of Escript here.

• More information including user documentation are available here.

• Contact: Lutz Gross

 

Overview

esys.escript is a software environment for the rapid development of 3D parallel simulations based on patial differential eqautions (PDEs). It has been developed as a module extension of the scripting language python. The finite element library finley has been specifically designed for solving large-scale problems on parallel computers with multi-core processors and has been incorporated as a differential equation solver into escript. In the escript programming model, python scripts orchestrate numerical algorithms that are parallelised behind-the-scenes in escript module calls, without low-level technical implementation by the escript user. In particular, the finley module extends the escript

functionality of data array abstraction with functions that assemble a system of linear equations from a linear boundary value problem defined over a mesh.

A direct application of the finite element method (FEM) is the solution of boundary value problems (BVPs) over a domain Ω. The escript module provides an environment to solve these problems through its core finite element library finley. Within escript, high-level numerical algorithms including specialised nonlinear solvers and time-differencing schemes can be rapidly developed. For instance, an unsteady initial BVP can be transformed into a sequence of steady BVPs to be solved at each time step. The coefficients of the steady, linear BVP can then be provided to finley to assemble a system matrix for the given domain and unstructured mesh.

The steady, linear second order BVP for an unknown function u is processed by finley in the following templated system of PDEs in tensorial notation:

-(A_ijkl u_k,l+B_ijku_k)_,j+C_ikl u_k,l+D_ik u_k=-X_ij,j+Y_i (1)

where u_k denote the k-th components of the solution u, and for any function Z, Z_,j denotes the derivative of Z with respect to the j-th spatial coordinate. The Einstein summation convention also applies where a duplicate subscript within a term implies summation over all possible values of that subscript. The coefficients A, B, C, D, X and Y are functions of their location in the physical domain that must be defined at the quadrature points within each finite element.

finley accepts a system of (natural) boundary conditions given by:

(A_ijkl u_k,l+B_ijku_k) n_j+d_ik u_k=X_ij n_j +y_i on Γ_i^N (2)

 

where n_j denotes the outer normal field of the domain and A, B and X are as for the PDE (1). The functions d and y are coefficients defined on the boundary, Γ. Here, Γ_i^N is a portion of the boundary where the boundary condition applies. Moreover, the Dirichlet boundary condition:

u_i=r_i on Γ_i^D

 

where Γ_i^D is subset of Γ such that Γ_i^D and Γ_i^N partion Γ, is also accepted. The right-hand side r is a function defined on the boundary. Within finley, the sets Γ_i^D are represented through a characteristic function q defined

q_i(x) = 1 if x is in Γ_i^D and q_i(x) = 0 otherwise (3).

The functions r and q are defined at the nodes of the finite element mesh. The BVP defined by equations (1), (2) and (3) will be referred to as the finley Boundary Value Template or finley BVT.

A special form of the Finley BVT for the case of single PDE for a scalar, unknown function u can be derived from (1) as:

-(A_jl u_,l+B_ju)_,j+ C_l u_,l+D u=-X_j,j+Y (4)

 

In this case, the accepted natural boundary condition (2) is given in the form:

(A_jl u_,l+B_ju) n_j+d_i u=X_j n_j +y on Γ^N (5)

 

and constraints are given by:

u=r on Γ^D .

 

finley obtains a discretisation of the finley BVT from the variational formulation. A solution u is sought which fulfills the constraints (3) and

∫_Ω v_i,j (A_ijkl u_k,l+B_ijku_k)+v_i (C_ikl u_k,l+D_ik u_k ) dΩ + ∫_Γ v_i d_ik u_k dΓ = ∫_Ω v_i,j X_ij+v_iX_i dΩ + ∫_Γ v_i y_i dΓ (6)

for all trial functions v with v_i=0 on Γ_i^D. The set Γ is a superset of all Γ_i^N where the functions d and y are set to zero on Γ-Γ_i^N.

finley uses isoparametric finite elements on unstructured meshes to discretise the variational problem (6). Available elements shapes are line, triangle, quadrilateral, tetrahedron, and hexahedron of orders one and two. The discretisation leads to a system of linear equations with a sparse coefficient system matrix. finley also provides a set of iterative Krylov subspace methods to solve the arising linear system.

A Simple Example

The current section demonstrates the use of finley and escript. Assume that a solution u to the diffusion equation

u_,t - u_,jj=f

for the scalar, time-dependent function u on the unit square is required. In this case, u_,t is the derivative of u with respect of time and we will assume f is a given source term identical to one. At an initial time, u is known to be a Gaussian profile with its peak at c=(½,½) and analytical form u(x)=0.1 e^-20|x-c|. The values of u at the bottom boundary x₁=0 are fixed in time.

Application of the backward Euler scheme

u^(t-Δt)=u^(t)-u_,t^(t) Δt

 

with step size Δt produces a steady linear differential equation

u^(t)-Δt u_,jj^(t)=Δt f +u^(t-Δt)

 

at each time step t. Equating like terms against the finley BVT for the single PDE (4) yields the template relations:

Ajl	=	Δtδjl
Bj	=	0
Cl	=	0
D	=	1
Xj	=	0
Y	=	Δt f +u(t-Δt)

 

 

As the implicit natural boundary condition can be written:

u_,jn_j=0

 

like terms are similarly equated against the finley BVT for boundary conditions (5) yielding further template relations:

d=y=0.

Finally, in implementation of the Dirichlet condition, the finley BVT equation (3) requires a characteristic function q which given by:

q(x₀,x₁) = 1 if x₁=0 and q(x₀,x₁) = 0 otherwise.

 

With the template parameters correctly identified, a python script is then written to solve this problem:

from esys.escript import *  from esys.escript.linearPDEs import LinearPDE from esys.finley import  Rectangle   # step size  dt=0.1  # create domain (20x20 mesh over [0,1]x[0,1])  dom=Rectangle(20,20,1,l0=1,l1=1)  # source term: constant 1.0  f=Scalar(1.0,Function(dom))  # get the locations in the domain x=dom.getX()  # initial value for u:  u=1./10.*exp(-20.*length(x-[0.5,0.5]))  # create the PDE object pde=LinearPDE(dom) pde.setValue(A=[[dt,0],[0,dt]],D=1,q=whereZero(x[1])) # start of iteration  t=dt  while t<=1.:   pde.setValue(Y=dt*f+u, r=u)    u=M.solve(b)     # next time step     t+=dt

f=Scalar(...) creates the source term as an escript data object describing a scalar value. The argument Function(dom) indicates the domain and the 'smoothness' of the function. In this case values are associates to the elements of the finite element mesh, dom. Values associated with the elements of a finite element mesh are actually associated with the element's quadrature nodes. The value of f is initialised to 1. escript avoids copying this constant value to the elements and stores the value efficiently in this case. The statement x=dom.getX() returns the coordinates of the nodes of the finite element mesh as an escript data object representing a vector value associated with the nodes. In the next statement, the initial value for u is set. escript automatically converts the list [0.5,0.5] into an escript data object associated with nodes in context with the subtraction from the data object x associated with nodes. The result is a data object representing a vector tied to nodes. The final result for u is a scalar object on nodes.

The LinearPDE object pde controls the calculation of the stiffness matrix and right hand side vector, in this case using the finite element discretisation defined by the mesh dom. It is expected that the characteristic functions q and the values of the constraints r are associated with nodes. The coefficients A, D and Y must be defined on elements. In this case, A and D are defined as a list and constant value, respectively. escript automatically converts the given values into matrix and scalar data objects associated with elements. In a manner similar to that described for f, the values are not copied to the elements but stored as single values used by all elements when referring to the corresponding coefficient. The situation for Y is slightly more complicated. In fact the value for Y is defined by an expression involving data tied to elements (f) and data tied to nodes (u). In this situation, pde object automatically invokes an interpolation of the data associated with the nodes to data associated with the elements before the assemblage is started. Notice, that f is set to a constant value (=1.0). As u will have an individual for each element the constant value is copied to each element before the addition can be performed. It is pointed out that in each iteration step only values of the PDE coefficients are specified that are modified over time. The pde object will try to reuse as much information – in this case the stiffness matrix.