 More revisions to cookbook, fixed cblib.py for Scons, new figures.

 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 19 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 infact warped by the structure of the model and is heavily dependant on the material properties and shape. 20 21 \begin{figure}[h!] 22 \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 23 \caption{Heat refraction model with point and line labels.} 24 \label{fig:anticlinehrmodel} 25 \end{figure} 26 27 We will start by defining our domain and material boundaries using \pycad. This example is located in the \esc directory under \exf and a labeled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have: 28 \begin{verbatim} 29 from esys.pycad import * #domain constructor 30 from esys.pycad.gmsh import Design #Finite Element meshing package 31 from esys.finley import MakeDomain #Converter for escript 32 import os #file path tool 33 import numpy as np #numerial python for arrays 34 from math import * # math package 35 #used to construct polygons for plotting from pycad 36 from cblib import getLoopCoords 37 \end{verbatim} 38 There are two modes available in our example. When \verb modal=1 we indicate to the script that we wish to model an anticline. Othewise when \verb modal=-1 we wish to modal a syncline. The modal operator simply changes the phase 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. 39 40 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 one we will use are: 41 \begin{verbatim} 42 Point() #Create a point in space. 43 Line() #Creates a line from a number of points. 44 CurveLoop() #Creates a closed loop from a number of lines. 45 Surface() #Creates a surface based on a CurveLoop 46 \end{verbatim} 47 We start by inputting the variables we need to construct the model. 48 \begin{verbatim} 49 #Model Parameters 50 width=5000.0 #width of model 51 depth=-6000.0 #depth of model 52 sspl=51 #number of discrete points in spline 53 dsp=width/(sspl-1) #dx of spline steps for width 54 dep_sp=2500.0 #avg depth of spline 55 h_sp=1500.0 #heigh of spline 56 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down 57 \end{verbatim} 58 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. 59 \begin{verbatim} 60 # Overall Domain 61 p0=Point(0.0, 0.0, 0.0) 62 p1=Point(0.0, depth, 0.0) 63 p2=Point(width, depth, 0.0) 64 p3=Point(width, 0.0, 0.0) 65 \end{verbatim} 66 Now lines are defined in an \textbf{anti-clockwise} direction using our points. This forms a rectangle around our domain. 67 \begin{verbatim} 68 l01=Line(p0, p1) 69 l12=Line(p1, p2) 70 l23=Line(p2, p3) 71 l30=Line(p3, p0) 72 \end{verbatim} 73 These lines form the basis for our domain boundary, which is a closed loop. 74 \begin{verbatim} 75 c=CurveLoop(l01, l12, l23, l30) 76 \end{verbatim} 77 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. 78 \begin{verbatim} 79 # Material Boundary 80 x=[ Point(i*dsp\ 81 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 82 for i in range(0,sspl)\ 83 ] 84 mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 85 \end{verbatim} 86 The start and end points of the spline can be returned to help define the material boundaries. 87 \begin{verbatim} 88 x1=Spline.getStartPoint(mysp) 89 x2=Spline.getEndPoint(mysp) 90 \end{verbatim} 91 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. 92 \begin{verbatim} 93 # TOP BLOCK 94 tbl1=Line(p0,x1) 95 tbl2=mysp 96 tbl3=Line(x2,p3) 97 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 98 tblock = PlaneSurface(tblockloop) 99 \end{verbatim} 100 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 retreive its point coordinates with the function \verb getLoopCoords() and then output it with \modnumpy for our solution script. 101 \begin{verbatim} 102 tpg = getLoopCoords(tblockloop) 103 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 104 \end{verbatim} 105 We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is the fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 106 \begin{verbatim} 107 bbl4=-mysp 108 \end{verbatim} 109 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. 110 \begin{verbatim} 111 #clockwise check 112 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 113 bpg = getLoopCoords(bblockloop2) 114 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 115 \end{verbatim} 116 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. 117 To initial 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 simlpe 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 scirpt. We then save the geometry, mesh and create our \esc domain. 118 \begin{verbatim} 119 # Create a Design which can make the mesh 120 d=Design(dim=2, element_size=200) 121 # Add the subdomains and flux boundaries. 122 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),PropertySet("linebottom",l12)) 123 # Create the geometry, mesh and Escript domain 124 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo")) 125 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 126 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True) 127 # Create a file that can be read back in to python with mesh=ReadMesh(fileName) 128 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 129 \end{verbatim} 130 The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem. 131 132 \subsection{Steady-state PDE solution} 133 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. 134 \begin{equation} 135 \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 136 \label{eqn:hd2} 137 \end{equation} 138 In the steady state PDE however, there is no time dependance. This means that \refEq{eqn:hd2} reduces to 139 \begin{equation} 140 - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 141 \label{eqn:sshd} 142 \end{equation} 143 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 144 \begin{equation} 145 u=r \textit{ where } q>0 146 \label{eqn:bdr} 147 \end{equation} 148 using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true. 149 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. 150 \begin{verbatim} 151 ##ESTABLISHING VARIABLES 152 qin=70*Milli*W/(m*m) #our heat source temperature is now zero 153 Ti=290.15*K # Kelvin #the starting temperature of our iron bar 154 width=5000.0*m 155 depth=-6000.0*m 156 \end{verbatim} 157 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. 158 \begin{verbatim} 159 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 160 tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 161 tpgx = tpg[:,0] 162 tpgy = tpg[:,1] 163 bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 164 bpgx = bpg[:,0] 165 bpgy = bpg[:,1] 166 167 # set up kappa (thermal conductivity across domain) using tags 168 kappa=Scalar(0,Function(mymesh)) 169 kappa.setTaggedValue("top",2.0) 170 kappa.setTaggedValue("bottom",4.0) 171 \end{verbatim} 172 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() . 173 \begin{verbatim} 174 #... generate functionspace... 175 #... open PDE ... 176 mypde=LinearPDE(mymesh) 177 #define first coefficient 178 mypde.setValue(A=kappa*kronecker(mymesh)) 179 180 # ... set initial temperature .... 181 x=mymesh.getX() 182 183 qH=Scalar(0,FunctionOnBoundary(mymesh)) 184 qH.setTaggedValue("linebottom",qin) 185 mypde.setValue(q=whereZero(x),r=Ti) 186 mypde.setValue(y=qH)#,r=17*Celsius) 187 188 # get steady state solution and export to vtk. 189 T=mypde.getSolution() 190 \end{verbatim} 191 192 \subsection{Quiver plots in \mpl} 193 \begin{verbatim} 194 # rearrage mymesh to suit solution function space 195 oldspacecoords=mymesh.getX() 196 coords=Data(oldspacecoords, T.getFunctionSpace()) 197 198 # calculate gradient of solution for quiver plot 199 qu=-kappa*grad(T) 200 201 # create quiver locations 202 quivshape = [20,20] #quivers x and quivers y 203 numquiv = quivshape*quivshape # total number of quivers 204 maxx = 5000. # maximum x displacement 205 dx = maxx/quivshape+1. # quiver x spacing 206 maxy = -6000. # maximum y displacement 207 dy = maxy/quivshape+1. # quiver y spacing 208 qulocs=np.zeros([numquiv,2],float) # memory for quiver locations 209 # fill qulocs 210 for i in range(0,quivshape-1): 211 for j in range(0,quivshape-1): 212 qulocs[i*quivshape+j,:] = [dx+dx*i,dy+dy*j] 213 # retreive values for quivers direction and shape from qu 214 quL = Locator(qu.getFunctionSpace(),qulocs.tolist()) 215 qu = quL(qu) #array of dx,dy for quivers 216 qu = np.array(qu) #turn into a numpy array 217 qulocs = quL.getX() #returns actual locations from data 218 qulocs = np.array(qulocs) #turn into a numpy array 219 220 kappaT = kappa.toListOfTuples(scalarastuple=False) 221 coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 222 coordsK = np.array(coordsK.toListOfTuples()) 223 coordKX = coordsK[:,0] 224 coordKY = coordsK[:,1] 225 226 coords = np.array(coords.toListOfTuples()) 227 tempT = T.toListOfTuples(scalarastuple=False) 228 229 coordX = coords[:,0] 230 coordY = coords[:,1] 231 232 xi = np.linspace(0.0,5000.0,100) 233 yi = np.linspace(-6000.0,0.0,100) 234 # grid the data. 235 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 236 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 237 # contour the gridded data, plotting dots at the randomly spaced data points. 238 239 pl.matplotlib.pyplot.autumn() 240 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000) 241 #~ CK = pl.contourf(xi,yi,ziK,2) 242 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 243 pl.clabel(CS, inline=1, fontsize=8) 244 pl.title("Heat Refraction across a clinal structure.") 245 pl.xlabel("Horizontal Displacement (m)") 246 pl.ylabel("Depth (m)") 247 #~ CB = pl.colorbar(CS, shrink=0.8, extend='both') 248 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png")) 249 250 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white") 251 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 252 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png")) 253 \end{verbatim} 254 255 256 \begin{figure} 257 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}} 258 \caption{Heat refraction model with gradient indicated by quivers.} 259 \label{fig:hr001qumodel} 260 \end{figure} 261