proof reading and review of cookbook, minor changes only + updated depth profiles to cb from last week

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2009 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 \section{Steady-state heat refraction} 15 Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in escript to define properties like material constants and source locations. 16 17 \subsection{Creating the model with \pycad} 18 \sslist{heatrefraction_mesher001.py and cblib.py} 19 20 Our first heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependant on the material properties and shape. 21 22 \begin{figure}[h!] 23 \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 24 \caption{Heat refraction model with point and line labels.} 25 \label{fig:anticlinehrmodel} 26 \end{figure} 27 28 We will start by defining our domain and material boundaries using the \pycad module. This example is located in the \esc directory under \exf and a labelled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have: 29 \begin{verbatim} 30 from esys.pycad import * #domain constructor 31 from esys.pycad.gmsh import Design #Finite Element meshing package 32 from esys.finley import MakeDomain #Converter for escript 33 import os #file path tool 34 import numpy as np #numerial python for arrays 35 from math import * # math package 36 #used to construct polygons for plotting from pycad 37 from cblib import getLoopCoords 38 \end{verbatim} 39 There are two modes available in our example. When \verb modal=1 we indicate to the script that we wish to model an anticline. Otherwise when \verb modal=-1 we wish to modal a syncline. The modal operator simply changes the orientation of the boundary function so that it is either upwards or downwards curving. A \verb save_path has also been defined so that we can easily separate our data from other examples and our scripts. 40 41 It is now possible to start defining our domain and boundaries. There are a few primary constructors in \pycad that build upon each other to define domains and boundaries; the ones we will use are: 42 \begin{verbatim} 43 Point() #Create a point in space. 44 Line() #Creates a line from a number of points. 45 CurveLoop() #Creates a closed loop from a number of lines. 46 Surface() #Creates a surface based on a CurveLoop 47 \end{verbatim} 48 We start by inputting the variables we need to construct the model. 49 \begin{verbatim} 50 #Model Parameters 51 width=5000.0 #width of model 52 depth=-6000.0 #depth of model 53 sspl=51 #number of discrete points in spline 54 dsp=width/(sspl-1) #dx of spline steps for width 55 dep_sp=2500.0 #avg depth of spline 56 h_sp=1500.0 #heigh of spline 57 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down 58 \end{verbatim} 59 The variables are then used to construct the four corners of our domain, which from the origin has the dimensions of 5000 meters width and -6000 meters depth. This is done with the \verb Point() function which accepts x, y and z coordinates. Our domain is in two dimensions so z should always be zero. Be careful to define all your constructs in an \textbf{anti-clockwise} manner otherwise the meshing algorithm may fail. 60 \begin{verbatim} 61 # Overall Domain 62 p0=Point(0.0, 0.0, 0.0) 63 p1=Point(0.0, depth, 0.0) 64 p2=Point(width, depth, 0.0) 65 p3=Point(width, 0.0, 0.0) 66 \end{verbatim} 67 Now lines are defined in an \textbf{anti-clockwise} direction using our points. This forms a rectangle around our domain. 68 \begin{verbatim} 69 l01=Line(p0, p1) 70 l12=Line(p1, p2) 71 l23=Line(p2, p3) 72 l30=Line(p3, p0) 73 \end{verbatim} 74 These lines form the basis for our domain boundary, which is a closed loop. 75 \begin{verbatim} 76 c=CurveLoop(l01, l12, l23, l30) 77 \end{verbatim} 78 The curve defining our clinal structure is located in approximately the middle of the domain and has a cosinusoidal shape. We define the curve by generating points at discrete intervals; 51 in this case, and then create a smooth curve through the points using the \verb Spline() function. 79 \begin{verbatim} 80 # Material Boundary 81 x=[ Point(i*dsp\ 82 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 83 for i in range(0,sspl)\ 84 ] 85 mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 86 \end{verbatim} 87 The start and end points of the spline can be returned to help define the material boundaries. 88 \begin{verbatim} 89 x1=Spline.getStartPoint(mysp) 90 x2=Spline.getEndPoint(mysp) 91 \end{verbatim} 92 Our top block or material above the clinal/spline boundary is defined in a \textbf{anti-clockwise} manner by creating lines and then a closed loop. As we will be meshing the subdomain we also need to assign it a planar surface. 93 \begin{verbatim} 94 # TOP BLOCK 95 tbl1=Line(p0,x1) 96 tbl2=mysp 97 tbl3=Line(x2,p3) 98 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 99 tblock = PlaneSurface(tblockloop) 100 \end{verbatim} 101 It is also necessary to create and export a polygon for \esc so that we can plot the boundaries of our subdomains. First we take the loop from our block and retrieve its point coordinates with the function \verb getLoopCoords() and then output it with \modnumpy for our solution script. 102 \begin{verbatim} 103 tpg = getLoopCoords(tblockloop) 104 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 105 \end{verbatim} 106 We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 107 \begin{verbatim} 108 bbl4=-mysp 109 \end{verbatim} 110 This reverse spline option unfortunately does not work for the getLoopCoords command, however, the polygon tool will accept clock-wise oriented points so we can define a new curve. 111 \begin{verbatim} 112 #clockwise check 113 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 114 bpg = getLoopCoords(bblockloop2) 115 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 116 \end{verbatim} 117 The last few steps in creating our model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 118 To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements we wish to apply to the domain \verb element_size . It then becomes a simple task of adding the subdomains and flux boundaries to the design. If we wish to give the elements names as in the example we can use the \verb PropertySet() function. This will allow us to easily define the subdomain properties in the solution script. We then save the geometry, mesh and create our \esc domain. 119 \begin{verbatim} 120 # Create a Design which can make the mesh 121 d=Design(dim=2, element_size=200) 122 # Add the subdomains and flux boundaries. 123 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 124 PropertySet("linebottom",l12)) 125 # Create the geometry, mesh and Escript domain 126 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo")) 127 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 128 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\ 129 optimizeLabeling=True) 130 # Create a file that can be read back in to python with mesh=ReadMesh(fileName) 131 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 132 \end{verbatim} 133 The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem. 134 135 \subsection{Steady-state PDE solution} 136 \sslist{heatrefraction.py and cblib.py} 137 While a steady-state solution will not give an indication as to the quantitative heat flow in a model, it is useful because it allows us to visualise the direction of heat transport and the equilibrium state of the model. Also, a steady-state model only requires one iteration to find a solution. We know from \refCh{CHAP HEAT DIFF} that the full form of the heat diffusion PDE can be expressed as follows. 138 \begin{equation} 139 \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 140 \label{eqn:hd2} 141 \end{equation} 142 In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to 143 \begin{equation} 144 - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 145 \label{eqn:sshd} 146 \end{equation} 147 where we see our only variables are $\kappa$ the thermal conductivity and the heat source term $q\hackscore H$. Our boundary conditions are similar to earlier problems with some exceptions. Our heat source $q\hackscore H$ is now located at the base of the model and is implemented using our \modpycad identity \textit{linebottom} and we have a constant temperature along the surface of the model where $z=0$. The \modLPDE module's general form as a provision for these cases of the form 148 \begin{equation} 149 u=r \textit{ where } q>0 150 \label{eqn:bdr} 151 \end{equation} 152 using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true. 153 The structure of our solution script is similar to those we have used already. Firstly, the dependencies must be imported. Then the variables must be defined. In this case there are four major variables. \verb qin is the temperature of our source, \verb Ti is the temperature on the surface at $z=0$ and our width and depth of the model. 154 \begin{verbatim} 155 ##ESTABLISHING VARIABLES 156 qin=70*Milli*W/(m*m) #our heat source temperature is now zero 157 Ti=290.15*K # Kelvin #the starting temperature of our iron bar 158 width=5000.0*m 159 depth=-6000.0*m 160 \end{verbatim} 161 The mesh is now imported via the \verb ReadMesh() function and loaded into a domain variable \verb mymesh . The structural polygons are next. With the mesh imported it is now possible to use our tagging property to set up our PDE coefficients. In this case $\kappa$ is set via the \verb setTaggedValue() function which takes two arguments, the name of the tagged points and the value to assign to them. 162 \begin{verbatim} 163 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 164 tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 165 tpgx = tpg[:,0] 166 tpgy = tpg[:,1] 167 bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 168 bpgx = bpg[:,0] 169 bpgy = bpg[:,1] 170 171 # set up kappa (thermal conductivity across domain) using tags 172 kappa=Scalar(0,Function(mymesh)) 173 kappa.setTaggedValue("top",2.0) 174 kappa.setTaggedValue("bottom",4.0) 175 \end{verbatim} 176 Setting up the PDE is a little more complicated. The linear PDE \verb mypde is again assigned the domain \verb mymesh and \verb A is set. \verb qH is set on the boundary by creating \verb qH as a function space of type \verb FunctionOnBoundary which is in turn based on \verb mymesh . We can then set the tagged value \textit{linebottom} to the condition \verb qin . We also need to use our condition from \refEq{eqn:bdr}. By extracting \verb x we can set the domain at zero to 1 and everywhere else 0 via the \verb whereZero() function. This is put into \verb q and \verb r is set to \verb Ti our temperature at the surface. Note \verb y=qH . As there is no time derivative in the equation we do not have any iterative steps and there is one solution to the problem. This is found in the step \verb T=mypde.getSolution() . 177 \begin{verbatim} 178 #... generate functionspace... 179 #... open PDE ... 180 mypde=LinearPDE(mymesh) 181 #define first coefficient 182 mypde.setValue(A=kappa*kronecker(mymesh)) 183 184 # ... set initial temperature .... 185 x=mymesh.getX() 186 187 qH=Scalar(0,FunctionOnBoundary(mymesh)) 188 qH.setTaggedValue("linebottom",qin) 189 mypde.setValue(q=whereZero(x[1]),r=Ti) 190 mypde.setValue(y=qH) 191 192 # get steady state solution. 193 T=mypde.getSolution() 194 \end{verbatim} 195 The problem is now solved and plotting is required to visualise the data. 196 197 \subsection{Line profiles of 2D data} 198 It is possible to take slices or profile sections of 2d data using \mpl. To acheive this there are three major steps, firstly the solution must be smoothed across the domain using the \verb projector() function. 199 \begin{verbatim} 200 proj=Projector(mymesh) 201 smthT=proj(T) 202 \end{verbatim} 203 The data must then be moved to a regular grid. 204 \begin{verbatim} 205 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth) 206 \end{verbatim} 207 The profile can then be plotted by taking a slice of the data output \verb zi . 208 \begin{verbatim} 209 cut=int(len(xi)/2) 210 pl.clf() 211 pl.plot(zi[:,cut],yi) 212 pl.title("Heat Refraction Temperature Depth Profile") 213 pl.xlabel("Temperature (K)") 214 pl.ylabel("Depth (m)") 215 \end{verbatim} 216 This process can be repeated for various variations of the solultion. In this case we have temperature, temperature gradient, thermal conductivity and heat flow \reffig{figs:dps}. 217 \begin{figure} 218 \centering 219 \subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tdp.png}\label{fig:tdp}} 220 \subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tgdp.png}\label{fig:tgdp}} 221 \subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefraction001tcdp.png}\label{fig:tcdp}} 222 \subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001hf.png}\label{fig:hf}} 223 \caption{Depth profiles down centre of model.} 224 \label{figs:dps} 225 \end{figure} 226 227 228 229 230 \subsection{Quiver plots in \mpl} 231 To visualise this data we are going to do three things. Firstly we will generate the gradient of the solution and create some quivers or arrows that describe the gradient. Then we need to generate some temperatur contours and thirdly, we will create some polygons from our meshing output files to show our material boundaries. 232 233 The first step is very easily provided by \ESCRIPT, where we can use the \verb grad() function to generate our gradient. It is then necessary to mold the solution to the domain node points. We do this by first extracting the nodes to \verb oldspacecoords and transforming them via interpolation to suit our solution which is in the function space form of \verb T . 234 \begin{verbatim} 235 # calculate gradient of solution for quiver plot 236 qu=-kappa*grad(T) 237 238 # rearrage mymesh to suit solution function space 239 oldspacecoords=mymesh.getX() 240 coords=Data(oldspacecoords, T.getFunctionSpace()) 241 \end{verbatim} 242 Now the number of quivers we want is specified as $20 \times 20 = 400$ and their locations calculated based on our coordinates in \verb qu . All our date is then transformed to tuples for \modmpl . 243 \begin{verbatim} 244 quivshape = [20,20] #quivers x and quivers y 245 # function to calculate quiver locations 246 qu,qulocs = toQuivLocs(quivshape,width,depth,qu) 247 248 kappaT = kappa.toListOfTuples(scalarastuple=False) 249 coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 250 coordKX, coordKY = toXYTuple(coordsK) 251 252 tempT = T.toListOfTuples(scalarastuple=False) 253 coordX, coordY = toXYTuple(coords) 254 \end{verbatim} 255 In a similar fashion to the contouring, our data needs to be transferred to a regular grid to contour the temperature. Our polygons are handled by the \verb fill() function making sure they are ordered to the lowest level of the plot beneath all other data. Labeling is then added and finally the Quivers using the \verb quiver() function. More detail on each \modmpl function is available from the \modmpl website. 256 \begin{verbatim} 257 xi = np.linspace(0.0,width,100) 258 yi = np.linspace(depth,0.0,100) 259 # grid the data. 260 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 261 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 262 # contour the gridded data, plotting dots at the randomly spaced data points. 263 pl.matplotlib.pyplot.autumn() 264 265 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000) 266 #~ CK = pl.contourf(xi,yi,ziK,2) 267 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 268 269 pl.clabel(CS, inline=1, fontsize=8) 270 pl.title("Heat Refraction across a clinal structure.") 271 pl.xlabel("Horizontal Displacement (m)") 272 pl.ylabel("Depth (m)") 273 #~ CB = pl.colorbar(CS, shrink=0.8, extend='both') 274 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png")) 275 276 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white") 277 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 278 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png")) 279 \end{verbatim} 280 The data has now been plotted. 281 \begin{figure}[ht] 282 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}} 283 \caption{Heat refraction model with gradient indicated by quivers.} 284 \label{fig:hr001qumodel} 285 \end{figure} 286 287 \newpage 288 \subsection{Fault and Overburden Model} 289 A slightly more complicated model can be found in the examples \textit{heatrefraction_mesher002.py} and \textit{heatrefractoin002.py} where three blocks are used within the model. 290 \begin{figure}[ht] 291 \centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}} 292 \caption{Heat refraction model with three blocks and gradient quivers.} 293 \end{figure} 294