/[escript]/trunk/doc/cookbook/example01.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 2 weeks ago) by jfenwick
File MIME type: application/x-tex
File size: 37830 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26