/[escript]/trunk/doc/cookbook/example01.tex
ViewVC logotype

Diff of /trunk/doc/cookbook/example01.tex

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

revision 2979 by ahallam, Tue Mar 9 02:54:32 2010 UTC revision 2982 by caltinay, Wed Mar 10 04:27:47 2010 UTC
# Line 32  cooler one until both Line 32  cooler one until both
32  blocks have the same temperature.  blocks have the same temperature.
33    
34  \subsection{1D Heat Diffusion Equation}  \subsection{1D Heat Diffusion Equation}
35  We can model the heat distribution of this problem over time using one  We can model the heat distribution of this problem over time using the one
36  dimensional heat diffusion equation\footnote{A detailed discussion on how the  dimensional heat diffusion equation\footnote{A detailed discussion on how the
37  heat diffusion equation is derived can be found at  heat diffusion equation is derived can be found at
38  \url{  \url{
# Line 58  derivatives in \refEq{eqn:hd}; $\frac{\p Line 58  derivatives in \refEq{eqn:hd}; $\frac{\p
58  change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is  change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is
59  the spatial change of temperature. As there is only a single spatial dimension  the spatial change of temperature. As there is only a single spatial dimension
60  to our problem, our temperature solution $T$ is only dependent on the time $t$  to our problem, our temperature solution $T$ is only dependent on the time $t$
61  and our signed distance from the the block-block interface $x$.  and our signed distance from the block-block interface $x$.
62    
63  \subsection{PDEs and the General Form}  \subsection{PDEs and the General Form}
64  It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact  It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact
# Line 75  Firstly, we discretise the PDE \refEq{eq Line 75  Firstly, we discretise the PDE \refEq{eq
75  steady linear PDE which involves spatial derivatives only and needs to be solved  steady linear PDE which involves spatial derivatives only and needs to be solved
76  in each time step to progress in time. \esc can help us here.  in each time step to progress in time. \esc can help us here.
77    
78  For time discretization we use the Backwards Euler approximation  For time discretisation we use the Backward Euler approximation
79  scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is  scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is
80  based on the  based on the approximation
 approximation  
81  \begin{equation}  \begin{equation}
82  \frac{\partial T(t)}{\partial t} \approx \frac{T(t)-T(t-h)}{h}  \frac{\partial T(t)}{\partial t} \approx \frac{T(t)-T(t-h)}{h}
83  \label{eqn:beuler}  \label{eqn:beuler}
# Line 107  T^{(n)}}{\partial x^{2}} = q\hackscore H Line 106  T^{(n)}}{\partial x^{2}} = q\hackscore H
106  Notice that we evaluate the spatial derivative term at the current time  Notice that we evaluate the spatial derivative term at the current time
107  $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively,  $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively,
108  one can evaluate the spatial derivative term at the previous time $t^{(n-1)}$.  one can evaluate the spatial derivative term at the previous time $t^{(n-1)}$.
109  This  This approach is called the \textbf{forward Euler} scheme. This scheme can
110  approach is called the \textbf{forward Euler} scheme. This scheme can provide  provide some computational advantages, which
 some computational advantages, which  
111  are not discussed here. However, the \textbf{forward Euler} scheme has a major  are not discussed here. However, the \textbf{forward Euler} scheme has a major
112  disadvantage. Namely, depending on the  disadvantage. Namely, depending on the
113  material parameters as well as the domain discretization of the spatial  material parameters as well as the domain discretization of the spatial
114  derivative term, the time step size $h$ needs to be chosen sufficiently small to  derivative term, the time step size $h$ needs to be chosen sufficiently small to
115  achieve a stable temperature when progressing in time. Stabiliy is achieved if  achieve a stable temperature when progressing in time. Stability is achieved if
116  the temperature does not grow beyond its initial bounds and become  the temperature does not grow beyond its initial bounds and becomes
117  non-physical.  non-physical.
118  The backward Euler scheme, which we use here, is unconditionally stable meaning  The backward Euler scheme, which we use here, is unconditionally stable meaning
119  that under the assumption of a  that under the assumption of a
120  physically correct problem set-up the temperature approximation remains physical  physically correct problem set-up the temperature approximation remains physical
121  for all time steps.  for all time steps.
122  The user needs to keep in mind that the discretization error introduced by  The user needs to keep in mind that the discretisation error introduced by
123  \refEq{eqn:beuler}  \refEq{eqn:beuler}
124  is sufficiently small, thus a good approximation of the true temperature is  is sufficiently small, thus a good approximation of the true temperature is
125  computed. It is  computed. It is therefore very important that any results are viewed with
126  therefore very important that any results are viewed with caution. For example,  caution. For example, one may compare the results for different time and
127  one may compare the results for different time and spatial step sizes.  spatial step sizes.
128    
129  To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear  To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear
130  differential equation \refEq{eqn:hddisc} which only includes spatial  differential equation \refEq{eqn:hddisc} which only includes spatial
131  derivatives. To solve this problem  derivatives. To solve this problem we want to use \esc.
 we want to use \esc.  
132    
133  In \esc any given PDE can be described by the general form. For the purpose of  In \esc any given PDE can be described by the general form. For the purpose of
134  this introduction we illustrate a simpler version of the general form for full  this introduction we illustrate a simpler version of the general form for full
# Line 188  respectively. Line 185  respectively.
185  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to
186  prescribe a known value to the unknown solution (in our example the temperature)  prescribe a known value to the unknown solution (in our example the temperature)
187  on parts of the boundary or on the entire boundary of the region of interest.  on parts of the boundary or on the entire boundary of the region of interest.
188  We discuss Dirichlet boundary condition in our second example presented in  We discuss the Dirichlet boundary condition in our second example presented in
189  Section~\ref{Sec:1DHDv0}.  Section~\ref{Sec:1DHDv0}.
190    
191  However, for this example we have made the model assumption that the system is  However, for this example we have made the model assumption that the system is
192  insulated, so we need  insulated, so we need to add an appropriate boundary condition to prevent
 to add an appropriate boundary condition to prevent  
193  any loss or inflow of energy at the boundary of our domain. Mathematically this  any loss or inflow of energy at the boundary of our domain. Mathematically this
194  is expressed by prescribing  is expressed by prescribing
195  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified
# Line 208  or in a more general case as Line 204  or in a more general case as
204  \end{equation}  \end{equation}
205  where $n$  is the outer normal field \index{outer normal field} at the surface  where $n$  is the outer normal field \index{outer normal field} at the surface
206  of the domain.  of the domain.
207  The $\cdot$ (dot) refers to the  dot product of the vectors $\nabla T$ and $n$.  The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$.
208  In fact, the term $\nabla T \cdot n$ is the normal derivative of  In fact, the term $\nabla T \cdot n$ is the normal derivative of
209  the temperature $T$. Other notations used here are\footnote{The \esc notation  the temperature $T$. Other notations used here are\footnote{The \esc notation
210  for the normal  for the normal
# Line 216  derivative is $T\hackscore{,i} n\hacksco Line 212  derivative is $T\hackscore{,i} n\hacksco
212  \begin{equation}  \begin{equation}
213  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .
214  \end{equation}  \end{equation}
215  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neumann boundary
216  condition} for the PDE.  condition} for the PDE.
217    
218  The PDE \refEq{eqn:hdgenf}  The PDE \refEq{eqn:hdgenf}
219  and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with  and the Neumann boundary condition~\ref{eqn:hdgenf} (potentially together with
220  the Dirichlet boundary conditions)  define a \textbf{boundary value problem}.  the Dirichlet boundary conditions)  define a \textbf{boundary value problem}.
221  It is the nature of a boundary value problem to allow making statements about  It is the nature of a boundary value problem to allow making statements about
222  the solution in the  the solution in the
223  interior of the domain from information known on the boundary only. In most  interior of the domain from information known on the boundary only. In most
224  cases  cases we use the term partial differential equation but in fact it is a
225  we use the term partial differential equation but in fact it is a boundary value  boundary value problem.
 problem.  
226  It is important to keep in mind that boundary conditions need to be complete and  It is important to keep in mind that boundary conditions need to be complete and
227  consistent in the sense that  consistent in the sense that
228  at any point on the boundary either a Dirichlet or a Neuman boundary condition  at any point on the boundary either a Dirichlet or a Neumann boundary condition
229  must be set.  must be set.
230    
231  Conveniently, \esc makes default assumption on the boundary conditions which the  Conveniently, \esc makes a default assumption on the boundary conditions which
232  user may modify where appropriate.  the user may modify where appropriate.
233  For a problem of the form in~\refEq{eqn:commonform nabla} the default  For a problem of the form in~\refEq{eqn:commonform nabla} the default
234  condition\footnote{In the form of the \esc users guide which is using the  condition\footnote{In the \esc user guide which uses the Einstein convention
235  Einstein convention is written as  this is written as
236  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;
237  \begin{equation}\label{NEUMAN 2}  \begin{equation}\label{NEUMAN 2}
238  -n\cdot A \cdot\nabla u = 0  -n\cdot A \cdot\nabla u = 0
# Line 252  condition translates into Line 247  condition translates into
247  \begin{equation}\label{NEUMAN 2b}  \begin{equation}\label{NEUMAN 2b}
248  \kappa \frac{\partial T}{\partial x} = 0  \kappa \frac{\partial T}{\partial x} = 0
249  \end{equation}  \end{equation}
250  for the boundary of the domain. This is identical to the Neuman boundary  for the boundary of the domain. This is identical to the Neumann boundary
251  condition we want to set. \esc will take care of this condition for us. We  condition we want to set. \esc will take care of this condition for us. We
252  discuss the Dirichlet boundary condition later.  discuss the Dirichlet boundary condition later.
253    
# Line 268  which has a rectangular shape. Secondly, Line 263  which has a rectangular shape. Secondly,
263  each time step to get the updated temperature. Thirdly, we need to define the  each time step to get the updated temperature. Thirdly, we need to define the
264  coefficients of the PDE and finally we need to solve the PDE. The last two steps  coefficients of the PDE and finally we need to solve the PDE. The last two steps
265  need to be repeated until the final time marker has been reached. The work flow  need to be repeated until the final time marker has been reached. The work flow
266  described in \reffig{fig:wf}.  is described in \reffig{fig:wf}.
267  % \begin{enumerate}  % \begin{enumerate}
268  %  \item create domain  %  \item create domain
269  %  \item create PDE  %  \item create PDE
# Line 284  described in \reffig{fig:wf}. Line 279  described in \reffig{fig:wf}.
279  \begin{figure}[h!]  \begin{figure}[h!]
280   \centering   \centering
281     \includegraphics[width=1in]{figures/workflow.png}     \includegraphics[width=1in]{figures/workflow.png}
282     \caption{Workflow for developing an \esc model and solution.}     \caption{Workflow for developing an \esc model and solution}
283     \label{fig:wf}     \label{fig:wf}
284  \end{figure}  \end{figure}
285    
# Line 310  advanced features, but these are not req Line 305  advanced features, but these are not req
305  \subsection{The Domain Constructor in \esc}  \subsection{The Domain Constructor in \esc}
306  \label{ss:domcon}  \label{ss:domcon}
307  Whilst it is not strictly relevant or necessary, a better understanding of  Whilst it is not strictly relevant or necessary, a better understanding of
308   how values are spatially distributed (\textit{e.g.} Temperature) and how PDE  how values are spatially distributed (\textit{e.g.} Temperature) and how PDE
309  coefficients are interpreted in \esc can be helpful.  coefficients are interpreted in \esc can be helpful.
310    
311  There are various ways to construct domain objects. The simplest form is a  There are various ways to construct domain objects. The simplest form is a
# Line 324  representation\footnote{We use the finit Line 319  representation\footnote{We use the finit
319  \url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}.  \url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}.
320  Therefore we will have access to an approximation of the true PDE solution  Therefore we will have access to an approximation of the true PDE solution
321  only.  only.
322  The quality of the approximation depends - besides other factors- mainly on the  The quality of the approximation depends - besides other factors - mainly on the
323  number of elements being used. In fact, the  number of elements being used. In fact, the
324  approximation becomes better when more elements are used. However, computational  approximation becomes better when more elements are used. However, computational
325  cost grows with the number of  cost grows with the number of
# Line 335  In general, one can think about a domain Line 330  In general, one can think about a domain
330  elements.  elements.
331  As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to  As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to
332  describe its vertices.  describe its vertices.
333  To represent spatial distributed values the user can use  To represent spatially distributed values the user can use
334  the values at the nodes, at the elements in the interior of the domain or at the  the values at the nodes, at the elements in the interior of the domain or at the
335  elements located at the surface of the domain.  elements located on the surface of the domain.
336  The different approach used to represent values is called \textbf{function  The different approach used to represent values is called \textbf{function
337  space} and is attached to all objects  space} and is attached to all objects
338  in \esc representing a spatial distributed value such as the solution of a PDE.  in \esc representing a spatially distributed value such as the solution of
339  The three  a PDE. The three function spaces we use at the moment are;
 function spaces we use at the moment are;  
340  \begin{enumerate}  \begin{enumerate}
341  \item the nodes, called by \verb|ContinuousFunction(domain)| ;  \item the nodes, called by \verb|ContinuousFunction(domain)| ;
342  \item the elements/cells, called by \verb|Function(domain)| ; and  \item the elements/cells, called by \verb|Function(domain)| ; and
343  \item the boundary, called by \verb|FunctionOnBoundary(domain)| .  \item the boundary, called by \verb|FunctionOnBoundary(domain)|.
344  \end{enumerate}  \end{enumerate}
345  A function space object such as \verb|ContinuousFunction(domain)| has the method  A function space object such as \verb|ContinuousFunction(domain)| has the method
346  \verb|getX| attached to it. This method returns the  \verb|getX| attached to it. This method returns the
# Line 354  location of the so-called \textbf{sample Line 348  location of the so-called \textbf{sample
348  particular function space. So the  particular function space. So the
349  call \verb|ContinuousFunction(domain).getX()| will return the coordinates of the  call \verb|ContinuousFunction(domain).getX()| will return the coordinates of the
350  nodes used to describe the domain while  nodes used to describe the domain while
351  the  \verb|Function(domain).getX()| returns the coordinates of numerical  \verb|Function(domain).getX()| returns the coordinates of numerical
352  integration points within elements, see  integration points within elements, see \reffig{fig:fs}.
 \reffig{fig:fs}.  
353    
354  This distinction between different representations of spatially distributed  This distinction between different representations of spatially distributed
355  values  values
# Line 366  The coefficients of a PDE do not need to Line 359  The coefficients of a PDE do not need to
359  \verb|Function()| type.  \verb|Function()| type.
360  On the other hand a temperature distribution must be continuous and needs to be  On the other hand a temperature distribution must be continuous and needs to be
361  represented with a \verb|ContinuousFunction()| function space.  represented with a \verb|ContinuousFunction()| function space.
362  An influx may only be defined at the boundary and is therefore a \verb  An influx may only be defined at the boundary and is therefore a
363  FunctionOnBoundary()  object.    \verb|FunctionOnBoundary()| object.  
364  \esc allows certain transformations of the function spaces. A \verb  \esc allows certain transformations of the function spaces. A
365  ContinuousFunction()  can be transformed into a \verb|FunctionOnBoundary()|  \verb|ContinuousFunction()| can be transformed into a
366  or \verb|Function()|. On the other hand there is not enough information in a  \verb|FunctionOnBoundary()| or \verb|Function()|. On the other hand there is
367  \verb FunctionOnBoundary()  to transform it to a \verb ContinuousFunction()  .  not enough information in a \verb|FunctionOnBoundary()| to transform it to a
368    \verb|ContinuousFunction()|.
369  These transformations, which are called \textbf{interpolation} are invoked  These transformations, which are called \textbf{interpolation} are invoked
370  automatically by \esc if needed.  automatically by \esc if needed.
371    
372  Later in this introduction we discuss how  Later in this introduction we discuss how
373  to define specific areas of geometry with different materials which are  to define specific areas of geometry with different materials which are
374  represented by different material coefficients such the  represented by different material coefficients such as the
375  thermal conductivities $\kappa$. A very powerful technique to define these types  thermal conductivities $\kappa$. A very powerful technique to define these types
376  of PDE  of PDE
377  coefficients is tagging. Blocks of materials and boundaries can be named and  coefficients is tagging. Blocks of materials and boundaries can be named and
# Line 390  Section~\ref{STEADY-STATE HEAT REFRACTIO Line 384  Section~\ref{STEADY-STATE HEAT REFRACTIO
384  \subsection{A Clarification for the 1D Case}  \subsection{A Clarification for the 1D Case}
385  \label{SEC: 1D CLARIFICATION}  \label{SEC: 1D CLARIFICATION}
386  It is necessary for clarification that we revisit our general PDE from  It is necessary for clarification that we revisit our general PDE from
387  \refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently  \refeq{eqn:commonform nabla} for a two dimensional domain. \esc is inherently
388  designed to solve problems that are greater than one dimension and so  designed to solve problems that are multi-dimensional and so
389  \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem.  \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem.
390  In the case of two spatial dimensions the \textit{Nabla operator} has in fact  In the case of two spatial dimensions the \textit{Nabla operator} has in fact
391  two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial  two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial
# Line 440  It is generally a good idea to import al Line 434  It is generally a good idea to import al
434  if the functions and classes required are known they can be specified  if the functions and classes required are known they can be specified
435  individually. The function \verb|LinearPDE| has been imported explicitly for  individually. The function \verb|LinearPDE| has been imported explicitly for
436  ease of use later in the script. \verb|Rectangle| is going to be our type of  ease of use later in the script. \verb|Rectangle| is going to be our type of
437  domain. The module \verb unitsSI  provides support for SI unit definitions with  domain. The module \verb|unitsSI| provides support for SI unit definitions with
438  our variables.  our variables.
439    
440  Once our library dependencies have been established, defining the problem  Once our library dependencies have been established, defining the problem
# Line 466  small we increase the accuracy of our so Line 460  small we increase the accuracy of our so
460  the number of calculations that must be solved per time step. This means more  the number of calculations that must be solved per time step. This means more
461  computational time is required to produce a solution. In this \textit{1D}  computational time is required to produce a solution. In this \textit{1D}
462  problem, the bar is defined as being 1 metre long. An appropriate step size  problem, the bar is defined as being 1 metre long. An appropriate step size
463  \verb|ndx| would be 1 to 10\% of the length. Our \verb|ndy| need only be 1, this  \verb|ndx| would be 1 to 10\% of the length. Our \verb|ndy| needs only be 1,
464  is because our problem stipulates no partial derivatives in the $y$ direction.  this is because our problem stipulates no partial derivatives in the $y$
465    direction.
466  Thus the temperature does not vary with $y$. Hence, the model parameters can be  Thus the temperature does not vary with $y$. Hence, the model parameters can be
467  defined as follows; note we have used the \verb unitsSI  convention to make sure  defined as follows; note we have used the \verb|unitsSI| convention to make sure
468  all our input units are converted to SI.  all our input units are converted to SI.
469  \begin{python}  \begin{python}
470  mx = 500.*m #meters - model length  mx = 500.*m #meters - model length
# Line 502  print "Expected Number of time outputs i Line 497  print "Expected Number of time outputs i
497  i=0 #loop counter  i=0 #loop counter
498  \end{python}  \end{python}
499  Now that we know our inputs we will build a domain using the  Now that we know our inputs we will build a domain using the
500   \verb Rectangle() function from \verb finley . The four arguments allow us to define our domain  \verb|Rectangle()| function from \FINLEY. The four arguments allow us to
501  \verb model  as:  define our domain \verb|model| as:
502  \begin{python}  \begin{python}
503  #generate domain using rectangle  #generate domain using rectangle
504  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
505  \end{python}  \end{python}
506  \verb blocks  now describes a domain in the manner of Section \ref{ss:domcon}.  \verb|blocks| now describes a domain in the manner of Section \ref{ss:domcon}.
507    
508  With a domain and all our required variables established, it is now possible to  With a domain and all the required variables established, it is now possible to
509  set up our PDE so that it can be solved by \esc. The first step is to define the  set up our PDE so that it can be solved by \esc. The first step is to define the
510  type of PDE that we are trying to solve in each time step. In this example it is  type of PDE that we are trying to solve in each time step. In this example it is
511  a single linear PDE\footnote{in comparison to a system of PDEs which we discuss  a single linear PDE\footnote{in contrast to a system of PDEs which we discuss
512  later.}. We also need to state the values of our general form variables.  later.}. We also need to state the values of our general form variables.
513  \begin{python}  \begin{python}
514  mypde=LinearPDE(blocks)  mypde=LinearPDE(blocks)
# Line 521  A=zeros((2,2))) Line 516  A=zeros((2,2)))
516  A[0,0]=kappa  A[0,0]=kappa
517  mypde.setValue(A=A, D=rhocp/h)  mypde.setValue(A=A, D=rhocp/h)
518  \end{python}  \end{python}
519  In a many cases it may be possible to decrease the computational time of the  In many cases it may be possible to decrease the computational time of the
520  solver if the PDE is symmetric.  solver if the PDE is symmetric.
521  Symmetry of a PDE is defined by;  Symmetry of a PDE is defined by;
522  \begin{equation}\label{eqn:symm}  \begin{equation}\label{eqn:symm}
# Line 529  A\hackscore{jl}=A\hackscore{lj} Line 524  A\hackscore{jl}=A\hackscore{lj}
524  \end{equation}  \end{equation}
525  Symmetry is only dependent on the $A$ coefficient in the general form and the  Symmetry is only dependent on the $A$ coefficient in the general form and the
526  other coefficients $D$ as well as the right hand side $Y$. From the above  other coefficients $D$ as well as the right hand side $Y$. From the above
527  definition we can see that our PDE is symmetric. The \verb LinearPDE  class  definition we can see that our PDE is symmetric. The \verb|LinearPDE| class
528  provides the method \method{checkSymmetry} to check if the given PDE is  provides the method \method{checkSymmetry} to check if the given PDE is
529  symmetric. As our PDE is symmetrical we enable symmetry via;  symmetric. As our PDE is symmetrical we enable symmetry via;
530  \begin{python}  \begin{python}
531   myPDE.setSymmetryOn()  myPDE.setSymmetryOn()
532  \end{python}  \end{python}
533  Next we need to establish the initial temperature distribution \verb|T|. We need  Next we need to establish the initial temperature distribution \verb|T|. We need
534  to  to
535  assign the value \verb|T1| to all sample points left to the contact interface at  assign the value \verb|T1| to all sample points left to the contact interface at
536  $x\hackscore{0}=\frac{mx}{2}$  $x\hackscore{0}=\frac{mx}{2}$
537  and the value \verb|T2| right to the contact interface. \esc  and the value \verb|T2| right to the contact interface. \esc
538  provides the \verb|whereNegative| function to construct this. In fact,  provides the \verb|whereNegative| function to construct this. More
539  \verb|whereNegative| returns the value $1$ at those sample points where the  specifically, \verb|whereNegative| returns the value $1$ at those sample points
540  argument  where the argument has a negative value. Otherwise zero is returned.
541  has a negative value. Otherwise zero is returned. If \verb|x| are the  If \verb|x| are the $x\hackscore{0}$
 $x\hackscore{0}$  
542  coordinates of the sample points used to represent the temperature distribution  coordinates of the sample points used to represent the temperature distribution
543  then \verb|x[0]-boundloc| gives us a negative value for  then \verb|x[0]-boundloc| gives us a negative value for
544  all sample points left to the interface and non-negative value to  all sample points left to the interface and non-negative value to
# Line 554  the right of the interface. So with; Line 548  the right of the interface. So with;
548  T= T1*whereNegative(x[0]-boundloc)+T2*(1-whereNegative(x[0]-boundloc))  T= T1*whereNegative(x[0]-boundloc)+T2*(1-whereNegative(x[0]-boundloc))
549  \end{python}  \end{python}
550  we get the desired temperature distribution. To get the actual sample points  we get the desired temperature distribution. To get the actual sample points
551  \verb|x| we use  \verb|x| we use the \verb|getX()| method of the function space
552  the  \verb|getX()| method of the function space \verb|Solution(blocks)|  \verb|Solution(blocks)| which is used to represent the solution of a PDE;
 which is used to represent the solution of a PDE;  
553  \begin{python}  \begin{python}
554  x=Solution(blocks).getX()  x=Solution(blocks).getX()
555  \end{python}  \end{python}
# Line 565  As \verb|x| are the sample points for th Line 558  As \verb|x| are the sample points for th
558  the initial temperature \verb|T| is using these sample points for  the initial temperature \verb|T| is using these sample points for
559  representation.  representation.
560  Although \esc is trying to be forgiving with the choice of sample points and to  Although \esc is trying to be forgiving with the choice of sample points and to
561  convert  convert
562  where necessary the adjustment of the function space is not always possible. So  where necessary the adjustment of the function space is not always possible. So
563  it is  it is advisable to make a careful choice on the function space used.  
 advisable to make a careful choice on the function space used.    
564    
565  Finally we initialise an iteration loop to solve our PDE for all the time steps  Finally we initialise an iteration loop to solve our PDE for all the time steps
566  we specified in the variable section. As the right hand side of the general form  we specified in the variable section. As the right hand side of the general form
# Line 584  while t < tend: Line 576  while t < tend:
576      totE = integrate(rhocp*T) #get the total heat (energy) in the system      totE = integrate(rhocp*T) #get the total heat (energy) in the system
577  \end{python}  \end{python}
578  The last statement in this script calculates the total energy in the system as  The last statement in this script calculates the total energy in the system as
579  volume integral  the volume integral of $\rho c\hackscore{p} T$ over the block.
580  of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy  As the blocks are insulated no energy should be lost or added.
 should be get lost or added.  
581  The total energy should stay constant for the example discussed here.  The total energy should stay constant for the example discussed here.
582    
583  \subsection{Running the Script}  \subsection{Running the Script}
584  The script presented so far is available under  The script presented so far is available under
585  \verb|example01a.py|. You can edit this file with your favourite text editor.  \verb|example01a.py|. You can edit this file with your favourite text editor.
586  On most operating systems\footnote{The you can use \texttt{run-escript} launcher  On most operating systems\footnote{The \texttt{run-escript} launcher is not
587  is not supported under {\it MS Windows} yet.} you can use the  supported under {\it MS Windows} yet.} you can use the
588  \program{run-escript} command  \program{run-escript} command
589  to launch {\it escript} scripts. For the example script use;  to launch {\it escript} scripts. For the example script use;
590  \begin{verbatim}  \begin{verbatim}
# Line 604  the python interpreter directly; Line 595  the python interpreter directly;
595  \begin{verbatim}  \begin{verbatim}
596  python example01a.py  python example01a.py
597  \end{verbatim}  \end{verbatim}
598  if the system is configured correctly (Please talk to your system  if the system is configured correctly (please talk to your system
599  administrator).  administrator).
600    
601  \begin{figure}  \begin{figure}
602  \begin{center}  \begin{center}
603  \includegraphics[width=4in]{figures/ttblockspyplot150}  \includegraphics[width=4in]{figures/ttblockspyplot150}
604  \caption{Example 1b: Total Energy in the Blocks over Time (in seconds).}  \caption{Example 1b: Total Energy in the Blocks over Time (in seconds)}
605  \label{fig:onedheatout1}  \label{fig:onedheatout1}
606  \end{center}  \end{center}
607  \end{figure}  \end{figure}
# Line 622  administrator). Line 613  administrator).
613  use a variety of free \pyt packages for visualisation.  use a variety of free \pyt packages for visualisation.
614  Two types will be demonstrated in this cookbook;  Two types will be demonstrated in this cookbook;
615  \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and  \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and
616  \verb VTK  \footnote{\url{http://www.vtk.org/}} visualisation.  \verb|VTK|\footnote{\url{http://www.vtk.org/}}.
617  The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}}  The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}}
618  and is good for basic graphs and plots.  and is good for basic graphs and plots.
619  For more complex visualisation tasks in particular, two and three dimensional  For more complex visualisation tasks, in particular two and three dimensional
620  problems we recommend the use of more advanced tools. For instance,  \mayavi  problems we recommend the use of more advanced tools. For instance,  \mayavi
621  \footnote{\url{http://code.enthought.com/projects/mayavi/}}  \footnote{\url{http://code.enthought.com/projects/mayavi/}}
622  which is based upon the \verb|VTK| toolkit. The usage of \verb|VTK| based  which is based upon the \verb|VTK| toolkit. The usage of \verb|VTK| based
623  visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two  visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two
624  dimensional PDE.  dimensional PDE.
625    
626  For our simple granite block problem, we have two plotting tasks. Firstly, we  For our simple granite block problem, we have two plotting tasks. Firstly, we
627  are interested in showing the  are interested in showing the
628  behaviour of the total energy over time and secondly, how the temperature  behaviour of the total energy over time and secondly, how the temperature
629  distribution within the block is  distribution within the block is developing over time.
630  developing over time. Lets start with the first task.  Let us start with the first task.
631    
632  The trick is to create a record of the time marks and the corresponding total  The idea is to create a record of the time marks and the corresponding total
633  energies observed.  energies observed.
634  \pyt provides the concept of lists for this. Before  \pyt provides the concept of lists for this. Before
635  the time loop is opened we create empty lists for the time marks \verb|t_list|  the time loop is opened we create empty lists for the time marks \verb|t_list|
# Line 660  while t<tend: Line 651  while t<tend:
651        E_list.append(totE) # add current total energy to record        E_list.append(totE) # add current total energy to record
652  \end{python}  \end{python}
653  To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs  To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs
654  to be loaded before used;  to be loaded before use;
655  \begin{python}  \begin{python}
656  import pylab as pl # plotting package.  import pylab as pl # plotting package.
657  \end{python}  \end{python}
658  Here we are not using the \verb|from pylab import *| in order to avoid name  Here we are not using \verb|from pylab import *| in order to avoid name
659  clashes for function names  clashes for function names within \esc.
 within \esc.  
660    
661  The following statements are added to the script after the time loop has been  The following statements are added to the script after the time loop has been
662  completed;  completed;
# Line 677  pl.axis([0,max(t_list),0,max(E_list)*1.1 Line 667  pl.axis([0,max(t_list),0,max(E_list)*1.1
667  pl.savefig("totE.png")  pl.savefig("totE.png")
668  \end{python}  \end{python}
669  The first statement hands over the time marks and corresponding total energies  The first statement hands over the time marks and corresponding total energies
670  to the plotter.  to the plotter.
671  The second statment is setting the title for the plot. The third statement  The second statement sets the title for the plot. The third statement
672  sets the axis ranges. In most cases these are set appropriately by the plotter.  sets the axis ranges. In most cases these are set appropriately by the plotter.
673    
674  The last statement renders the plot and writes the  The last statement generates the plot and writes the result into the file
675  result into the file \verb|totE.png| which can be displayed by (almost) any  \verb|totE.png| which can be displayed by (almost) any image viewer.
 image viewer.  
676  As expected the total energy is constant over time, see  As expected the total energy is constant over time, see
677  \reffig{fig:onedheatout1}.  \reffig{fig:onedheatout1}.
678    
# Line 691  As expected the total energy is constant Line 680  As expected the total energy is constant
680  \label{sec: plot T}  \label{sec: plot T}
681  \sslist{example01c.py}  \sslist{example01c.py}
682  For plotting the spatial distribution of the temperature we need to modify the  For plotting the spatial distribution of the temperature we need to modify the
683  strategy we have used  strategy we have used for the total energy.
684  for the total energy. Instead of producing a final plot at the end we will  Instead of producing a final plot at the end we will generate a
 generate a  
685  picture at each time step which can be browsed as a slide show or composed into  picture at each time step which can be browsed as a slide show or composed into
686  a movie.  a movie.
687  The first problem we encounter is that if we produce an image at each time step  The first problem we encounter is that if we produce an image at each time step
688  we need  we need to make sure that the images previously generated are not overwritten.
 to make sure that the images previously generated are not overwritten.  
689    
690  To develop an incrementing file name we can use the following convention. It is  To develop an incrementing file name we can use the following convention. It is
691  convenient to  convenient to put all image files showing the same variable - in our case the
692  put all image file showing the same variable - in our case the temperature  temperature distribution - into a separate directory.
693  distribution -  As part of the \verb|os| module\footnote{The \texttt{os} module provides
 into a separate directory. As part of the \verb|os| module\footnote{The  
 \texttt{os} module provides  
694  a powerful interface to interact with the operating system, see  a powerful interface to interact with the operating system, see
695  \url{http://docs.python.org/library/os.html}.} \pyt  \url{http://docs.python.org/library/os.html}.} \pyt
696  provides the  \verb|os.path.join| command to build file and  provides the \verb|os.path.join| command to build file and directory names in a
697  directory names in a platform independent way. Assuming that  platform independent way. Assuming that
698  \verb|save_path| is name of directory we want to put the results the command  \verb|save_path| is the name of the directory we want to put the results in the
699  is;  command is;
700  \begin{python}  \begin{python}
701  import os  import os
702  os.path.join(save_path, "tempT%03d.png"%i )  os.path.join(save_path, "tempT%03d.png"%i )
703  \end{python}  \end{python}
704  where \verb|i| is the time step counter.  where \verb|i| is the time step counter.
705  There are two arguments to the \verb join  command. The  There are two arguments to the \verb|join| command. The \verb|save_path|
706  \verb save_path  variable is a predefined string pointing to the directory we want to save our  variable is a predefined string pointing to the directory we want to save our
707  data, for example a single sub-folder called \verb data  would be defined by;  data, for example a single sub-folder called \verb|data| would be defined by;
708  \begin{verbatim}  \begin{verbatim}
709  save_path = "data"  save_path = "data"
710  \end{verbatim}  \end{verbatim}
711  while a sub-folder of \verb data  called \verb example01  would be defined by;  while a sub-folder of \verb|data| called \verb|example01| would be defined by;
712  \begin{verbatim}  \begin{verbatim}
713  save_path = os.path.join("data","example01")  save_path = os.path.join("data","example01")
714  \end{verbatim}  \end{verbatim}
715  The second argument of \verb join \xspace contains a string which is the file  The second argument of \verb|join| contains a string which is the file
716  name or subdirectory name. We can use the operator \verb|%| to use the value of  name or subdirectory name. We can use the operator \verb|%| to use the value of
717  \verb|i| as part of our filename. The sub-string \verb %03d  indicates that we  \verb|i| as part of our filename. The sub-string \verb|%03d| indicates that we
718  want to substitude a value into the name;  want to substitute a value into the name;
719  \begin{itemize}  \begin{itemize}
720   \item \verb 0  means that small numbers should have leading zeroes;   \item \verb 0  means that small numbers should have leading zeroes;
721   \item \verb 3  means that numbers should be written using at least 3 digits;   \item \verb 3  means that numbers should be written using at least 3 digits;
722  and  and
723   \item \verb d  means that the value to substitute will be an integer.   \item \verb d  means that the value to substitute will be a decimal integer.
724  \end{itemize}  \end{itemize}
725  To actually substitute the value of \verb|i| into the name write \verb %i  after  To actually substitute the value of \verb|i| into the name write \verb|%i| after
726  the string.  the string.
727  When done correctly, the output files from this command would be place in the  When done correctly, the output files from this command will be placed in the
728  directory defined by \verb save_path  as;  directory defined by \verb save_path  as;
729  \begin{verbatim}  \begin{verbatim}
730  blockspyplot001.png  blockspyplot001.png
# Line 753  A sub-folder check/constructor is availa Line 738  A sub-folder check/constructor is availa
738  \begin{verbatim}  \begin{verbatim}
739  mkDir(save_path)  mkDir(save_path)
740  \end{verbatim}  \end{verbatim}
741  will check for the existence of \verb save_path  and if missing, make the  will check for the existence of \verb save_path  and if missing, create the
742  required directories.  required directories.
743    
744  We start by modifying our solution script.  We start by modifying our solution script.
745  Prior to the \verb|while| loop we will need to extract our finite solution  Prior to the \verb|while| loop we need to extract our finite solution
746  points to a data object that is compatible with \mpl. First we create the node  points to a data object that is compatible with \mpl. First we create the node
747  coordinates of the sample points used to represent  coordinates of the sample points used to represent
748  the temperature as a \pyt list of tuples or a \numpy array as requested by the  the temperature as a \pyt list of tuples or a \numpy array as requested by the
# Line 780  plx = plx[:,0] # extract x locations Line 765  plx = plx[:,0] # extract x locations
765  \includegraphics[width=4in]{figures/blockspyplot050}  \includegraphics[width=4in]{figures/blockspyplot050}
766  \includegraphics[width=4in]{figures/blockspyplot200}  \includegraphics[width=4in]{figures/blockspyplot200}
767  \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps  \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps
768  $1$, $50$ and $200$.}  $1$, $50$ and $200$}
769  \label{fig:onedheatout}  \label{fig:onedheatout}
770  \end{center}  \end{center}
771  \end{figure}  \end{figure}
# Line 789  We use the same techniques provided by \ Line 774  We use the same techniques provided by \
774  energy over time.  energy over time.
775  For each time step we generate a plot of the temperature distribution and save  For each time step we generate a plot of the temperature distribution and save
776  each to a file.  each to a file.
777  The following is appended to the end of the \verb while  loop and creates one  The following is appended to the end of the \verb|while| loop and creates one
778  figure of the temperature distribution. We start by converting the solution to a  figure of the temperature distribution. We start by converting the solution to a
779  tuple and then plotting this against our \textit{x coordinates} \verb plx  we  tuple and then plotting this against our \textit{x coordinates} \verb|plx| we
780  have generated before. We add a title to the diagram before it is rendered into  have generated before. We add a title to the diagram before it is rendered into
781  a file.  a file.
782  Finally, the figure is saved to a \verb|*.png| file and cleared for the  Finally, the figure is saved to a \verb|*.png| file and cleared for the
# Line 812  while t<tend: Line 797  while t<tend:
797  \end{python}  \end{python}
798  Some results are shown in \reffig{fig:onedheatout}.  Some results are shown in \reffig{fig:onedheatout}.
799    
800  \subsection{Make a video}  \subsection{Making a Video}
801  Our saved plots from the previous section can be cast into a video using the  Our saved plots from the previous section can be cast into a video using the
802  following command appended to the end of the script. The \verb mencoder command  following command appended to the end of the script. The \verb mencoder command
803  is Linux only, so other platform users need to use an alternative video encoder.  is not available on every platform, so some users need to use an alternative
804    video encoder.
805  \begin{python}  \begin{python}
806  # compile the *.png files to create a *.avi videos that show T change  # compile the *.png files to create a *.avi video that shows T change
807  # with time. This operation uses Linux mencoder. For other operating  # with time. This operation uses Linux mencoder. For other operating
808  # systems it is possible to use your favourite video compiler to  # systems it is possible to use your favourite video compiler to
809  # convert image files to videos.  # convert image files to videos.

Legend:
Removed from v.2979  
changed lines
  Added in v.2982

  ViewVC Help
Powered by ViewVC 1.1.26