# Diff of /trunk/doc/cookbook/onedheatdiff001.tex

revision 2877 by gross, Mon Jan 25 08:34:15 2010 UTC revision 2878 by gross, Thu Jan 28 00:59:30 2010 UTC
# Line 14  Line 14
14  We will start by examining a simple one dimensional heat diffusion example. This problem will provide a good launch pad to build our knowledge of \esc and demonstrate how to solve simple partial differential equations (PDEs)\footnote{Wikipedia provides an excellent and comprehensive introduction to \textit{Partial Differential Equations} \url{http://en.wikipedia.org/wiki/Partial_differential_equation}, however their relevance to \esc and implementation should become a clearer as we develop our understanding further into the cookbook.}  We will start by examining a simple one dimensional heat diffusion example. This problem will provide a good launch pad to build our knowledge of \esc and demonstrate how to solve simple partial differential equations (PDEs)\footnote{Wikipedia provides an excellent and comprehensive introduction to \textit{Partial Differential Equations} \url{http://en.wikipedia.org/wiki/Partial_differential_equation}, however their relevance to \esc and implementation should become a clearer as we develop our understanding further into the cookbook.}
15
16  \section{One Dimensional Heat Diffusion in an Iron Rod}  \section{One Dimensional Heat Diffusion in an Iron Rod}
17  \sslist{onedheatdiff001.py and cblib.py}
18  %\label{Sec:1DHDv0}  %\label{Sec:1DHDv0}
19  The first model consists of a simple cold iron bar at a constant temperature of zero \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end. Intuition tells us that as heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source.  The first model consists of a simple cold iron bar at a constant temperature of zero \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end. Intuition tells us that as heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source.
20  \begin{figure}[h!]  \begin{figure}[h!]
# Line 119  A = \kappa; D = \frac{\rho c \hackscore{ Line 119  A = \kappa; D = \frac{\rho c \hackscore{
119  % for all $x$ in the domain.  % for all $x$ in the domain.
120
121  \subsection{Boundary Conditions}  \subsection{Boundary Conditions}
122    \label{SEC BOUNDARY COND}
123  With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary conditions\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, respectively.  With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary conditions\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, respectively.
124  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown - in our example the temperature - on parts of or the entire boundary of the region of interest.  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown - in our example the temperature - on parts of or the entire boundary of the region of interest.
125  For our model problem we want to keep the initial temperature setting on the left side of the  For our model problem we want to keep the initial temperature setting on the left side of the
126  iron bar over time. This defines a Dirichlet boundary condition for the PDE \refEq{eqn:hddisc} to be solved at each time step.  iron bar over time. This defines a Dirichlet boundary condition for the PDE \refEq{eqn:hddisc} to be solved at each time step.
127
128  On the other end of the iron rod we want to add an appropriate boundary condition to define insolation to prevent  On the other end of the iron rod we want to add an appropriate boundary condition to define insulation to prevent
129  any loss or inflow of energy at the right end of the rod. Mathematically this is expressed by prescribing  any loss or inflow of energy at the right end of the rod. Mathematically this is expressed by prescribing
130  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero on the right end of the rod  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero on the right end of the rod
131  In our simplified one dimensional model this is expressed  In our simplified one dimensional model this is expressed
# Line 154  we use the term partial differential equ Line 155  we use the term partial differential equ
155  It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that  It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that
156  at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set.  at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set.
157
158  Conviniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate.  Conveniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate.
159  For a problem of the form in~\refEq{eqn:commonform nabla} the default condition\footnote{In the form of the \esc users guide which is using the Einstein convention is written as  For a problem of the form in~\refEq{eqn:commonform nabla} the default condition\footnote{In the form of the \esc users guide which is using the Einstein convention is written as
160  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;
161  \label{NEUMAN 2}  \label{NEUMAN 2}
# Line 173  for the right hand side of the rod. This Line 174  for the right hand side of the rod. This
174  \label{sec:outline}  \label{sec:outline}
175  To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will write a simple \pyt script. At this point we assume that you have some basic understanding of the \pyt programming language. If not there are some pointers and links available in Section \ref{sec:escpybas}. The script we will discuss later in details will have four major steps. Firstly we need to define the domain where we want to  To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will write a simple \pyt script. At this point we assume that you have some basic understanding of the \pyt programming language. If not there are some pointers and links available in Section \ref{sec:escpybas}. The script we will discuss later in details will have four major steps. Firstly we need to define the domain where we want to
176  calculate the temperature. For our problem this is the iron rod which has a rectangular shape. Secondly we need to define the PDE  calculate the temperature. For our problem this is the iron rod which has a rectangular shape. Secondly we need to define the PDE
177  we need to solve in each time step to get the updated temperature. Thirdly we need to define the the coefficients of the PDE and finally we need to solve the PDE. The last two steps need to be repeated until the final time marker has been reached.  we need to solve in each time step to get the updated temperature. Thirdly we need to define the the coefficients of the PDE and finally we need to solve the PDE. The last two steps need to be repeated until the final time marker has been reached. As a work flow this takes the form;
178    \begin{enumerate}
179     \item create domain
180     \item create PDE
181     \item while end time not reached:
182    \begin{enumerate}
183     \item set PDE coefficients
184     \item solve PDE
185     \item update time marker
186    \end{enumerate}
187    \item end of calculation
188    \end{enumerate}
189  In the terminology of \pyt the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it defined by it usage and features  In the terminology of \pyt the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it defined by it usage and features
190  rather than its actual representation. So we will create a domain object to decribe the geometry of our iron rod. The main feature  rather than its actual representation. So we will create a domain object to describe the geometry of our iron rod. The main feature
191  of the object we will use is the fact that we can define PDEs and spatially distributed values such as the temperature  of the object we will use is the fact that we can define PDEs and spatially distributed values such as the temperature
192  on a domain. In fact the domain object has many more features - most of them you will  on a domain. In fact the domain object has many more features - most of them you will
193  never use and do not need to understand. Similar a PDE object is defined by the fact that we can define the coefficients of the PDE and solve the PDE. At a  never use and do not need to understand. Similar a PDE object is defined by the fact that we can define the coefficients of the PDE and solve the PDE. At a
194  later stage you may use more advanved features of the PDE class but you need to worry about them only at the point when you use them.  later stage you may use more advanced features of the PDE class but you need to worry about them only at the point when you use them.
195
196
197  \begin{figure}[t]  \begin{figure}[t]
# Line 197  from the user's point of view the repres Line 208  from the user's point of view the repres
208
209  There are various ways to construct domain objects. The simplest form is as rectangular shaped region with a length and height. There is  There are various ways to construct domain objects. The simplest form is as rectangular shaped region with a length and height. There is
210  a ready to use function call for this. Besides the spatial dimensions the function call will require you to specify the number  a ready to use function call for this. Besides the spatial dimensions the function call will require you to specify the number
211  elements or cells to be used along the length and height, see Figure~\ref{fig:fs}. Any spacially distributed value  elements or cells to be used along the length and height, see Figure~\ref{fig:fs}. Any specially distributed value
212  and the PDE is represented in discrete form using this element representation\footnote{We will use the finite element method (FEM), see \url{http://en.wikipedia.org/wiki/Finite_element_method}i for details.}. Therefore we will have access to an approximation of the true PDE solution only.  and the PDE is represented in discrete form using this element representation\footnote{We will use the finite element method (FEM), see \url{http://en.wikipedia.org/wiki/Finite_element_method}i for details.}. Therefore we will have access to an approximation of the true PDE solution only.
213  The quality of the approximation depends - besides other factors- mainly on the number of elements being used. In fact, the  The quality of the approximation depends - besides other factors- mainly on the number of elements being used. In fact, the
214  approximation becomes better the more elements are used. However, computational costs and compute time grow with the number of  approximation becomes better the more elements are used. However, computational costs and compute time grow with the number of
215  elements being used. It therefore important that you find the right balance between the demand in accuarcy and acceptable resource usage.  elements being used. It therefore important that you find the right balance between the demand in accuracy and acceptable resource usage.
216
217  In general, one can thinks about a domain object as a composition of nodes and elements.  In general, one can thinks about a domain object as a composition of nodes and elements.
218  As shown in Figure~\ref{fig:fs}, an element is defined by the nodes used to describe its vertices.  As shown in Figure~\ref{fig:fs}, an element is defined by the nodes used to describe its vertices.
# Line 221  call \verb|ContinuousFunction(domain).ge Line 232  call \verb|ContinuousFunction(domain).ge
232  the  \verb|Function(domain).getX()| returns the coordinates of numerical integration points within elements, see  the  \verb|Function(domain).getX()| returns the coordinates of numerical integration points within elements, see
233  Figure~\ref{fig:fs}.  Figure~\ref{fig:fs}.
234
235  This distiction between different representations of spatial distributed values  This distinction between different representations of spatial distributed values
236  is important in order to be able to vary the degrees of smoothness in a PDE problem.  is important in order to be able to vary the degrees of smoothness in a PDE problem.
237  The coefficients of a PDE need not be continuous thus this qualifies as a \verb|Function()| type.  The coefficients of a PDE need not be continuous thus this qualifies as a \verb|Function()| type.
238  On the other hand a temperature distribution must be continuous and needs to be represented with a \verb|ContinuousFunction()| function space.  On the other hand a temperature distribution must be continuous and needs to be represented with a \verb|ContinuousFunction()| function space.
239  An influx may only be defined at the boundary and is therefore a \verb FunctionOnBoundary()  object.    An influx may only be defined at the boundary and is therefore a \verb FunctionOnBoundary()  object.
240  \esc allows certain transformations of the function spaces. A \verb ContinuousFunction()  can be transformed into a \verb|FunctionOnBoundary()|  \esc allows certain transformations of the function spaces. A \verb ContinuousFunction()  can be transformed into a \verb|FunctionOnBoundary()|
241  or \verb|Function()|. On the other hand there is not enough information in a \verb FunctionOnBoundary()  to transform it to a \verb ContinuousFunction()  .  or \verb|Function()|. On the other hand there is not enough information in a \verb FunctionOnBoundary()  to transform it to a \verb ContinuousFunction()  .
242  These transformations, which ar ecalled \textbf{interpolation} are invoked autmatically by \esc if needed.  These transformations, which are called \textbf{interpolation} are invoked automatically by \esc if needed.
243
244  Later in this introduction we will discuss how  Later in this introduction we will discuss how
245  to define specific areas of geometry with different materials which are represented by different material coefficients such the  to define specific areas of geometry with different materials which are represented by different material coefficients such the
246  thermal conductivities $kappa$. A very powerful techique to define these types of PDE  thermal conductivities $kappa$. A very powerful technique to define these types of PDE
247  coefficients is tagging. Blocks of materials and boundaries can be named and values can be defined on subregions based on their names.  coefficients is tagging. Blocks of materials and boundaries can be named and values can be defined on subregions based on their names.
248  This is simplifing PDE coefficient and flux definitions. It makes for much easier scripting. We will discuss this technique in Section~\ref{STEADY-STATE HEAT REFRACTION}.  This is simplifying PDE coefficient and flux definitions. It makes for much easier scripting. We will discuss this technique in Section~\ref{STEADY-STATE HEAT REFRACTION}.
249
250
251  \subsection{A Clarification for the 1D Case}  \subsection{A Clarification for the 1D Case}
# Line 258  A\hackscore{00}=A; A\hackscore{01}=A\hac Line 269  A\hackscore{00}=A; A\hackscore{01}=A\hac
269
270  \subsection{Developing a PDE Solution Script}  \subsection{Developing a PDE Solution Script}
271  \label{sec:key}  \label{sec:key}
272    \sslist{onedheatdiffbase.py}
273  We will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules.  We will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules.
274  By developing a script for \esc, the heat diffusion equation can be solved at successive time steps for a predefined period using our general form \refEq{eqn:hdgenf}. Firstly it is necessary to import all the libraries\footnote{The libraries contain predefined scripts that are required to solve certain problems, these can be simple like $sine$ and $cosine$ functions or more complicated like those from our \esc library.}  By developing a script for \esc, the heat diffusion equation can be solved at successive time steps for a predefined period using our general form \refEq{eqn:hdgenf}. Firstly it is necessary to import all the libraries\footnote{The libraries contain predefined scripts that are required to solve certain problems, these can be simple like $sine$ and $cosine$ functions or more complicated like those from our \esc library.}
275  that we will require.  that we will require.
# Line 270  from esys.finley import Rectangle Line 282  from esys.finley import Rectangle
282  # A useful unit handling package which will make sure all our units  # A useful unit handling package which will make sure all our units
283  # match up in the equations under SI.  # match up in the equations under SI.
284  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
import pylab as pl #Plotting package.
import numpy as np #Array package.
import os #This package is necessary to handle saving our data.
285  \end{python}  \end{python}
286  It is generally a good idea to import all of the \modescript library, although if the functions and classes required are known they can be specified individually. The function \verb|LinearPDE| has been imported explicitly for ease of use later in the script. \verb|Rectangle| is going to be our type of model. The module \verb unitsSI  provides support for SI unit definitions with our variables; and the \verb|os| module is needed to handle file outputs once our PDE has been solved. \verb pylab  and \verb numpy  are modules developed independently of \esc. They are used because they have efficient plotting and array handling capabilities.  It is generally a good idea to import all of the \modescript library, although if the functions and classes required are known they can be specified individually. The function \verb|LinearPDE| has been imported explicitly for ease of use later in the script. \verb|Rectangle| is going to be our type of model. The module \verb unitsSI  provides support for SI unit definitions with our variables.
287
288  Once our library dependencies have been established, defining the problem specific variables is the next step. In general the number of variables needed will vary between problems. These variables belong to two categories. They are either directly related to the PDE and can be used as inputs into the \esc solver, or they are script variables used to control internal functions and iterations in our problem. For this PDE there are a number of constants which will need values. Firstly, the model upon which we wish to solve our problem needs to be defined. There are many different types of models in \modescript which we will demonstrate in later tutorials but for our iron rod, we will simply use a rectangular model.  Once our library dependencies have been established, defining the problem specific variables is the next step. In general the number of variables needed will vary between problems. These variables belong to two categories. They are either directly related to the PDE and can be used as inputs into the \esc solver, or they are script variables used to control internal functions and iterations in our problem. For this PDE there are a number of constants which will need values. Firstly, the model upon which we wish to solve our problem needs to be defined. There are many different types of models in \modescript which we will demonstrate in later tutorials but for our iron rod, we will simply use a rectangular model.
289
# Line 289  ndy = 1 # mesh steps in y direction - on Line 298  ndy = 1 # mesh steps in y direction - on
298  The material constants and the temperature variables must also be defined. For the iron rod in the model they are defined as:  The material constants and the temperature variables must also be defined. For the iron rod in the model they are defined as:
299  \begin{python}  \begin{python}
300  #PDE related  #PDE related
q=200. * Celsius #Kelvin - our heat source temperature
Tref = 0. * Celsius #Kelvin - starting temp of iron bar
301  rho = 7874. *kg/m**3 #kg/m^{3} density of iron  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
302  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity  cp = 449.*J/(kg*K) # J/Kg.K thermal capacity
303  rhocp = rho*cp  rhocp = rho*cp
304  kappa = 80.*W/m/K #watts/m.Kthermal conductivity  kappa = 80.*W/m/K   # watts/m.Kthermal conductivity
305    qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source
306    T0=100 * Celsius # initial temperature at left end of rod
307    Tref=20 * Celsius # base temperature
308  \end{python}  \end{python}
309  Finally, to control our script we will have to specify our timing controls and where we would like to save the output from the solver. This is simple enough:  Finally, to control our script we will have to specify our timing controls and where we would like to save the output from the solver. This is simple enough:
310  \begin{python}  \begin{python}
311  t=0 #our start time, usually zero  t=0 * day  #our start time, usually zero
312  tend=5.*minute #seconds - time to end simulation  tend=1. * day # - time to end simulation
313  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
314  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
315  #user warning statement  #user warning statement
316  print "Expected Number of time outputs is: ", (tend-t)/h  print "Expected Number of time outputs is: ", (tend-t)/h
317  i=0 #loop counter  i=0 #loop counter
#the folder to put our outputs in, leave blank "" for script path
save_path="data/onedheatdiff001"
318  \end{python}  \end{python}
319  Now that we know our inputs we will build a domain using the \verb Rectangle() function from \verb finley . The four arguments allow us to define our domain \verb rod  as:  Now that we know our inputs we will build a domain using the \verb Rectangle() function from \verb finley . The four arguments allow us to define our domain \verb rod  as:
320  \begin{python}  \begin{python}
321  #generate domain using rectangle  #generate domain using rectangle
322  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
323  \end{python}  \end{python}
324  \verb rod  now describes a domain in the manner of Section \ref{ss:domcon}. As we define our variables, various function spaces will be created to accommodate them. There is an easy way to extract finite points from the domain \verb|rod| using the domain property function \verb|getX()| . This function sets the vertices of each cell as finite points to solve in the solution. If we let \verb|x| be these finite points, then;  \verb rod  now describes a domain in the manner of Section \ref{ss:domcon}. There is an easy way to extract
325    the coordinates of the nodes used to describe the domain \verb|rod| using the
326    domain property function \verb|getX()| . This function sets the vertices of each cell as finite points to solve in the solution. If we let \verb|x| be these finite points, then;
327  \begin{python}  \begin{python}
328  #extract finite points - the solution points  #extract data points - the solution points
329  x=rod.getX()  x=rod.getX()
330  \end{python}  \end{python}
331  The data locations of specific function spaces can be returned in a similar manner by extracting the relevant function space from the domain followed by the \verb .getX()  operator.  The data locations of specific function spaces can be returned in a similar manner by extracting the relevant function space from the domain followed by the \verb|.getX()| method.
332
333  With a domain and all our required variables established, it is now possible to set up our PDE so that it can be solved by \esc. The first step is to define the type of PDE that we are trying to solve in each time step. In this example it is a single linear PDE\footnote{in comparison to a system of PDEs which will be discussed later.}. We also need to state the values of our general form variables.  With a domain and all our required variables established, it is now possible to set up our PDE so that it can be solved by \esc. The first step is to define the type of PDE that we are trying to solve in each time step. In this example it is a single linear PDE\footnote{in comparison to a system of PDEs which will be discussed later.}. We also need to state the values of our general form variables.
334  \begin{python}  \begin{python}
335  mypde=LinearSinglePDE(rod)  mypde=LinearSinglePDE(rod)
336  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)  A=zeros((2,2)))
337  \end{python}  A[0,0]=kappa
338    q=whereZero(x[0])
339    mypde.setValue(A=A, D=rhocp/h, q=q, r=T0)
340    \end{python}
341    The argument \verb|q| has not been discussed yet: In fact the arguments \verb|q| and \verb|r| are used to define
342    Dirichlet boundary condition as discussed in Section~\ref{SEC BOUNDARY COND}. In the \esc
343    PDE from the argument \verb|q| indicates by a positive value for which nodes we want to apply a
344    Dirichlet boundary condition, ie. where we want to prescribe the value of the PDE solution
345    rather then using the PDE. The actually value for the solution to be taken is set by the argument \verb|r|.
346    In our case we want to keep the initial temperature $T0$ on the left face of the rode for all times. Notice,
347    that as set to a constant value \verb|r| is assumed to have the same value
348    at all nodes, however only the value at those nodes marked by a positive value by \verb|q| are actually used.
349
350    In order to set \verb|q| we use
351    \verb|whereZero| function. The function returns the value (positive) one for those data points (=nodes) where the argument is equal to zero and otherwise returns (non-positive) value zero.
352    As \verb|x[0]| given the $x$-coordinates of the nodes for the domain,
353    \verb|whereZero(x[0])| gives the value $1$ for the nodes at the left end of the rod $x=x_0=0$ and
354    zero elsewhere which is exactly what we need.
355
356  In a few special cases it may be possible to decrease the computational time of the solver if our PDE is symmetric. Symmetry of a PDE is defined by;  In a many cases it may be possible to decrease the computational time of the solver if the PDE is symmetric.
357    Symmetry of a PDE is defined by;
358  \label{eqn:symm}  \label{eqn:symm}
359  A\hackscore{jl}=A\hackscore{lj}  A\hackscore{jl}=A\hackscore{lj}
360
361  Symmetry is only dependent on the $A$ coefficient in the general form and the other coefficients $D$ and $d$ as well as the RHS $Y$ and $y$ may take any value. From the above definition we can see that our PDE is symmetric. The \verb LinearPDE  class provides the method \method{checkSymmetry} to check if the given PDE is symmetric. As our PDE is symmetrical we will enable symmetry via;  Symmetry is only dependent on the $A$ coefficient in the general form and the other coefficients $D$ as well as the right hand side $Y$ may take any value. From the above definition we can see that our PDE is symmetric. The \verb LinearPDE  class provides the method \method{checkSymmetry} to check if the given PDE is symmetric. As our PDE is symmetrical we will enable symmetry via;
362  \begin{python}  \begin{python}
363   myPDE.setSymmetryOn()   myPDE.setSymmetryOn()
364  \end{python}  \end{python}
365    Next we need to establish the initial temperature distribution \verb|T|. We want to have this initial
366  We now need to specify our boundary conditions and initial values. The initial values required to solve this PDE are temperatures for each discrete point in our model. We will set our bar to:  value to be \verb|Tref| except at the left end of the rod $x=0$ where we have the temperature \verb|T0|. We use;
\begin{python}
T = Tref
\end{python}
Boundary conditions are a little more difficult. Fortunately the \esc solver will handle our insulated boundary conditions by default with a zero flux operator. However, we will need to apply our heat source $q_{H}$ to the end of the bar at $x=0$ . \esc makes this easy by letting us define areas in our model. The finite points in the model were previously defined as \verb x  and it is possible to set all of points that satisfy $x=0$ to \verb q  via the \verb whereZero()  function. There are a few \verb where  functions available in \esc. They will return a value \verb 1  where they are satisfied and \verb 0  where they are not. In this case our \verb qH  is only applied to the far LHS of our model as required.
367  \begin{python}  \begin{python}
368  # ... set heat source: ....  # ... set initial temperature ....
369  qH=q*whereZero(x[0])  T = T0*whereZero(x[0])+Tref*(1-whereZero(x[0]))
370  \end{python}  \end{python}
371    Finally we will initialize an iteration loop to solve our PDE for all the time steps we specified in the variable section. As the right hand side of the general form is dependent on the previous values for temperature \verb T  across the bar this must be updated in the loop. Our output at each time step is \verb T  the heat distribution and \verb totT  the total heat in the system.
Finally we will initialise an iteration loop to solve our PDE for all the time steps we specified in the variable section. As the RHS of the general form is dependent on the previous values for temperature \verb T  across the bar this must be updated in the loop. Our output at each time step is \verb T  the heat distribution and \verb totT  the total heat in the system.
372  \begin{python}  \begin{python}
373  while t<=tend:  while t < tend:
374      i+=1 #increment the counter      i+=1 #increment the counter
375      t+=h #increment the current time      t+=h #increment the current time
376      mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients      mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients
377      T=mypde.getSolution() #get the PDE solution      T=mypde.getSolution() #get the PDE solution
378      totT = rhocp*T #get the total heat solution in the system      totE = integrate(rhocp*T) #get the total heat (energy) in the system
379  \end{python}  \end{python}
380    The last statement in this script calculates the total energy in the system as volume integral
381    of $\rho \c_p T$ over the rod.
382
383    \subsection{Plotting the Total Energy}
384    \sslist{onedheatdiff001.py}
385
386  \subsection{Plotting the heat solutions}  \esc does not include its own plotting capabilities. However, it is possible to use a variety of free \pyt packages for visualization.
387  Visualisation of the solution can be achieved using \mpl a module contained within \pylab. We start by modifying our solution script from before. Prior to the \verb while  loop we will need to extract our finite solution points to a data object that is compatible with \mpl. First it is necessary to convert \verb x  to a list of tuples. These are then converted to a \numpy array and the $x$ locations extracted via an array slice to the variable \verb plx  .  Two types will be demonstrated in this cookbook; \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and \verb VTK \footnote{\url{http://www.vtk.org/}} visualisation.
388    The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots.
389    For more complex visualisation tasks in particular when it comes to  two and three dimensional problems it is recommended to us more advanced tools for instance  \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}}
390    which bases the \verb|VTK| toolkit. We will discuss the usage of \verb|VTK| based
391    visualization in Chapter~\ref{Sec:2DHD} where will discuss a two dimensional PDE.
392
393    For our simple problem we have two plotting tasks: Firstly we are interested in showing the
394    behavior of the total energy over time and secondly in how the temperature distribution within the rod is
396
397    \begin{figure}
398    \begin{center}
399    \includegraphics[width=4in]{figures/ttrodpyplot150}
400    \caption{Total Energy in Rod over Time (in seconds).}
401    \label{fig:onedheatout1}
402    \end{center}
403    \end{figure}
404
405    The trick is to create a record of the time marks and the corresponding total energies observed.
406    \pyt provides the concept of lists for this. Before
407    the time loop is opened we create empty lists for the time marks \verb|t_list| and the total energies \verb|E_list|.
408    After the new temperature as been calculated by solving the PDE we append the new time marker and total energy
409    to the corresponding list using the \verb|append| method. With these modifications the script looks as follows:
410    \begin{python}
411    t_list=[]
412    E_list=[]
413    # ... start iteration:
414    while t<tend:
415          t+=h
416          mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients
417          T=mypde.getSolution() #get the PDE solution
418          totE=integrate(rhocp*T)
419          t_list.append(t)   # add current time mark to record
420          E_list.append(totE) # add current total energy to record
421    \end{python}
422    To plot $t$ over $totE$ we use the \mpl a module contained within \pylab which needs to be loaded before used;
423    \begin{python}
424    import pylab as pl # plotting package.
425    \end{python}
426    Here we are not using the \verb|from pylab import *| in order to name clashes for function names
427    with \esc.
428
429    The following statements are added to the script after the time loop has been completed;
430    \begin{python}
431    pl.plot(t_list,E_list)
432    pl.title("Total Energy")
433    pl.savefig("totE.png")
434    \end{python}
435    The first statement hands over the time marks and corresponding total energies to the plotter.
436    The second statement is setting the title for the plot. The last statement renders the plot and writes the
437    result into the file \verb|totE.png| which can be displayed by (almost) any image viewer. Your result should look
438    similar to Figure~\ref{fig:onedheatout1}.
439
440    \subsection{Plotting the Temperature Distribution}
441    \sslist{onedheatdiff001b.py}
442    For plotting the spatial distribution of the temperature we need to modify the strategy we have used
443    for the total energy. Instead of producing a final plot at the end we will generate a
444    picture at each time step which can be browsed as slide show or composed to a movie.
445    The first problem we encounter is that if we produce an image in each time step we need
446    to make sure that the images previously generated are not overwritten.
447
448    To develop an incrementing file name we can use the following convention. It is convenient to
449    put all image file showing the same variable - in our case the temperature distribution -
450    into a separate directory. As part of the \verb|os| module\footnote{The \texttt{os} module provides
451    a powerful interface to interact with the operating system, see \url{http://docs.python.org/library/os.html}.} \pyt
452    provides the  \verb|os.path.join| command to build file and
453    directory names in a platform independent way. Assuming that
454    \verb|save_path| is name of directory we want to put the results the command is;
455    \begin{python}
456    import os
457    os.path.join(save_path, "tempT%03d.png"%i )
458    \end{python}
459    where \verb|i| is the time step counter.
460    There are two arguments to the \verb join  command. The \verb save_path  variable is a predefined string pointing to the directory we want to save our data in, for example a single sub-folder called \verb data  would be defined by;
461    \begin{verbatim}
462    save_path = "data"
463    \end{verbatim}
464    while a sub-folder of \verb data  called \verb onedheatdiff001  would be defined by;
465    \begin{verbatim}
466    save_path = os.path.join("data","onedheatdiff001")
467    \end{verbatim}
468    The second argument of \verb join \xspace contains a string which is the filename or subdirectory name. We can use the operator \verb|%| to increment our file names with the value \verb|i| denoting a incrementing counter. The substring \verb %03d  does this by defining the following parameters;
469    \begin{itemize}
470     \item \verb 0  becomes the padding number;
471     \item \verb 3  tells us the amount of padding numbers that are required; and
472     \item \verb d  indicates the end of the \verb %  operator.
473    \end{itemize}
474    To increment the file name a \verb %i  is required directly after the operation the string is involved in. When correctly implemented the output files from this command would be place in the directory defined by \verb save_path  as;
475    \begin{verbatim}
476    rodpyplot.png
477    rodpyplot.png
478    rodpyplot.png
479    ...
480    \end{verbatim}
481    and so on.
482
483    A sub-folder check/constructor is available in \esc. The command;
484    \begin{verbatim}
485    mkDir(save_path)
486    \end{verbatim}
487    will check for the existence of \verb save_path  and if missing, make the required directories.
488
489    We start by modifying our solution script from before.
490    Prior to the \verb|while| loop we will need to extract our finite solution points to a data object that is compatible with \mpl. First we create the node coordinates of the data points used to represent
491    the temperature as a \pyt list of tuples. As a solution of a PDE
492    the temperature has the \verb|Solution(rod)| function space attribute. We use
493    the \verb|getX()| method to get the coordinates of the data points as an \esc object
494    which is then converted to a \numpy array. The $x$ component is then extracted via an array slice to the variable \verb|plx|;
495  \begin{python}  \begin{python}
496    import numpy as np # array package.
497  #convert solution points for plotting  #convert solution points for plotting
498  plx = x.toListOfTuples()  plx = Solution(rod).getX().toListOfTuples()
499  plx = np.array(plx) #convert to tuple to numpy array  plx = np.array(plx) # convert to tuple to numpy array
500  plx = plx[:,0] #extract x locations  plx = plx[:,0] # extract x locations
\end{python}
As there are two solution outputs, we will generate two plots and save each to a file for every time step in the solution. The following is appended to the end of the \verb while  loop and creates two figures. The first figure is for the temperature distribution, and the second the total temperature in the bar. Both cases are similar with a few minor changes for scale and labelling. We start by converting the solution to a tuple and then plotting this against our \textit{x coordinates} \verb plx  from before. The axis is then standardised and a title applied. Finally, the figure is saved to a *.png file and cleared for the following iteration.
\begin{python}
#establish figure 1 for temperature vs x plots
tempT = T.toListOfTuples(scalarastuple=False)
pl.figure(1) #current figure
pl.plot(plx,tempT) #plot solution
#define axis extents and title
pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
pl.title("Temperature accross Rod")
#save figure to file
pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i)
pl.clf() #clear figure

#establish figure 2 for total temperature vs x plots and repeat
tottempT = totT.toListOfTuples(scalarastuple=False)
pl.figure(2)
pl.plot(plx,tottempT)
pl.axis([0,1.0,9.657E08,12000+9.657E08])
pl.title("Total temperature across Rod")
pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i)
pl.clf()
501  \end{python}  \end{python}
502
503  \begin{figure}  \begin{figure}
504  \begin{center}  \begin{center}
505  \includegraphics[width=4in]{figures/ttrodpyplot150}  \includegraphics[width=4in]{figures/rodpyplot001}
506  \caption{Total temperature ($T$) distribution in rod at $t=150$}  \includegraphics[width=4in]{figures/rodpyplot050}
507    \includegraphics[width=4in]{figures/rodpyplot200}
508    \caption{Temperature ($T$) distribution in rod at time steps $1$, $50$ and $200$.}
509  \label{fig:onedheatout}  \label{fig:onedheatout}
510  \end{center}  \end{center}
511  \end{figure}  \end{figure}
512
513  \subsubsection{Parallel scripts (MPI)}  For each time step we will generate a plot of the temperature distribution and save each to a file. We use the same
514  In some of the example files for this cookbook the plotting commands are a little different.  techniques provided by \mpl as we have used to plot the total energy over time.
515  For example,  The following is appended to the end of the \verb while  loop and creates one figure of the temperature distribution. We start by converting the solution to a tuple and then plotting this against our \textit{x coordinates} \verb plx  we have generated before. We add a title to the diagram before it is rendered into a file.
516  \begin{python}  Finally, the figure is saved to a \verb|*.png| file and cleared for the following iteration.
517      if getMPIRankWorld() == 0:  \begin{python}
518          pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i)  # ... start iteration:
519      pl.clf()      while t<tend:
520            ....
521        T=mypde.getSolution() #get the PDE solution
522            tempT = T.toListOfTuples() # convert to a tuple
523            pl.plot(plx,tempT) # plot solution
524        # set scale (Temperature should be between Tref and T0)
525            pl.axis([0,mx,Tref*.9,T0*1.1])
527        pl.title("Temperature across rod at time %e minutes"%(t/minutes))
528        #save figure to file
529        pl.savefig(os.path.join(save_path,"tempT","rodpyplot%03d.png") %i)
530  \end{python}  \end{python}
531    Some results are shown in Figure~\ref{fig:onedheatout}.
The additional \verb if  statement is not necessary for normal desktop use.
It becomes important for scripts run on parallel computers.
Its purpose is to ensure that only one copy of the file is written.
For more details on writing scripts for parallel computing please consult the \emph{user's guide}.
532
533  \subsection{Make a video}  \subsection{Make a video}
534  Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. \verb mencoder  is linux only however, and other platform users will need to use an alternative video encoder.  Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. \verb mencoder  is Linux only however, and other platform users will need to use an alternative video encoder.
535  \begin{python}  \begin{python}
536  # compile the *.png files to create two *.avi videos that show T change  # compile the *.png files to create a *.avi videos that show T change
537  # with time. This operation uses linux mencoder. For other operating  # with time. This operation uses Linux mencoder. For other operating
538  # systems it is possible to use your favourite video compiler to  # systems it is possible to use your favorite video compiler to
539  # convert image files to videos.  # convert image files to videos.
540
541  os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\  os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
542  w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \             w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
543  onedheatdiff001tempT.avi")             onedheatdiff001tempT.avi")

os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
onedheatdiff001totT.avi")
544  \end{python}  \end{python}
545

Legend:
 Removed from v.2877 changed lines Added in v.2878