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