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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2982 - (hide annotations)
Wed Mar 10 04:27:47 2010 UTC (9 years, 8 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 ahallam 2401
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3     %
4 jfenwick 2881 % Copyright (c) 2003-2010 by University of Queensland
5 ahallam 2401 % Earth Systems Science Computational Center (ESSCC)
6     % http://www.uq.edu.au/esscc
7     %
8     % Primary Business: Queensland, Australia
9     % Licensed under the Open Software License version 3.0
10     % http://www.opensource.org/licenses/osl-3.0.php
11     %
12     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13    
14 gross 2931 \begin{figure}[t]
15 ahallam 2606 \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}
16 caltinay 2982 \caption{Example 3: 2D model: granitic intrusion of sandstone country rock}
17 ahallam 2606 \label{fig:twodhdmodel}
18     \end{figure}
19    
20 gross 2950 \sslist{example03a.py and cblib.py}
21 gross 2931
22 ahallam 2979 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 caltinay 2982 Assuming that the cylinder is very long
28 ahallam 2979 we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement
29     the same
30 caltinay 2982 diffusion model as we have used for the granite blocks in \refSec{Sec:1DHDv00}
31 gross 2931 but will add the second spatial dimension and show how to define
32     variables depending on the location of the domain.
33 caltinay 2982 We use \file{onedheatdiff001b.py} as the starting point to develop this model.
34 gross 2931
35 ahallam 2979 \section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic
36     Intrusion}
37 gross 2931 \label{Sec:2DHD}
38    
39 ahallam 2979 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 caltinay 2982 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 ahallam 2801 \begin{python}
49 ahallam 2645 mx = 600*m #meters - model length
50     my = 600*m #meters - model width
51 gross 2931 ndx = 150 #mesh steps in x direction
52     ndy = 150 #mesh steps in y direction
53 ahallam 2645 r = 200*m #meters - radius of intrusion
54 gross 2931 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 ahallam 2801 \end{python}
57 gross 2931 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 ahallam 2979 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 caltinay 2982 need to deal with the second spatial dimension. We can investigate these two
69     modifications at the same time.
70 ahallam 2979 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 ahallam 2401 \begin{equation}
75 ahallam 2979 \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 gross 2931 \label{eqn:hd2}
79 ahallam 2401 \end{equation}
80 caltinay 2982 Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and
81 ahallam 2979 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 gross 2931 \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 ahallam 2979 \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 gross 2931 \label{eqn:hdgenf2}
95     \end{equation}
96 ahallam 2979 which is very similar to \refEq{eqn:hdgenf} used to define the temperature
97     update in the simple 1D case.
98 caltinay 2982 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 gross 2931 \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 ahallam 2979 $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 gross 2931 \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 ahallam 2979 The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set
116     to zero but $\kappa$,
117 caltinay 2982 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 gross 2931
123 caltinay 2982 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 ahallam 2979 \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The
126 gross 2931 \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 ahallam 2979 As the argument \verb|model| represents a two dimensional domain the matrix is
136 caltinay 2982 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 ahallam 2801 \begin{python}
140 gross 2931 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 caltinay 2982 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 gross 2931 \begin{equation}
148 gross 2948 -n \cdot \kappa \nabla T^{(n)} =
149 ahallam 2979 -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 gross 2948 \label{eq:hom flux}
152 gross 2931 \end{equation}
153 ahallam 2979 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 caltinay 2982 that the region is insulated which is what we want.
156 ahallam 2979 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 gross 2931
163 caltinay 2982 \section{Setting variable PDE Coefficients}
164 ahallam 2975 Now we need to look into the problem of how we define the material coefficients
165 artak 2969 $\kappa$ and $\rho c\hackscore p$ depending on their location in the domain.
166 ahallam 2979 We can make use of the technique used in the granite block example in
167     \refSec{Sec:1DHDv00}
168 artak 2969 to set up the initial temperature. However,
169 gross 2931 the situation is more complicated here as we have to deal with a
170 caltinay 2982 curved interface between the two sub-domains.
171 gross 2931
172 ahallam 2979 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 gross 2931 \begin{equation}
177     (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}
178     \end{equation}
179 ahallam 2979 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 caltinay 2982 intrusion points and $s-r > 0$ for all country rock points.
187 ahallam 2979 Defining these conditions within the script is quite simple and is done using
188     the following command:
189 gross 2931 \begin{python}
190 caltinay 2982 bound = length(x-ic)-r #where the boundary will be located
191 ahallam 2801 \end{python}
192 caltinay 2982 This definition of the boundary can now be used with the \verb|whereNegative|
193 ahallam 2979 command again to help define the material constants and temperatures in our
194     domain.
195 gross 2931 If \verb|kappai| and \verb|kappac| are the
196 ahallam 2979 thermal conductivities for the intrusion material granite and for the
197     surrounding sandstone, then we set;
198 ahallam 2801 \begin{python}
199 gross 2931 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 ahallam 2979 Notice that we are using the sample points of the \verb|Function| function space
205 caltinay 2982 as expected for the PDE coefficient \verb|A|\footnote{For the experienced user: use
206     \texttt{x=mypde.getFunctionSpace("A").getX()}.}.
207 gross 2931 The corresponding statements are used to set $\rho c\hackscore p$.
208 ahallam 2401
209 ahallam 2979 Our PDE has now been properly established. The initial conditions for
210     temperature are set out in a similar manner:
211 ahallam 2801 \begin{python}
212 ahallam 2658 #defining the initial temperatures.
213 gross 2931 x=Solution(model).getX()
214     bound = length(x-ic)-r
215     T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
216 ahallam 2801 \end{python}
217 caltinay 2982 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 artak 2969 notice that we reset \verb|x| and \verb|bound| to refer to the appropriate
220 ahallam 2979 sample points of a PDE solution\footnote{For the experienced user: use
221     \texttt{x=mypde.getFunctionSpace("r").getX()}.}.
222 ahallam 2606
223 gross 2950 \begin{figure}[ht]
224 gross 2931 \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 ahallam 2979 \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step
228     $1$, $20$ and $200$. Contour lines show temperature.}
229 gross 2931 \label{fig:twodhdans}
230     \end{figure}
231    
232 caltinay 2982 \section{Contouring \esc Data using \modmpl.}
233 gross 2948 \label{Sec:2DHD plot}
234 gross 2931 To complete our transition from a 1D to a 2D model we also need to modify the
235 caltinay 2982 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 ahallam 2979 used to represent the solution. Therefore it is necessary to interpolate our
239     solution onto a regular grid.
240 gross 2931
241 ahallam 2979 In \refSec{sec: plot T} we have already learned how to extract the $x$
242     coordinates of sample points as
243 caltinay 2982 \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 ahallam 2801 \begin{python}
246 gross 2931 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 caltinay 2982 For convenience we have put this function into \file{clib.py} so we can use
254     this function in other scripts more easily.
255 gross 2931
256     We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$
257 ahallam 2979 are the dimensions in the $x$ and $y$ directions) which is done using the
258 caltinay 2982 \modnumpy function \verb|linspace|.
259 gross 2931 \begin{python}
260     from clib import toXYTuple
261     # get sample points for temperature as for contouring
262     coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
263 ahallam 2606 # create regular grid
264 gross 2931 xi = np.linspace(0.0,mx,75)
265     yi = np.linspace(0.0,my,75)
266 ahallam 2801 \end{python}
267 gross 2931 The values \verb|[xi[k], yi[l]]| are the grid points.
268 ahallam 2606
269 caltinay 2982 The remainder of our contouring commands resides within a \verb|while| loop so
270 ahallam 2979 that a new contour is generated for each time step. For each time step the
271 caltinay 2982 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 ahallam 2979 rectangular grid defined by
275 caltinay 2982 \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 ahallam 2979 gradient of our plots is set via the command
278 caltinay 2982 \verb|pl.matplotlib.pyplot.autumn()|, other colours are listed on the \modmpl web page\footnote{see
279 ahallam 2979 \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured,
280 caltinay 2982 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 ahallam 2801 \begin{python}
283 ahallam 2645 #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 gross 2952 pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i)
295 ahallam 2645 pl.clf()
296 ahallam 2801 \end{python}
297 caltinay 2982 The function \verb|pl.contour| is used to add labelled contour lines to the
298 ahallam 2979 plot.
299 gross 2931 The results for selected time steps are shown in \reffig{fig:twodhdans}.
300 ahallam 2645
301 gross 2878
302 caltinay 2982 \section{Advanced Visualisation using VTK}
303 gross 2878
304 gross 2950 \sslist{example03b.py}
305 caltinay 2982 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 ahallam 2979 \url{http://code.enthought.com/projects/mayavi/}} as an example.
309 gross 2931
310     \mayavi is an interactive, GUI driven tool which is
311 caltinay 2982 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 ahallam 2979 problems.
314 caltinay 2982 The decision of which tool is the best can be subjective and users should
315 ahallam 2979 determine which package they require and are most comfortable with. The main
316 caltinay 2982 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 gross 2931
324 ahallam 2979 To write the temperatures at each time step to data files in the VTK file format
325 caltinay 2982 one needs to insert a \verb|saveVTK| call into the code;
326 gross 2878 \begin{python}
327 gross 2931 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 gross 2878 \end{python}
334 caltinay 2982 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 gross 2878
339 gross 2948 \begin{figure}[ht]
340 gross 2931 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}}
341 caltinay 2982 \caption{Example 3b: \mayavi start up window}
342 gross 2931 \label{fig:mayavi window}
343     \end{figure}
344    
345 gross 2948 \begin{figure}[ht]
346 gross 2931 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}}
347 caltinay 2982 \caption{Example 3b: \mayavi data control window}
348 gross 2931 \label{fig:mayavi window2}
349     \end{figure}
350 artak 2969 After you run the script you will find the
351 ahallam 2979 result files \file{data.*.vtu} in the result directory \file{data/example03}.
352 caltinay 2982 Run the command
353 gross 2931 \begin{python}
354     >> mayavi2 -d data.001.vtu -m Surface &
355     \end{python}
356 ahallam 2979 from the result directory. \mayavi will start up a window similar to
357     \reffig{fig:mayavi window}.
358 gross 2931 The right hand side shows the temperature at the first time step. To show
359 ahallam 2979 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 caltinay 2982 shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top
363 ahallam 2979 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 caltinay 2982 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 gross 2931
370 caltinay 2982 \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 ahallam 2979 For more details on writing scripts for parallel computing please consult the
375 caltinay 2982 \emph{user's guide}.
376    

  ViewVC Help
Powered by ViewVC 1.1.26