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

Revision 2975 - (show annotations)
Thu Mar 4 05:44:12 2010 UTC (10 years, 4 months ago) by ahallam
File MIME type: application/x-tex
File size: 17680 byte(s)
Cookbook modifications Review -> Not quite finished yet.
New figure.

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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 domain by another dimension. 23 For this example we use 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{Example 3: 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 essentially ignored the second dimension. In our definition phase we create a square domain in $x$ and $y$\footnote{In \esc the notation 34 $x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the number of discrete spatial cells to 150 in both direction and the radius of the intrusion to $200$ meters with the centre located at the $300$ meter mark on the $x$-axis. Thus, the domain variables are; 35 \begin{python} 36 mx = 600*m #meters - model length 37 my = 600*m #meters - model width 38 ndx = 150 #mesh steps in x direction 39 ndy = 150 #mesh steps in y direction 40 r = 200*m #meters - radius of intrusion 41 ic = [300*m, 0] #coordinates of the centre of intrusion (meters) 42 qH=0.*J/(sec*m**3) #our heat source temperature is zero 43 \end{python} 44 As before we use 45 \begin{python} 46 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 47 \end{python} 48 to generate the domain. 49 50 There are two fundamental changes to the PDE that we have discussed in \refSec{Sec:1DHDv00}. Firstly, 51 because the material coefficients for granite and sandstone are different, we need to deal with 52 PDE coefficients which vary with their location in the domain. Secondly, we need 53 to deal with the second spatial dimension. We can investigate these two modification at the same time. 54 In fact, the temperature model \refEq{eqn:hd} we utilised in \refSec{Sec:1DHDv00} applied for the 55 1D case with a constant material parameter only. For the more general case examined in this chapter, the correct model equation is 56 \begin{equation} 57 \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 58 \label{eqn:hd2} 59 \end{equation} 60 Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 61 for the case of constant material parameters $\frac{\partial }{\partial x} \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides with a simplified model equation for this case. It is more convenient 62 to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 63 \begin{equation}\label{eqn:Tform nabla} 64 \rho c\hackscore p \frac{\partial T}{\partial t} 65 -\nabla \cdot \kappa \nabla T = q\hackscore H 66 \end{equation} 67 We can easily apply the backward Euler scheme as before to end up with 68 \begin{equation} 69 \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)} 70 \label{eqn:hdgenf2} 71 \end{equation} 72 which is very similar to \refEq{eqn:hdgenf} used to define the temperature solution in the simple 1D case. 73 The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 74 the light of the more general case we need to revisit the \esc PDE form as 75 shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 76 the form needs to be read as 77 \begin{equation}\label{eqn:commonform2D 2} 78 -\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x} 79 -\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y} 80 -\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x} 81 -\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y} 82 + Du = f 83 \end{equation} 84 So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 85 $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 86 \begin{equation}\label{eqn: kappa general} 87 A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 88 \end{equation} 89 The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$, 90 which brings in the second dimension. Important to notice that the coefficients 91 of the PDE may depend on their location in the domain and 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. 92 93 A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can be carried out using the 94 Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 95 \esc function \verb|kronecker| returns this matrix; 96 \begin{equation} 97 \verb|kronecker(model)| = \left[ 98 \begin{array}{cc} 99 1 & 0 \\ 100 0 & 1 \\ 101 \end{array} 102 \right] 103 \end{equation} 104 As the argument \verb|model| represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix 105 (In case of a three-dimensional domain a $3 \times 3$ matrix is returned). The call 106 \begin{python} 107 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 108 \end{python} 109 sets the PDE coefficients according to \refEq{eqn: kappa general}. 110 111 We need to check the boundary conditions before we turn to the question: how we set $\kappa$. As 112 pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 113 this assumptions translates to; 114 \begin{equation} 115 -n \cdot \kappa \nabla T^{(n)} = 116 -n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} - n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 117 \label{eq:hom flux} 118 \end{equation} 119 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. 120 On the left and right face of the domain where we have $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 121 this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain 122 where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 123 124 \section{Setting Variable PDE Coefficients} 125 Now we need to look into the problem of how we define the material coefficients 126 $\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 127 We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 128 to set up the initial temperature. However, 129 the situation is more complicated here as we have to deal with a 130 curved interface between the two sub-domain. 131 132 Prior to setting up the PDE, the interface between the two materials must be established. 133 The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 134 \begin{equation} 135 (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2} 136 \end{equation} 137 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$. 138 All the points that fall within the given radius $r$ of our intrusion will have a corresponding 139 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$ 140 for all country rock points. 141 Defining these conditions within the script is quite simple and is done using the following command: 142 \begin{python} 143 bound = length(x-ic)-r #where the boundary will be located 144 \end{python} 145 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. 146 If \verb|kappai| and \verb|kappac| are the 147 thermal conductivities for the intrusion material granite and for the surrounding sandstone, then we set; 148 \begin{python} 149 x=Function(model).getX() 150 bound = length(x-ic)-r 151 kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound)) 152 mypde.setValue(A=kappa*kronecker(model)) 153 \end{python} 154 Notice that we are using the sample points of the \verb|Function| function space as expected for the 155 PDE coefficient \verb|A|\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 156 The corresponding statements are used to set $\rho c\hackscore p$. 157 158 Our PDE has now been properly established. The initial conditions for temperature are set out in a similar manner: 159 \begin{python} 160 #defining the initial temperatures. 161 x=Solution(model).getX() 162 bound = length(x-ic)-r 163 T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound)) 164 \end{python} 165 where \verb|Ti| and \verb|Tc| are the initial temperature 166 in the regions of the granite and surrounding sandstone, respectively. It is important to 167 notice that we reset \verb|x| and \verb|bound| to refer to the appropriate 168 sample points of a PDE solution\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 169 170 \begin{figure}[ht] 171 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 172 \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 173 \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 174 \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.} 175 \label{fig:twodhdans} 176 \end{figure} 177 178 \section{Contouring \esc data using \modmpl.} 179 \label{Sec:2DHD plot} 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 over 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 \begin{python} 190 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 the $x$ and $y$ directions) 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 # create regular grid 208 xi = np.linspace(0.0,mx,75) 209 yi = np.linspace(0.0,my,75) 210 \end{python} 211 The values \verb|[xi[k], yi[l]]| are the grid points. 212 213 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 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 \begin{python} 216 #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,"Tcontour%03d.png") %i) 228 pl.clf() 229 \end{python} 230 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 233 234 235 \section{Advanced Visualization using VTK} 236 237 \sslist{example03b.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 are variety 240 of packages available. Here we 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 of which tool is the best can be subjective and the user should determine which package they require and are most comfortable with. The main difference between using \mayavi (and other VTK based tools) 246 or \modmpl is that the actual visualization is detached from the 247 calculation by writing the results to external files 248 and importing them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 249 before the results are available. Therefore the user may need to try 250 different position before the best is found. Moreover, in many cases a 3D interactive 251 visualization is the only way to really understand the results (e.g. using stereographic projection). 252 253 To write the temperatures at each time step to data files in the VTK file format one 254 needs to insert a \verb|saveVTK| call into the code; 255 \begin{python} 256 while t<=tend: 257 i+=1 #counter 258 t+=h #current time 259 mypde.setValue(Y=qH+T*rhocp/h) 260 T=mypde.getSolution() 261 saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 262 \end{python} 263 The data files, eg. \file{data.001.vtu}, contains all necessary information to 264 visualize the temperature and can directly processed by \mayavi. Notice that there is no 265 regridding required. It is recommended to use the file extension \file{.vtu} for files 266 created by \verb|saveVTK|. 267 268 \begin{figure}[ht] 269 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}} 270 \caption{Example 3b: \mayavi start up Window.} 271 \label{fig:mayavi window} 272 \end{figure} 273 274 \begin{figure}[ht] 275 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}} 276 \caption{Example 3b: \mayavi data control window.} 277 \label{fig:mayavi window2} 278 \end{figure} 279 After you run the script you will find the 280 result files \file{data.*.vtu} in the result directory \file{data/example03}. Run the 281 command 282 \begin{python} 283 >> mayavi2 -d data.001.vtu -m Surface & 284 \end{python} 285 from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}. 286 The right hand side shows the temperature at the first time step. To show 287 the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)} 288 at the top left hand side. This will bring up a new window similar to the 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. 289 You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The 290 name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used 291 in the \verb|saveVTK| call. In this menu item you can select other results 292 you may have written to the output file for visualization. 293 294 \textbf{For the advanced user}: Using the \modmpl to visualize spatially distributed data 295 is not MPI compatible. However, the \verb|saveVTK| function can be used with MPI. In fact, 296 the function collects the values of the sample points spread across processor ranks into a single. 297 For more details on writing scripts for parallel computing please consult the \emph{user's guide}.