ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 2982 - (show annotations)
Wed Mar 10 04:27:47 2010 UTC (9 years, 10 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).

2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % 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 \begin{figure}[h!]
15 \centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}}
16 \caption{Example 1: Temperature differential along a single interface between
17 two granite blocks.}
18 \label{fig:onedgbmodel}
19 \end{figure}
21 \section{Example 1: One Dimensional Heat Diffusion in Granite}
22 \label{Sec:1DHDv00}
24 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 We assume that the system is insulated.
29 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
31 cooler one until both
32 blocks have the same temperature.
34 \subsection{1D Heat Diffusion Equation}
35 We can model the heat distribution of this problem over time using the one
36 dimensional heat diffusion equation\footnote{A detailed discussion on how the
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:
41 \begin{equation}
42 \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}
45 \end{equation}
46 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 parameters are \textbf{constant}.
53 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 block-block interface $x$.
63 \subsection{PDEs and the General Form}
64 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.
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.
78 For time discretisation we use the Backward Euler approximation
79 scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is
80 based on the approximation
81 \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 where the upper index $n$ denotes the n\textsuperscript{th} time step. So one
92 has
93 \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 \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 \label{eqn:hddisc}
105 \end{equation}
106 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 This approach is called the \textbf{forward Euler} scheme. This scheme can
110 provide some computational advantages, which
111 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 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 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 The user needs to keep in mind that the discretisation error introduced by
123 \refEq{eqn:beuler}
124 is sufficiently small, thus a good approximation of the true temperature is
125 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.
129 To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear
130 differential equation \refEq{eqn:hddisc} which only includes spatial
131 derivatives. To solve this problem we want to use \esc.
133 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 $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}
139 is described by;
140 \begin{equation}\label{eqn:commonform nabla}
141 -\nabla\cdot(A\cdot\nabla u) + Du = f
142 \end{equation}
143 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 \begin{equation}
149 \nabla = \frac{\partial}{\partial x}
150 \end{equation}
151 and we can write \refEq{eqn:commonform nabla} as;
152 \begin{equation}\label{eqn:commonform}
153 -A\frac{\partial^{2}u}{\partial x^{2}} + Du = f
154 \end{equation}
155 if $A$ is constant. To match this simplified general form to our problem
156 \refEq{eqn:hddisc}
157 we rearrange \refEq{eqn:hddisc};
158 \begin{equation}
159 \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 \label{eqn:hdgenf}
162 \end{equation}
163 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 \begin{equation}\label{ESCRIPT SET}
172 u=T^{(n)};
173 A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho
174 c\hackscore p}{h} T^{(n-1)}
175 \end{equation}
177 \subsection{Boundary Conditions}
179 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 We discuss the Dirichlet boundary condition in our second example presented in
189 Section~\ref{Sec:1DHDv0}.
191 However, for this example we have made the model assumption that the system is
192 insulated, so we need to add an appropriate boundary condition to prevent
193 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 in the form;
198 \begin{equation}
199 \kappa \frac{\partial T}{\partial x} = 0
200 \end{equation}
201 or in a more general case as
202 \begin{equation}\label{NEUMAN 1}
203 \kappa \nabla T \cdot n = 0
204 \end{equation}
205 where $n$ is the outer normal field \index{outer normal field} at the surface
206 of the domain.
207 The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$.
208 In fact, the term $\nabla T \cdot n$ is the normal derivative of
209 the temperature $T$. Other notations used here are\footnote{The \esc notation
210 for the normal
211 derivative is $T\hackscore{,i} n\hackscore i$.};
212 \begin{equation}
213 \nabla T \cdot n = \frac{\partial T}{\partial n} \; .
214 \end{equation}
215 A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neumann boundary
216 condition} for the PDE.
218 The PDE \refEq{eqn:hdgenf}
219 and the Neumann boundary condition~\ref{eqn:hdgenf} (potentially together with
220 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 cases we use the term partial differential equation but in fact it is a
225 boundary value problem.
226 It is important to keep in mind that boundary conditions need to be complete and
227 consistent in the sense that
228 at any point on the boundary either a Dirichlet or a Neumann boundary condition
229 must be set.
231 Conveniently, \esc makes a default assumption on the boundary conditions which
232 the user may modify where appropriate.
233 For a problem of the form in~\refEq{eqn:commonform nabla} the default
234 condition\footnote{In the \esc user guide which uses the Einstein convention
235 this is written as
236 $n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is;
237 \begin{equation}\label{NEUMAN 2}
238 -n\cdot A \cdot\nabla u = 0
239 \end{equation}
240 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 condition translates into
247 \begin{equation}\label{NEUMAN 2b}
248 \kappa \frac{\partial T}{\partial x} = 0
249 \end{equation}
250 for the boundary of the domain. This is identical to the Neumann boundary
251 condition we want to set. \esc will take care of this condition for us. We
252 discuss the Dirichlet boundary condition later.
254 \subsection{Outline of the Implementation}
255 \label{sec:outline}
256 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 is described in \reffig{fig:wf}.
267 % \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 \caption{Workflow for developing an \esc model and solution}
283 \label{fig:wf}
284 \end{figure}
286 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.
298 \begin{figure}[t]
299 \centering
300 \includegraphics[width=6in]{figures/functionspace.pdf}
301 \caption{\esc domain construction overview}
302 \label{fig:fs}
303 \end{figure}
305 \subsection{The Domain Constructor in \esc}
306 \label{ss:domcon}
307 Whilst it is not strictly relevant or necessary, a better understanding of
308 how values are spatially distributed (\textit{e.g.} Temperature) and how PDE
309 coefficients are interpreted in \esc can be helpful.
311 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 The quality of the approximation depends - besides other factors - mainly on the
323 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.
329 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 To represent spatially distributed values the user can use
334 the values at the nodes, at the elements in the interior of the domain or at the
335 elements located on the surface of the domain.
336 The different approach used to represent values is called \textbf{function
337 space} and is attached to all objects
338 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 \begin{enumerate}
341 \item the nodes, called by \verb|ContinuousFunction(domain)| ;
342 \item the elements/cells, called by \verb|Function(domain)| ; and
343 \item the boundary, called by \verb|FunctionOnBoundary(domain)|.
344 \end{enumerate}
345 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 \verb|Function(domain).getX()| returns the coordinates of numerical
352 integration points within elements, see \reffig{fig:fs}.
354 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 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 These transformations, which are called \textbf{interpolation} are invoked
370 automatically by \esc if needed.
372 Later in this introduction we discuss how
373 to define specific areas of geometry with different materials which are
374 represented by different material coefficients such as the
375 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
384 \subsection{A Clarification for the 1D Case}
386 It is necessary for clarification that we revisit our general PDE from
387 \refeq{eqn:commonform nabla} for a two dimensional domain. \esc is inherently
388 designed to solve problems that are multi-dimensional and so
389 \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 \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 Notice that for the higher dimensional case $A$ becomes a matrix. It is also
402 important to notice that the usage of the Nabla operator creates
403 a compact formulation which is also independent from the spatial dimension.
404 To make the general PDE \refEq{eqn:commonform2D} one dimensional as
405 shown in \refEq{eqn:commonform} we need to set
406 \begin{equation}
407 A\hackscore{00}=A; A\hackscore{01}=A\hackscore{10}=A\hackscore{11}=0
408 \end{equation}
411 \subsection{Developing a PDE Solution Script}
412 \label{sec:key}
413 \sslist{example01a.py}
414 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 that we will require.
423 \begin{python}
424 from esys.escript import *
425 # 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 \end{python}
433 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 domain. The module \verb|unitsSI| provides support for SI unit definitions with
438 our variables.
440 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.
451 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 \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 Thus the temperature does not vary with $y$. Hence, the model parameters can be
467 defined as follows; note we have used the \verb|unitsSI| convention to make sure
468 all our input units are converted to SI.
469 \begin{python}
470 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 \end{python}
476 The material constants and the temperature variables must also be defined. For
477 the granite in the model they are defined as:
478 \begin{python}
479 #PDE related
480 rho = 2750. *kg/m**3 #kg/m^{3} density of iron
481 cp = 790.*J/(kg*K) # J/Kg.K thermal capacity
482 rhocp = rho*cp
483 kappa = 2.2*W/m/K # watts/m.Kthermal conductivity
484 qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source
485 T1=20 * Celsius # initial temperature at Block 1
486 T2=2273. * Celsius # base temperature at Block 2
487 \end{python}
488 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 \begin{python}
491 t=0 * day #our start time, usually zero
492 tend=1. * day # - time to end simulation
493 outputs = 200 # number of time steps required.
494 h=(tend-t)/outputs #size of time step
495 #user warning statement
496 print "Expected Number of time outputs is: ", (tend-t)/h
497 i=0 #loop counter
498 \end{python}
499 Now that we know our inputs we will build a domain using the
500 \verb|Rectangle()| function from \FINLEY. The four arguments allow us to
501 define our domain \verb|model| as:
502 \begin{python}
503 #generate domain using rectangle
504 blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
505 \end{python}
506 \verb|blocks| now describes a domain in the manner of Section \ref{ss:domcon}.
508 With a domain and all the required variables established, it is now possible to
509 set up our PDE so that it can be solved by \esc. The first step is to define the
510 type of PDE that we are trying to solve in each time step. In this example it is
511 a single linear PDE\footnote{in contrast to a system of PDEs which we discuss
512 later.}. We also need to state the values of our general form variables.
513 \begin{python}
514 mypde=LinearPDE(blocks)
515 A=zeros((2,2)))
516 A[0,0]=kappa
517 mypde.setValue(A=A, D=rhocp/h)
518 \end{python}
519 In many cases it may be possible to decrease the computational time of the
520 solver if the PDE is symmetric.
521 Symmetry of a PDE is defined by;
522 \begin{equation}\label{eqn:symm}
523 A\hackscore{jl}=A\hackscore{lj}
524 \end{equation}
525 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 definition we can see that our PDE is symmetric. The \verb|LinearPDE| class
528 provides the method \method{checkSymmetry} to check if the given PDE is
529 symmetric. As our PDE is symmetrical we enable symmetry via;
530 \begin{python}
531 myPDE.setSymmetryOn()
532 \end{python}
533 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 and the value \verb|T2| right to the contact interface. \esc
538 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 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 \begin{python}
547 # ... set initial temperature ....
548 T= T1*whereNegative(x[0]-boundloc)+T2*(1-whereNegative(x[0]-boundloc))
549 \end{python}
550 we get the desired temperature distribution. To get the actual sample points
551 \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 \begin{python}
554 x=Solution(blocks).getX()
555 \end{python}
556 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 convert
562 where necessary the adjustment of the function space is not always possible. So
563 it is advisable to make a careful choice on the function space used.
565 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 \begin{python}
571 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 The last statement in this script calculates the total energy in the system as
579 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 The total energy should stay constant for the example discussed here.
583 \subsection{Running the Script}
584 The script presented so far is available under
585 \verb|example01a.py|. You can edit this file with your favourite text editor.
586 On most operating systems\footnote{The \texttt{run-escript} launcher is not
587 supported under {\it MS Windows} yet.} you can use the
588 \program{run-escript} command
589 to launch {\it escript} scripts. For the example script use;
590 \begin{verbatim}
591 run-escript example01a.py
592 \end{verbatim}
593 The program will print a progress report. Alternatively, you can use
594 the python interpreter directly;
595 \begin{verbatim}
596 python example01a.py
597 \end{verbatim}
598 if the system is configured correctly (please talk to your system
599 administrator).
601 \begin{figure}
602 \begin{center}
603 \includegraphics[width=4in]{figures/ttblockspyplot150}
604 \caption{Example 1b: Total Energy in the Blocks over Time (in seconds)}
605 \label{fig:onedheatout1}
606 \end{center}
607 \end{figure}
609 \subsection{Plotting the Total Energy}
610 \sslist{example01b.py}
612 \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 \verb|VTK|\footnote{\url{http://www.vtk.org/}}.
617 The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}}
618 and is good for basic graphs and plots.
619 For more complex visualisation tasks, in particular two and three dimensional
620 problems we recommend the use of more advanced tools. For instance, \mayavi
621 \footnote{\url{http://code.enthought.com/projects/mayavi/}}
622 which is based upon the \verb|VTK| toolkit. The usage of \verb|VTK| based
623 visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two
624 dimensional PDE.
626 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 distribution within the block is developing over time.
630 Let us start with the first task.
632 The idea is to create a record of the time marks and the corresponding total
633 energies observed.
634 \pyt provides the concept of lists for this. Before
635 the time loop is opened we create empty lists for the time marks \verb|t_list|
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 \begin{python}
642 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 \end{python}
653 To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs
654 to be loaded before use;
655 \begin{python}
656 import pylab as pl # plotting package.
657 \end{python}
658 Here we are not using \verb|from pylab import *| in order to avoid name
659 clashes for function names within \esc.
661 The following statements are added to the script after the time loop has been
662 completed;
663 \begin{python}
664 pl.plot(t_list,E_list)
665 pl.title("Total Energy")
666 pl.axis([0,max(t_list),0,max(E_list)*1.1])
667 pl.savefig("totE.png")
668 \end{python}
669 The first statement hands over the time marks and corresponding total energies
670 to the plotter.
671 The second statement sets the title for the plot. The third statement
672 sets the axis ranges. In most cases these are set appropriately by the plotter.
674 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 As expected the total energy is constant over time, see
677 \reffig{fig:onedheatout1}.
679 \subsection{Plotting the Temperature Distribution}
680 \label{sec: plot T}
681 \sslist{example01c.py}
682 For plotting the spatial distribution of the temperature we need to modify the
683 strategy we have used for the total energy.
684 Instead of producing a final plot at the end we will generate a
685 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 we need to make sure that the images previously generated are not overwritten.
690 To develop an incrementing file name we can use the following convention. It is
691 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 a powerful interface to interact with the operating system, see
695 \url{http://docs.python.org/library/os.html}.} \pyt
696 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 \begin{python}
701 import os
702 os.path.join(save_path, "tempT%03d.png"%i )
703 \end{python}
704 where \verb|i| is the time step counter.
705 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 \begin{verbatim}
709 save_path = "data"
710 \end{verbatim}
711 while a sub-folder of \verb|data| called \verb|example01| would be defined by;
712 \begin{verbatim}
713 save_path = os.path.join("data","example01")
714 \end{verbatim}
715 The second argument of \verb|join| contains a string which is the file
716 name or subdirectory name. We can use the operator \verb|%| to use the value of
717 \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 \begin{itemize}
720 \item \verb 0 means that small numbers should have leading zeroes;
721 \item \verb 3 means that numbers should be written using at least 3 digits;
722 and
723 \item \verb d means that the value to substitute will be a decimal integer.
724 \end{itemize}
725 To actually substitute the value of \verb|i| into the name write \verb|%i| after
726 the string.
727 When done correctly, the output files from this command will be placed in the
728 directory defined by \verb save_path as;
729 \begin{verbatim}
730 blockspyplot001.png
731 blockspyplot002.png
732 blockspyplot003.png
733 ...
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 will check for the existence of \verb save_path and if missing, create the
742 required directories.
744 We start by modifying our solution script.
745 Prior to the \verb|while| loop we need to extract our finite solution
746 points to a data object that is compatible with \mpl. First we create the node
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 \begin{python}
755 import numpy as np # array package.
756 #convert solution points for plotting
757 plx = x.toListOfTuples()
758 plx = np.array(plx) # convert to tuple to numpy array
759 plx = plx[:,0] # extract x locations
760 \end{python}
762 \begin{figure}
763 \begin{center}
764 \includegraphics[width=4in]{figures/blockspyplot001}
765 \includegraphics[width=4in]{figures/blockspyplot050}
766 \includegraphics[width=4in]{figures/blockspyplot200}
767 \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps
768 $1$, $50$ and $200$}
769 \label{fig:onedheatout}
770 \end{center}
771 \end{figure}
773 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 The following is appended to the end of the \verb|while| loop and creates one
778 figure of the temperature distribution. We start by converting the solution to a
779 tuple and then plotting this against our \textit{x coordinates} \verb|plx| we
780 have generated before. We add a title to the diagram before it is rendered into
781 a file.
782 Finally, the figure is saved to a \verb|*.png| file and cleared for the
783 following iteration.
784 \begin{python}
785 # ... 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 pl.title("Temperature across the blocks at time %e minutes"%(t/day))
795 #save figure to file
796 pl.savefig(os.path.join(save_path,"tempT","blockspyplot%03d.png") %i)
797 \end{python}
798 Some results are shown in \reffig{fig:onedheatout}.
800 \subsection{Making a Video}
801 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 is not available on every platform, so some users need to use an alternative
804 video encoder.
805 \begin{python}
806 # compile the *.png files to create a *.avi video that shows T change
807 # with time. This operation uses Linux mencoder. For other operating
808 # systems it is possible to use your favourite video compiler to
809 # convert image files to videos.
811 os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
812 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
813 example01tempT.avi")
814 \end{python}

  ViewVC Help
Powered by ViewVC 1.1.26