ViewVC logotype

Contents of /branches/stage3.1/doc/cookbook/steadystateheatdiff.tex

Parent Directory Parent Directory | Revision Log Revision Log

Revision 2945 - (show annotations)
Wed Feb 24 00:17:46 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 17865 byte(s)
Bringing release stage up to trunk version 2944

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{Steady-state Heat Refraction}
18 New concepts will include non-linear boundaries and the ability to prescribe location specific variables.
20 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 models. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modeling 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.
22 \subsection{Creating the model with \pycad}
23 \sslist{heatrefraction_mesher001.py and cblib.py}
25 Our first 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 dependant upon the material properties and shape.
27 \begin{figure}[h!]
28 \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}
29 \caption{Heat refraction model with point and line labels.}
30 \label{fig:anticlinehrmodel}
31 \end{figure}
33 We will start by defining our domain and material boundaries using the \pycad module. This example is located in the \esc directory under \exf and a labelled diagram is available in \reffig{fig:anticlinehrmodel}. The dependencies come first and for this example we have:
34 \begin{python}
35 from esys.pycad import * #domain constructor
36 from esys.pycad.gmsh import Design #Finite Element meshing package
37 from esys.finley import MakeDomain #Converter for escript
38 import os #file path tool
39 import numpy as np #numerial python for arrays
40 from math import * # math package
41 #used to construct polygons for plotting from pycad
42 from cblib import getLoopCoords
43 \end{python}
44 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.
46 It is now possible to start defining our domain and boundaries. There are a few primary constructors in \pycad that build upon each other to define domains and boundaries; the ones we will use are:
47 \begin{python}
48 Point() #Create a point in space.
49 Line() #Creates a line from a number of points.
50 CurveLoop() #Creates a closed loop from a number of lines.
51 Surface() #Creates a surface based on a CurveLoop
52 \end{python}
53 We start by inputting the variables we need to construct the model.
54 \begin{python}
55 #Model Parameters
56 width=5000.0 #width of model
57 depth=-6000.0 #depth of model
58 sspl=51 #number of discrete points in spline
59 dsp=width/(sspl-1) #dx of spline steps for width
60 dep_sp=2500.0 #avg depth of spline
61 h_sp=1500.0 #heigh of spline
62 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
63 \end{python}
64 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. Be careful to define all your constructs in an \textbf{anti-clockwise} manner otherwise the meshing algorithm may fail.
65 \begin{python}
66 # Overall Domain
67 p0=Point(0.0, 0.0, 0.0)
68 p1=Point(0.0, depth, 0.0)
69 p2=Point(width, depth, 0.0)
70 p3=Point(width, 0.0, 0.0)
71 \end{python}
72 Now lines are defined in an \textbf{anti-clockwise} direction using our points. This forms a rectangle around our domain.
73 \begin{python}
74 l01=Line(p0, p1)
75 l12=Line(p1, p2)
76 l23=Line(p2, p3)
77 l30=Line(p3, p0)
78 \end{python}
79 These lines form the basis for our domain boundary, which is a closed loop.
80 \begin{python}
81 c=CurveLoop(l01, l12, l23, l30)
82 \end{python}
83 The curve defining our clinal structure is located in approximately 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.
84 \begin{python}
85 # Material Boundary
86 x=[ Point(i*dsp\
87 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
88 for i in range(0,sspl)\
89 ]
90 mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple
91 \end{python}
92 The start and end points of the spline can be returned to help define the material boundaries.
93 \begin{python}
94 x1=Spline.getStartPoint(mysp)
95 x2=Spline.getEndPoint(mysp)
96 \end{python}
97 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.
98 \begin{python}
100 tbl1=Line(p0,x1)
101 tbl2=mysp
102 tbl3=Line(x2,p3)
103 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
104 tblock = PlaneSurface(tblockloop)
105 \end{python}
106 It is also necessary to create and export a polygon for \esc so that we can plot the boundaries of our sub-domains. First we take the loop from our block and retrieve its point coordinates with the function \verb getLoopCoords() and then output it with \modnumpy for our solution script.
107 \begin{python}
108 tpg = getLoopCoords(tblockloop)
109 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
110 \end{python}
111 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.
112 \begin{python}
113 bbl4=-mysp
114 \end{python}
115 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.
116 \begin{python}
117 #clockwise check
118 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
119 bpg = getLoopCoords(bblockloop2)
120 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
121 \end{python}
122 The last few steps in creating the model take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc.
123 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.
124 \begin{python}
125 # Create a Design which can make the mesh
126 d=Design(dim=2, element_size=200)
127 # Add the sub-domains and flux boundaries.
128 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
129 PropertySet("linebottom",l12))
130 # Create the geometry, mesh and \esc domain
131 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))
132 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))
133 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\
134 optimizeLabeling=True)
135 # Create a file that can be read back in to python with mesh=ReadMesh(fileName)
136 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
137 \end{python}
138 The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.
140 \subsection{Steady-state PDE solution}
141 \sslist{heatrefraction.py and cblib.py}
142 While a steady-state solution will not give an indication as to the quantitative heat flow in a model, it is useful because it allows us to visualise the direction of heat transport and the equilibrium state of the model. Also, a steady-state model only requires one iteration to find a solution. We know from \refCh{CHAP HEAT DIFF} that the full form of the heat diffusion PDE can be expressed as follows.
143 \begin{equation}
144 \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
145 \label{eqn:hd2}
146 \end{equation}
147 In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to
148 \begin{equation}
149 - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
150 \label{eqn:sshd}
151 \end{equation}
152 where we see our only variables are $\kappa$ the thermal conductivity and the heat source term $q\hackscore H$. Our boundary conditions are similar to earlier problems with some exceptions. Our heat source $q\hackscore H$ is now located at the base of the model and is implemented using our \modpycad identity \textit{linebottom} and we have a constant temperature along the surface of the model where $z=0$. The \modLPDE module's general form as a provision for these cases of the form
153 \begin{equation}
154 u=r \textit{ where } q>0
155 \label{eqn:bdr}
156 \end{equation}
157 using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true.
158 The structure of our solution script is similar to those we have used already. We begin by importing the dependencies and defining the PDE and script variables. In this case there are four major variables; \verb qin is the temperature of our source, \verb Ti is the temperature on the surface at $z=0$ and our width and depth of the model.
159 \begin{python}
161 qin=70*Milli*W/(m*m) #our heat source temperature is now zero
162 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
163 width=5000.0*m
164 depth=-6000.0*m
165 \end{python}
166 The mesh is now imported via the \verb ReadMesh() function and loaded into a domain variable \verb mymesh . This is followed by the structural polygons. 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.
167 \begin{python}
168 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
169 tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
170 tpgx = tpg[:,0]
171 tpgy = tpg[:,1]
172 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
173 bpgx = bpg[:,0]
174 bpgy = bpg[:,1]
176 # set up kappa (thermal conductivity across domain) using tags
177 kappa=Scalar(0,Function(mymesh))
178 kappa.setTaggedValue("top",2.0)
179 kappa.setTaggedValue("bottom",4.0)
180 \end{python}
181 Setting up the PDE is a little more complicated. The linear PDE \verb mypde is again assigned the domain \verb mymesh and \verb A is set. \verb qH is set on the boundary by creating \verb qH as a function space of type \verb FunctionOnBoundary which is in turn based on \verb mymesh . We can then set the tagged value \textit{linebottom} to the condition \verb qin . We also need to use our condition from \refEq{eqn:bdr}. By extracting \verb x we can set the domain at zero to 1 and everywhere else 0 via the \verb whereZero() function. This is put into \verb q and \verb r is set to \verb Ti our temperature at the surface. Note \verb y=qH . As there is no time derivative in the equation we do not have any iterative steps and there is one solution to the problem. This is found in the step \verb T=mypde.getSolution() .
182 \begin{python}
183 #... generate functionspace...
184 #... open PDE ...
185 mypde=LinearPDE(mymesh)
186 #define first coefficient
187 mypde.setValue(A=kappa*kronecker(mymesh))
189 # ... set initial temperature ....
190 x=mymesh.getX()
192 qH=Scalar(0,FunctionOnBoundary(mymesh))
193 qH.setTaggedValue("linebottom",qin)
194 mypde.setValue(q=whereZero(x[1]),r=Ti)
195 mypde.setValue(y=qH)
197 # get steady state solution.
198 T=mypde.getSolution()
199 \end{python}
200 The problem is now solved and plotting is required to visualise the data.
202 \subsection{Line profiles of 2D data}
203 It is possible to take slices or profile sections of 2d data using \mpl. To achieve this there are three major steps, firstly the solution must be smoothed across the domain using the \verb projector() function.
204 \begin{python}
205 proj=Projector(mymesh)
206 smthT=proj(T)
207 \end{python}
208 The data must then be moved to a regular grid.
209 \begin{python}
210 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
211 \end{python}
212 The profile can then be plotted by taking a slice of the data output \verb zi .
213 \begin{python}
214 cut=int(len(xi)/2)
215 pl.clf()
216 pl.plot(zi[:,cut],yi)
217 pl.title("Heat Refraction Temperature Depth Profile")
218 pl.xlabel("Temperature (K)")
219 pl.ylabel("Depth (m)")
220 \end{python}
221 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}.
222 \begin{figure}
223 \centering
224 \subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tdp.png}\label{fig:tdp}}
225 \subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tgdp.png}\label{fig:tgdp}}
226 \subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefraction001tcdp.png}\label{fig:tcdp}}
227 \subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001hf.png}\label{fig:hf}}
228 \caption{Depth profiles down centre of model.}
229 \label{figs:dps}
230 \end{figure}
235 \subsection{Quiver plots in \mpl}
236 To visualise this data we are going to do three things. Firstly we will generate the gradient of the solution and create some quivers or arrows that describe the gradient. Then we need to generate some temperature contours and thirdly, we will create some polygons from our meshing output files to show our material boundaries.
238 The first step is very easily provided by \esc, where we can use the \verb grad() function to generate our gradient. It is then necessary to mould the solution to the domain node points. We do this by first extracting the nodes to \verb oldspacecoords and transforming them via interpolation to suit our solution which is in the function space form of \verb T .
239 \begin{python}
240 # calculate gradient of solution for quiver plot
241 qu=-kappa*grad(T)
243 # rearrange mymesh to suit solution function space
244 oldspacecoords=mymesh.getX()
245 coords=Data(oldspacecoords, T.getFunctionSpace())
246 \end{python}
247 Now the number of quivers we want is specified as $20 \times 20 = 400$ and their locations calculated based on our coordinates in \verb qu . All our data is then transformed to tuples for \modmpl .
248 \begin{python}
249 quivshape = [20,20] #quivers x and quivers y
250 # function to calculate quiver locations
251 qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
253 kappaT = kappa.toListOfTuples(scalarastuple=False)
254 coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
255 coordKX, coordKY = toXYTuple(coordsK)
257 tempT = T.toListOfTuples(scalarastuple=False)
258 coordX, coordY = toXYTuple(coords)
259 \end{python}
260 In a similar fashion to the contouring, our data needs to be transferred to a regular grid to contour the temperature. Our polygons are handled by the \verb fill() function making sure they are ordered to the lowest level of the plot beneath all other data. Labelling is then added and finally the Quivers using the \verb quiver() function. More detail on each \modmpl function is available from the \modmpl website.
261 \begin{python}
262 xi = np.linspace(0.0,width,100)
263 yi = np.linspace(depth,0.0,100)
264 # grid the data.
265 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
266 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
267 # contour the gridded data, plotting dots at the randomly spaced data points.
268 pl.matplotlib.pyplot.autumn()
270 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
271 #~ CK = pl.contourf(xi,yi,ziK,2)
272 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
274 pl.clabel(CS, inline=1, fontsize=8)
275 pl.title("Heat Refraction across a clinal structure.")
276 pl.xlabel("Horizontal Displacement (m)")
277 pl.ylabel("Depth (m)")
278 #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
279 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
281 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
282 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
283 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
284 \end{python}
285 The data has now been plotted.
286 \begin{figure}[ht]
287 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}}
288 \caption{Heat refraction model with gradient indicated by quivers.}
289 \label{fig:hr001qumodel}
290 \end{figure}
292 \newpage
293 \subsection{Fault and Overburden Model}
294 A slightly more complicated model can be found in the examples \textit{heatrefraction_mesher002.py} and \textit{heatrefractoin002.py} where three blocks are used within the model.
295 \begin{figure}[ht]
296 \centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}}
297 \caption{Heat refraction model with three blocks and gradient quivers.}
298 \end{figure}

  ViewVC Help
Powered by ViewVC 1.1.26