 # Diff of /trunk/esys2/doc/user/introduction.tex

revision 99 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 1  Line 1
1  % $Id$  % $Id$
2
3    \chapter{The First Steps}
4
5  \section{Introduction}  \section{Introduction}
6
7  \subsection{Getting the software}  \subsection{Getting the software}
8

9  \escript, \ESyS, all freely available.  Where do people get \finley from?  \escript, \ESyS, all freely available.  Where do people get \finley from?
10
11    \begin{enumerate}
12    \item how to get the software
13    \item a few words about the general structure
14    \item installation
15    \end{enumerate}
16
17    \section{How to solve a linear PDE}
18
19    \begin{figure}
20    \centerline{\includegraphics[width=\figwidth]{FirstStepDomain}}
21    \caption{Domain $\Omega$ with outer normal field $n$.}
22    \label{fig:FirstSteps.1}
23    \end{figure}
24
25    \begin{figure}
26    \centerline{\includegraphics[width=\figwidth]{FirstStepMesh}}
27    \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here
28    each element is a quadrilateral and described by four nodes, namely
29    the corner points. The solution is interpolated by a bi-linear
30    polynomial.}
31    \label{fig:FirstSteps.2}
32    \end{figure}
33
34    We want to solve the \index{partial differential equation}(\index{PDE})
35    \begin{equation}
36    -\Delta u + \alpha u =f
37    \label{eq:FirstSteps.1}
38    \end{equation}
39    for the solution $u$ on the domain $\Omega$.  Here we assume that the
40    domain is the rectangle of length $1$ and height $2$, see
41    \fig{fig:FirstSteps.1}.  $\Delta$ denotes the \index{Laplace
42    operator} which is defined by
43    \begin{equation}
44    \Delta u = (u\hackscore {,1})\hackscore{,1}+(u\hackscore{,2})\hackscore{,2}
45    \label{eq:FirstSteps.1.1}
46    \end{equation}
47    where for any function $w$ and any direction $i$ $u\hackscore{,i}$
48    denotes the derivative of $w$ with respect to $i$.  $\alpha$ is a
49    known constant (we will set $\alpha=10$) and $f$ is a given function
50    which may depend on the location in the domain.  On the boundary of
51    the domain $\Omega$ the solution $u$ shall fulfill the so-called
52    homogeneous \index{Neumann boundary condition}
53    \begin{equation}
54    \frac{\partial u}{\partial n} = 0
55    \label{eq:FirstSteps.2}
56    \end{equation}
57    where $n=(n\hackscore1,n\hackscore2)$ denotes the outer normal field
58    of the domain, see Figure~\ref{fig:FirstSteps.1} and
59    \begin{equation}
60    \frac{\partial u}{\partial n} = n\hackscore1 u\hackscore{,1} +
61    n\hackscore2 u\hackscore{,2}
62    \label{eq:FirstSteps.2.1}
63    \end{equation}
64    denotes the normal derivative on the boundary.  The partial
65    differential \eqn{eq:FirstSteps.1}) together with the
66    boundary condition~(\eqn{eq:FirstSteps.2}) forms a boundary value
67    problem (\index{BVP}) for unknown function $u$.
68
69    In most cases, the BVP cannot be solved analytically and numerical
70    methods have to be used construct an approximation of the solution
71    $u$. Here we will use the \index{finite element method}
72    (\index{FEM}). The basic idea is to fill the domain with a set of
73    points, so called nodes. The solution is approximated by its values on
74    the nodes. Moreover, the domain is subdivide into small subdomain,
75    so-called elements. On each element the solution is represented by a
76    polynomial of a certain degree through its values at the nodes located
77    in the element. The nodes and its connection through elements is
78    called a \index{mesh}.  \fig{fig:FirstSteps.2} shows an example
79    of a FEM mesh with four elements in the $x_0$ and for elements in the
80    $x_1$ direction over a rectangular domain. On more details we referring
81    to the literature, for instance \cite{X1,X2,X3}.
82
83    \escript provides the class \linearPDE to define a
84    general linear, steady differential partial differential equation of
85    second order. We will discuss the most general form that can be
86    defined through \class{linearPDE} later. The components which are
87    relevant for us here is as follows:
88    \begin{equation}
89    -\sum\hackscore{i,j=0}^k (A\hackscore{ij}
90     u\hackscore{,j})\hackscore{,i} + D u = Y
91    \label{eq:FirstSteps.3}
92    \end{equation}
93    In this form $D$ and $Y$ are scalars and $A$ is a $k \times k$ matrix
94    where $k$ denotes the spatial dimension (in our example we have
95    $k=2$).  By comparing the template~(\ref{eq:FirstSteps.3}) with the
96    differential equation~(\ref{eq:FirstSteps.1}) we want to solve we can
97    immediately identify the appropriate values for $A$, $D$ and $Y$:
98    \begin{equation}
99    \begin{array}{lcccc}
100    A\hackscore{ij} & = & \delta\hackscore{ij}&  =&
101                     \left[
102                                      \begin{array}{cc}
103                                      1 & 0\\
104                                      0 & 1
105                                      \end{array}
106                     \right] \\
107    D & =& \alpha \\
108    Y & = & f \\
109    \end{array}
110    \label{eq:FirstSteps.3.1}
111    \end{equation}
112    When the PDE is defined via~(\eqn{eq:FirstSteps.3}), \class{linearPDE}
113    makes a particular assumptions about the boundary conditions:
114    \begin{equation}
115    -\sum\hackscore{i,j=0}^k n\hackscore{i} A\hackscore{ij}
116     u\hackscore{,j} = 0
117    \label{eq:FirstSteps.4}
118    \end{equation}
119    Note that this boundary condition does not require extra information
120    as it only refers to the coefficient $A$ which already appears in the
121    PDE and so natural for it. Therefore this boundary condition is called
122    a \index{natural boundary condition}.
123
124    We will set $\alpha=10$ and $f=10$ such that $u=1$ becomes the exact
125    solution of the boundary value problem.  We make this very simple
126    choice to be able to test our program as we can compare the result
127    with the exact solution. Later we will set $\alpha$ and $f$ to
128    functions of their locations in the domain in which case we will not
129    be able to give an analytic solution. However, after testing our
130    program on this very simple case, we can be confident that it working
131    correctly before we apply it is a more complicated situation.
132
133    This is the program to solve the boundary value problem: (Remember that
134    lines starting with '\#' are commend lines in Python)
135    %\verbatiminput{examples/FirstSteps1.py}
136    \begin{python}
137    # import ESyS and finley
138    from ESyS import *
139    import finley
140    # set a value for alpha:
141    alpha=10
142    # generate mesh:
143    mydomain=finley.Rectangle(n0=40,n1=20,l0=2.,l1=1.)
144    # generate a system:
145    mypde=linearPDE(A=[[1,0],[0,1]],D=alpha,Y=10,domain=mydomain)
146    # generate a test solution:
147    u=mypde.getSolution()
148    # calculate the error of the solution
149    error=u-1.
150    print "norm of the approximation error is ",Lsup(error)
151    \end{python}
152    Line 2 import \escript and a few other tools for \ESyS. In Line 3 is
153    importing \finley which is used to solve the partial differential
154    equation. In line 7, a rectangular domain of length $l\hackscore 0=2$
155    and height $l\hackscore 1=1$ is generated and subdivided in
156    $n\hackscore 0=40$ and $n\hackscore 1=20$ elements in $x\hackscore 0$
157    and $x\hackscore 1$ direction, respectively.
158
159    We are using a function of \finley. This determines later in the code
160    which solver for the PDE is actually being used to solve the PDE.
161
162    The solution is done in three steps:
163  \begin{enumerate}  \begin{enumerate}
164   \item general structure  \item generate a finite element mesh subdividing the domain into elements.
165   \item how to get the software  \item assemble the system of linear equations $Mu=b$ from the BVP
166   \item a few words about the general structure  \item solve the linear system to get $u$
\item installation
167  \end{enumerate}  \end{enumerate}
168    The returned $u$ is given an approximation of the solution of the BVP
169    at the nodes of the finite element mesh. The quality of the
170    approximation depends on the size of the elements of the finite
171    element mesh: As smaller the element size the better the
172    approximation.  In our example we know the solution of the BVP so we
173    can compare the returned approximation with the true solution.  In
174    fact, as the true solution is simple, we can expect that the
175    approximation is exact.
176
177    The first step imports the package \ESyS which includes among
178    others the module \finley.
179
180    \section{A Time Dependent Problem}
181
182    % \verbatiminput{exam/finley\hackscoretime.py}
183
184    \section{With Dirichlet Conditions}
185
186    % \verbatiminput{exam/finley\hackscoredirichlet.py}
187
188    \section{Systems of PDEs}
189
190    % \verbatiminput{exam/system\hackscoretime.py}
191
192    \section{Explicit Schemes}
193    % \verbatiminput{exam/explicit\hackscoretime.py}

Legend:
 Removed from v.99 changed lines Added in v.100