# Annotation of /trunk/doc/cookbook/example05.tex

Revision 3153 - (hide annotations)
Sun Sep 5 01:31:49 2010 UTC (10 years, 3 months ago) by ahallam
File MIME type: application/x-tex
File size: 13104 byte(s)
Updates to cookbook. Fixes to SWP and new Gravitational Potential

 1 ahallam 3153 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 \section{Example 5: A Heat Refraction Model} 15 \label{example5} 16 \sslist{example05a.py and cblib.py} 17 Our heat refraction model will be a large anticlinal structure that is subject 18 to a constant temperature at the surface and experiencing a steady heat flux at 19 its base. Our aim is to show that the temperature flux across the surface is 20 not linear from bottom to top, but is in fact warped by the structure of the 21 model. The heat flow pattern demonstrates the dependence upon the material 22 properties and the shape of the interface. 23 24 The script of \refSec{example4} is modified by subdividing the block into two 25 parts. The curve separating the two blocks is given as a spline, see 26 \reffig{fig:anticlinehrmodel}. The data points 27 used to define the curve may be imported from a database of measurements 28 (\textit{e.g.} borehole depth data), but for simplicity it is assumed here that 29 the coordinates are known in an analytic form. 30 31 \begin{figure}[ht] 32 \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 33 \caption{Example 5a: Heat refraction model with point and line labels} 34 \label{fig:anticlinehrmodel} 35 \end{figure} 36 37 There are two modes available in this example. When \verb|modal=1|, this 38 indicates to the script that the model should be an anticline. Otherwise, when 39 \verb|modal=-1|, the model is a syncline. The modal operator simply changes the 40 orientation of the boundary function so that it is either upwards or downwards 41 curving. A \verb|save_path| has also been defined so that we can easily separate 42 our data from other examples and our scripts. 43 44 It is now possible to start defining our domain and boundaries. 45 46 The curve defining our clinal structure is located approximately in the middle 47 of the domain and has a sinusoidal shape. We define the curve by generating 48 points at discrete intervals; $51$ in this case, and then create a smooth curve 49 through the points using the \verb|Spline()| function. 50 \begin{python} 51 # Material Boundary 52 x=[ Point(i*dsp\ 53 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 54 for i in range(0,sspl)\ 55 ] 56 mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 57 \end{python} 58 The start and end points of the spline can be returned to help define the 59 material boundaries. 60 \begin{python} 61 x1=mysp.getStartPoint() 62 x2=mysp.getEndPoint() 63 \end{python} 64 The top block or material above the clinal/spline boundary is defined in an 65 \textbf{anti-clockwise} manner by creating lines and then a closed loop. By 66 meshing the sub-domain we also need to assign it a planar surface. 67 \begin{python} 68 # TOP BLOCK 69 tbl1=Line(p0,x1) 70 tbl2=mysp 71 tbl3=Line(x2,p3) 72 l30=Line(p3, p0) 73 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 74 tblock = PlaneSurface(tblockloop) 75 \end{python} 76 This process is repeated for every other sub-domain. In this example there is 77 only one other, the bottom block. The process is similar to the top block but 78 with a few differences. The spline points must be reversed by setting the spline 79 as negative. 80 \begin{python} 81 bbl4=-mysp 82 \end{python} 83 This reverse spline option unfortunately does not work for the 84 \verb|getLoopCoords| command, however, the \modmpl polygon tool will accept 85 clock-wise oriented points so we can define a new curve. 86 \begin{python} 87 #clockwise check 88 bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 89 \end{python} 90 The last few steps in creating the domain require that the previously defined 91 domain and sub-domain points are submitted to generate a mesh that can be 92 imported into \esc. 93 To initialise the mesh it first needs some design parameters. In this case we 94 have 2 dimensions \verb|dim| and a specified number of finite elements that need 95 to be applied to the domain \verb|element_size|. It then becomes a simple task 96 of adding the sub-domains and flux boundaries to the design. Each element of our 97 model can be given an identifier which makes it easier to define the sub-domain 98 properties in the solution script. This is done using the 99 \verb|PropertySet()| function. The geometry and mesh are then saved so the 100 \esc domain can be created. 101 \begin{python} 102 # Create a Design which can make the mesh 103 d=Design(dim=2, element_size=200) 104 # Add the sub-domains and flux boundaries. 105 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 106 PropertySet("linebottom",l12)) 107 # Create the geometry, mesh and escript domain 108 d.setScriptFileName(os.path.join(save_path,"example05.geo")) 109 d.setMeshFileName(os.path.join(save_path,"example05.msh")) 110 domain=MakeDomain(d,optimizeLabeling=True) 111 \end{python} 112 The creation of our domain and its mesh is now complete. 113 114 With the mesh imported it is now possible to use our tagging property to set up 115 our PDE coefficients. In this case $\kappa$ is set via the 116 \verb|setTaggedValue()| function which takes two arguments, the name of the 117 tagged points and the value to assign to them. 118 \begin{python} 119 # set up kappa (thermal conductivity across domain) using tags 120 kappa=Scalar(0,Function(domain)) 121 kappa.setTaggedValue("top",2.0*W/m/K) 122 kappa.setTaggedValue("bottom",4.0*W/m/K) 123 \end{python} 124 No further changes are required to the PDE solution step, see 125 \reffig{fig:anticlinetemp} for the result. 126 127 \begin{figure}[ht] 128 \centerline{\includegraphics[width=4.in]{figures/heatrefraction}} 129 \caption{Example 5a: Temperature Distribution in the Heat Refraction Model} 130 \label{fig:anticlinetemp} 131 \end{figure} 132 \clearpage 133 134 \section{Line Profiles of 2D Data} 135 \sslist{example05b.py and cblib.py} 136 We want to investigate the profile of the data of the last example. 137 Of particular interest is the depth profile of the heat flux which is the 138 second component of $-\kappa \nabla T$. The script from the previous section 139 is extended to show how a vertical profile can be plotted. 140 141 The first important piece of information, is that \esc assumes that $-\kappa 142 \nabla T$ is not smooth and that the point values of this solution are defined 143 at numerical interpolation points. This assumption is reasonable as 144 the flux is the product of the piecewise constant function $\kappa$ and 145 the gradient of the temperature $T$ which has a discontinuity at the rock 146 interface. 147 Before plotting this function we need to smooth the solution using the 148 \verb|Projector()| class; 149 \begin{python} 150 from esys.escript.pdetools import Projector 151 proj=Projector(domain) 152 qu=proj(-kappa*grad(T)) 153 \end{python} 154 The \verb|proj| object provides a mechanism to distribute values given at the 155 numerical interpolation points to the nodes of the FEM mesh - the heat flux in 156 this example. \verb|qu| has the same function space as the temperature 157 \verb|T|. The smoothed flux is interpolated to a regular $200\times 200$ grid 158 via; 159 \begin{python} 160 xiq,yiq,ziq = toRegGrid(qu[1],200,200) 161 \end{python} 162 using the \verb|toRegGrid| function from the cookbook library which we are 163 using for the contour plot. 164 At return \verb|ziq[j,i]| is the value of vertical heat flux at point 165 (\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by 166 plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script 167 creates a deep profile at $x_{0}=\frac{width}{2}$; 168 \begin{python} 169 cut=int(len(xiq)/2) 170 pl.plot(ziq[:,cut]*1000.,yiq) 171 pl.title("Vertical Heat Flow Depth Profile") 172 pl.xlabel("Heat Flow (mW/m^2)") 173 pl.ylabel("Depth (m)") 174 pl.savefig(os.path.join(save_path,"hf.png")) 175 \end{python} 176 This process can be repeated for various variations of the solution. 177 \reffig{figs:dps} shows temperature, temperature gradient, thermal conductivity 178 and heat flow. 179 180 \begin{figure}[htp] 181 \centering 182 \subfigure[Temperature Depth 183 Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{ 184 fig:tdp}} 185 \subfigure[Temperature Gradient Depth 186 Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{ 187 fig:tgdp}} 188 \subfigure[Thermal Conductivity 189 Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{ 190 fig:tcdp}} 191 \subfigure[Heat Flow Depth 192 Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf} 193 } 194 \caption{Example 5b: Depth profiles down centre of model} 195 \label{figs:dps} 196 \end{figure} 197 \clearpage 198 199 \section{Arrow Plots in \mpl} 200 \sslist{example05c.py and cblib.py} 201 The distribution of the flux $-\kappa \nabla T$ is now visualised over the 202 domain and the results are plotted in \reffig{fig:hr001qumodel}. 203 The plot puts together three components. A contour plot of the temperature, 204 a coloured representation of the two sub-domains where colour represents the 205 thermal conductivity in the particular region, and finally the arrows 206 representing the local direction of the steepest gradient of the flux. 207 208 Contours have already been discussed in \refSec{Sec:2DHD plot}. To show 209 sub-domains, we need to go back to \pycad data to get the points used to 210 describe the boundary of the sub-domains. We have created the \verb|CurveLoop| 211 class object \verb|tblockloop| to define the boundary of the upper sub-domain. 212 We use the \verb|getPolygon()| method of \verb|CurveLoop| to get 213 access to the \verb|Point|s used to define the boundary. The statement 214 \begin{python} 215 [ p.getCoordinates() for p in tblockloop.getPolygon() ] 216 \end{python} 217 creates a list of the node coordinates of all the points in question. In order 218 to simplify the selection of the $x$ and $y$ coordinates the list is converted 219 into \modnumpy array. To add the area coloured in brown to the plot we use; 220 \begin{python} 221 import pylab as pl 222 import numarray as np 223 tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ]) 224 pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000) 225 \end{python} 226 The same code is applied to \verb|bblockloop| to create the red area for this 227 sub-domain. 228 229 To plot vectors representing the flux orientation we use the 230 \verb|quiver| function in \pylab. The function places vectors at locations in 231 the domain. 232 For instance one can plot vectors at the locations of the sample points used by 233 \esc to represent the flux \verb|-kappa*grad(T)|. As a vector is plotted at 234 each sample point one typically ends up with too many vectors. So one needs to 235 select a subset of points as follows: 236 237 First we create a coarse grid of points on a rectangular mesh, e.g. $20 \times 238 20$ points. Here we choose a grid of points which are located at the centre of 239 a \verb|nx| $\times$ \verb|ny| grid; 240 \begin{python} 241 dx = width/nx # x spacing 242 dy = depth/ny # y spacing 243 grid = [ ] # the grid points 244 for j in xrange(0,ny-1): 245 for i in xrange(0,nx-1): 246 grid.append([dx/2+dx*i,dy/2+dy*j]) 247 \end{python} 248 With the \verb|Locator| function \esc provides a mechanism to identify sample 249 points that are closest to the grid points we have selected and to retrieve the 250 data at these points; 251 \begin{python} 252 from esys.escript.pdetools import Locator 253 flux=-kappa*grad(T) 254 fluxLoc = Locator(flux.getFunctionSpace(),grid) 255 subflux= fluxLoc(flux) 256 \end{python} 257 \verb|subflux| now contains a list of flux components at certain sample points. 258 To get the list of the sample point coordinates one can use the \verb|getX()| 259 method of the \verb|Locator|; 260 \begin{python} 261 subfluxloc = fluxLoc.getX() 262 \end{python} 263 To simplify the selection of $x$ and $y$ components it is convenient to 264 transform \verb|subflux| and \verb|subfluxloc| to \numpy arrays 265 \verb|xflux|, \verb|flux|. 266 This function is implemented in the \verb|subsample| function within the 267 \file{clib.py} file so we can use it in other examples. One can easily use this 268 function to create a vector plot of the flux; 269 \begin{python} 270 from cblib import subsample 271 xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20) 272 pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white") 273 \end{python} 274 Finally, we add a title and labels; 275 \begin{python} 276 pl.title("Heat Refraction across a clinal structure.") 277 pl.xlabel("Horizontal Displacement (m)") 278 pl.ylabel("Depth (m)") 279 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 280 pl.savefig(os.path.join(saved_path,"flux.png")) 281 \end{python} 282 to get the desired result, see \reffig{fig:hr001qumodel}. 283 284 \begin{figure}[ht] 285 \centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}} 286 \caption{Example 5c: Heat refraction model with gradient indicated by vectors} 287 \label{fig:hr001qumodel} 288 \end{figure} 289 \clearpage 290 291 \section{Example 6:Fault and Overburden Model} 292 \sslist{example06.py and cblib.py} 293 A slightly more complicated model can be found in the example file 294 \textit{heatrefraction2_solver.py} where three blocks are used within the 295 model, see~\reffig{fig:hr002qumodel}. It is left to the reader to work through 296 this example. 297 298 \begin{figure}[ht] 299 \centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}} 300 \caption{Example 6: Heat refraction model with three blocks and heat flux} 301 \label{fig:hr002qumodel} 302 \end{figure} 303 \clearpage

 ViewVC Help Powered by ViewVC 1.1.26