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

revision 2974 by artak, Wed Mar 3 01:21:09 2010 UTC revision 2975 by ahallam, Thu Mar 4 05:44:12 2010 UTC
# Line 37  which is defined as: Line 37  which is defined as:
37  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 $\kappa$ is the thermal
38  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  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
39  parameters are \textbf{constant}.  parameters are \textbf{constant}.
40  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 position along the iron bar $x$.  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$.
41
42  \subsection{PDEs and the General Form}  \subsection{PDEs and the General Form}
43  Potentially, 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 a problem this way. Alternatively, computers can be used to solve these kinds of problems. To do this, a numerical approach is required to discretise  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
44  the PDE \refEq{eqn:hd} in time and space, then we left with a finite number of equations for a finite number of spatial points and time steps in the model. While discretisation introduces approximations and a degree of error, a sufficiently sampled model is generally accurate enough to satisfy the requirements of the modeller.  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.
45
46  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.  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.
47
# Line 72  Substituting \refEq{eqn:Tbeuler} into \r Line 72  Substituting \refEq{eqn:Tbeuler} into \r
72  \end{equation}  \end{equation}
73  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 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
74  approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which  approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which
75  are not discussed here. However, this scheme has a major disadvantage, namely depending on the  are not discussed here. However, the \textbf{forward Euler} scheme has a major disadvantage. Depending on the
76  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. The term \textit{stable} means  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.
that the approximation of the temperature will not grow beyond its initial bounds and become non-physical.
77  The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of  The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of
78  physically correct problem set-up the temperature approximation remains physical for all time steps.  physically correct problem set-up the temperature approximation remains physical for all time steps.
79  The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}  The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}
80  is sufficiently small, thus a good approximation of the true temperature is computed. It is  is sufficiently small, thus a good approximation of the true temperature is computed. It is
81  therefore crucial that the user remains critical about his/her results and for instance compares  therefore crucial that the user remains sceptical about their results and for instance compares
82  the results for different time and spatial step sizes.  the results for different time and spatial step sizes for correlation.
83
84  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
85  differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem  differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem
86  we want to use \esc.  we want to use \esc.
87
88  In \esc any given PDE described by 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 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
89  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}
90  is described by;  is described by;
91  \begin{equation}\label{eqn:commonform nabla}  \begin{equation}\label{eqn:commonform nabla}
# Line 121  With the PDE sufficiently modified, cons Line 120  With the PDE sufficiently modified, cons
120  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.  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.
121  We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}.  We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}.
122
123  We make 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 insulated, so we need
124  to add an appropriate boundary condition to prevent  to add an appropriate boundary condition to prevent
125  any loss or inflow of energy at boundary of our domain. Mathematically this is expressed by prescribing  any loss or inflow of energy at the boundary of our domain. Mathematically this is expressed by prescribing
126  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified one dimensional model this is expressed  the heat flux $\kappa \frac{\partial T}{\partial x}$  to zero. In our simplified one dimensional model this is expressed
127  in the form;  in the form;
128  \begin{equation}  \begin{equation}
# Line 167  for the boundary of the domain. This is Line 166  for the boundary of the domain. This is
166
167  \subsection{Outline of the Implementation}  \subsection{Outline of the Implementation}
168  \label{sec:outline}  \label{sec:outline}
169  To solve the heat diffusion equation (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 we discuss later in details have 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 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
170  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 is as follows:  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}.
171  \begin{enumerate}  % \begin{enumerate}
172   \item create domain  %  \item create domain
173   \item create PDE  %  \item create PDE
174   \item while end time not reached:  %  \item while end time not reached:
175  \begin{enumerate}  % \begin{enumerate}
176   \item set PDE coefficients  %  \item set PDE coefficients
177   \item solve PDE  %  \item solve PDE
178   \item update time marker  %  \item update time marker
179  \end{enumerate}  % \end{enumerate}
180  \item end of calculation  % \item end of calculation
181  \end{enumerate}  % \end{enumerate}
182
183    \begin{figure}[h!]
184     \centering
185       \includegraphics[width=1in]{figures/workflow.png}
186       \caption{Workflow for developing an \esc model and solution.}
187       \label{fig:wf}
188    \end{figure}
189
190  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 \textbf{objects}. The nice feature of an object is that it is defined by its usage and features
191  rather than its actual representation. So we will create a domain object to describe the geometry of the two  rather than its actual representation. So we will create a domain object to describe the geometry of the two
192  granite blocks. Then we define PDEs and spatially distributed values such as the temperature  granite blocks. Then we define PDEs and spatially distributed values such as the temperature
# Line 189  on this domain. Similarly, to define a P Line 196  on this domain. Similarly, to define a P
196  \begin{figure}[t]  \begin{figure}[t]
197   \centering   \centering
198     \includegraphics[width=6in]{figures/functionspace.pdf}     \includegraphics[width=6in]{figures/functionspace.pdf}
\label{fig:fs}
199     \caption{\esc domain construction overview}     \caption{\esc domain construction overview}
200       \label{fig:fs}
201  \end{figure}  \end{figure}
202
203  \subsection{The Domain Constructor in \esc}  \subsection{The Domain Constructor in \esc}
204  \label{ss:domcon}  \label{ss:domcon}
205  It is helpful to have a better understanding how spatially distributed value such as temperature or PDE coefficients are interpreted in \esc. Again  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.
from the user's point of view the representation of these spatially distributed values is not relevant.
206
207  There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is  There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is
208  a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of  a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of
209  elements or cells to be used along the length and height, see \reffig{fig:fs}. Any spatially distributed value  elements or cells to be used along the length and height, see \reffig{fig:fs}. Any spatially distributed value
210  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.  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.
211  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
212  approximation becomes better when more elements are used. However, computational costs and compute time grow with the number of  approximation becomes better when more elements are used. However, computational costs and time grow with the number of
213  elements being used. It is therefore important that you find the right balance between the demand in accuracy and acceptable resource usage.  elements being used. It is therefore important that you find the right balance between the demand in accuracy and acceptable resource usage.
214
215  In general, one can think about a domain object as a composition of nodes and elements.  In general, one can think about a domain object as a composition of nodes and elements.
# Line 226  the  \verb|Function(domain).getX()| retu Line 232  the  \verb|Function(domain).getX()| retu
232
233  This distinction between different representations of spatially distributed values  This distinction between different representations of spatially distributed values
234  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.
235  The coefficients of a PDE do not need be continuous, thus this qualifies as a \verb|Function()| type.  The coefficients of a PDE do not need to be continuous, thus this qualifies as a \verb|Function()| type.
236  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.
237  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.
238  \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()|
# Line 235  These transformations, which are called Line 241  These transformations, which are called
241
242  Later in this introduction we discuss how  Later in this introduction we discuss how
243  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
244  thermal conductivities $kappa$. A very powerful technique to define these types of PDE  thermal conductivities $\kappa$. A very powerful technique to define these types of PDE
245  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.
246  This is simplifying PDE coefficient and flux definitions. It makes scripting much easier. We will discuss this technique in Section~\ref{STEADY-STATE HEAT REFRACTION}.  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}.
247
248
249  \subsection{A Clarification for the 1D Case}  \subsection{A Clarification for the 1D Case}
# Line 280  It is generally a good idea to import al Line 286  It is generally a good idea to import al
286
287  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.  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.
288
289  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. 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.  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.
290  \begin{python}  \begin{python}
291  mx = 500.*m #meters - model length  mx = 500.*m #meters - model length
292  my = 100.*m #meters - model width  my = 100.*m #meters - model width
# Line 372  of $\rho c\hackscore{p} T$ over the bloc Line 378  of $\rho c\hackscore{p} T$ over the bloc
378  The total energy should stay constant for the example discussed here.  The total energy should stay constant for the example discussed here.
379
380  \subsection{Running the Script}  \subsection{Running the Script}
381  The script presented so for is available under  The script presented so far is available under
382  \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.
383  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 is not supported under {\it MS Windows} yet.} you can use the \program{run-escript} command
384  to launch {\it escript} scripts. For the example script use;  to launch {\it escript} scripts. For the example script use;
# Line 400  if the system is configured correctly (P Line 406  if the system is configured correctly (P
406  \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 use a variety of free \pyt packages for visualisation.
407  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.  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.
408  The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots.  The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots.
409  For more complex visualisation tasks in particular when it comes to  two and three dimensional problems we recommend to use more advanced tools for instance  \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}}  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/}}
410  which bases on the \verb|VTK| toolkit. We discuss the usage of \verb|VTK| based  which is based upon the \verb|VTK| toolkit. The usage of \verb|VTK| based
411  visualization in Chapter~\ref{Sec:2DHD} where focus on a two dimensional PDE.  visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two dimensional PDE.
412
413  For our simple problem we have two plotting tasks: we are interested in showing the  For our simple granite block problem, we have two plotting tasks. Firstly, we are interested in showing the
414  behaviour of the total energy over time and how the temperature distribution within the block is  behaviour of the total energy over time and secondly, how the temperature distribution within the block is
416
417  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 energies observed.
418  \pyt provides the concept of lists for this. Before  \pyt provides the concept of lists for this. Before
419  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| and the total energies \verb|E_list|.
420  After the new temperature as been calculated by solving the PDE we append the new time marker and total energy  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
421  to the corresponding list using the \verb|append| method. With these modifications our script looks as follows:  to the corresponding list using the \verb|append| method. With these modifications our script looks as follows:
422  \begin{python}  \begin{python}
423  t_list=[]  t_list=[]
# Line 425  while t<tend: Line 431  while t<tend:
431        t_list.append(t)   # add current time mark to record        t_list.append(t)   # add current time mark to record
432        E_list.append(totE) # add current total energy to record        E_list.append(totE) # add current total energy to record
433  \end{python}  \end{python}
434  To plot $t$ over $totE$ we use the \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 to be loaded before used;
435  \begin{python}  \begin{python}
436  import pylab as pl # plotting package.  import pylab as pl # plotting package.
437  \end{python}  \end{python}
# Line 451  As expected the total energy is constant Line 457  As expected the total energy is constant
457  \sslist{example01c.py}  \sslist{example01c.py}
458  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 strategy we have used
459  for the total energy. Instead of producing a final plot at the end we will generate a  for the total energy. Instead of producing a final plot at the end we will generate a
460  picture at each time step which can be browsed as slide show or composed to a movie.  picture at each time step which can be browsed as a slide show or composed into a movie.
461  The first problem we encounter is that if we produce an image in each time step we need  The first problem we encounter is that if we produce an image at each time step we need
462  to make sure that the images previously generated are not overwritten.  to make sure that the images previously generated are not overwritten.
463
464  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 convenient to

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