ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 3373 - (show annotations)
Tue Nov 23 00:29:07 2010 UTC (10 years, 2 months ago) by ahallam
File MIME type: application/x-tex
File size: 37654 byte(s)
New figures.
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 \section{Example 1: One Dimensional Heat Diffusion in Granite}
15 \label{Sec:1DHDv00}
17 The first model consists of two blocks of isotropic material, for instance
18 granite, sitting next to each other (\autoref{fig:onedgbmodel}).
19 Initial temperature in \textit{Block 1} is \verb|T1| and in \textit{Block 2} is
20 \verb|T2|.
21 We assume that the system is insulated.
22 What would happen to the temperature distribution in each block over time?
23 Intuition tells us that heat will be transported from the hotter block to the
24 cooler one until both
25 blocks have the same temperature.
27 \begin{figure}[ht]
28 \centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}}
29 \caption{Example 1: Temperature differential along a single interface between
30 two granite blocks.}
31 \label{fig:onedgbmodel}
32 \end{figure}
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_p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2}
43 T}{\partial x^{2}} = q_H
44 \label{eqn:hd}
45 \end{equation}
46 where $\rho$ is the material density, $c_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_{H}$; this can take the form of a constant or a function of time and
55 space. For example $q_{H} = q_{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_p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2}
103 T^{(n)}}{\partial x^{2}} = q_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_{jl} u_{,l})_{,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_p}{h} T^{(n)} - \kappa \frac{\partial^2 T^{(n)}}{\partial
160 x^2} = q_H + \frac{\rho c_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 _{p}}{h}; f = q _{H} + \frac{\rho
174 c_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_{,i} n_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_{j}A_{jl} u_{,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}[htp]
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_{00}\frac{\partial^{2}u}{\partial x^{2}}
396 -A_{01}\frac{\partial^{2}u}{\partial x\partial y}
397 -A_{10}\frac{\partial^{2}u}{\partial y\partial x}
398 -A_{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_{00}=A; A_{01}=A_{10}=A_{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=50 * yr # - 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_{jl}=A_{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_{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_{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_{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 \subsection{Plotting the Total Energy}
602 \sslist{example01b.py}
604 \esc does not include its own plotting capabilities. However, it is possible to
605 use a variety of free \pyt packages for visualisation.
606 Two types will be demonstrated in this cookbook;
607 \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and
608 \verb|VTK|\footnote{\url{http://www.vtk.org/}}.
609 The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}}
610 and is good for basic graphs and plots.
611 For more complex visualisation tasks, in particular two and three dimensional
612 problems we recommend the use of more advanced tools. For instance, \mayavi
613 \footnote{\url{http://code.enthought.com/projects/mayavi/}}
614 which is based upon the \verb|VTK| toolkit. The usage of \verb|VTK| based
615 visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two
616 dimensional PDE.
618 For our simple granite block problem, we have two plotting tasks. Firstly, we
619 are interested in showing the
620 behaviour of the total energy over time and secondly, how the temperature
621 distribution within the block is developing over time.
622 Let us start with the first task.
624 The idea is to create a record of the time marks and the corresponding total
625 energies observed.
626 \pyt provides the concept of lists for this. Before
627 the time loop is opened we create empty lists for the time marks \verb|t_list|
628 and the total energies \verb|E_list|.
629 After the new temperature has been calculated by solving the PDE we append the
630 new time marker and the total energy value for that time
631 to the corresponding list using the \verb|append| method. With these
632 modifications our script looks as follows:
633 \begin{python}
634 t_list=[]
635 E_list=[]
636 # ... start iteration:
637 while t<tend:
638 t+=h
639 mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients
640 T=mypde.getSolution() #get the PDE solution
641 totE=integrate(rhocp*T)
642 t_list.append(t) # add current time mark to record
643 E_list.append(totE) # add current total energy to record
644 \end{python}
645 To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs
646 to be loaded before use;
647 \begin{python}
648 import pylab as pl # plotting package.
649 \end{python}
650 Here we are not using \verb|from pylab import *| in order to avoid name
651 clashes for function names within \esc.
653 The following statements are added to the script after the time loop has been
654 completed;
655 \begin{python}
656 pl.plot(t_list,E_list)
657 pl.title("Total Energy")
658 pl.axis([0,max(t_list),0,max(E_list)*1.1])
659 pl.savefig("totE.png")
660 \end{python}
661 The first statement hands over the time marks and corresponding total energies
662 to the plotter.
663 The second statement sets the title for the plot. The third statement
664 sets the axis ranges. In most cases these are set appropriately by the plotter.
666 The last statement generates the plot and writes the result into the file
667 \verb|totE.png| which can be displayed by (almost) any image viewer.
668 As expected the total energy is constant over time, see
669 \reffig{fig:onedheatout1}.
671 \begin{figure}[ht]
672 \begin{center}
673 \includegraphics[width=4in]{figures/ttblockspyplot150}
674 \caption{Example 1b: Total Energy in the Blocks over Time (in seconds)}
675 \label{fig:onedheatout1}
676 \end{center}
677 \end{figure}
678 \clearpage
680 \subsection{Plotting the Temperature Distribution}
681 \label{sec: plot T}
682 \sslist{example01c.py}
683 For plotting the spatial distribution of the temperature we need to modify the
684 strategy we have used for the total energy.
685 Instead of producing a final plot at the end we will generate a
686 picture at each time step which can be browsed as a slide show or composed into
687 a movie.
688 The first problem we encounter is that if we produce an image at each time step
689 we need to make sure that the images previously generated are not overwritten.
691 To develop an incrementing file name we can use the following convention. It is
692 convenient to put all image files showing the same variable - in our case the
693 temperature distribution - into a separate directory.
694 As part of the \verb|os| module\footnote{The \texttt{os} module provides
695 a powerful interface to interact with the operating system, see
696 \url{http://docs.python.org/library/os.html}.} \pyt
697 provides the \verb|os.path.join| command to build file and directory names in a
698 platform independent way. Assuming that
699 \verb|save_path| is the name of the directory we want to put the results in the
700 command is;
701 \begin{python}
702 import os
703 os.path.join(save_path, "tempT%03d.png"%i )
704 \end{python}
705 where \verb|i| is the time step counter.
706 There are two arguments to the \verb|join| command. The \verb|save_path|
707 variable is a predefined string pointing to the directory we want to save our
708 data, for example a single sub-folder called \verb|data| would be defined by;
709 \begin{verbatim}
710 save_path = "data"
711 \end{verbatim}
712 while a sub-folder of \verb|data| called \verb|example01| would be defined by;
713 \begin{verbatim}
714 save_path = os.path.join("data","example01")
715 \end{verbatim}
716 The second argument of \verb|join| contains a string which is the file
717 name or subdirectory name. We can use the operator \verb|%| to use the value of
718 \verb|i| as part of our filename. The sub-string \verb|%03d| indicates that we
719 want to substitute a value into the name;
720 \begin{itemize}
721 \item \verb 0 means that small numbers should have leading zeroes;
722 \item \verb 3 means that numbers should be written using at least 3 digits;
723 and
724 \item \verb d means that the value to substitute will be a decimal integer.
725 \end{itemize}
727 To actually substitute the value of \verb|i| into the name write \verb|%i| after
728 the string.
729 When done correctly, the output files from this command will be placed in the
730 directory defined by \verb save_path as;
731 \begin{verbatim}
732 blockspyplot001.png
733 blockspyplot002.png
734 blockspyplot003.png
735 ...
736 \end{verbatim}
737 and so on.
739 A sub-folder check/constructor is available in \esc. The command;
740 \begin{verbatim}
741 mkDir(save_path)
742 \end{verbatim}
743 will check for the existence of \verb save_path and if missing, create the
744 required directories.
746 We start by modifying our solution script.
747 Prior to the \verb|while| loop we need to extract our finite solution
748 points to a data object that is compatible with \mpl. First we create the node
749 coordinates of the sample points used to represent
750 the temperature as a \pyt list of tuples or a \numpy array as requested by the
751 plotting function.
752 We need to convert the array \verb|x| previously set as
753 \verb|Solution(blocks).getX()| into a \pyt list
754 and then to a \numpy array. The $x_{0}$ component is then extracted via
755 an array slice to the variable \verb|plx|;
756 \begin{python}
757 import numpy as np # array package.
758 #convert solution points for plotting
759 plx = x.toListOfTuples()
760 plx = np.array(plx) # convert to tuple to numpy array
761 plx = plx[:,0] # extract x locations
762 \end{python}
764 \begin{figure}
765 \begin{center}
766 \includegraphics[width=4in]{figures/blockspyplot001}
767 \includegraphics[width=4in]{figures/blockspyplot050}
768 \includegraphics[width=4in]{figures/blockspyplot200}
769 \caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps
770 $1$, $50$ and $200$}
771 \label{fig:onedheatout}
772 \end{center}
773 \end{figure}
774 \clearpage
776 We use the same techniques provided by \mpl as we have used to plot the total
777 energy over time.
778 For each time step we generate a plot of the temperature distribution and save
779 each to a file.
780 The following is appended to the end of the \verb|while| loop and creates one
781 figure of the temperature distribution. We start by converting the solution to a
782 tuple and then plotting this against our \textit{x coordinates} \verb|plx| we
783 have generated before. We add a title to the diagram before it is rendered into
784 a file.
785 Finally, the figure is saved to a \verb|*.png| file and cleared for the
786 following iteration.
787 \begin{python}
788 # ... start iteration:
789 while t<tend:
790 i+=1
791 t+=h
792 mypde.setValue(Y=qH+rhocp/h*T)
793 T=mypde.getSolution()
794 totE=integrate(rhocp*T)
795 print "time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE)
796 t_list.append(t)
797 E_list.append(totE)
799 #establish figure 1 for temperature vs x plots
800 tempT = T.toListOfTuples()
801 pl.figure(1) #current figure
802 pl.plot(plx,tempT) #plot solution
803 # add title
804 pl.axis([0,mx,T1*.9,T2*1.1])
805 pl.title("Temperature across blocks at time %d days"%(t/day))
806 #save figure to file
807 pl.savefig(os.path.join(save_path,"tempT", "blockspyplot%03d.png"%i))
808 pl.clf() #clear figure
809 \end{python}
810 Some results are shown in \reffig{fig:onedheatout}.
812 \subsection{Making a Video}
813 Our saved plots from the previous section can be cast into a video using the
814 following command appended to the end of the script. The \verb mencoder command
815 is not available on every platform, so some users need to use an alternative
816 video encoder.
817 \begin{python}
818 # compile the *.png files to create a *.avi video that shows T change
819 # with time. This operation uses Linux mencoder. For other operating
820 # systems it is possible to use your favourite video compiler to
821 # convert image files to videos.
823 os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
824 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
825 example01tempT.avi")
826 \end{python}

  ViewVC Help
Powered by ViewVC 1.1.26