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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26