ViewVC logotype

Annotation of /trunk/esys2/doc/user/firststep.tex

Parent Directory Parent Directory | Revision Log Revision Log

Revision 107 - (hide annotations)
Thu Jan 27 06:21:48 2005 UTC (18 years, 2 months ago) by jgs
File MIME type: application/x-tex
File size: 12361 byte(s)
commit of branch dev-01 back to main trunk on 2005-01-27

1 jgs 102 % $Id$
3     \chapter{The First Steps}
4     \label{FirstSteps}
6     \begin{figure}
7     \centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}
8     \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
9     \label{fig:FirstSteps.1}
10     \end{figure}
12 jgs 107 In this chapter we will give an introduction how to use \escript to solve
13     a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). The reader should be familiar with Python. The knowledge presented at the Python tutorial at \url{http://docs.python.org/tut/tut.html}
14     is sufficient. It is helpful if the reader has some basic knowledge of PDEs \index{partial differential equation}.
15 jgs 102
16 jgs 107 The PDE \index{partial differential equation} we wish to solve is the Poisson equation \index{Poisson equation}
17 jgs 102 \begin{equation}
18     -\Delta u =f
19     \label{eq:FirstSteps.1}
20     \end{equation}
21 jgs 107 for the solution $u$. The function $f$ is the given right hand side. The domain of interest, denoted by $\Omega$
22 jgs 102 is the unit square
23     \begin{equation}
24     \Omega=[0,1]^2=\{ (x\hackscore 0;x\hackscore 1) | 0\le x\hackscore{0} \le 1 \mbox{ and } 0\le x\hackscore{1} \le 1 \}
25     \label{eq:FirstSteps.1b}
26     \end{equation}
27 jgs 107 The domain is shown in \fig{fig:FirstSteps.1}.
28 jgs 102
29     $\Delta$ denotes the Laplace operator\index{Laplace operator} which is defined by
30     \begin{equation}
31     \Delta u = (u\hackscore {,0})\hackscore{,0}+(u\hackscore{,1})\hackscore{,1}
32     \label{eq:FirstSteps.1.1}
33     \end{equation}
34 jgs 107 where, for any function $w$ and any direction $i$, $u\hackscore{,i}$
35     denotes the partial derivative \index{partial derivative} of $u$ with respect to $i$.
36 jgs 102 \footnote{Some readers
37     may be more familiar with the Laplace operator\index{Laplace operator} being written
38     as $\nabla^2$, and written in the form
39     \begin{equation*}
40 jgs 107 \nabla^2 u = \nabla^t \cot \nabla u = \frac{\partial^2 u}{\partial x\hackscore 0^2}
41 jgs 102 + \frac{\partial^2 u}{\partial x\hackscore 1^2}
42     \end{equation*}
43     and \eqn{eq:FirstSteps.1} as
44     \begin{equation*}
45     -\nabla^2 u = f
46     \end{equation*}
47     }
48 jgs 107 Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect
49 jgs 102 to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
50     which leads to
51     \begin{equation}
52     \Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
53     \label{eq:FirstSteps.1.1b}
54     \end{equation}
55     In some cases, and we will see examples for this in the next chapter,
56     the usage of the nested $\sum$ symbols blows up the formulas and therefore
57 jgs 107 it is convenient to use the Einstein summation convention \index{summation convention}. This
58     drops the $\sum$ sign and assumes that a summation over a repeated index is performed
59 jgs 102 ("repeated index means summation"). For instance we write
60     \begin{eqnarray}
61     x\hackscore{i}y\hackscore{i}=\sum\hackscore{i=0}^2 x\hackscore{i}y\hackscore{i} \\
62     x\hackscore{i}u\hackscore{,i}=\sum\hackscore{i=0}^2 x\hackscore{i}u\hackscore{,i} \\
63     u\hackscore{,ii}=\sum\hackscore{i=0}^2 u\hackscore{,ii} \\
64 jgs 107 x\hackscore{ij}u\hackscore{i,j}=\sum\hackscore{j=0}^2\sum\hackscore{i=0}^2 x\hackscore{ij}u\hackscore{i,j} \\
65 jgs 102 \label{eq:FirstSteps.1.1c}
66     \end{eqnarray}
67     With the summation convention we can write the Poisson equation \index{Poisson equation} as
68     \begin{equation}
69     - u\hackscore{,ii} =1
70     \label{eq:FirstSteps.1.sum}
71     \end{equation}
72     On the boundary of the domain $\Omega$ the normal derivative $n\hackscore{i} u\hackscore{,i}$
73     of the solution $u$ shall be zero, ie. $u$ shall fulfill
74     the homogeneous Neumann boundary condition\index{Neumann
75     boundary condition!homogeneous}
76     \begin{equation}
77     n\hackscore{i} u\hackscore{,i}= 0 \;.
78     \label{eq:FirstSteps.2}
79     \end{equation}
80     $n=(n\hackscore{i})$ denotes the outer normal field
81     of the domain, see \fig{fig:FirstSteps.1}. Remember that we
82     are applying the Einstein summation convention \index{summation convention}, i.e
83 jgs 107 $n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
84     n\hackscore{1} u\hackscore{,1}$.
85 jgs 102 \footnote{Some readers may familiar with the notation
86     \begin{equation*}
87     \frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
88     \end{equation*}
89     for the normal derivative.}
90     The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
91     set $\Gamma^N$ which is the top and right edge of the domain:
92     \begin{equation}
93     \Gamma^N=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=1 \mbox{ or } x\hackscore{1}=1 \}
94     \label{eq:FirstSteps.2b}
95     \end{equation}
96 jgs 107 On the bottom and the left edge of the domain which is defined
97 jgs 102 as
98     \begin{equation}
99     \Gamma^D=\{(x\hackscore 0;x\hackscore 1) \in \Omega | x\hackscore{0}=0 \mbox{ or } x\hackscore{1}=0 \}
100     \label{eq:FirstSteps.2c}
101     \end{equation}
102     the solution shall be identically zero:
103     \begin{equation}
104     u=0 \; .
105     \label{eq:FirstSteps.2d}
106     \end{equation}
107 jgs 107 This kind of boundary condition is called a homogeneous Dirichlet boundary condition
108 jgs 102 \index{Dirichlet boundary condition!homogeneous}. The partial differential equation in \eqn{eq:FirstSteps.1.sum} together
109     with the Neumann boundary condition \eqn{eq:FirstSteps.2} and
110     Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
111     called boundary value
112 jgs 107 problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
113     the unknown
114 jgs 102 function $u$.
117     \begin{figure}
118     \centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}
119     \caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here
120     each element is a quadrilateral and described by four nodes, namely
121     the corner points. The solution is interpolated by a bi-linear
122     polynomial.}
123     \label{fig:FirstSteps.2}
124     \end{figure}
126     In general the BVP\index{boundary value problem!BVP} cannot be solved analytically and numerical
127     methods have to be used construct an approximation of the solution
128     $u$. Here we will use the finite element method\index{finite element
129     method} (FEM\index{finite element
130     method!FEM}). The basic idea is to fill the domain with a
131 jgs 107 set of points called nodes. The solution is approximated by its
132 jgs 102 values on the nodes\index{finite element
133 jgs 107 method!nodes}. Moreover, the domain is subdivided into small,
134     sub-domain called elements \index{finite element
135 jgs 102 method!element}. On each element the solution is
136     represented by a polynomial of a certain degree through its values at
137     the nodes located in the element. The nodes and its connection through
138     elements is called a mesh\index{finite element
139 jgs 107 method!mesh}. \fig{fig:FirstSteps.2} shows an
140 jgs 102 example of a FEM mesh with four elements in the $x_0$ and four elements
141     in the $x_1$ direction over the unit square.
142     For more details we refer the reader to the literature, for instance
143     \Ref{Zienc,NumHand}.
145     \escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.
146     (We will discuss a more general form of a PDE \index{partial differential equation!PDE}
147     that can be defined through the \LinearPDE class later). The instantiation of
148     a \Poisson class object requires the specification of the domain $\Omega$. In \escript
149     the \Domain class objects are used to describe the geometry of a domain but it also
150     contains information about the discretization methods and the actual solver which is used
151     to solve the PDE. Here we are using the FEM library \finley \index{finite element
152     method}. The following statements create the \Domain object \var{mydomain} from the
153     \finley method \method{Rectangle}
154     \begin{python}
155 jgs 107 from esys.finley import Rectangle
156     mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
157 jgs 102 \end{python}
158     In this case the domain is a rectangle with the lower, left corner at point $(0,0)$ and
159     the right, upper corner at $(\var{l0},\var{l1})=(1,1)$.
160 jgs 107 The arguments \var{n0} and \var{n1} define the number of elements in $x\hackscore{0}$ and
161 jgs 102 $x\hackscore{1}$-direction respectively. For more details on \method{Rectangle} and
162     other \Domain generators within the \finley module,
163     see \Chap{CHAPTER ON FINLEY}.
165 jgs 107 The following statements define the \Poisson class object \var{mypde} with domain \var{mydomain} and
166 jgs 102 the right hand side $f$ of the PDE to constant $1$:
167     \begin{python}
168 jgs 107 from esys.escript import Poisson
169     mypde = Poisson(mydomain)
170     mypde.setValue(f=1)
171 jgs 102 \end{python}
172     We have not specified any boundary condition but the
173     \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann
174     boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary
175     condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$
176     and any constant $C$ the function $u+C$ becomes a solution as well. We have to add
177     a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done
178 jgs 107 by defining a characteristic function \index{characteristic function}
179     which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set
180 jgs 102 and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c},
181 jgs 107 we need a function which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$:
182 jgs 102 \begin{python}
183     x=mydomain.getX()
184     gammaD=x[0].whereZero()+x[1].whereZero()
185     \end{python}
186 jgs 107 In the first statement, the method \method{getX} of the \Domain \var{mydomain}
187     gives access to locations
188     in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object
189     which we will learn more about later. \code{x[0]} returns the $x\hackscore{0}$ coordinates of the locations and
190 jgs 102 \code{x[0].whereZero()} creates function which equals $1$ where \code{x[0]} is (nearly) equal to zero
191 jgs 107 and $0$ elsewhere.
192     Similarly, \code{x[1].whereZero()} creates function which equals $1$ where \code{x[1]} is
193     equal to zero and $0$ elsewhere.
194     The sum of the results of \code{x[0].whereZero()} and \code{x[1].whereZero()} gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
195 jgs 102
196 jgs 107 The additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the
197 jgs 102 characteristic function \index{characteristic function} of the locations
198     of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}
199     are set. The complete definition of our example is now:
200     \begin{python}
201 jgs 107 from esys.linearPDEs import Poisson
202 jgs 102 x = mydomain.getX()
203     gammaD = x[0].whereZero()+x[1].whereZero()
204 jgs 107 mypde = Poisson(domain=mydomain)
205     mypde = setValue(f=1,q=gammaD)
206 jgs 102 \end{python}
207 jgs 107 The first statement imports the \Poisson class definition form the \linearPDEsPack module which is part of the \ESyS package.
208     To get the solution of the Poisson equation defined by \var{mypde} we just have to call its
209 jgs 102 \method{getSolution}.
211     Now we can write the script to solve our test problem (Remember that
212 jgs 107 lines starting with '\#' are comment lines in Python) (available as \file{mypoisson.py}
213 jgs 102 in the \ExampleDirectory):
214     \begin{python}
215 jgs 107 from esys.finley import Rectangle
216 jgs 102 from esys.linearPDEs import Poisson
217     # generate domain:
218 jgs 107 mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
219 jgs 102 # define characteristic function of Gamma^D
220     x = mydomain.getX()
221     gammaD = x[0].whereZero()+x[1].whereZero()
222     # define PDE and get its solution u
223     mypde = Poisson(domain=mydomain,f=1,q=gammaD)
224     u = mypde.getSolution()
225     # write u to an external file
226     u.saveDX("u.dx")
227     \end{python}
228     The last statement writes the solution to the external file \file{u.dx} in
229     \OpenDX file format. \OpenDX is a software package
230     for the visualization of scientific, engineering and analytical data and is freely available
231     from \url{http://www.opendx.org}.
233     \begin{figure}
234     \centerline{\includegraphics[width=\figwidth]{FirstStepResult.eps}}
235 jgs 107 \caption{\OpenDX Visualization of the Possion Equation Solution for $f=1$}
236 jgs 102 \label{fig:FirstSteps.3}
237     \end{figure}
239 jgs 107 You can edit the script file using your favourite text editor (or the Integrated DeveLopment Environment IDLE
240     for Python, see \url{http://idlefork.sourceforge.net}). If the script file has the name \file{mypoisson.py} \index{scripts!\file{mypoisson.py}} you can run the
241 jgs 102 script from any shell using the command:
242     \begin{verbatim}
243     python mypoisson.py
244     \end{verbatim}
245     After the script has (hopefully successfully) been completed you will find the file \file{u.dx} in the current
246     directory. An easy way to visualize the results is the command
247     \begin{verbatim}
248 jgs 107 dx -prompter &
249 jgs 102 \end{verbatim}
250 jgs 107 to start the generic data visualization interface of \OpenDX. \fig{fig:FirstSteps.3} shows the result.


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

  ViewVC Help
Powered by ViewVC 1.1.26