/[escript]/trunk/doc/user/diffusion.tex
ViewVC logotype

Diff of /trunk/doc/user/diffusion.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 567 by gross, Mon Feb 27 00:04:17 2006 UTC revision 568 by gross, Tue Feb 28 03:59:06 2006 UTC
# Line 23  In Section~\ref{DIFFUSION TRANS SEC} the Line 23  In Section~\ref{DIFFUSION TRANS SEC} the
23  solver for the temperature diffusion problem.  solver for the temperature diffusion problem.
24    
25  \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion}  \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
   
26  The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation  The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
27  in the interior of the domain is given by  in the interior of the domain is given by
28  \begin{equation}  \begin{equation}
# Line 117  takes the form Line 116  takes the form
116  \kappa u\hackscore{,i} n\hackscore{i} =  g - \eta u\mbox{ on } \Gamma  \kappa u\hackscore{,i} n\hackscore{i} =  g - \eta u\mbox{ on } \Gamma
117  \label{DIFFUSION HELM EQ 2}  \label{DIFFUSION HELM EQ 2}
118  \end{equation}  \end{equation}
119  The partial differential  The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
 \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}  
120  is called the Helmholtz equation \index{Helmholtz equation}.  is called the Helmholtz equation \index{Helmholtz equation}.
121    
122  We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the  We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the
123  Helmholtz equation. We have used a special case of the \LinearPDE class, namely the  Helmholtz equation. For a single PDE the \LinearPDE class supports the following form:
124  \Poisson class already in \Chap{FirstSteps}.  \begin{equation}\label{LINEARPDE.SINGLE.1}
125  Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation  -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+(B\hackscore{j} u)\hackscore{,j}+C\hackscore{l} u\hackscore{,l}+D u =-X\hackscore{j,j}+Y \; .
126  and use the \method{getSolution} method of parent class to actually solve the problem.  \end{equation}
127    The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the
128  The form of a single PDE that can be handled by the \LinearPDE class is  \Function on the PDE or objects that can be converted into such \Data objects.
129  \begin{equation}\label{EQU.FEM.1}  $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar.
130  -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .  The following natural
131  \end{equation}  boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
132  We show here the terms which are relevant for the Helmholtz problem.  \begin{equation}\label{LINEARPDE.SINGLE.2}
133  The general form and systems is discussed in \Sec{SEC LinearPDE}.    n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y  \;.
134  $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients}  \end{equation}
135  Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar.  Notice that the coefficients $A$, $B$ and $X$ are the same like in the PDE~\eqn{LINEARPDE.SINGLE.1}. The coefficients $d$ and $y$ are  
136  They may be constant or may depend on their  each a \Scalar in the \FunctionOnBoundary.  Constraints \index{constraint} for the solution prescribing the value of the
137  location in the domain but must not depend on the unknown solution $u$.  solution at certain locations in the domain. They have the form
138  The following natural boundary conditions \index{boundary condition!natural} that  \begin{equation}\label{LINEARPDE.SINGLE.3}
139  are used in the \LinearPDE class have the form  u=r \mbox{ where } q>0
140  \begin{equation}\label{EQU.FEM.2}  \end{equation}
141  n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y  \;.  $r$ and $q$ are each \Scalar where $q$ is the characteristic function
142  \end{equation}  \index{characteristic function} defining where the constraint is applied.
143  where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that  The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1}
144  the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients.  or \eqn{LINEARPDE.SINGLE.2}.
145    The \Poisson class of the \linearPDEs module,
146    which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
147    \LinearPDE class. The \linearPDEs module provides a \Helmholtz class but
148    we will make direct use of the general \LinearPDE class.
149    
150  By inspecting the Helmholtz equation \index{Helmholtz equation}  By inspecting the Helmholtz equation \index{Helmholtz equation}
151  we can easily assign values to the coefficients in the  we can easily assign values to the coefficients in the
# Line 155  d=\eta & y= g &  \\ Line 157  d=\eta & y= g &  \\
157  \end{array}  \end{array}
158  \end{equation}  \end{equation}
159  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
160  $i=j$ and $0$ otherwise.  $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present\footnote{There is a difference
161    in \escript of being not present and set to zero. As not present coefficients are not processed,
162    it is more efficient to leave a coefficient undefined insted assigning zero to it.}
163    
164    Defining and solving the Helmholtz equation is very easy now:
165    \begin{python}
166    from esys.escript import *
167    from linearPDEs import LinearPDE
168    mypde=LinearPDE(mydomain)
169    mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g)
170    u=mypde.getSolution()
171    \end{python}
172    where we assume that \code{mydomain} is a \Domain object and
173    \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values
174    typically \code{float} or \Data objects. The \method{setValue} method
175    assigns values to the coefficients of the general PDE. The \method{getSolution} method solves
176    the PDE and returns the solution \code{u} of the PDE. \code{kronecker} is \escript function
177    returning the Kronecker symbol.
178    
179    The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily.
180    If a value is assigned to a coefficint several times, the last assigned value is used when
181    the solution is calculated:
182    \begin{python}
183    mypde=LinearPDE(mydomain)
184    mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
185    mypde.setValue(D=omega,Y=f,y=g)
186    mypde.setValue(d=2*eta) # overwrites d=eta
187    u=mypde.getSolution()
188    \end{python}
189    In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the
190    PDE is called symmetric if
191    \begin{equation}\label{LINEARPDE.SINGLE.4}
192    A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j} \;.
193    \end{equation}
194    Note that $D$ and $d$ may have any value and the right hand side $X$, $Y$, $y$ as well as the constraints
195    are not relevant. The Helmholtz problem is symmetric.
196    The \LinearPDE class provides the method \code{checkSymmetry} method to check if the given PDE is symmetric.
197    \begin{python}
198    mypde=LinearPDE(mydomain)
199    mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
200    print mypde.checkSymmetry() # returns True
201    mypde.setValue(B=kronecker(mydomain)[0])
202    print mypde.checkSymmetry() # returns False
203    mypde.setValue(C=kronecker(mydomain)[0])
204    print mypde.checkSymmetry() # returns True
205    \end{python}
206    Unfortunately, a \code{checkSymmetry} is very expensive and is recommendable to use for
207    testing and debugging purposes only. The \code{setSymmetryOn} method is used to
208    declare a PDE symmetric:
209    \begin{python}
210    mypde = LinearPDE(mydomain)
211    mypde.setValue(A=kappa*kronecker(mydomain))
212    mypde.setSymmetryOn()
213    \end{python}
214    
215    
216    
217    
218    
219    
220    
221    =====
222    
223    Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
224    and use the \method{getSolution} method of parent class to actually solve the problem.
225    
226  We want to implement a  We want to implement a
227  new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but  new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
228  is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,  is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,
229  $g$ rather than the general form given by \eqn{EQU.FEM.1}.  $g$ rather than the general form given by \eqn{EQU.FEM.1}.
230  Python's mechanism of subclasses allows us to do this in a very simple way.  Python's mechanism of subclasses allows us to do this in a very simple way.
231  The \Poisson class of the \linearPDEs module,  That means that all methods (such as the \method{getSolution})
 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general  
 \LinearPDE class. That means that all methods (such as the \method{getSolution})  
232  from the parent class \LinearPDE are available for any \Poisson object. However, new  from the parent class \LinearPDE are available for any \Poisson object. However, new
233  methods can be added and methods of the parent class can be redefined. In fact,  methods can be added and methods of the parent class can be redefined. In fact,
234  the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign  the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
235  values to the coefficients of the PDE. This is exactly what we will do when we define  values to the coefficients of the PDE. This is exactly what we will do when we define
236  our new \class{Helmholtz} class:  our new \class{Helmholtz} class:
 \begin{python}  
 from esys.linearPDEs import LinearPDE  
 import numarray  
 class Helmholtz(LinearPDE):  
    def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)  
         ndim=self.getDim()  
         kronecker=numarray.identity(ndim)  
         self._LinearPDE_setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)  
 \end{python}  
 \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass  
 of \LinearPDE which we have imported in the first line of the script.  
 We add the method \method{setValue} to the class which overwrites the  
 \method{setValue} method of the \LinearPDE class. The new method which has the  
 parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments  
 maps the parameters of the coefficients of the general PDE defined  
 in \eqn{EQU.FEM.1}. We are actually using the \method{_LinearPDE__setValue} of  
 the \LinearPDE class.  
 The coefficient \var{A} involves the Kronecker symbol. We use the  
 \numarray function \function{identity} which returns a square matrix with ones on the  
 main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix.  
 The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the  
 PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to  
 put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.  
 You can use your favourite editor to create and edit the file.    
   
 An object of the \class{Helmholtz} class is created through the statments:  
 \begin{python}  
 from mytools import Helmholtz  
 mypde = Helmholtz(mydomain)  
 mypde.setValue(kappa=10.,omega=0.1,f=12)  
 u = mypde.getSolution()  
 \end{python}  
 In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure  
 that \file{mytools.py} is in the directory from where you started Python.  
 \var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are  
 assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are  
 specified, the default values $0$ are used. \footnote{It would be better to use the default value  
 \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and  
 would not be processed when the PDE is evaluated}. In the fourth statement the solution of the  
 PDE is returned.  
237    
238  To test our \class{Helmholtz} class on a rectangular domain  To test our \class{Helmholtz} class on a rectangular domain
239  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
# Line 218  the test solution becomes the solution o Line 242  the test solution becomes the solution o
242  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
243  $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.  
244  So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}  So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}
245  \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
246  implements this test problem using the \finley PDE solver:  implements this test problem using the \finley PDE solver:
247  \begin{python}  \begin{python}
248  from mytools import Helmholtz  from esys.escript import *
249  from esys.escript import Lsup  from linearPDEs import LinearPDE
250  from esys.finley import Rectangle  from esys.finley import Rectangle
251  #... set some parameters ...  #... set some parameters ...
252  kappa=1.  kappa=1.
# Line 231  eta=10. Line 255  eta=10.
255  #... generate domain ...  #... generate domain ...
256  mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)  mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
257  #... open PDE and set coefficients ...  #... open PDE and set coefficients ...
258  mypde=Helmholtz(mydomain)  mypde=LinearPDE(mydomain)
259    mypde.setSymmetryOn()
260  n=mydomain.getNormal()  n=mydomain.getNormal()
261  x=mydomain.getX()  x=mydomain.getX()
262  mypde.setValue(kappa,omega,omega*x[0],eta,kappa*n[0]+eta*x[0])  mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0],d=eta,y=kappa*n[0]+eta*x[0])
263  #... calculate error of the PDE solution ...  #... calculate error of the PDE solution ...
264  u=mypde.getSolution()  u=mypde.getSolution()
265  print "error is ",Lsup(u-x[0])  print "error is ",Lsup(u-x[0])

Legend:
Removed from v.567  
changed lines
  Added in v.568

  ViewVC Help
Powered by ViewVC 1.1.26