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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2972 - (show annotations)
Thu Mar 4 01:33:38 2010 UTC (10 years, 10 months ago) by artak
File MIME type: application/x-tex
File size: 24137 byte(s)
DONE. Einstain Summation Convention (ESC, what an irony, reminds me ESC button) is removed from cookbook as it is not used anywhere in the cookbook.
1
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14 \section{Steady-state Heat Refraction}
15 \label{STEADY-STATE HEAT REFRACTION}
16
17 In this chapter we show how to handle more complex geometries.
18
19 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 generate finite element meshes 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 \esc to define properties like material constants and source locations.
20
21 We proceed in this chapter by first looking at a very simple geometry. In fact, we solve
22 the steady-state heat equation over a rectangular domain. This is not a very interesting problem but then we extend our script by introducing an internal, curved interface.
23
24 \begin{figure}[ht]
25 \centerline{\includegraphics[width=4.in]{figures/pycadrec}}
26 \caption{Example 4: Rectangular Domain for \pycad.}
27 \label{fig:pycad rec}
28 \end{figure}
29
30 \section{Example 4: Creating the domain with \pycad}
31 \sslist{example04a.py}
32
33 We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look the steady state
34 case with slightly modified boundary conditions and use a more flexible tool
35 to generate to generate the geometry. Lets look at the geometry first.
36
37 We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and
38 insulation to both sides we want to calculate the steady-state temperature distribution.
39
40 In \pycad there are a few primary constructors that build upon each other to define domains and boundaries;
41 the ones we use are:
42 \begin{python}
43 from esys.pycad import *
44 Point() #Create a point in space.
45 Line() #Creates a line from a number of points.
46 CurveLoop() #Creates a closed loop from a number of lines.
47 PlaneSurface() #Creates a surface based on a CurveLoop
48 \end{python}
49 So to build-up our domain as shown in \reffig{fig:pycad rec} we first need to create
50 the corner points. From the corner points we build the four edges of the rectangle. The four edges
51 form then a closed loop which defines our domain as a surface.
52 We start by inputting the variables we need to construct the model.
53 \begin{python}
54 width=5000.0*m #width of model
55 depth=-6000.0*m #depth of model
56 \end{python}
57 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.
58 \begin{python}
59 # Overall Domain
60 p0=Point(0.0, 0.0, 0.0)
61 p1=Point(0.0, depth, 0.0)
62 p2=Point(width, depth, 0.0)
63 p3=Point(width, 0.0, 0.0)
64 \end{python}
65 Now lines are defined using our points. This forms a rectangle around our domain;
66 \begin{python}
67 l01=Line(p0, p1)
68 l12=Line(p1, p2)
69 l23=Line(p2, p3)
70 l30=Line(p3, p0)
71 \end{python}
72 Notice that lines have a direction. These lines form the basis for our domain boundary, which is a closed loop.
73 \begin{python}
74 c=CurveLoop(l01, l12, l23, l30)
75 \end{python}
76 Be careful to define the curved loop in an \textbf{anti-clockwise} manner otherwise the meshing algorithm may fail.
77 Finally we can define the domain as
78 \begin{python}
79 rec = PlaneSurface(c)
80 \end{python}
81 At this point the introduction of the curved loop seems to be unnecessary but this concept plays an important role
82 if holes are introduced.
83
84 Now we are ready to handover the domain \verb|rec| to a mesher which subdivides the domain into triangles (or tetrahedron in 3D). In our case we use \gmsh. We create
85 an instance of the \verb|Design| class which will handle the interface to \gmsh:
86 \begin{python}
87 from esys.pycad.gmsh import Design
88 d=Design(dim=2, element_size=200*m)
89 \end{python}
90 The argument \verb|dim| defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verb|element_size| defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be chosen with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m
91 and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So we end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily.
92 We can easily add our domain \verb|rec| to the \verb|Design|;
93 \begin{python}
94 d.addItem(rec)
95 \end{python}
96 We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this
97 but \pycad offers a more convenient technique called tagging. With this technique items in the domain are
98 named using the \verb|PropertySet| class. We can then later use this name to set values secifically for
99 those sample points located on the named items. Here we name the bottom face of the
100 domain where we will set the heat influx:
101 \begin{python}
102 ps=PropertySet("linebottom",l12))
103 d.addItem(ps)
104 \end{python}
105 Now we are ready to hand over the \verb|Design| to \FINLEY:
106 \begin{python}
107 from esys.finley import MakeDomain
108 domain=MakeDomain(d)
109 \end{python}
110 The \verb|domain| object can now be used in the same way like the return object of the \verb|Rectangle|
111 object we have used previously to generate a mesh. It is common practice to separate the
112 mesh generation from the PDE solution. The main reason for this is that mesh generation can be computationally very expensive in particular in 3D. So it is more efficient to generate the mesh once and write it to a file. The mesh
113 can then be read in every time a new simulation is run. \FINLEY supports this in the following
114 way~\footnote{An alternative are the \texttt{dump} and \texttt{load} functions. They using a binary format and tend to be much smaller.}
115 \begin{python}
116 # write domain to an text file
117 domain.write("example04.fly")
118 \end{python}
119 and then for reading in another script;
120 \begin{python}
121 # read domain from text file
122 from esys.finley import ReadMesh
123 domain =ReadMesh("example04.fly")
124 \end{python}
125
126 \begin{figure}[ht]
127 \centerline{\includegraphics[width=4.in]{figures/simplemesh}}
128 \caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}.}
129 \label{fig:pycad rec mesh}
130 \end{figure}
131
132 Before we discuss how we solve the PDE for the
133 problem it is useful to present two additional options of the \verb|Design|.
134 They allow the user accessing the script which is used by \gmsh to generate the mesh as well as
135 the mesh as it has been generated by \gmsh. This is done by setting specific names for these files:
136 \begin{python}
137 d.setScriptFileName("example04.geo")
138 d.setMeshFileName("example04.msh")
139 \end{python}
140 Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and
141 the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage.
142 Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using
143 \begin{verbatim}
144 gmsh example04.geo # show geometry
145 gmsh example04.msh # show mesh
146 \end{verbatim}
147 The mesh is shown in \reffig{fig:pycad rec mesh}.
148 \begin{figure}[ht]
149 \centerline{\includegraphics[width=4.in]{figures/simpleheat}}
150 \caption{Example 4b: Result of simple steady state heat problem.}
151 \label{fig:steady state heat}
152 \end{figure}
153
154
155 \section{The Steady-state Heat Equation}
156 \sslist{example04b.py, cblib}
157
158 Steady state is reached in the temperature when it is not changing in time. So to calculate the
159 the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero;
160 \begin{equation}\label{eqn:Tform nabla steady}
161 -\nabla \cdot \kappa \nabla T = q\hackscore H
162 \end{equation}
163 This PDE is easier to solve, than the PDE in \refEq{eqn:hdgenf2} in each time step. We can just drop
164 the \verb|D| term;
165 \begin{python}
166 mypde=LinearPDE(domain)
167 mypde.setValue(A=kappa*kronecker(model), Y=qH)
168 \end{python}
169 The temperature at the top face of the domain is known as \verb|Ttop| (=$20 C$). In \refSec{Sec:1DHDv0} we have
170 already discussed how this constraint is added to the PDE:
171 \begin{python}
172 x=Solution(domain).getX()
173 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
174 \end{python}
175 Notice that we use the \verb|sup| function to calculate the maximum of $y$ coordinates of the relevant sample points.
176
177 In all cases so far we have assumed that the domain is insulated which translates
178 into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling
179 set-up of this chapter we want to set the normal heat flux at the bottom to \verb|qin| while we still
180 maintain insulation at the left and right face. Mathematically we can express this as the format
181 \begin{equation}
182 -n \cdot \kappa \nabla T = q\hackscore{S}
183 \label{eq:inhom flux}
184 \end{equation}
185 where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero
186 for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face.
187 Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here.
188 We could define $q\hackscore{S}$ by using the masking techniques demonstrated earlier. The tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise
189 constant functions such as $q\hackscore{S}$. You need to remember now that we
190 marked the bottom face with the name \verb|linebottom| when we defined the domain.
191 We can use this now to create $q\hackscore{S}$;
192 \begin{python}
193 qS=Scalar(0,FunctionOnBoundary(domain))
194 qS.setTaggedValue("linebottom",qin)
195 \end{python}
196 In the first line \verb|qS| is defined as a scalar value over the sample points on the boundary of the domain. It is
197 initialized to zero for all sample points. In the second statement the values for those sample points which
198 on the line marked by \verb|linebottom| are set to \verb|qin|.
199
200 The Neuman boundary condition assumed by \esc has in fact the form
201 \begin{equation}\label{NEUMAN 2b}
202 n\cdot A \cdot\nabla u = y
203 \end{equation}
204 In comparison to the version in \refEq{NEUMAN 2} we have used so far the right hand side is now
205 the new PDE coefficient $y$. As we have not specific $y$ in our previous examples \esc has assumed
206 the value zero for $y$. A comparison of \refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one need to
207 choose $y=-q\hackscore{S}$;
208 \begin{python}
209 qS=Scalar(0,FunctionOnBoundary(domain))
210 qS.setTaggedValue("linebottom",qin)
211 mypde.setValue(y=-qS)
212 \end{python}
213 To plot the results we are using the \modmpl library as shown \refSec{Sec:2DHD plot}. For convenience
214 the interpolation of the temperature to a rectangular grid for contour plotting is made available
215 via the \verb|toRegGrid| function in the \verb|cblib| module. Your result should look similar to
216 \reffig{fig:steady state heat}.
217
218 \begin{figure}[ht]
219 \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}
220 \caption{Example 5a: Heat refraction model with point and line labels.}
221 \label{fig:anticlinehrmodel}
222 \end{figure}
223
224 \begin{figure}[ht]
225 \centerline{\includegraphics[width=4.in]{figures/heatrefraction}}
226 \caption{Example 5a: Temperature Distribution in the Heat Refraction Model.}
227 \label{fig:anticlinetemp}
228 \end{figure}
229
230 \section{Example 5: A Heat Refraction model}
231 \sslist{example05a.py and cblib.py}
232
233 Our 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 in fact warped by the structure of the model and is heavily dependent upon the material properties and shape.
234
235
236 We modify the example of the previous section by subdividing the block into two parts. The curve
237 separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points
238 used to define the curve may be measurement but for simplicity we assume here that there coordinates are
239 known in analytic form.
240
241 There are two modes available in our example. When \verb modal=1 this indicates to the script that we wish to model an anticline. Otherwise when \verb modal=-1 this will model a syncline. The modal operator simply changes the orientation 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.
242
243 It is now possible to start defining our domain and boundaries.
244
245 The curve defining our clinal structure is located approximately in the middle of the domain and has a sinusoidal 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.
246 \begin{python}
247 # Material Boundary
248 x=[ Point(i*dsp\
249 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
250 for i in range(0,sspl)\
251 ]
252 mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple
253 \end{python}
254 The start and end points of the spline can be returned to help define the material boundaries.
255 \begin{python}
256 x1=mysp.getStartPoint()
257 x2=mysp.getEndPoint()
258 \end{python}
259 The 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 sub-domain we also need to assign it a planar surface.
260 \begin{python}
261 # TOP BLOCK
262 tbl1=Line(p0,x1)
263 tbl2=mysp
264 tbl3=Line(x2,p3)
265 l30=Line(p3, p0)
266 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
267 tblock = PlaneSurface(tblockloop)
268 \end{python}
269 We must repeat the above for every other sub-domain. In this example there is only one other, the bottom block. The process is fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative.
270 \begin{python}
271 bbl4=-mysp
272 \end{python}
273 This reverse spline option unfortunately does not work for the getLoopCoords command, however, the \modmpl polygon tool will accept clock-wise oriented points so we can define a new curve.
274 \begin{python}
275 #clockwise check
276 bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
277 \end{python}
278 The last few steps in creating the model are: take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc.
279 To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements that need to be applied to the domain \verb element_size . It then becomes a simple task of adding the sub-domains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the sub-domain properties in the solution script. This is done using the \verb PropertySet() function. The geometry and mesh are then saved so the \esc domain can be created.
280 \begin{python}
281 # Create a Design which can make the mesh
282 d=Design(dim=2, element_size=200)
283 # Add the sub-domains and flux boundaries.
284 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
285 PropertySet("linebottom",l12))
286 # Create the geometry, mesh and \esc domain
287 d.setScriptFileName(os.path.join(save_path,"example05.geo"))
288 d.setMeshFileName(os.path.join(save_path,"example05.msh"))
289 domain=MakeDomain(d,optimizeLabeling=True)
290 \end{python}
291 The creation of our domain and its mesh is complete.
292
293 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.
294 \begin{python}
295 # set up kappa (thermal conductivity across domain) using tags
296 kappa=Scalar(0,Function(domain))
297 kappa.setTaggedValue("top",2.0*W/m/K)
298 kappa.setTaggedValue("bottom",4.0*W/m/K)
299 \end{python}
300 No further changes are required to the PDE solution step, see \reffig{fig:anticlinetemp} for the result.
301
302 \begin{figure}
303 \centering
304 \subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{fig:tdp}}
305 \subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{fig:tgdp}}
306 \subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{fig:tcdp}}
307 \subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf}}
308 \caption{Example 5b: Depth profiles down centre of model.}
309 \label{figs:dps}
310 \end{figure}
311
312 \section{Line profiles of 2D data}
313 \sslist{example05b.py and cblib.py}
314
315 We want to investigate the profile of the data of the last example.
316 We are particularly interested in the depth profile of the heat flux which is
317 the second component of $-\kappa \nabla T$. We extend the script developed in the
318 previous section to show how for instance vertical profile can be plotted.
319
320 The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and
321 will in fact represent the values at numerical interpolation points. This assumption is reasonable as
322 the flux is the product of the piecewise constant function $\kappa$ and
323 the gradient of the temperature $T$ which has a kink across the rock interface.
324 Before plotting this function we need to smooth it using the
325 \verb|Projector()| class;
326 \begin{python}
327 from esys.escript.pdetools import Projector
328 proj=Projector(domain)
329 qu=proj(-kappa*grad(T))
330 \end{python}
331 The \verb|proj| object provides a mechanism to distribute values given at the numerical interpolation points
332 - in this case the heat flux - to the nodes of the FEM mesh. \verb|qu| has the same attached function space
333 like the temperature \verb|T|. The smoothed flux is interpolated
334 to a regular $200\times 200$ grid;
335 \begin{python}
336 xiq,yiq,ziq = toRegGrid(qu[1],200,200)
337 \end{python}
338 using the \verb|toRegGrid| function from the cookbook library which we are using for the contour plot.
339 At return \verb|ziq[j,i]| is the value of vertical heat flux at point
340 (\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by
341 plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script
342 creates a deep profile at $x_{0}=\frac{width}{2}$;
343 \begin{python}
344 cut=int(len(xiq)/2)
345 pl.plot(ziq[:,cut]*1000.,yiq)
346 pl.title("Vertical Heat Flow Depth Profile")
347 pl.xlabel("Heat Flow (mW/m^2)")
348 pl.ylabel("Depth (m)")
349 pl.savefig(os.path.join(save_path,"hf.png"))
350 \end{python}
351 This process can be repeated for various variations of the solution. In this case we have temperature, temperature gradient, thermal conductivity and heat flow \reffig{figs:dps}.
352
353 \begin{figure}[ht]
354 \centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}}
355 \caption{Example 5c: Heat refraction model with gradient indicated by vectors.}
356 \label{fig:hr001qumodel}
357 \end{figure}
358
359 \section{Arrow plots in \mpl}
360 \sslist{example05c.py and cblib.py}
361 We would like to visualise the distribution of the flux $-\kappa \nabla T$ over the domain
362 and produce a plot like shown in \reffig{fig:hr001qumodel}.
363 The plot puts together three components. A contour plot of the temperature,
364 a colored representation of the two sub-domains where colour represents the thermal conductivity
365 in the particular region and finally the arrows representing the direction of flux.
366
367 Contours have already been discussed in \refSec{Sec:2DHD plot}. To show sub-domains,
368 we need to go back to \pycad data to get the points used to describe the boundary of the
369 sub-domains. We have created the \verb|CurveLoop| class object
370 \verb|tblockloop| to define the boundary of the upper sub-domain.
371 We use the \verb|getPolygon()| method of \verb|CurveLoop| to get
372 access to the \verb|Point|s used top define the boundary. The statement
373 \begin{python}
374 [ p.getCoordinates() for p in tblockloop.getPolygon() ])
375 \end{python}
376 creates a list of the node coordinates of all the points in question. In order
377 to simplify the selection of the $x$ and $y$ coordinates the list is converted
378 into \modnumpy array. To add the area colored in brown to the plot we use;
379 \begin{python}
380 import pylab as pl
381 import numarray as np
382 tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
383 pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
384 \end{python}
385 We can apply the same code to \verb|bblockloop| to a red area for this sub-domain to the block.
386
387 To plot vectors representing the flux orientation we use the
388 \verb|quiver| function in \pylab. The function places vectors at locations in the domain.
389 For instance one can plot vectors at the locations of the sample points used by \esc
390 to represent the flux \verb|-kappa*grad(T)|. As a vector is plotted at each sample point one typically ends
391 up with two many vectors. So one needs to select a subset of points:
392 First we create a coarse grid of points on a rectangular mesh, e.g. $20 \times 20$ points. Here we choose a grid of points which are located at the center of a \verb|nx| $\times$ \verb|ny| grid;
393 \begin{python}
394 dx = width/nx # x spacing
395 dy = depth/ny # y spacing
396 grid = [ ] # the grid points
397 for j in xrange(0,ny-1):
398 for i in xrange(0,nx-1):
399 grid.append([dx/2+dx*i,dy/2+dy*j])
400 \end{python}
401 With the \verb|Locator| \esc provides a mechanism to identify sample points that are closest
402 to the the grid points we have selected and to get the data at these data points;
403 \begin{python}
404 from esys.escript.pdetools import Locator
405 flux=-kappa*grad(T)
406 fluxLoc = Locator(flux.getFunctionSpace(),grid)
407 subflux= fluxLoc(flux)
408 \end{python}
409 \verb|subflux| gives now a list of flux component at certain sample points. To get the
410 list of the sample point coordinates one can use the \verb|getX()| method of the
411 \verb|Locator|;
412 \begin{python}
413 subfluxloc = fluxLoc.getX()
414 \end{python}
415 To simplify the selection of $x$ and $y$ components it is convenient
416 to transform \verb|subflux| and \verb|subfluxloc| to \numpy arrays
417 \verb|xflux|, \verb|flux|.
418 This function is implemented in the \verb|subsample|
419 in the \file{clib.py} file so we can use it in other examples. One can easily use this function
420 to create a vector plot of the flux;
421 \begin{python}
422 from cblib import subsample
423 xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
424 pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
425 \end{python}
426 We add title and labels;
427 \begin{python}
428 pl.title("Heat Refraction across a clinal structure.")
429 pl.xlabel("Horizontal Displacement (m)")
430 pl.ylabel("Depth (m)")
431 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
432 pl.savefig(os.path.join(saved_path,"flux.png"))
433 \end{python}
434 to get the desired result, see \reffig{fig:hr001qumodel}.
435
436 \begin{figure}[ht]
437 \centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}}
438 \caption{Example 6: Heat refraction model with three blocks and heat flux.}
439 \label{fig:hr002qumodel}
440 \end{figure}
441
442 \section{Example 6:Fault and Overburden Model}
443 \sslist{example06.py and cblib.py}
444 A slightly more complicated model can be found in the examples \textit{heatrefraction2_solver.py} where three blocks are used within the model, see~\reffig{fig:hr002qumodel}. It is left to the reader to work through this example.
445
446

  ViewVC Help
Powered by ViewVC 1.1.26