In this chapter we give an introduction how to use esys.escript to solve a partial differential equation (PDE ). We assume you are at least a little familiar with Python. The knowledge presented at the Python tutorial at http://docs.python.org/tut/tut.htmlis more than sufficient.
The PDE we wish to solve is the Poisson equation
denotes the Laplace operator, which is defined by
On the boundary of the domain the normal derivative
of the solution
shall be zero, ie.
shall fulfill
the homogeneous Neumann boundary condition
![]() |
In general the BVP cannot be solved analytically and numerical
methods have to be used construct an approximation of the solution
. Here we will use the finite element method (FEM). The basic idea is to fill the domain with a
set of points called nodes. The solution is approximated by its
values on the nodes. Moreover, the domain is subdivided into smaller
sub-domains called elements . On each element the solution is
represented by a polynomial of a certain degree through its values at
the nodes located in the element. The nodes and its connection through
elements is called a mesh. Figure 1.2 shows an
example of a FEM mesh with four elements in the
and four elements
in the
direction over the unit square.
For more details we refer the reader to the literature, for instance Reference [28,5].
The esys.escript solver we want to use to solve this problem is embedded into the python interpreter language. So you can solve the problem interactively but you will learn quickly
that is more efficient to use scripts which you can edit with your favorite editor.
To enter the escript environment you use escript command:
escriptwhich will pass you on to the python prompt
Python 2.5.2 (r252:60911, Oct 5 2008, 19:29:17) [GCC 4.3.2] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>>Here you can use all available python commands and language features, for instance
esys.escript provides the class Poisson to define a Poisson equation .
(We will discuss a more general form of a PDE
that can be defined through the LinearPDE class later). The instantiation of
a Poisson class object requires the specification of the domain . In esys.escript
the Domain class objects are used to describe the geometry of a domain but it also
contains information about the discretization methods and the actual solver which is used
to solve the PDE. Here we are using the FEM library esys.finley . The following statements create the Domain object mydomain from the
esys.finley method Rectangle
In this case the domain is a rectangle with the lower, left corner at point and
the right, upper corner at
.
The arguments n0 and n1 define the number of elements in
and
-direction respectively. For more details on Rectangle and
other Domain generators within the esys.finley module,
see Chapter 7.
The following statements define the Poisson class object mypde with domain mydomain and
the right hand side of the PDE to constant
:
We have not specified any boundary condition but the
Poisson class implicitly assumes homogeneous Neuman boundary conditions defined by Equation (1.11). With this boundary
condition the BVP we have defined has no unique solution. In fact, with any solution
and any constant
the function
becomes a solution as well. We have to add
a Dirichlet boundary condition . This is done
by defining a characteristic function
which has positive values at locations
where Dirichlet boundary condition is set
and 0 elsewhere. In our case of
defined by Equation (1.13),
we need to construct a function gammaD which is positive for the cases
or
. To get
an object x which contains the coordinates of the nodes in the domain use
The method getX of the Domain mydomain
gives access to locations
in the domain defined by mydomain. The object x is actually a Data object which will be
discussed in Chapter 3 in more detail. What we need to know here is that
x has rank (number of dimensions) and a shape (list of dimensions) which can be viewed by
calling the getRank and getShape methods:
This will print something like
The Data object also maintains type information which is represented by the
FunctionSpace of the object. For instance
will print
which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a esys.finley mesh.
To get the
coordinates of the locations we use the
statement
Object x0
is again a Data object now with rank 0 and
shape . It inherits the FunctionSpace from x:
will print
We can now construct a function gammaD which is only non-zero on the bottom and left edges
of the domain with
whereZero(x[0]) creates function which equals where x[0] is (almost) equal to zero
and 0 elsewhere.
Similarly, whereZero(x[1]) creates function which equals
where x[1] is
equal to zero and 0 elsewhere.
The sum of the results of whereZero(x[0]) and whereZero(x[1])
gives a function on the domain mydomain which is strictly positive where
or
is equal to zero.
Note that gammaD has the same rank , shape and FunctionSpace like x0 used to define it. So from
one gets
An additional parameter q of the setValue method of the Poisson class defines the
characteristic function of the locations
of the domain where homogeneous Dirichlet boundary condition
are set. The complete definition of our example is now:
The first statement imports the Poisson class definition from the esys.escript.linearPDEs module esys.escript package.
To get the solution of the Poisson equation defined by mypde we just have to call its
getSolution.
Now we can write the script to solve our Poisson problem
The question is what we do with the calculated solution u. Besides postprocessing, eg. calculating the gradient or the average value, which will be discussed later, plotting the solution is one one things you might want to do. esys.escript offers two ways to do this, both base on external modules or packages and so data need to converted
to hand over the solution. The first option is using the matplotlib module which allows plotting 2D results relatively quickly, see [12]. However, there are limitations when using this tool, eg. in problem size and when solving 3D problems. Therefore esys.escript provides a second options based on VTK files which is especially
designed for large scale and 3D problem and which can be read by a variety of software packages such as mayavi [2], VisIt [24].