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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3153 - (hide annotations)
Sun Sep 5 01:31:49 2010 UTC (10 years, 3 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
1 ahallam 3153
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{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.
23    
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.
30    
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}
36    
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.
43    
44     It is now possible to start defining our domain and boundaries.
45    
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}
68     # TOP BLOCK
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.
113    
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.
126    
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
133    
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.
140    
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.
179    
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
198    
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.
207    
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.
228    
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:
236    
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}.
283    
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
290    
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.
297    
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