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

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


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