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

revision 2975 by ahallam, Thu Mar 4 05:44:12 2010 UTC revision 2979 by ahallam, Tue Mar 9 02:54:32 2010 UTC
# Line 13  Line 13
13
14  \begin{figure}[h!]  \begin{figure}[h!]
15  \centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}}  \centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}}
16  \caption{Example 1: Temperature differential along a single interface between two granite blocks.}  \caption{Example 1: Temperature differential along a single interface between
17    two granite blocks.}
18  \label{fig:onedgbmodel}  \label{fig:onedgbmodel}
19  \end{figure}  \end{figure}
20
21  \section{Example 1: One Dimensional Heat Diffusion in Granite}  \section{Example 1: One Dimensional Heat Diffusion in Granite}
22  \label{Sec:1DHDv00}  \label{Sec:1DHDv00}
23
24  The first model consists of two blocks of isotropic material, for instance granite, sitting next to each other.  The first model consists of two blocks of isotropic material, for instance
25  Initial temperature in \textit{Block 1} is \verb|T1| and in \textit{Block 2} is \verb|T2|.  granite, sitting next to each other.
26    Initial temperature in \textit{Block 1} is \verb|T1| and in \textit{Block 2} is
27    \verb|T2|.
28  We assume that the system is insulated.  We assume that the system is insulated.
29  What would happen to the temperature distribution in each block over time?  What would happen to the temperature distribution in each block over time?
30  Intuition tells us that heat will be transported from the hotter block to the cooler one until both  Intuition tells us that heat will be transported from the hotter block to the
31    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 dimensional heat diffusion equation\footnote{A detailed discussion on how the heat diffusion equation is derived can be found at \url{http://online.redwoods.edu/instruct/darnold/DEProj/sp02/AbeRichards/paper.pdf}};  We can model the heat distribution of this problem over time using one
36    dimensional heat diffusion equation\footnote{A detailed discussion on how the
37    heat diffusion equation is derived can be found at
38    \url{
39    http://online.redwoods.edu/instruct/darnold/DEProj/sp02/AbeRichards/paper.pdf}};
40  which is defined as:  which is defined as:
41  \begin{equation}  \begin{equation}
42  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2}
43    T}{\partial x^{2}} = q\hackscore H
44  \label{eqn:hd}  \label{eqn:hd}
45  \end{equation}  \end{equation}
46  where $\rho$ is the material density, $c\hackscore p$ is the specific heat and $\kappa$ is the thermal  where $\rho$ is the material density, $c\hackscore p$ is the specific heat and
47  conductivity\footnote{A list of some common thermal conductivities is available from Wikipedia \url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}. Here we assume that these material  $\kappa$ is the thermal
48    conductivity\footnote{A list of some common thermal conductivities is available
49    from Wikipedia
50    \url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}. Here we
51    assume that these material
52  parameters are \textbf{constant}.  parameters are \textbf{constant}.
53  The heat source is defined by the right hand side of \refEq{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = q\hackscore{0}e^{-\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our signed distance from the the block-block interface $x$.  The heat source is defined by the right hand side of \refEq{eqn:hd} as
54    $q\hackscore{H}$; this can take the form of a constant or a function of time and
55    space. For example $q\hackscore{H} = q\hackscore{0}e^{-\gamma t}$ where we have
56    the output of our heat source decaying with time. There are also two partial
57    derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the
58    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
60    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$.
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 solution to our problem. However, it is not always possible or practical to solve the problem this way. Alternatively, computers can be used to find the solution. To do this, a numerical approach is required to discretise  It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact
65  the PDE \refEq{eqn:hd} across time and space, this reduces the problem to a finite number of equations for a finite number of spatial points and time steps. These parameters together define the model. While discretisation introduces approximations and a degree of error, a sufficiently sampled model is generally accurate enough to satisfy the accuracy requirements for the final solution.  solution to our problem. However, it is not always practical to solve the
66    problem this way. Alternatively, computers can be used to find the solution. To
67  Firstly, we discretise the PDE \refEq{eqn:hd} in time. This leaves us with a steady linear PDE which involves spatial derivatives only and needs to be solved in each time step to progress in time. \esc can help us here.  do this, a numerical approach is required to discretise
68    the PDE \refEq{eqn:hd} across time and space, this reduces the problem to a
69  For time discretization we use the Backwards Euler approximation scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is based on the  finite number of equations for a finite number of spatial points and time steps.
70    These parameters together define the model. While discretisation introduces
71    approximations and a degree of error, a sufficiently sampled model is generally
72    accurate enough to satisfy the accuracy requirements for the final solution.
73
74    Firstly, we discretise the PDE \refEq{eqn:hd} in time. This leaves us with a
75    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.
77
78    For time discretization we use the Backwards Euler approximation
79    scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is
80    based on the
81  approximation  approximation
82  \begin{equation}  \begin{equation}
83  \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}
# Line 57  where $h$ is the time step size. This ca Line 89  where $h$ is the time step size. This ca
89  \frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)} - T^{(n-1)}}{h}  \frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)} - T^{(n-1)}}{h}
90  \label{eqn:Tbeuler}  \label{eqn:Tbeuler}
91  \end{equation}  \end{equation}
92  where the upper index $n$ denotes the n\textsuperscript{th} time step. So one has  where the upper index $n$ denotes the n\textsuperscript{th} time step. So one
93    has
94  \begin{equation}  \begin{equation}
95  \begin{array}{rcl}  \begin{array}{rcl}
96  t^{(n)} & = & t^{(n-1)}+h \\  t^{(n)} & = & t^{(n-1)}+h \\
# Line 67  T^{(n)} & = & T(t^{(n-1)}) \\ Line 100  T^{(n)} & = & T(t^{(n-1)}) \\
100  \end{equation}  \end{equation}
101  Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get;  Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get;
102  \begin{equation}  \begin{equation}
103  \frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2} T^{(n)}}{\partial x^{2}} = q\hackscore H  \frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2}
104    T^{(n)}}{\partial x^{2}} = q\hackscore H
105  \label{eqn:hddisc}  \label{eqn:hddisc}
106  \end{equation}  \end{equation}
107  Notice that we evaluate the spatial derivative term at current time $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively, one can evaluate the spatial derivative term at the previous time $t^{(n-1)}$. This  Notice that we evaluate the spatial derivative term at the current time
108  approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which  $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively,
109  are not discussed here. However, the \textbf{forward Euler} scheme has a major disadvantage. Depending on the  one can evaluate the spatial derivative term at the previous time $t^{(n-1)}$.
110  material parameter as well as the discretization of the spatial derivative term, the time step size $h$ needs to be chosen sufficiently small to achieve a stable temperature when progressing in time. Stabiliy is achieved if the temperature does not grow beyond its initial bounds and become non-physical.  This
111  The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of  approach is called the \textbf{forward Euler} scheme. This scheme can provide
112  physically correct problem set-up the temperature approximation remains physical for all time steps.  some computational advantages, which
113  The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}  are not discussed here. However, the \textbf{forward Euler} scheme has a major
114  is sufficiently small, thus a good approximation of the true temperature is computed. It is  disadvantage. Namely, depending on the
115  therefore crucial that the user remains sceptical about their results and for instance compares  material parameters as well as the domain discretization of the spatial
116  the results for different time and spatial step sizes for correlation.  derivative term, the time step size $h$ needs to be chosen sufficiently small to
117    achieve a stable temperature when progressing in time. Stabiliy is achieved if
118    the temperature does not grow beyond its initial bounds and become
119    non-physical.
120    The backward Euler scheme, which we use here, is unconditionally stable meaning
121    that under the assumption of a
122    physically correct problem set-up the temperature approximation remains physical
123    for all time steps.
124    The user needs to keep in mind that the discretization error introduced by
125    \refEq{eqn:beuler}
126    is sufficiently small, thus a good approximation of the true temperature is
127    computed. It is
128    therefore very important that any results are viewed with caution. For example,
129    one may compare the results for different time and spatial step sizes.
130
131  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
132  differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem  differential equation \refEq{eqn:hddisc} which only includes spatial
133    derivatives. To solve this problem
134  we want to use \esc.  we want to use \esc.
135
136  In \esc any given PDE can be described by the general form. For the purpose of this introduction we illustrate a simpler version of the general form for full linear PDEs which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{The form in the \esc users guide which uses the Einstein convention is written as  In \esc any given PDE can be described by the general form. For the purpose of
137    this introduction we illustrate a simpler version of the general form for full
138    linear PDEs which is available in the \esc user's guide. A simplified form that
139    suits our heat diffusion problem\footnote{The form in the \esc users guide which
140    uses the Einstein convention is written as
141  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}
142  is described by;  is described by;
143  \begin{equation}\label{eqn:commonform nabla}  \begin{equation}\label{eqn:commonform nabla}
144  -\nabla\cdot(A\cdot\nabla u) + Du = f  -\nabla\cdot(A\cdot\nabla u) + Du = f
145  \end{equation}  \end{equation}
146  where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del operator} represents  where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The
147  the spatial derivative of its subject - in this case $u$. Lets assume for a moment that we deal with a one-dimensional problem then ;  symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del
148    operator} represents
149    the spatial derivative of its subject - in this case $u$. Lets assume for a
150    moment that we deal with a one-dimensional problem then ;
151  \begin{equation}  \begin{equation}
152  \nabla = \frac{\partial}{\partial x}  \nabla = \frac{\partial}{\partial x}
153  \end{equation}  \end{equation}
# Line 100  and we can write \refEq{eqn:commonform n Line 155  and we can write \refEq{eqn:commonform n
155  \begin{equation}\label{eqn:commonform}  \begin{equation}\label{eqn:commonform}
156  -A\frac{\partial^{2}u}{\partial x^{2}} + Du = f  -A\frac{\partial^{2}u}{\partial x^{2}} + Du = f
157  \end{equation}  \end{equation}
158  if $A$ is constant. To match this simplified general form to our problem \refEq{eqn:hddisc}  if $A$ is constant. To match this simplified general form to our problem
159    \refEq{eqn:hddisc}
160  we rearrange \refEq{eqn:hddisc};  we rearrange \refEq{eqn:hddisc};
161  \begin{equation}  \begin{equation}
162  \frac{\rho c\hackscore p}{h} T^{(n)} - \kappa \frac{\partial^2 T^{(n)}}{\partial x^2} = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}  \frac{\rho c\hackscore p}{h} T^{(n)} - \kappa \frac{\partial^2 T^{(n)}}{\partial
163    x^2} = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}
164  \label{eqn:hdgenf}  \label{eqn:hdgenf}
165  \end{equation}  \end{equation}
166  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is
167  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed to be constant.  required for \esc to solve our PDE. This can be done by generating a solution
168  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see that;  for successive increments in the time nodes $t^{(n)}$ where
169    $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed
170    to be constant.
171    In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$.
172    Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see
173    that;
174  \begin{equation}\label{ESCRIPT SET}  \begin{equation}\label{ESCRIPT SET}
175  u=T^{(n)};  u=T^{(n)};
176  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho
177    c\hackscore p}{h} T^{(n-1)}
178  \end{equation}  \end{equation}
179
180  \subsection{Boundary Conditions}  \subsection{Boundary Conditions}
181  \label{SEC BOUNDARY COND}  \label{SEC BOUNDARY COND}
182  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
183  A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown solution (in our example the temperature) on parts of the boundary or on the entire boundary of the region of interest.  boundary conditions of our model. Typically there are two main types of boundary
184  We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}.  conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary
186    Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}},
187    respectively.
188    A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to
189    prescribe a known value to the unknown solution (in our example the temperature)
190    on parts of the boundary or on the entire boundary of the region of interest.
191    We discuss Dirichlet boundary condition in our second example presented in
192    Section~\ref{Sec:1DHDv0}.
193
194  However, for this example we have made the model assumption that the system is insulated, so we need  However, for this example we have made the model assumption that the system is
195    insulated, so we need
196  to add an appropriate boundary condition to prevent  to add an appropriate boundary condition to prevent
197  any loss or inflow of energy at the boundary of our domain. Mathematically this is expressed by prescribing  any loss or inflow of energy at the boundary of our domain. Mathematically this
198  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified one dimensional model this is expressed  is expressed by prescribing
199    the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified
200    one dimensional model this is expressed
201  in the form;  in the form;
202  \begin{equation}  \begin{equation}
203  \kappa \frac{\partial T}{\partial x}  = 0  \kappa \frac{\partial T}{\partial x}  = 0
# Line 132  or in a more general case as Line 206  or in a more general case as
206  \begin{equation}\label{NEUMAN 1}  \begin{equation}\label{NEUMAN 1}
207  \kappa \nabla T \cdot n  = 0  \kappa \nabla T \cdot n  = 0
208  \end{equation}  \end{equation}
209  where $n$  is the outer normal field \index{outer normal field} at the surface of the domain.  where $n$  is the outer normal field \index{outer normal field} at the surface
210  The $\cdot$ (dot) refers to the  dot product of the vectors $\nabla T$ and $n$. In fact, the term $\nabla T \cdot n$ is the normal derivative of  of the domain.
211  the temperature $T$. Other notations used here are\footnote{The \esc notation for the normal  The $\cdot$ (dot) refers to the  dot product of the vectors $\nabla T$ and $n$.
212    In fact, the term $\nabla T \cdot n$ is the normal derivative of
213    the temperature $T$. Other notations used here are\footnote{The \esc notation
214    for the normal
215  derivative is $T\hackscore{,i} n\hackscore i$.};  derivative is $T\hackscore{,i} n\hackscore i$.};
216  \begin{equation}  \begin{equation}
217  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .  \nabla T \cdot n  = \frac{\partial T}{\partial n} \; .
218  \end{equation}  \end{equation}
219  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary condition} for the PDE.  A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary
220    condition} for the PDE.
221
222  The PDE \refEq{eqn:hdgenf}  The PDE \refEq{eqn:hdgenf}
223  and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with the Dirichlet boundary conditions)  define a \textbf{boundary value problem}.  and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with
224  It is the nature of a boundary value problem to allow making statements about the solution in the  the Dirichlet boundary conditions)  define a \textbf{boundary value problem}.
225  interior of the domain from information known on the boundary only. In most cases  It is the nature of a boundary value problem to allow making statements about
226  we use the term partial differential equation but in fact it is a boundary value problem.  the solution in the
227  It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that  interior of the domain from information known on the boundary only. In most
228  at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set.  cases
229    we use the term partial differential equation but in fact it is a boundary value
230  Conveniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate.  problem.
231  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  It is important to keep in mind that boundary conditions need to be complete and
232    consistent in the sense that
233    at any point on the boundary either a Dirichlet or a Neuman boundary condition
234    must be set.
235
236    Conveniently, \esc makes default assumption on the boundary conditions which the
237    user may modify where appropriate.
238    For a problem of the form in~\refEq{eqn:commonform nabla} the default
239    condition\footnote{In the form of the \esc users guide which is using the
240    Einstein convention is written as
241  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;  $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;
242  \begin{equation}\label{NEUMAN 2}  \begin{equation}\label{NEUMAN 2}
243  -n\cdot A \cdot\nabla u = 0  -n\cdot A \cdot\nabla u = 0
244  \end{equation}  \end{equation}
245  which is used everywhere on the boundary. Again $n$ denotes the outer normal field.  which is used everywhere on the boundary. Again $n$ denotes the outer normal
246  Notice that the coefficient $A$ is the same as in the \esc PDE~\ref{eqn:commonform nabla}.  field.
247  With the settings for the coefficients we have already identified in \refEq{ESCRIPT SET} this  Notice that the coefficient $A$ is the same as in the \esc
248    PDE~\ref{eqn:commonform nabla}.
249    With the settings for the coefficients we have already identified in
250    \refEq{ESCRIPT SET} this
251  condition translates into  condition translates into
252  \begin{equation}\label{NEUMAN 2b}  \begin{equation}\label{NEUMAN 2b}
253  \kappa \frac{\partial T}{\partial x} = 0  \kappa \frac{\partial T}{\partial x} = 0
254  \end{equation}  \end{equation}
255  for the boundary of the domain. This is identical to the Neuman boundary condition we want to set. \esc will take care of this condition for us. We discuss the Dirichlet boundary condition later.  for the boundary of the domain. This is identical to the Neuman boundary
256    condition we want to set. \esc will take care of this condition for us. We
257    discuss the Dirichlet boundary condition later.
258
259  \subsection{Outline of the Implementation}  \subsection{Outline of the Implementation}
260  \label{sec:outline}  \label{sec:outline}
261  To solve the heat diffusion equation (\refEq{eqn:hd}) we 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 (discussed in \refSec{sec:key}) has four major steps. Firstly, we need to define the domain where we want to  To solve the heat diffusion equation (\refEq{eqn:hd}) we write a simple \pyt
262  calculate the temperature. For our problem this is the joint blocks of granite which has a rectangular shape. Secondly, we need to define the PDE to solve in each time step to get the updated temperature. Thirdly, we need to define 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. The work flow described in \reffig{fig:wf}.  script. At this point we assume that you have some basic understanding of the
263    \pyt programming language. If not, there are some pointers and links available
264    in Section \ref{sec:escpybas}. The script (discussed in \refSec{sec:key}) has
265    four major steps. Firstly, we need to define the domain where we want to
266    calculate the temperature. For our problem this is the joint blocks of granite
267    which has a rectangular shape. Secondly, we need to define the PDE to solve in
268    each time step to get the updated temperature. Thirdly, we need to define the
269    coefficients of the PDE and finally we need to solve the PDE. The last two steps
270    need to be repeated until the final time marker has been reached. The work flow
271    described in \reffig{fig:wf}.
272  % \begin{enumerate}  % \begin{enumerate}
273  %  \item create domain  %  \item create domain
274  %  \item create PDE  %  \item create PDE
# Line 187  calculate the temperature. For our probl Line 288  calculate the temperature. For our probl
288     \label{fig:wf}     \label{fig:wf}
289  \end{figure}  \end{figure}
290
291  In the terminology of \pyt, the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it is defined by its usage and features  In the terminology of \pyt, the domain and PDE are represented by
292  rather than its actual representation. So we will create a domain object to describe the geometry of the two  \textbf{objects}. The nice feature of an object is that it is defined by its
293  granite blocks. Then we define PDEs and spatially distributed values such as the temperature  usage and features
294  on this domain. Similarly, to define a PDE object we use the fact that one needs only to define the coefficients of the PDE and solve the PDE. The PDE object has advanced features, but these are not required in simple cases.  rather than its actual representation. So we will create a domain object to
295    describe the geometry of the two
296    granite blocks. Then we define PDEs and spatially distributed values such as the
297    temperature
298    on this domain. Similarly, to define a PDE object we use the fact that one needs
299    only to define the coefficients of the PDE and solve the PDE. The PDE object has
300    advanced features, but these are not required in simple cases.
301
302
303  \begin{figure}[t]  \begin{figure}[t]
# Line 202  on this domain. Similarly, to define a P Line 309  on this domain. Similarly, to define a P
309
310  \subsection{The Domain Constructor in \esc}  \subsection{The Domain Constructor in \esc}
311  \label{ss:domcon}  \label{ss:domcon}
312  While it is not strictly relevant or necessary, it can be helpful to have a better understanding of how values are spatially distributed and how PDE coefficients are interpreted in \esc. The example in this case would be temperature.  Whilst it is not strictly relevant or necessary, a better understanding of
313     how values are spatially distributed (\textit{e.g.} Temperature) and how PDE
314  There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is  coefficients are interpreted in \esc can be helpful.
315  a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of
316  elements or cells to be used along the length and height, see \reffig{fig:fs}. Any spatially distributed value  There are various ways to construct domain objects. The simplest form is a
317  and the PDE is represented in discrete form using this element representation\footnote{We use the finite element method (FEM), see \url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}. Therefore we will have access to an approximation of the true PDE solution only.  rectangular shaped region with a length and height. There is
318  The quality of the approximation depends - besides other factors- mainly on the number of elements being used. In fact, the  a ready to use function for this named \verb rectangle(). Besides the spatial
319  approximation becomes better when more elements are used. However, computational costs and time grow with the number of  dimensions this function requires to specify the number of
320  elements being used. It is therefore important that you find the right balance between the demand in accuracy and acceptable resource usage.  elements or cells to be used along the length and height, see \reffig{fig:fs}.
321    Any spatially distributed value
322  In general, one can think about a domain object as a composition of nodes and elements.  and the PDE is represented in discrete form using this element
323  As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to describe its vertices.  representation\footnote{We use the finite element method (FEM), see
324    \url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}.
325    Therefore we will have access to an approximation of the true PDE solution
326    only.
327    The quality of the approximation depends - besides other factors- mainly on the
328    number of elements being used. In fact, the
329    approximation becomes better when more elements are used. However, computational
330    cost grows with the number of
331    elements being used. It is therefore important that you find the right balance
332    between the demand in accuracy and acceptable resource usage.
333
334    In general, one can think about a domain object as a composition of nodes and
335    elements.
336    As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to
337    describe its vertices.
338  To represent spatial distributed values the user can use  To represent spatial distributed values the user can use
339  the values at the nodes, at the elements in the interior of the domain or at the elements located at the surface of the domain.  the values at the nodes, at the elements in the interior of the domain or at the
340  The different approach used to represent values is called \textbf{function space} and is attached to all objects  elements located at the surface of the domain.
341  in \esc representing a spatial distributed value such as the solution of a PDE. The three  The different approach used to represent values is called \textbf{function
342    space} and is attached to all objects
343    in \esc representing a spatial distributed value such as the solution of a PDE.
344    The three
345  function spaces we use at the moment are;  function spaces we use at the moment are;
346  \begin{enumerate}  \begin{enumerate}
347  \item the nodes, called by \verb|ContinuousFunction(domain)| ;  \item the nodes, called by \verb|ContinuousFunction(domain)| ;
348  \item the elements/cells, called by \verb|Function(domain)| ; and  \item the elements/cells, called by \verb|Function(domain)| ; and
349  \item the boundary, called by \verb|FunctionOnBoundary(domain)| .  \item the boundary, called by \verb|FunctionOnBoundary(domain)| .
350  \end{enumerate}  \end{enumerate}
351  A function space object such as \verb|ContinuousFunction(domain)| has the method \verb|getX| attached to it. This method returns the  A function space object such as \verb|ContinuousFunction(domain)| has the method
352  location of the so-called \textbf{sample points} used to represent values of the particular function space. So the  \verb|getX| attached to it. This method returns the
353  call \verb|ContinuousFunction(domain).getX()| will return the coordinates of the nodes used to describe the domain while  location of the so-called \textbf{sample points} used to represent values of the
354  the  \verb|Function(domain).getX()| returns the coordinates of numerical integration points within elements, see  particular function space. So the
355    call \verb|ContinuousFunction(domain).getX()| will return the coordinates of the
356    nodes used to describe the domain while
357    the  \verb|Function(domain).getX()| returns the coordinates of numerical
358    integration points within elements, see
359  \reffig{fig:fs}.  \reffig{fig:fs}.
360
361  This distinction between different representations of spatially distributed values  This distinction between different representations of spatially distributed
362  is important in order to be able to vary the degrees of smoothness in a PDE problem.  values
363  The coefficients of a PDE do not need to be continuous, thus this qualifies as a \verb|Function()| type.  is important in order to be able to vary the degrees of smoothness in a PDE
364  On the other hand a temperature distribution must be continuous and needs to be represented with a \verb|ContinuousFunction()| function space.  problem.
365  An influx may only be defined at the boundary and is therefore a \verb FunctionOnBoundary()  object.    The coefficients of a PDE do not need to be continuous, thus this qualifies as a
366  \esc allows certain transformations of the function spaces. A \verb ContinuousFunction()  can be transformed into a \verb|FunctionOnBoundary()|  \verb|Function()| type.
367  or \verb|Function()|. On the other hand there is not enough information in a \verb FunctionOnBoundary()  to transform it to a \verb ContinuousFunction()  .  On the other hand a temperature distribution must be continuous and needs to be
368  These transformations, which are called \textbf{interpolation} are invoked automatically by \esc if needed.  represented with a \verb|ContinuousFunction()| function space.
369    An influx may only be defined at the boundary and is therefore a \verb
370    FunctionOnBoundary()  object.
371    \esc allows certain transformations of the function spaces. A \verb
372    ContinuousFunction()  can be transformed into a \verb|FunctionOnBoundary()|
373    or \verb|Function()|. On the other hand there is not enough information in a
374    \verb FunctionOnBoundary()  to transform it to a \verb ContinuousFunction()  .
375    These transformations, which are called \textbf{interpolation} are invoked
376    automatically by \esc if needed.
377
378  Later in this introduction we discuss how  Later in this introduction we discuss how
379  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
380  thermal conductivities $\kappa$. A very powerful technique to define these types of PDE  represented by different material coefficients such the
381  coefficients is tagging. Blocks of materials and boundaries can be named and values can be defined on subregions based on their names.  thermal conductivities $\kappa$. A very powerful technique to define these types
382  This is a method for simplifying PDE coefficient and flux definitions. It makes scripting much easier and we will discuss this technique in Section~\ref{STEADY-STATE HEAT REFRACTION}.  of PDE
383    coefficients is tagging. Blocks of materials and boundaries can be named and
384    values can be defined on subregions based on their names.
385    This is a method for simplifying PDE coefficient and flux definitions. It makes
386    scripting much easier and we will discuss this technique in
388
389
390  \subsection{A Clarification for the 1D Case}  \subsection{A Clarification for the 1D Case}
391  \label{SEC: 1D CLARIFICATION}  \label{SEC: 1D CLARIFICATION}
392  It is necessary for clarification that we revisit our general PDE from \refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently designed to solve problems that are greater than one dimension and so \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. In the case of two spatial dimensions the \textit{Nabla operator} has in fact two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$. Assuming the coefficient $A$ is constant, the \refEq{eqn:commonform nabla} takes the following form;  It is necessary for clarification that we revisit our general PDE from
393    \refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently
394    designed to solve problems that are greater than one dimension and so
395    \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem.
396    In the case of two spatial dimensions the \textit{Nabla operator} has in fact
397    two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial 398 y})$. Assuming the coefficient $A$ is constant, the \refEq{eqn:commonform nabla}
399    takes the following form;
400  \begin{equation}\label{eqn:commonform2D}  \begin{equation}\label{eqn:commonform2D}
401  -A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}}  -A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}}
402  -A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y}  -A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y}
# Line 269  A\hackscore{00}=A; A\hackscore{01}=A\hac Line 417  A\hackscore{00}=A; A\hackscore{01}=A\hac
417  \subsection{Developing a PDE Solution Script}  \subsection{Developing a PDE Solution Script}
418  \label{sec:key}  \label{sec:key}
419  \sslist{example01a.py}  \sslist{example01a.py}
420  We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules.  We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl
421  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.}  modules.
422    By developing a script for \esc, the heat diffusion equation can be solved at
423    successive time steps for a predefined period using our general form
424    \refEq{eqn:hdgenf}. Firstly it is necessary to import all the
425    libraries\footnote{The libraries contain predefined scripts that are required to
426    solve certain problems, these can be simple like sine and cosine functions or
427    more complicated like those from our \esc library.}
428  that we will require.  that we will require.
429  \begin{python}  \begin{python}
430  from esys.escript import *  from esys.escript import *
# Line 282  from esys.finley import Rectangle Line 436  from esys.finley import Rectangle
436  # match up in the equations under SI.  # match up in the equations under SI.
437  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
438  \end{python}  \end{python}
439  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 domain. The module \verb unitsSI  provides support for SI unit definitions with our variables.  It is generally a good idea to import all of the \modescript library, although
440    if the functions and classes required are known they can be specified
441  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 need values. Firstly, the domain upon which we wish to solve our problem needs to be defined. There are different types of domains in \modescript which we demonstrate in later tutorials but for our granite blocks, we simply use a rectangular domain.  individually. The function \verb|LinearPDE| has been imported explicitly for
442    ease of use later in the script. \verb|Rectangle| is going to be our type of
443  Using a rectangular domain simplifies our granite blocks (which would in reality be a \textit{3D} object) into a single dimension. The granite blocks will have a lengthways cross section that looks like a rectangle.  As a result we do not need to model the volume of the block due to symmetry. There are four arguments we must consider when we decide to create a rectangular domain, the domain \textit{length}, \textit{width} and \textit{step size} in each direction. When defining the size of our problem it will help us determine appropriate values for our model arguments. If we make our dimensions large but our step sizes very small we increase the accuracy of our solution. Unfortunately we also increase the number of calculations that must be solved per time step. This means more computational time is required to produce a solution. In this \textit{1D} problem, the bar is defined as being 1 metre long. An appropriate step size \verb|ndx| would be 1 to 10\% of the length. Our \verb|ndy| need only be 1, this is because our problem stipulates no partial derivatives in the $y$ direction. Thus the temperature does not vary with $y$. Hence, the model parameters can be defined as follows; note we have used the \verb unitsSI  convention to make sure all our input units are converted to SI.  domain. The module \verb unitsSI  provides support for SI unit definitions with
444    our variables.
445
446    Once our library dependencies have been established, defining the problem
447    specific variables is the next step. In general the number of variables needed
448    will vary between problems. These variables belong to two categories. They are
449    either directly related to the PDE and can be used as inputs into the \esc
450    solver, or they are script variables used to control internal functions and
451    iterations in our problem. For this PDE there are a number of constants which
452    need values. Firstly, the domain upon which we wish to solve our problem needs
453    to be defined. There are different types of domains in \modescript which we
454    demonstrate in later tutorials but for our granite blocks, we simply use a
455    rectangular domain.
456
457    Using a rectangular domain simplifies our granite blocks (which would in reality
458    be a \textit{3D} object) into a single dimension. The granite blocks will have a
459    lengthways cross section that looks like a rectangle.  As a result we do not
460    need to model the volume of the block due to symmetry. There are four arguments
461    we must consider when we decide to create a rectangular domain, the domain
462    \textit{length}, \textit{width} and \textit{step size} in each direction. When
463    defining the size of our problem it will help us determine appropriate values
464    for our model arguments. If we make our dimensions large but our step sizes very
465    small we increase the accuracy of our solution. Unfortunately we also increase
466    the number of calculations that must be solved per time step. This means more
467    computational time is required to produce a solution. In this \textit{1D}
468    problem, the bar is defined as being 1 metre long. An appropriate step size
469    \verb|ndx| would be 1 to 10\% of the length. Our \verb|ndy| need only be 1, this
470    is because our problem stipulates no partial derivatives in the $y$ direction.
471    Thus the temperature does not vary with $y$. Hence, the model parameters can be
472    defined as follows; note we have used the \verb unitsSI  convention to make sure
473    all our input units are converted to SI.
474  \begin{python}  \begin{python}
475  mx = 500.*m #meters - model length  mx = 500.*m #meters - model length
476  my = 100.*m #meters - model width  my = 100.*m #meters - model width
# Line 294  ndx = 50 # mesh steps in x direction Line 478  ndx = 50 # mesh steps in x direction
478  ndy = 1 # mesh steps in y direction  ndy = 1 # mesh steps in y direction
479  boundloc = mx/2 # location of boundary between the two blocks  boundloc = mx/2 # location of boundary between the two blocks
480  \end{python}  \end{python}
481  The material constants and the temperature variables must also be defined. For the granite in the model they are defined as:  The material constants and the temperature variables must also be defined. For
482    the granite in the model they are defined as:
483  \begin{python}  \begin{python}
484  #PDE related  #PDE related
485  rho = 2750. *kg/m**3 #kg/m^{3} density of iron  rho = 2750. *kg/m**3 #kg/m^{3} density of iron
# Line 305  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h Line 490  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h
490  T1=20 * Celsius # initial temperature at Block 1  T1=20 * Celsius # initial temperature at Block 1
491  T2=2273. * Celsius # base temperature at Block 2  T2=2273. * Celsius # base temperature at Block 2
492  \end{python}  \end{python}
493  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
494    where we would like to save the output from the solver. This is simple enough:
495  \begin{python}  \begin{python}
496  t=0 * day  #our start time, usually zero  t=0 * day  #our start time, usually zero
497  tend=1. * day # - time to end simulation  tend=1. * day # - time to end simulation
# Line 315  h=(tend-t)/outputs #size of time step Line 501  h=(tend-t)/outputs #size of time step
501  print "Expected Number of time outputs is: ", (tend-t)/h  print "Expected Number of time outputs is: ", (tend-t)/h
502  i=0 #loop counter  i=0 #loop counter
503  \end{python}  \end{python}
504  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 model  as:  Now that we know our inputs we will build a domain using the
505     \verb Rectangle() function from \verb finley . The four arguments allow us to define our domain
506    \verb model  as:
507  \begin{python}  \begin{python}
508  #generate domain using rectangle  #generate domain using rectangle
509  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
510  \end{python}  \end{python}
511  \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}.
512
513  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 we discuss 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
514    set up our PDE so that it can be solved by \esc. The first step is to define the
515    type of PDE that we are trying to solve in each time step. In this example it is
516    a single linear PDE\footnote{in comparison to a system of PDEs which we discuss
517    later.}. We also need to state the values of our general form variables.
518  \begin{python}  \begin{python}
519  mypde=LinearPDE(blocks)  mypde=LinearPDE(blocks)
520  A=zeros((2,2)))  A=zeros((2,2)))
521  A[0,0]=kappa  A[0,0]=kappa
522  mypde.setValue(A=A, D=rhocp/h)  mypde.setValue(A=A, D=rhocp/h)
523  \end{python}  \end{python}
524  In a many cases it may be possible to decrease the computational time of the solver if the PDE is symmetric.  In a many cases it may be possible to decrease the computational time of the
525    solver if the PDE is symmetric.
526  Symmetry of a PDE is defined by;  Symmetry of a PDE is defined by;
527  \begin{equation}\label{eqn:symm}  \begin{equation}\label{eqn:symm}
528  A\hackscore{jl}=A\hackscore{lj}  A\hackscore{jl}=A\hackscore{lj}
529  \end{equation}  \end{equation}
530  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$. 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 enable symmetry via;  Symmetry is only dependent on the $A$ coefficient in the general form and the
531    other coefficients $D$ as well as the right hand side $Y$. From the above
532    definition we can see that our PDE is symmetric. The \verb LinearPDE  class
533    provides the method \method{checkSymmetry} to check if the given PDE is
534    symmetric. As our PDE is symmetrical we enable symmetry via;
535  \begin{python}  \begin{python}
536   myPDE.setSymmetryOn()   myPDE.setSymmetryOn()
537  \end{python}  \end{python}
538  Next we need to establish the initial temperature distribution \verb|T|. We need to  Next we need to establish the initial temperature distribution \verb|T|. We need
539  assign the value \verb|T1| to all sample points left to the contact interface at $x\hackscore{0}=\frac{mx}{2}$  to
540    assign the value \verb|T1| to all sample points left to the contact interface at
541    $x\hackscore{0}=\frac{mx}{2}$
542  and the value \verb|T2| right to the contact interface. \esc  and the value \verb|T2| right to the contact interface. \esc
543  provides the \verb|whereNegative| function to construct this. In fact,  provides the \verb|whereNegative| function to construct this. In fact,
544  \verb|whereNegative| returns the value $1$ at those sample points where the argument  \verb|whereNegative| returns the value $1$ at those sample points where the
545  has a negative value. Otherwise zero is returned. If \verb|x| are the $x\hackscore{0}$  argument
546    has a negative value. Otherwise zero is returned. If \verb|x| are the
547    $x\hackscore{0}$
548  coordinates of the sample points used to represent the temperature distribution  coordinates of the sample points used to represent the temperature distribution
549  then \verb|x-boundloc| gives us a negative value for  then \verb|x-boundloc| gives us a negative value for
550  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 352  the right of the interface. So with; Line 553  the right of the interface. So with;
553  # ... set initial temperature ....  # ... set initial temperature ....
554  T= T1*whereNegative(x-boundloc)+T2*(1-whereNegative(x-boundloc))  T= T1*whereNegative(x-boundloc)+T2*(1-whereNegative(x-boundloc))
555  \end{python}  \end{python}
556  we get the desired temperature distribution. To get the actual sample points \verb|x| we use  we get the desired temperature distribution. To get the actual sample points
557    \verb|x| we use
558  the  \verb|getX()| method of the function space \verb|Solution(blocks)|  the  \verb|getX()| method of the function space \verb|Solution(blocks)|
559  which is used to represent the solution of a PDE;  which is used to represent the solution of a PDE;
560  \begin{python}  \begin{python}
561  x=Solution(blocks).getX()  x=Solution(blocks).getX()
562  \end{python}  \end{python}
563  As \verb|x| are the sample points for the function space \verb|Solution(blocks)|  As \verb|x| are the sample points for the function space
564  the initial temperature \verb|T| is using these sample points for representation.  \verb|Solution(blocks)|
565  Although \esc is trying to be forgiving with the choice of sample points and to convert  the initial temperature \verb|T| is using these sample points for
566  where necessary the adjustment of the function space is not always possible. So it is  representation.
567    Although \esc is trying to be forgiving with the choice of sample points and to
568    convert
569    where necessary the adjustment of the function space is not always possible. So
570    it is
571  advisable to make a careful choice on the function space used.    advisable to make a careful choice on the function space used.
572
573  Finally we initialise 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 initialise an iteration loop to solve our PDE for all the time steps
574    we specified in the variable section. As the right hand side of the general form
575    is dependent on the previous values for temperature \verb T  across the bar this
576    must be updated in the loop. Our output at each time step is \verb T  the heat
577    distribution and \verb totT  the total heat in the system.
578  \begin{python}  \begin{python}
579  while t < tend:  while t < tend:
580      i+=1 #increment the counter      i+=1 #increment the counter
# Line 373  while t < tend: Line 583  while t < tend:
583      T=mypde.getSolution() #get the PDE solution      T=mypde.getSolution() #get the PDE solution
584      totE = integrate(rhocp*T) #get the total heat (energy) in the system      totE = integrate(rhocp*T) #get the total heat (energy) in the system
585  \end{python}  \end{python}
586  The last statement in this script calculates the total energy in the system as volume integral  The last statement in this script calculates the total energy in the system as
587  of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy should be get lost or added.  volume integral
588    of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy
589    should be get lost or added.
590  The total energy should stay constant for the example discussed here.  The total energy should stay constant for the example discussed here.
591
592  \subsection{Running the Script}  \subsection{Running the Script}
593  The script presented so far is available under  The script presented so far is available under
594  \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.
595  On most operating systems\footnote{The you can use \texttt{run-escript} launcher is not supported under {\it MS Windows} yet.} you can use the \program{run-escript} command  On most operating systems\footnote{The you can use \texttt{run-escript} launcher
596    is not supported under {\it MS Windows} yet.} you can use the
597    \program{run-escript} command
598  to launch {\it escript} scripts. For the example script use;  to launch {\it escript} scripts. For the example script use;
599  \begin{verbatim}  \begin{verbatim}
600  run-escript example01a.py  run-escript example01a.py
# Line 390  the python interpreter directly; Line 604  the python interpreter directly;
604  \begin{verbatim}  \begin{verbatim}
605  python example01a.py  python example01a.py
606  \end{verbatim}  \end{verbatim}
609
610  \begin{figure}  \begin{figure}
611  \begin{center}  \begin{center}
# Line 403  if the system is configured correctly (P Line 618  if the system is configured correctly (P
618  \subsection{Plotting the Total Energy}  \subsection{Plotting the Total Energy}
619  \sslist{example01b.py}  \sslist{example01b.py}
620
621  \esc does not include its own plotting capabilities. However, it is possible to use a variety of free \pyt packages for visualisation.  \esc does not include its own plotting capabilities. However, it is possible to
622  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.  use a variety of free \pyt packages for visualisation.
623  The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots.  Two types will be demonstrated in this cookbook;
624  For more complex visualisation tasks in particular, two and three dimensional problems we recommend the use of more advanced tools. For instance,  \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}}  \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and
625    \verb VTK  \footnote{\url{http://www.vtk.org/}} visualisation.
626    The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}}
627    and is good for basic graphs and plots.
628    For more complex visualisation tasks in particular, two and three dimensional
629    problems we recommend the use of more advanced tools. For instance,  \mayavi
630    \footnote{\url{http://code.enthought.com/projects/mayavi/}}
631  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
632  visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two dimensional PDE.  visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two
633    dimensional PDE.
634
635  For our simple granite block problem, we have two plotting tasks. Firstly, we are interested in showing the  For our simple granite block problem, we have two plotting tasks. Firstly, we
636  behaviour of the total energy over time and secondly, how the temperature distribution within the block is  are interested in showing the
637    behaviour of the total energy over time and secondly, how the temperature
638    distribution within the block is
640
641  The trick is to create a record of the time marks and the corresponding total energies observed.  The trick is to create a record of the time marks and the corresponding total
642    energies observed.
643  \pyt provides the concept of lists for this. Before  \pyt provides the concept of lists for this. Before
644  the time loop is opened we create empty lists for the time marks \verb|t_list| and the total energies \verb|E_list|.  the time loop is opened we create empty lists for the time marks \verb|t_list|
645  After the new temperature has been calculated by solving the PDE we append the new time marker and the total energy value for that time  and the total energies \verb|E_list|.
646  to the corresponding list using the \verb|append| method. With these modifications our script looks as follows:  After the new temperature has been calculated by solving the PDE we append the
647    new time marker and the total energy value for that time
648    to the corresponding list using the \verb|append| method. With these
649    modifications our script looks as follows:
650  \begin{python}  \begin{python}
651  t_list=[]  t_list=[]
652  E_list=[]  E_list=[]
# Line 431  while t<tend: Line 659  while t<tend:
659        t_list.append(t)   # add current time mark to record        t_list.append(t)   # add current time mark to record
660        E_list.append(totE) # add current total energy to record        E_list.append(totE) # add current total energy to record
661  \end{python}  \end{python}
662  To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs to be loaded before used;  To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs
663    to be loaded before used;
664  \begin{python}  \begin{python}
665  import pylab as pl # plotting package.  import pylab as pl # plotting package.
666  \end{python}  \end{python}
667  Here we are not using the \verb|from pylab import *| in order to avoid name clashes for function names  Here we are not using the \verb|from pylab import *| in order to avoid name
668    clashes for function names
669  within \esc.  within \esc.
670
671  The following statements are added to the script after the time loop has been completed;  The following statements are added to the script after the time loop has been
672    completed;
673  \begin{python}  \begin{python}
674  pl.plot(t_list,E_list)  pl.plot(t_list,E_list)
675  pl.title("Total Energy")  pl.title("Total Energy")
676  pl.axis([0,max(t_list),0,max(E_list)*1.1])  pl.axis([0,max(t_list),0,max(E_list)*1.1])
677  pl.savefig("totE.png")  pl.savefig("totE.png")
678  \end{python}  \end{python}
679  The first statement hands over the time marks and corresponding total energies to the plotter.  The first statement hands over the time marks and corresponding total energies
680    to the plotter.
681  The second statment is setting the title for the plot. The third statement  The second statment is setting the title for the plot. The third statement
682  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.
683
684  The last statement renders the plot and writes the  The last statement renders the plot and writes the
685  result into the file \verb|totE.png| which can be displayed by (almost) any image viewer.  result into the file \verb|totE.png| which can be displayed by (almost) any
686  As expected the total energy is constant over time, see \reffig{fig:onedheatout1}.  image viewer.
687    As expected the total energy is constant over time, see
688    \reffig{fig:onedheatout1}.
689
690  \subsection{Plotting the Temperature Distribution}  \subsection{Plotting the Temperature Distribution}
691  \label{sec: plot T}  \label{sec: plot T}
692  \sslist{example01c.py}  \sslist{example01c.py}
693  For plotting the spatial distribution of the temperature we need to modify the strategy we have used  For plotting the spatial distribution of the temperature we need to modify the
694  for the total energy. Instead of producing a final plot at the end we will generate a  strategy we have used
695  picture at each time step which can be browsed as a slide show or composed into a movie.  for the total energy. Instead of producing a final plot at the end we will
696  The first problem we encounter is that if we produce an image at each time step we need  generate a
697    picture at each time step which can be browsed as a slide show or composed into
698    a movie.
699    The first problem we encounter is that if we produce an image at each time step
700    we need
701  to make sure that the images previously generated are not overwritten.  to make sure that the images previously generated are not overwritten.
702
703  To develop an incrementing file name we can use the following convention. It is convenient to  To develop an incrementing file name we can use the following convention. It is
704  put all image file showing the same variable - in our case the temperature distribution -  convenient to
705  into a separate directory. As part of the \verb|os| module\footnote{The \texttt{os} module provides  put all image file showing the same variable - in our case the temperature
706  a powerful interface to interact with the operating system, see \url{http://docs.python.org/library/os.html}.} \pyt  distribution -
707    into a separate directory. As part of the \verb|os| module\footnote{The
708    \texttt{os} module provides
709    a powerful interface to interact with the operating system, see
710    \url{http://docs.python.org/library/os.html}.} \pyt
711  provides the  \verb|os.path.join| command to build file and  provides the  \verb|os.path.join| command to build file and
712  directory names in a platform independent way. Assuming that  directory names in a platform independent way. Assuming that
713  \verb|save_path| is name of directory we want to put the results the command is;  \verb|save_path| is name of directory we want to put the results the command
714    is;
715  \begin{python}  \begin{python}
716  import os  import os
717  os.path.join(save_path, "tempT%03d.png"%i )  os.path.join(save_path, "tempT%03d.png"%i )
718  \end{python}  \end{python}
719  where \verb|i| is the time step counter.  where \verb|i| is the time step counter.
720  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, for example a single sub-folder called \verb data  would be defined by;  There are two arguments to the \verb join  command. The
721    \verb save_path  variable is a predefined string pointing to the directory we want to save our
722    data, for example a single sub-folder called \verb data  would be defined by;
723  \begin{verbatim}  \begin{verbatim}
724  save_path = "data"  save_path = "data"
725  \end{verbatim}  \end{verbatim}
# Line 481  while a sub-folder of \verb data  called Line 727  while a sub-folder of \verb data  called
727  \begin{verbatim}  \begin{verbatim}
728  save_path = os.path.join("data","example01")  save_path = os.path.join("data","example01")
729  \end{verbatim}  \end{verbatim}
730  The second argument of \verb join \xspace contains a string which is the file name or subdirectory name. We can use the operator \verb|%| to use the value of \verb|i| as part of our filename. The sub-string \verb %03d  indicates that we want to substitude a value into the name;  The second argument of \verb join \xspace contains a string which is the file
731    name or subdirectory name. We can use the operator \verb|%| to use the value of
732    \verb|i| as part of our filename. The sub-string \verb %03d  indicates that we
733    want to substitude a value into the name;
734  \begin{itemize}  \begin{itemize}
735   \item \verb 0  means that small numbers should have leading zeroes;   \item \verb 0  means that small numbers should have leading zeroes;
736   \item \verb 3  means that numbers should be written using at least 3 digits; and   \item \verb 3  means that numbers should be written using at least 3 digits;
737    and
738   \item \verb d  means that the value to substitute will be an integer.   \item \verb d  means that the value to substitute will be an integer.
739  \end{itemize}  \end{itemize}
740  To actually substitute the value of \verb|i| into the name write \verb %i  after the string.  To actually substitute the value of \verb|i| into the name write \verb %i  after
741  When done correctly, the output files from this command would be place in the directory defined by \verb save_path  as;  the string.
742    When done correctly, the output files from this command would be place in the
743    directory defined by \verb save_path  as;
744  \begin{verbatim}  \begin{verbatim}
745  blockspyplot001.png  blockspyplot001.png
746  blockspyplot002.png  blockspyplot002.png
# Line 501  A sub-folder check/constructor is availa Line 753  A sub-folder check/constructor is availa
753  \begin{verbatim}  \begin{verbatim}
754  mkDir(save_path)  mkDir(save_path)
755  \end{verbatim}  \end{verbatim}
756  will check for the existence of \verb save_path  and if missing, make the required directories.  will check for the existence of \verb save_path  and if missing, make the
757    required directories.
758
759  We start by modifying our solution script.  We start by modifying our solution script.
760  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 sample points used to represent  Prior to the \verb|while| loop we will need to extract our finite solution
761  the temperature as a \pyt list of tuples or a \numpy array as requested by the plotting function.  points to a data object that is compatible with \mpl. First we create the node
762  We need to convert the array \verb|x| previously set as \verb|Solution(blocks).getX()| into a \pyt list  coordinates of the sample points used to represent
763  and then to a \numpy array. The $x\hackscore{0}$ component is then extracted via an array slice to the variable \verb|plx|;  the temperature as a \pyt list of tuples or a \numpy array as requested by the
764    plotting function.
765    We need to convert the array \verb|x| previously set as
766    \verb|Solution(blocks).getX()| into a \pyt list
767    and then to a \numpy array. The $x\hackscore{0}$ component is then extracted via
768    an array slice to the variable \verb|plx|;
769  \begin{python}  \begin{python}
770  import numpy as np # array package.  import numpy as np # array package.
771  #convert solution points for plotting  #convert solution points for plotting
# Line 521  plx = plx[:,0] # extract x locations Line 779  plx = plx[:,0] # extract x locations
779  \includegraphics[width=4in]{figures/blockspyplot001}  \includegraphics[width=4in]{figures/blockspyplot001}
780  \includegraphics[width=4in]{figures/blockspyplot050}  \includegraphics[width=4in]{figures/blockspyplot050}
781  \includegraphics[width=4in]{figures/blockspyplot200}  \includegraphics[width=4in]{figures/blockspyplot200}
782  \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps $1$, $50$ and $200$.}  \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps
783    $1$, $50$ and $200$.}
784  \label{fig:onedheatout}  \label{fig:onedheatout}
785  \end{center}  \end{center}
786  \end{figure}  \end{figure}
787
788  We use the same techniques provided by \mpl as we have used to plot the total energy over time.  We use the same techniques provided by \mpl as we have used to plot the total
789  For each time step we generate a plot of the temperature distribution and save each to a file.  energy over time.
790  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.  For each time step we generate a plot of the temperature distribution and save
791  Finally, the figure is saved to a \verb|*.png| file and cleared for the following iteration.  each to a file.
792    The following is appended to the end of the \verb while  loop and creates one
793    figure of the temperature distribution. We start by converting the solution to a
794    tuple and then plotting this against our \textit{x coordinates} \verb plx  we
795    have generated before. We add a title to the diagram before it is rendered into
796    a file.
797    Finally, the figure is saved to a \verb|*.png| file and cleared for the
798    following iteration.
799  \begin{python}  \begin{python}
800  # ... start iteration:  # ... start iteration:
801  while t<tend:  while t<tend:
# Line 547  while t<tend: Line 813  while t<tend:
813  Some results are shown in \reffig{fig:onedheatout}.  Some results are shown in \reffig{fig:onedheatout}.
814
815  \subsection{Make a video}  \subsection{Make a video}
816  Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. The \verb mencoder command is Linux only, so other platform users need to use an alternative video encoder.  Our saved plots from the previous section can be cast into a video using the
817    following command appended to the end of the script. The \verb mencoder command
818    is Linux only, so other platform users need to use an alternative video encoder.
819  \begin{python}  \begin{python}
820  # compile the *.png files to create a *.avi videos that show T change  # compile the *.png files to create a *.avi videos that show T change
821  # with time. This operation uses Linux mencoder. For other operating  # with time. This operation uses Linux mencoder. For other operating

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