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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26