 # Diff of /trunk/doc/user/beam.tex

revision 567 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 568 by gross, Tue Feb 28 03:59:06 2006 UTC
# Line 1  Line 1
1  % $Id$  % $Id$
2
3    The \LinearPDE class is used to define a general linear, steady, second order PDE
4    for an unknown function $u$ on a given $\Omega$ defined through a \Domain object.
5    In the following $\Gamma$ denotes the boundary of the domain $\Omega$. $n$ denotes
6    the outer normal field on $\Gamma$.
7
8    For a single PDE with a solution with a single component the linear PDE is defined in the
9    following form:
10    \begin{equation}\label{LINEARPDE.SINGLE.1}
11    -(A\hackscore{jl} u\hackscore{,l}){,j}+(B\hackscore{j} u)\hackscore{,j}+C\hackscore{l} u\hackscore{,l}+D u =-X\hackscore{j,j}+Y \; .
12    \end{equation}
13    $u_{,j}$ denotes the derivative of $u$ with respect to the $j$-th spatial direction. Einstein's summation convention, ie. summation over indexes appearing twice in a term of a sum is performed, is used.
14    The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the
15    \Function on the PDE or objects that can be converted into such \Data objects.
16    $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar.
17    The following natural
18    boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
19    \begin{equation}\label{LINEARPDE.SINGLE.2}
20    n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y  \;.
21    \end{equation}
22    Notice that the coefficients $A$, $B$ and $X$ are defined in the PDE. The coefficients $d$ and $y$ are
23    each a \Scalar in the \FunctionOnBoundary.  Constraints \index{constraint} for the solution prescribing the value of the
24    solution at certain locations in the domain. They have the form
25    \begin{equation}\label{LINEARPDE.SINGLE.3}
26    u=r \mbox{ where } q>0
27    \end{equation}
28    $r$ and $q$ are each \Scalar where $q$ is the characteristic function
29    \index{characteristic function} defining where the constraint is applied.
30    The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1}
31    or \eqn{LINEARPDE.SINGLE.2}. The PDE is symmetrical \index{symmetrical} if
32    \begin{equation}\label{LINEARPDE.SINGLE.4}
33    A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j}
34    \end{equation}
35    For a system of PDEs and a solution with several components the PDE has the form
36    \begin{equation}\label{LINEARPDE.SYSTEM.1}
37    -(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u_k)\hackscore{,j}+C\hackscore{ikl} u\hackscore{k,l}+D\hackscore{ik} u_k =-X\hackscore{ij,j}+Y\hackscore{i} \; .
38    \end{equation}
39    $A$ is a \RankFour, $B$ and $C$ are each a \RankThree, $D$ and $X$ are each a \RankTwo and $Y$ is a \RankOne.
40    The natural boundary conditions \index{boundary condition!natural} take the form:
41    \begin{equation}\label{LINEARPDE.SYSTEM.2}
42    n\hackscore{j}(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u_k)+d\hackscore{ik} u_k=n\hackscore{j}-X\hackscore{ij}+y\hackscore{i}  \;.
43    \end{equation}
44    The coefficient $d$ is a \RankTwo and $y$ is a
45    \RankOne both in the \FunctionOnBoundary. Constraints \index{constraint} take the form
46    \begin{equation}\label{LINEARPDE.SYSTEM.3}
47    u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0
48    \end{equation}
49    $r$ and $q$ are each \RankOne. Notice that at some locations not necessarily all components must
50    have a constraint. The system of PDEs is symmetrical \index{symmetrical} if
51    \begin{eqnarray}\label{LINEARPDE.SYSTEM.4}
52    A\hackscore{ijkl}=A\hackscore{klij} \\
53    B\hackscore{ijk}=C\hackscore{kij} \\
54    D\hackscore{ik}=D\hackscore{ki} \\
55    d\hackscore{ik}=d\hackscore{ki} \
56    \end{eqnarray}
57    \LinearPDE also supports solution discontinuities \index{discontinuity} over contact region $\Gamma^{contact}$
58    in the domain $\Omega$. To specify the conditions across the discontinuity we are using the
59    generalised flux $J$ which is in the case of a systems of PDEs and several components of the solution
60    defined as
61    \begin{equation}\label{LINEARPDE.SYSTEM.5}
62    J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}-X\hackscore{ij}
63    \end{equation}
64    For the case of single solution component and single PDE $J$ is defined
65    \begin{equation}\label{LINEARPDE.SINGLE.5}
66    J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}-X\hackscore{j}
67    \end{equation}
68    In the context of discontinuities \index{discontinuity} $n$ denotes the normal on the
69    discontinuity pointing from side 0 towards side 1. For a system of PDEs
70    the contact condition takes the form
71    \begin{equation}\label{LINEARPDE.SYSTEM.6}
72    n\hackscore{j} J^{0}\hackscore{ij}=n\hackscore{j} J^{1}\hackscore{ij}=y^{contact}\hackscore{i} - d^{contact}\hackscore{ik} [u]\hackscore{k} \; .
73    \end{equation}
74    where $J^{0}$ and $J^{1}$ are the fluxes on side $0$ and side $1$ of the
75    discontinuity $\Gamma^{contact}$, respectively. $[u]$, which is the difference
76    of the solution at side 1 and at side 0, denotes the jump of $u$ across $\Gamma^{contact}$.
77    The coefficient $d^{contact}$ is a \RankTwo and $y^{contact}$ is a
78    \RankOne both in the \FunctionOnContactZero or \FunctionOnContactOne.
79    In case of a single PDE and a single component solution the contact condition takes the form
80    \begin{equation}\label{LINEARPDE.SINGLE.6}
81    n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact} - d^{contact}[u]
82    \end{equation}
83    In this case the the coefficient $d^{contact}$ and $y^{contact}$ are eaach \Scalar
84    both in the \FunctionOnContactZero or \FunctionOnContactOne.
85    ======================
86
87
88    We have used a special case of the \LinearPDE class, namely the
89    \Poisson class already in \Chap{FirstSteps}.
90    Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
91    and use the \method{getSolution} method of parent class to actually solve the problem.
92
93    The form of a single PDE that can be handled by the \LinearPDE class is
94    \begin{equation}\label{EQU.FEM.1}
95    -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .
96    \end{equation}
97    We show here the terms which are relevant for the Helmholtz problem.
98    The general form and systems is discussed in \Sec{SEC LinearPDE}.
99    $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients}
100    Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar.
101    They may be constant or may depend on their
102    location in the domain but must not depend on the unknown solution $u$.
103    The following natural boundary conditions \index{boundary condition!natural} that
104    are used in the \LinearPDE class have the form
105    \begin{equation}\label{EQU.FEM.2}
106    n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y  \;.
107    \end{equation}
108    where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that
109    the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients.
110
111    By inspecting the Helmholtz equation \index{Helmholtz equation}
112    we can easily assign values to the coefficients in the
113    general PDE of the \LinearPDE class:
114    \begin{equation}\label{DIFFUSION HELM EQ 3}
115    \begin{array}{llllll}
116    A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\
117    d=\eta & y= g &  \\
118    \end{array}
119    \end{equation}
120    $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
121    $i=j$ and $0$ otherwise.
122
123    We want to implement a
124    new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
125    is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,
126    $g$ rather than the general form given by \eqn{EQU.FEM.1}.
127    Python's mechanism of subclasses allows us to do this in a very simple way.
128    The \Poisson class of the \linearPDEs module,
129    which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
130    \LinearPDE class. That means that all methods (such as the \method{getSolution})
131    from the parent class \LinearPDE are available for any \Poisson object. However, new
132    methods can be added and methods of the parent class can be redefined. In fact,
133    the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
134    values to the coefficients of the PDE. This is exactly what we will do when we define
135    our new \class{Helmholtz} class:
136    \begin{python}
137    from esys.linearPDEs import LinearPDE
138    import numarray
139    class Helmholtz(LinearPDE):
140       def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
141            ndim=self.getDim()
142            kronecker=numarray.identity(ndim)
143            self.setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)
144    \end{python}
145    \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
146    of \LinearPDE which we have imported in the first line of the script.
147    We add the method \method{setValue} to the class which overwrites the
148    \method{setValue} method of the \LinearPDE class. The new method which has the
149    parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
150    maps the parameters of the coefficients of the general PDE defined
151    in \eqn{EQU.FEM.1}. We are actually using the \method{_LinearPDE__setValue} of
152    the \LinearPDE class.
153    The coefficient \var{A} involves the Kronecker symbol. We use the
154    \numarray function \function{identity} which returns a square matrix with ones on the
155    main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix.
156    The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the
157    PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to
158    put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.
159    You can use your favourite editor to create and edit the file.
160
161    An object of the \class{Helmholtz} class is created through the statments:
162    \begin{python}
163    from mytools import Helmholtz
164    mypde = Helmholtz(mydomain)
165    mypde.setValue(kappa=10.,omega=0.1,f=12)
166    u = mypde.getSolution()
167    \end{python}
168    In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure
169    that \file{mytools.py} is in the directory from where you started Python.
170    \var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are
171    assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
172    specified, the default values $0$ are used. \footnote{It would be better to use the default value
173    \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
174    would not be processed when the PDE is evaluated}. In the fourth statement the solution of the
175    PDE is returned.
176
177    To test our \class{Helmholtz} class on a rectangular domain
178    of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
179    we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that
180    the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
181    an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
182    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
183    So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}
184    \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory
185    implements this test problem using the \finley PDE solver:
186    \begin{python}
187    from mytools import Helmholtz
188    from esys.escript import Lsup
189    from esys.finley import Rectangle
190    #... set some parameters ...
191    kappa=1.
192    omega=0.1
193    eta=10.
194    #... generate domain ...
195    mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
196    #... open PDE and set coefficients ...
197    mypde=Helmholtz(mydomain)
198    n=mydomain.getNormal()
199    x=mydomain.getX()
200    mypde.setValue(kappa,omega,omega*x,eta,kappa*n+eta*x)
201    #... calculate error of the PDE solution ...
202    u=mypde.getSolution()
203    print "error is ",Lsup(u-x)
204    \end{python}
205    The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.
206    \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
207    is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument.
208    The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
209    used to approximate the solution and our solution is a linear function of the spatial coordinates one might
210    expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
211    However most PDE packages use an iterative solver which is terminated
212    when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
213    \method{setTolerance} of the \LinearPDE class.
214
215    \subsection{The Transition Problem}
216    \label{DIFFUSION TRANS SEC}
217    Now we are ready to solve the original time dependent problem. The main
218    part of the script is the loop over time $t$ which takes the following form:
219    \begin{python}
220    t=0
221    T=Tref
222    mypde=Helmholtz(mydomain)
223    while t<t_end:
224          mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)
225          T=mypde.getSolution()
226          t+=h
227    \end{python}
228    \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
229    in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde}
230    is created before the loop over time is entered while in each time step only the coefficients
231    are reset in each time step. This way some information about the representation of the PDE can be reused
232    \footnote{The efficience can be improved further by setting the coefficients in the operator
233    \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only updating the coefficients
234    in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
235    method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator.}. The variable \var{T}
236    holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
237    Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
238    temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
239    initial temperature distribution, which equal to $T\hackscore{ref}$.
240
241    The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
242    in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
243    \var{q0} is a fixed constant. The following script defines \var{q} as desired:
244    \begin{python}
245    from esys.escript import length
246    xc=[0.02,0.002]
247    r=0.001
248    x=mydomain.getX()
249    q=q0*(length(x-xc)-r).whereNegative()
250    \end{python}
251    \var{x} is a \Data class object of
252    the \escript module defining locations in the \Domain \var{mydomain}.
253    The \function{length()} imported from the \escript module returns the
254    Euclidean norm:
255    \begin{equation}\label{DIFFUSION HELM EQ 3aba}
256    \|y\|=\sqrt{
257    y\hackscore{i}
258    y\hackscore{i}
259    } = \function{esys.escript.length}(y)
260    \end{equation}
261    So \code{length(x-xc)} calculates the distances
262    of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
263    Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
264    converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of
265    a \Data class object, in this case the result of the expression
266    \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
267    zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
268
269    Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
270    \index{scripts!\file{diffusion.py}}:
271    \begin{python}
272    from mytools import Helmholtz
273    from esys.escript import Lsup
274    from esys.finley import Rectangle
275    #... set some parameters ...
276    xc=[0.02,0.002]
277    r=0.001
278    qc=50.e6
279    Tref=0.
280    rhocp=2.6e6
281    eta=75.
282    kappa=240.
283    tend=5.
284    # ... time, time step size and counter ...
285    t=0
286    h=0.1
287    i=0
288    #... generate domain ...
289    mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
290    #... open PDE ...
291    mypde=Helmholtz(mydomain)
292    # ... set heat source: ....
293    x=mydomain.getX()
294    q=qc*(length(x-xc)-r).whereNegative()
295    # ... set initial temperature ....
296    T=Tref
297    # ... start iteration:
298    while t<tend:
299          i+=1
300          t+=h
301          print "time step :",t
302          mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
303          T=mypde.getSolution()
304          T.saveDX("T%d.dx"%i)
305    \end{python}
306    The script will create the files \file{T.1.dx},
307     \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the
308    temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.
309
310    \begin{figure}
311    \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
312    \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
313    \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
314    \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
315    \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
316    \label{DIFFUSION FIG 2}
317    \end{figure}
318
319    An easy way to visualize the results is the command
320    \begin{verbatim}
321    dx -edit diffusion.net &
322    \end{verbatim}
323    where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.
324    Use the \texttt{Sequencer} to move forward and and backwards in time.
325    \fig{DIFFUSION FIG 2} shows the result for some selected time steps.
326    ====