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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14 \begin{figure}[t]
15 \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}
16 \caption{Example 3: 2D model: granitic intrusion of sandstone country rock}
17 \label{fig:twodhdmodel}
18 \end{figure}
19
20 \sslist{example03a.py and cblib.py}
21
22 Building upon our success from the 1D models, it is now prudent to expand our
23 domain by another dimension.
24 For this example we use a very simple magmatic intrusion as the basis for our
25 model. The simulation will be a single event where some molten granite has
26 formed a cylindrical dome at the base of some cold sandstone country rock.
27 Assuming that the cylinder is very long
28 we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement
29 the same
30 diffusion model as we have used for the granite blocks in \refSec{Sec:1DHDv00}
31 but will add the second spatial dimension and show how to define
32 variables depending on the location of the domain.
33 We use \file{onedheatdiff001b.py} as the starting point to develop this model.
34
35 \section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic
36 Intrusion}
37 \label{Sec:2DHD}
38
39 To expand upon our 1D problem, the domain must first be expanded. In fact, we
40 have solved a two dimensional problem already but essentially ignored the second
41 dimension. In our definition phase we create a square domain in $x$ and
42 $y$\footnote{In \esc the notation
43 $x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.}
44 that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the
45 number of discrete spatial cells to $150$ in both directions and the radius of
46 the intrusion to $200$ meters with the centre located at the $300$ meter mark
47 on the $x$-axis. Thus, the domain variables are;
48 \begin{python}
49 mx = 600*m #meters - model length
50 my = 600*m #meters - model width
51 ndx = 150 #mesh steps in x direction
52 ndy = 150 #mesh steps in y direction
53 r = 200*m #meters - radius of intrusion
54 ic = [300*m, 0] #coordinates of the centre of intrusion (meters)
55 qH=0.*J/(sec*m**3) #our heat source temperature is zero
56 \end{python}
57 As before we use
58 \begin{python}
59 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
60 \end{python}
61 to generate the domain.
62
63 There are two fundamental changes to the PDE that we have discussed in
64 \refSec{Sec:1DHDv00}. Firstly,
65 because the material coefficients for granite and sandstone are different, we
66 need to deal with
67 PDE coefficients which vary with their location in the domain. Secondly, we
68 need to deal with the second spatial dimension. We can investigate these two
69 modifications at the same time.
70 In fact, the temperature model \refEq{eqn:hd} we utilised in
71 \refSec{Sec:1DHDv00} applied for the
72 1D case with a constant material parameter only. For the more general case
73 examined in this chapter, the correct model equation is
74 \begin{equation}
75 \rho c\hackscore p \frac{\partial T}{\partial t} - \frac{\partial }{\partial x}
76 \kappa \frac{\partial T}{\partial x} - \frac{\partial }{\partial y} \kappa
77 \frac{\partial T}{\partial y} = q\hackscore H
78 \label{eqn:hd2}
79 \end{equation}
80 Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and
81 for the case of constant material parameters $\frac{\partial }{\partial x}
82 \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides
83 with a simplified model equation for this case. It is more convenient
84 to write this equation using the $\nabla$ notation as we have already seen in
85 \refEq{eqn:commonform nabla};
86 \begin{equation}\label{eqn:Tform nabla}
87 \rho c\hackscore p \frac{\partial T}{\partial t}
88 -\nabla \cdot \kappa \nabla T = q\hackscore H
89 \end{equation}
90 We can easily apply the backward Euler scheme as before to end up with
91 \begin{equation}
92 \frac{\rho c\hackscore p}{h} T^{(n)} -\nabla \cdot \kappa \nabla T^{(n)} =
93 q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n-1)}
94 \label{eqn:hdgenf2}
95 \end{equation}
96 which is very similar to \refEq{eqn:hdgenf} used to define the temperature
97 update in the simple 1D case.
98 The difference is in the second order derivative term
99 $\nabla \cdot \kappa \nabla T^{(n)}$. Under the light of the more general case
100 we need to revisit the \esc PDE form as shown in \ref{eqn:commonform2D}.
101 For the 2D case with variable PDE coefficients the form needs to be read as
102 \begin{equation}\label{eqn:commonform2D 2}
103 -\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x}
104 -\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y}
105 -\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x}
106 -\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y}
107 + Du = f
108 \end{equation}
109 So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and
110 $f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}$ as we have used
111 before (see \refEq{ESCRIPT SET}) we need to set
112 \begin{equation}\label{eqn: kappa general}
113 A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0
114 \end{equation}
115 The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set
116 to zero but $\kappa$,
117 which brings in the second dimension. It is important to note that the
118 coefficients of the PDE may depend on their location in the domain which does
119 not influence the usage of the PDE form. This is very convenient as we can
120 introduce spatial dependence to the PDE coefficients without modification to
121 the way we call the PDE solver.
122
123 A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general}
124 can be carried out using the Kronecker $\delta$ symbol\footnote{see
125 \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The
126 \esc function \verb|kronecker| returns this matrix;
127 \begin{equation}
128 \verb|kronecker(model)| = \left[
129 \begin{array}{cc}
130 1 & 0 \\
131 0 & 1 \\
132 \end{array}
133 \right]
134 \end{equation}
135 As the argument \verb|model| represents a two dimensional domain the matrix is
136 returned as a $2 \times 2$ matrix
137 (in the case of a three-dimensional domain a $3 \times 3$ matrix is returned).
138 The call
139 \begin{python}
140 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
141 \end{python}
142 sets the PDE coefficients according to \refEq{eqn: kappa general}.
143
144 We need to check the boundary conditions before we turn to the question: how do
145 we set $\kappa$. As pointed out \refEq{NEUMAN 2} makes certain assumptions on
146 the boundary conditions. In our case these assumptions translate to;
147 \begin{equation}
148 -n \cdot \kappa \nabla T^{(n)} =
149 -n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} -
150 n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0
151 \label{eq:hom flux}
152 \end{equation}
153 which sets the normal component of the heat flux $- \kappa \cdot (\frac{\partial
154 T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means
155 that the region is insulated which is what we want.
156 On the left and right face of the domain where we have
157 $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$
158 this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on
159 the domain
160 where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is
161 $\frac{\partial T^{(n)}}{\partial y}=0$.
162
163 \section{Setting variable PDE Coefficients}
164 Now we need to look into the problem of how we define the material coefficients
165 $\kappa$ and $\rho c\hackscore p$ depending on their location in the domain.
166 We can make use of the technique used in the granite block example in
167 \refSec{Sec:1DHDv00}
168 to set up the initial temperature. However,
169 the situation is more complicated here as we have to deal with a
170 curved interface between the two sub-domains.
171
172 Prior to setting up the PDE, the interface between the two materials must be
173 established.
174 The distance $s\ge 0$ between two points $[x,y]$ and
175 $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as:
176 \begin{equation}
177 (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}
178 \end{equation}
179 If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes
180 the centre of the semi-circle of our intrusion, then for all the points $[x,y]$
181 in our model we can calculate a distance to $ic$.
182 All the points that fall within the given radius $r$ of our intrusion will have
183 a corresponding
184 value $s < r$ and all those belonging to the country rock will have a value $s >
185 r$. By subtracting $r$ from both of these conditions we find $s-r < 0$ for all
186 intrusion points and $s-r > 0$ for all country rock points.
187 Defining these conditions within the script is quite simple and is done using
188 the following command:
189 \begin{python}
190 bound = length(x-ic)-r #where the boundary will be located
191 \end{python}
192 This definition of the boundary can now be used with the \verb|whereNegative|
193 command again to help define the material constants and temperatures in our
194 domain.
195 If \verb|kappai| and \verb|kappac| are the
196 thermal conductivities for the intrusion material granite and for the
197 surrounding sandstone, then we set;
198 \begin{python}
199 x=Function(model).getX()
200 bound = length(x-ic)-r
201 kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound))
202 mypde.setValue(A=kappa*kronecker(model))
203 \end{python}
204 Notice that we are using the sample points of the \verb|Function| function space
205 as expected for the PDE coefficient \verb|A|\footnote{For the experienced user: use
206 \texttt{x=mypde.getFunctionSpace("A").getX()}.}.
207 The corresponding statements are used to set $\rho c\hackscore p$.
208
209 Our PDE has now been properly established. The initial conditions for
210 temperature are set out in a similar manner:
211 \begin{python}
212 #defining the initial temperatures.
213 x=Solution(model).getX()
214 bound = length(x-ic)-r
215 T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
216 \end{python}
217 where \verb|Ti| and \verb|Tc| are the initial temperature in the regions of the
218 granite and surrounding sandstone, respectively. It is important to
219 notice that we reset \verb|x| and \verb|bound| to refer to the appropriate
220 sample points of a PDE solution\footnote{For the experienced user: use
221 \texttt{x=mypde.getFunctionSpace("r").getX()}.}.
222
223 \begin{figure}[ht]
224 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}}
225 \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}}
226 \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}}
227 \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step
228 $1$, $20$ and $200$. Contour lines show temperature.}
229 \label{fig:twodhdans}
230 \end{figure}
231
232 \section{Contouring \esc Data using \modmpl.}
233 \label{Sec:2DHD plot}
234 To complete our transition from a 1D to a 2D model we also need to modify the
235 plotting procedure. As before we use \modmpl to do the plotting in this case a
236 contour plot. For 2D contour plots \modmpl expects that the data are regularly
237 gridded. We have no control over the location and ordering of the sample points
238 used to represent the solution. Therefore it is necessary to interpolate our
239 solution onto a regular grid.
240
241 In \refSec{sec: plot T} we have already learned how to extract the $x$
242 coordinates of sample points as
243 \verb|numpy| array to hand the values to \modmpl. This can easily be extended
244 to extract both the $x$ and the $y$ coordinates;
245 \begin{python}
246 import numpy as np
247 def toXYTuple(coords):
248 coords = np.array(coords.toListOfTuples()) #convert to Tuple
249 coordX = coords[:,0] #X components.
250 coordY = coords[:,1] #Y components.
251 return coordX,coordY
252 \end{python}
253 For convenience we have put this function into \file{clib.py} so we can use
254 this function in other scripts more easily.
255
256 We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$
257 are the dimensions in the $x$ and $y$ directions) which is done using the
258 \modnumpy function \verb|linspace|.
259 \begin{python}
260 from clib import toXYTuple
261 # get sample points for temperature as for contouring
262 coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
263 # create regular grid
264 xi = np.linspace(0.0,mx,75)
265 yi = np.linspace(0.0,my,75)
266 \end{python}
267 The values \verb|[xi[k], yi[l]]| are the grid points.
268
269 The remainder of our contouring commands resides within a \verb|while| loop so
270 that a new contour is generated for each time step. For each time step the
271 solution must be re-gridded for \modmpl using the \verb|griddata| function. This
272 function interprets irregularly located values \verb|tempT| at locations defined
273 by \verb|coordX| and \verb|coordY| as values at the new coordinates of a
274 rectangular grid defined by
275 \verb|xi| and \verb|yi|. The output is \verb|zi|. It is now possible to use the
276 \verb|contourf| function which generates colour filled contours. The colour
277 gradient of our plots is set via the command
278 \verb|pl.matplotlib.pyplot.autumn()|, other colours are listed on the \modmpl web page\footnote{see
279 \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured,
280 visually adjusted using the \modmpl functions and then saved to a file.
281 \verb|pl.clf()| clears the figure in readiness for the next time iteration.
282 \begin{python}
283 #grid the data.
284 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
285 # contour the gridded data, plotting dots at the randomly spaced data points.
286 pl.matplotlib.pyplot.autumn()
287 pl.contourf(xi,yi,zi,10)
288 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
289 pl.clabel(CS, inline=1, fontsize=8)
290 pl.axis([0,600,0,600])
291 pl.title("Heat diffusion from an intrusion.")
292 pl.xlabel("Horizontal Displacement (m)")
293 pl.ylabel("Depth (m)")
294 pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i)
295 pl.clf()
296 \end{python}
297 The function \verb|pl.contour| is used to add labelled contour lines to the
298 plot.
299 The results for selected time steps are shown in \reffig{fig:twodhdans}.
300
301
302 \section{Advanced Visualisation using VTK}
303
304 \sslist{example03b.py}
305 An alternative approach to \modmpl for visualisation is the usage of a package
306 which is based on the Visualization Toolkit (VTK) library\footnote{see \url{http://www.vtk.org/}}.
307 There is a variety of packages available. Here we use the package \mayavi\footnote{see
308 \url{http://code.enthought.com/projects/mayavi/}} as an example.
309
310 \mayavi is an interactive, GUI driven tool which is
311 really designed to visualise large three dimensional data sets where \modmpl
312 is not suitable. But it is also very useful when it comes to two dimensional
313 problems.
314 The decision of which tool is the best can be subjective and users should
315 determine which package they require and are most comfortable with. The main
316 difference between using \mayavi (or other VTK based tools) and \modmpl is that
317 the actual visualisation is detached from the calculation by writing the
318 results to external files and importing them into \mayavi. In 3D the best
319 camera position for rendering a scene is not obvious before the results are
320 available. Therefore the user may need to try different settings before the
321 best is found. Moreover, in many cases a 3D interactive visualisation is the
322 only way to really understand the results (e.g. using stereographic projection).
323
324 To write the temperatures at each time step to data files in the VTK file format
325 one needs to insert a \verb|saveVTK| call into the code;
326 \begin{python}
327 while t<=tend:
328 i+=1 #counter
329 t+=h #current time
330 mypde.setValue(Y=qH+T*rhocp/h)
331 T=mypde.getSolution()
332 saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T)
333 \end{python}
334 The data files, e.g. \file{data.001.vtu}, contain all necessary information to
335 visualise the temperature and can be directly processed by \mayavi. Note that
336 there is no re-gridding required. It is recommended to use the file extension
337 \file{.vtu} for files created by \verb|saveVTK|.
338
339 \begin{figure}[ht]
340 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}}
341 \caption{Example 3b: \mayavi start up window}
342 \label{fig:mayavi window}
343 \end{figure}
344
345 \begin{figure}[ht]
346 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}}
347 \caption{Example 3b: \mayavi data control window}
348 \label{fig:mayavi window2}
349 \end{figure}
350 After you run the script you will find the
351 result files \file{data.*.vtu} in the result directory \file{data/example03}.
352 Run the command
353 \begin{python}
354 >> mayavi2 -d data.001.vtu -m Surface &
355 \end{python}
356 from the result directory. \mayavi will start up a window similar to
357 \reffig{fig:mayavi window}.
358 The right hand side shows the temperature at the first time step. To show
359 the results at other time steps you can click at the item \texttt{VTK XML file
360 (data.001.vtu) (timeseries)}
361 at the top left hand side. This will bring up a new window similar to the window
362 shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top
363 right corner you can move forwards and backwards in time.
364 You will also notice the text \textbf{T} next to the item \texttt{Point scalars
365 name}. The
366 name \textbf{T} corresponds to the keyword argument name \texttt{T} we have
367 used in the \verb|saveVTK| call. In this menu item you can select other results
368 you may have written to the output file for visualisation.
369
370 \textbf{For the advanced user}: Using the \modmpl to visualise spatially
371 distributed data is not MPI compatible. However, the \verb|saveVTK| function
372 can be used with MPI. In fact, the function collects the values of the sample
373 points spread across processor ranks into a single file.
374 For more details on writing scripts for parallel computing please consult the
375 \emph{user's guide}.
376

  ViewVC Help
Powered by ViewVC 1.1.26