ViewVC logotype

Contents of /trunk/doc/cookbook/example05.tex

Parent Directory Parent Directory | Revision Log Revision Log

Revision 3153 - (show annotations)
Sun Sep 5 01:31:49 2010 UTC (10 years, 2 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
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % Earth Systems Science Computational Center (ESSCC)
6 % http://www.uq.edu.au/esscc
7 %
8 % Primary Business: Queensland, Australia
9 % Licensed under the Open Software License version 3.0
10 % http://www.opensource.org/licenses/osl-3.0.php
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
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.
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}
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.
44 It is now possible to start defining our domain and boundaries.
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}
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.
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.
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
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.
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.
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
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.
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.
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:
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}.
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
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.
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