/[escript]/trunk/doc/cookbook/example05.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26