ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 2982 - (hide annotations)
Wed Mar 10 04:27:47 2010 UTC (9 years, 6 months ago) by caltinay
File MIME type: application/x-tex
File size: 37699 byte(s)
Some minor corrections to the cookbook (mainly word repetitions and typos).

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

  ViewVC Help
Powered by ViewVC 1.1.26