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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 13285 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26