 # Diff of /trunk/doc/cookbook/steadystateheatdiff.tex

revision 2947 by gross, Mon Feb 8 06:45:41 2010 UTC revision 2948 by gross, Thu Feb 25 04:54:30 2010 UTC
# Line 14  Line 14
14  \section{Steady-state Heat Refraction}  \section{Steady-state Heat Refraction}
15  \label{STEADY-STATE HEAT REFRACTION}  \label{STEADY-STATE HEAT REFRACTION}
16
17    In this chapter we will learn how to handle more complex geometries.
18
19  New concepts will include non-linear boundaries and the ability to prescribe location specific variables.  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 generate finite element meshes 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.
20
21  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.  We will introduce in two stages. First we will look at a very simple geometry. In fact,  we will solve
22    the steady-state heat equation over a rectangular domain. This is not a very interesting problem but in a second
23    step we will extend our script by introducing an internal, curved interface.
24
25  \subsection{Creating the model with \pycad}  \begin{figure}[ht]
26  \sslist{heatrefraction_mesher001.py and cblib.py}  \centerline{\includegraphics[width=4.in]{figures/pycadrec}}
27    \caption{Rectangular Domain for \pycad.}
28    \label{fig:pycad rec}
29    \end{figure}
30
31  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.  \section{Creating the model with \pycad}
32    \sslist{steadyheat_mesher001.py}
33
34  \begin{figure}[h!]  We will modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: firstly we will look the steady state
35  \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}  case with slightly modified boundary conditions and secondly we will use a more flexible tool
36  \caption{Heat refraction model with point and line labels.}  to generate to generate the geometry. Lets look at the geometry first.
\label{fig:anticlinehrmodel}
\end{figure}
37
38  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:  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
39  \begin{python}  insulation to both sides we want to calculate the steady-state temperature distribution.
from esys.pycad import * #domain constructor
from esys.pycad.gmsh import Design #Finite Element meshing package
from esys.finley import MakeDomain #Converter for escript
import os #file path tool
import numpy as np #numerial python for arrays
from math import * # math package
#used to construct polygons for plotting from pycad
from cblib import getLoopCoords
\end{python}
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.
40
41  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:  In \pycad there are a few primary constructors that build upon each other to define domains and boundaries;
42    the ones we will use are:
43  \begin{python}  \begin{python}
44    from esys.pycad import *
45  Point() #Create a point in space.  Point() #Create a point in space.
46  Line() #Creates a line from a number of points.  Line() #Creates a line from a number of points.
47  CurveLoop() #Creates a closed loop from a number of lines.  CurveLoop() #Creates a closed loop from a number of lines.
48  Surface() #Creates a surface based on a CurveLoop  PlaneSurface() #Creates a surface based on a CurveLoop
49  \end{python}  \end{python}
50    So to build-up our domain as shown in \reffig{fig:pycad rec} we first need to create
51    the corner points. From the corner points we build the four edges of the rectangle. The four edges
52    form then a closed loop which defines our domain as a surface.
53  We start by inputting the variables we need to construct the model.  We start by inputting the variables we need to construct the model.
54  \begin{python}  \begin{python}
55  #Model Parameters  width=5000.0*m   #width of model
56  width=5000.0   #width of model  depth=-6000.0*m  #depth of model
depth=-6000.0  #depth of model
sspl=51 #number of discrete points in spline
dsp=width/(sspl-1) #dx of spline steps for width
dep_sp=2500.0 #avg depth of spline
h_sp=1500.0 #heigh of spline
orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
57  \end{python}  \end{python}
58  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.  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.
59  \begin{python}  \begin{python}
60  # Overall Domain  # Overall Domain
61  p0=Point(0.0,      0.0, 0.0)  p0=Point(0.0,      0.0, 0.0)
# Line 69  p1=Point(0.0,    depth, 0.0) Line 63  p1=Point(0.0,    depth, 0.0)
63  p2=Point(width, depth, 0.0)  p2=Point(width, depth, 0.0)
64  p3=Point(width,   0.0, 0.0)  p3=Point(width,   0.0, 0.0)
65  \end{python}  \end{python}
66  Now lines are defined in an \textbf{anti-clockwise} direction using our points. This forms a rectangle around our domain.  Now lines are defined using our points. This forms a rectangle around our domain;
67  \begin{python}  \begin{python}
68  l01=Line(p0, p1)  l01=Line(p0, p1)
69  l12=Line(p1, p2)  l12=Line(p1, p2)
70  l23=Line(p2, p3)  l23=Line(p2, p3)
71  l30=Line(p3, p0)  l30=Line(p3, p0)
72  \end{python}  \end{python}
73  These lines form the basis for our domain boundary, which is a closed loop.  Notice that lines have a direction. These lines form the basis for our domain boundary, which is a closed loop.
74  \begin{python}  \begin{python}
75  c=CurveLoop(l01, l12, l23, l30)  c=CurveLoop(l01, l12, l23, l30)
76  \end{python}  \end{python}
77  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.  Be careful to define the curved loop in an \textbf{anti-clockwise} manner otherwise the meshing algorithm may fail.
78    Finally we can define the domain as
79    \begin{python}
80    rec = PlaneSurface(c)
81    \end{python}
82    At this point the introduction of the curved loop seems to be unnecessary but this concept plays an important role
83    if holes are introduced.
84
85    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
86    an instance of the \verb|Design| class which will handle the interface to \gmsh:
87    \begin{python}
88    from esys.pycad.gmsh import Design
89    d=Design(dim=2, element_size=200*m)
90    \end{python}
91    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 choose 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
92    and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So will end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily.
93    We can easily add our domain \verb|rec| to the \verb|Design|;
94    \begin{python}
95    d.addItem(rec)
96    \end{python}
97    We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this
98    but \pycad offers a more convenient technique called tagging. With this techniques items in the domain are
99    named using the \verb|PropertySet| class. We can then later use this name to set values secifically for
100    those sample points located on the named items. Here we name the bottom face of the
101    domain where we will set the heat influx:
102    \begin{python}
103    ps=PropertySet("linebottom",l12))
104    d.addItem(ps)
105    \end{python}
106    Now we are ready to hand over the \verb|Design| to \FINLEY:
107    \begin{python}
108    from esys.finley import MakeDomain
109    domain=MakeDomain(d)
110    \end{python}
111    The \verb|domain| object can now be used in the same way like the return object of the \verb|Rectangle|
112    object we have used previously to generate a mesh. It is common practice to separate the
113    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
114    can then be read in every time a new simulation is run. \FINLEY supports this in the following
115    way~\footnote{An alternative are the \texttt{dump} and \texttt{load} functions. They using a binary format and tend to be much smaller.}
116    \begin{python}
117    # write domain to an text file
118    domain.write("simple.fly")
119    \end{python}
120    and then for recovery in another script;
121    \begin{python}
122    # read domain from text file
123    from esys.finley import ReadMesh
124    domain =ReadMesh("simple.fly")
125    \end{python}
126
127    \begin{figure}[ht]
128    \centerline{\includegraphics[width=4.in]{figures/simplemesh}}
129    \caption{Mesh over rectangular domain, see \reffig{fig:pycad rec}.}
130    \label{fig:pycad rec mesh}
131    \end{figure}
132
133    Before we discuss how we solve the PDE for the
134    problem it is useful to present two additional options of the \verb|Design|.
135    They allow the user accessing the script which is used to by \gmsh to generate the mesh as well as
136    the mesh as it has been generated by \gmsh. This is done by setting specific names for these files:
137    \begin{python}
138    d.setScriptFileName("simple.geo")
139    d.setMeshFileName("simple.msh")
140    \end{python}
141    Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and
142    the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage.
143    Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesher can be visualise with \gmsh from the command line using
144    \begin{verbatim}
145    gmsh simple.geo  # show geometry
146    gmsh simple.msh  # show mesh
147    \end{verbatim}
148    The mesh is shown in \reffig{fig:pycad rec mesh}.
149    \begin{figure}[ht]
150    \centerline{\includegraphics[width=4.in]{figures/simpleheat}}
151    \caption{REsult of simple steady state heat problem.}
152    \label{fig:steady state heat}
153    \end{figure}
154
155
156    \section{The Steady-state Heat Equation}
157    \sslist{simple_solver001.py, cblib}
158
159    Steady state is reached in the temperature is not changing in time. So to calculate the
160    the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero;
161    \begin{equation}\label{eqn:Tform nabla steady}
162    -\nabla \cdot \kappa \nabla T = q\hackscore H
163    \end{equation}
164    This PDE is easier then the PDE in \refEq{eqn:hdgenf2} we needed to solve in each time step. We can just drop
165    the \verb|D| term;
166    \begin{python}
167    mypde=LinearPDE(domain)
168    mypde.setValue(A=kappa*kronecker(model), Y=qH)
169    \end{python}
170    The temperature at the top face of the domain is known as \verb|Ttop| (=$20 C$). In \refSec{Sec:1DHDv0} we have
171    already discussed how this constraint is added to the PDE:
172    \begin{python}
173    x=Solution(domain).getX()
174    mypde.setValue(q=whereZero(x-sup(x)),r=Ttop)
175    \end{python}
176    Notice that we have use the \verb|sup| function to calculate the maximum $y$ coordinator of the relevant sample points.
177
178    In all cases so far we have assumed that the domain is insulated which translates
179    into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling
180    set-up of this chapter we want to set the normal heat flux at the bottom to \verb|qin| while we still
181    maintain insulation at the left and right face. Mathematically we can express this as the format
182    \begin{equation}
183    -n \cdot \kappa \nabla T = q\hackscore{S}
184    \label{eq:inhom flux}
185    \end{equation}
186    where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero
187    for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face.
188    Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here.
189    One could use the masking techniques we have already extensively used to define $q\hackscore{S}$ in \esc
190    but the tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise
191    constant functions such as $q\hackscore{S}$. You need to remember now that we
192    marked the bottom face with the name \verb|linebottom| when we defined the domain.
193    We can use this now to create $q\hackscore{S}$;
194    \begin{python}
195    qS=Scalar(0,FunctionOnBoundary(domain))
196    qS.setTaggedValue("linebottom",qin)
197    \end{python}
198    In the first line \verb|qS| is defined as a scalar value over the sample points on the boundary of the domain. It is
199    initialized to zero for all sample points. In the second statement the values for those sample points which
200    on the line marked by \verb|linebottom| are set to \verb|qin|.
201
202    The Neuman boundary condition assumed by \esc has in fact the form
203    \begin{equation}\label{NEUMAN 2b}
204    n\cdot A \cdot\nabla u = y
205    \end{equation}
206    In comparison to the version in \refEq{NEUMAN 2} we have used so far the right hand side is now
207    the new PDE coefficient $y$. As we have not specific $y$ in our previous examples \esc has assumed
208    the value zero for $y$. A comparison of \refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one need to
209    choose $y=-q\hackscore{S}$;
210    \begin{python}
211    qS=Scalar(0,FunctionOnBoundary(domain))
212    qS.setTaggedValue("linebottom",qin)
213    mypde.setValue(y=-qS)
214    \end{python}
215    To plot the results we are using the \modmpl library as shown \refSec{Sec:2DHD plot}. For convenience
216    the interpolation of the temperature to a rectangular grid for contour plotting is made available
217    via the \verb|toRegGrid| function in the \verb|cblib| module. Your result should look similar to
218    \reffig{fig:steady state heat}.
219
220    \begin{figure}[ht]
221    \centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}}
222    \caption{Heat refraction model with point and line labels.}
223    \label{fig:anticlinehrmodel}
224    \end{figure}
225
226    \begin{figure}[ht]
227    \centerline{\includegraphics[width=4.in]{figures/heatrefraction}}
228    \caption{Temperature Distribution in the Heat Refraction Model.}
229    \label{fig:anticlinetemp}
230    \end{figure}
231
232    \section{A Heat Refraction model}
233    \sslist{heatrefraction_solver001.py and cblib.py}
234
235    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.
236
237
238    We will modify the example of the previous section by subdividing the block into two parts. The curve
239    separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points
240    used to define the curve may be measurement but for simplicity we assume here that there coordinates are
241    known in analytic form.
242
243    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.
244
245    It is now possible to start defining our domain and boundaries.
246
247    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.
248  \begin{python}  \begin{python}
249  # Material Boundary  # Material Boundary
250  x=[ Point(i*dsp\  x=[ Point(i*dsp\
# Line 91  mysp = Spline(*tuple(x)) #*tuple() force Line 255  mysp = Spline(*tuple(x)) #*tuple() force
255  \end{python}  \end{python}
256  The start and end points of the spline can be returned to help define the material boundaries.  The start and end points of the spline can be returned to help define the material boundaries.
257  \begin{python}  \begin{python}
258  x1=Spline.getStartPoint(mysp)  x1=mysp.getStartPoint()
259  x2=Spline.getEndPoint(mysp)  x2=mysp.getEndPoint()
260  \end{python}  \end{python}
261  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.  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.
262  \begin{python}  \begin{python}
# Line 100  The top block or material above the clin Line 264  The top block or material above the clin
264  tbl1=Line(p0,x1)  tbl1=Line(p0,x1)
265  tbl2=mysp  tbl2=mysp
266  tbl3=Line(x2,p3)  tbl3=Line(x2,p3)
267    l30=Line(p3, p0)
268  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
269  tblock = PlaneSurface(tblockloop)  tblock = PlaneSurface(tblockloop)
270  \end{python}  \end{python}
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.
\begin{python}
tpg = getLoopCoords(tblockloop)
np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
\end{python}
271  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.  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.
272  \begin{python}  \begin{python}
273  bbl4=-mysp  bbl4=-mysp
# Line 115  bbl4=-mysp Line 275  bbl4=-mysp
275  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.  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.
276  \begin{python}  \begin{python}
277  #clockwise check  #clockwise check
278  bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))  bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
bpg = getLoopCoords(bblockloop2)
np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
279  \end{python}  \end{python}
280  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.  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.
281  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.  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.
# Line 128  d=Design(dim=2, element_size=200) Line 286  d=Design(dim=2, element_size=200)
286  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
287                     PropertySet("linebottom",l12))                     PropertySet("linebottom",l12))
288  # Create the geometry, mesh and \esc domain  # Create the geometry, mesh and \esc domain
289  d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))  d.setScriptFileName(os.path.join(save_path,"heatrefraction.geo"))
290  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))  d.setMeshFileName(os.path.join(save_path,"heatrefraction.msh"))
291  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\  domain=MakeDomain(d,optimizeLabeling=True)
optimizeLabeling=True)
# Create a file that can be read back in to python with mesh=ReadMesh(fileName)
domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
\end{python}
The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.

\subsection{Steady-state PDE solution}
\sslist{heatrefraction.py and cblib.py}
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.
\begin{equation}
\rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
\label{eqn:hd2}
\end{equation}
In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to
\begin{equation}
- \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
\label{eqn:sshd}
\end{equation}
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
\begin{equation}
u=r \textit{    where    } q>0
\label{eqn:bdr}
\end{equation}
using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true.
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.
\begin{python}
##ESTABLISHING VARIABLES
qin=70*Milli*W/(m*m) #our heat source temperature is now zero
Ti=290.15*K # Kelvin #the starting temperature of our iron bar
width=5000.0*m
depth=-6000.0*m
\end{python}
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.
\begin{python}
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
tpgx = tpg[:,0]
tpgy = tpg[:,1]
bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
bpgx = bpg[:,0]
bpgy = bpg[:,1]

# set up kappa (thermal conductivity across domain) using tags
kappa=Scalar(0,Function(mymesh))
kappa.setTaggedValue("top",2.0)
kappa.setTaggedValue("bottom",4.0)
\end{python}
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()  .
\begin{python}
#... generate functionspace...
#... open PDE ...
mypde=LinearPDE(mymesh)
#define first coefficient
mypde.setValue(A=kappa*kronecker(mymesh))

# ... set initial temperature ....
x=mymesh.getX()

qH=Scalar(0,FunctionOnBoundary(mymesh))
qH.setTaggedValue("linebottom",qin)
mypde.setValue(q=whereZero(x),r=Ti)
mypde.setValue(y=qH)

# get steady state solution.
T=mypde.getSolution()
292  \end{python}  \end{python}
293  The problem is now solved and plotting is required to visualise the data.  The creation of our domain and its mesh is complete.
294
295  \subsection{Line profiles of 2D data}  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.
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.
\begin{python}
proj=Projector(mymesh)
smthT=proj(T)
\end{python}
The data must then be moved to a regular grid.
296  \begin{python}  \begin{python}
297  xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)  # set up kappa (thermal conductivity across domain) using tags
298  \end{python}  kappa=Scalar(0,Function(domain))
299  The profile can then be plotted by taking a slice of the data output \verb zi  .  kappa.setTaggedValue("top",2.0*W/m/K)
300  \begin{python}  kappa.setTaggedValue("bottom",4.0*W/m/K)
cut=int(len(xi)/2)
pl.clf()
pl.plot(zi[:,cut],yi)
pl.title("Heat Refraction Temperature Depth Profile")
pl.xlabel("Temperature (K)")
pl.ylabel("Depth (m)")
301  \end{python}  \end{python}
302  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}.  No further changes are required to the PDE solution step, see \reffig{fig:anticlinetemp} for the result.
303
304  \begin{figure}  \begin{figure}
305  \centering  \centering
306      \subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tdp.png}\label{fig:tdp}}      \subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{fig:tdp}}
307      \subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tgdp.png}\label{fig:tgdp}}      \subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{fig:tgdp}}
308      \subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefraction001tcdp.png}\label{fig:tcdp}}      \subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{fig:tcdp}}
309      \subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001hf.png}\label{fig:hf}}      \subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf}}
310    \caption{Depth profiles down centre of model.}    \caption{Depth profiles down centre of model.}
311    \label{figs:dps}    \label{figs:dps}
312  \end{figure}  \end{figure}
313
314    \section{Line profiles of 2D data}
315    \sslist{heatrefraction_solver002.py and cblib.py}
316
317    We want to investigate the profile of the data of the last example.
318    We are particularly interested in the deepth profile of the heat flux which is
319    the second component of $-\kappa \nabla T$. We will extend the script developed in the
320    previous section to show how for instance vertical profile can be plotted.
321
322    The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and
323    will in fact represent the values at numerical interpolation points. This assumption is reasonable as
324    the flux is the product of the piecewise constant function $\kappa$ and
325    the gradient of the temperature $T$ which has a kink across the rock interface.
326    Before plotting this function we need to smooth it using the
327    \verb|Projector()| factory;
328    \begin{python}
329    from esys.escript.pdetools import Projector
330    proj=Projector(domain)
331    qu=proj(-kappa*grad(T))
332    \end{python}
333    The \verb|proj| object provides a mechanism to distribute values given at the numerical interpolation points
334    - in this case the heat flux - to the nodes of the FEM mesh. \verb|qu| has the same attached function space
335    like the temperature \verb|T|. The smoothed flux is interpolated
336    to a regular $200\times 200$ grid;
337    \begin{python}
338    xiq,yiq,ziq = toRegGrid(qu,200,200)
339    \end{python}
340    using the \verb|toRegGrid| function from the cookbook library which we are already using for the contour plot.
341    At return \verb|ziq[j,i]| is the value of vertical heat flux at point
342    (\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by
343    plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script
344    creates a deep profile at $x_{0}=\frac{width}{2}$;
345    \begin{python}
346    cut=int(len(xiq)/2)
347    pl.plot(ziq[:,cut]*1000.,yiq)
348    pl.title("Vertical Heat Flow Depth Profile")
349    pl.xlabel("Heat Flow (mW/m^2)")
350    pl.ylabel("Depth (m)")
351    pl.savefig(os.path.join(save_path,"heatrefraction_hf.png"))
352    \end{python}
353    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}.
354
355    \begin{figure}[ht]
356    \centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}}
357    \caption{Heat refraction model with gradient indicated by vectors.}
358    \label{fig:hr001qumodel}
359    \end{figure}
360
361  \subsection{Quiver plots in \mpl}  \section{Arrow plots in \mpl}
362  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.  \sslist{heatrefraction_solver003.py and cblib.py}
363    We would like to visualize the distribution of the flux $-\kappa \nabla T$ over the domain
364  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  .  and produce a plot like shown in \reffig{fig:hr001qumodel}.
365    The plot puts together three components. A contour plot of the temperature,
366    a colored representation of the two sub-domains where colour represents the thermal conductivity
367    in the particular region and finally the arrows representing the direction of flux.
368
369    Contours have already been discussed in \refSec{Sec:2DHD plot}. To show sub-domains,
370    we need to go back to \pycad data to get the points used to describe the boundary of the
371    sub-domains. We have created the \verb|CurveLoop| class object
372    \verb|tblockloop| to define the boundary of the upper sub-domain.
373    We use the \verb|getPolygon()| method of \verb|CurveLoop| to get
374    access to the \verb|Point|s used top define the boundary. The statement
375    \begin{python}
376    [ p.getCoordinates() for p in tblockloop.getPolygon() ])
377    \end{python}
378    creates a list of the node coordinates of all the points in question. In order
379    to simplify the selection of the $x$ and $y$ coordinates the list is converted
380    into \modnumpy array. To add the area colored in brown to the plot we use;
381    \begin{python}
382    import pylab as pl
383    import numarray as np
384    tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
385    pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
386    \end{python}
387    We can apply the same code to \verb|bblockloop| to a red area for this sub-domain to the block.
388
389    To plot vectors representing the flux orientation we use the
390    \verb|quiver| function in \pylab. The function places vectors at locations in the domain.
391    For instance one can plot vectors at the locations of the sample points used by \esc
392    to represent the flux \verb|-kappa*grad(T)|. As a vector is plotted at each sample point one typically ends
393    up with two many vectors. So one needs to select a subset of points:
394    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;
395    \begin{python}
396    dx = width/nx # x spacing
397    dy = depth/ny # y spacing
398    grid = [ ] # the grid points
399    for j in xrange(0,ny-1):
400        for i in xrange(0,nx-1):
401               grid.append([dx/2+dx*i,dy/2+dy*j])
402    \end{python}
403    With the \verb|Locator| \esc provides a mechanism to identify sample points that are closest
404    to the the grid points we have selected and to get the data at these data points;
405    \begin{python}
406    from esys.escript.pdetools import Locator
407    flux=-kappa*grad(T)
408    fluxLoc = Locator(flux.getFunctionSpace(),grid)
409    subflux= fluxLoc(flux)
410    \end{python}
411    \verb|subflux| gives now a list of flux component at certain sample points. To get the
412    list of the sample point coordinates one can use the \verb|getX()| method of the
413    \verb|Locator|;
414    \begin{python}
415    subfluxloc = fluxLoc.getX()
416    \end{python}
417    To simplify the selection of $x$ and $y$ components it is convenient
418    to transform \verb|subflux| and \verb|subfluxloc| to \numpy arrays
419    \verb|xflux|, \verb|flux|.
420    This function is implemented in the \verb|subsample|
421    in the  \file{clib.py} file so we can use it in other examples. One can easily use this function
422    to create a vector plot of the flux;
423    \begin{python}
424    from cblib import subsample
425    xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
426    pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
427    \end{python}
428    We add title and labels;
429  \begin{python}  \begin{python}
# calculate gradient of solution for quiver plot
qu=-kappa*grad(T)

# rearrange mymesh to suit solution function space
oldspacecoords=mymesh.getX()
coords=Data(oldspacecoords, T.getFunctionSpace())
\end{python}
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 .
\begin{python}
quivshape = [20,20] #quivers x and quivers y
# function to calculate quiver locations
qu,qulocs = toQuivLocs(quivshape,width,depth,qu)

kappaT = kappa.toListOfTuples(scalarastuple=False)
coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
coordKX, coordKY = toXYTuple(coordsK)

tempT = T.toListOfTuples(scalarastuple=False)
coordX, coordY = toXYTuple(coords)
\end{python}
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.
\begin{python}
xi = np.linspace(0.0,width,100)
yi = np.linspace(depth,0.0,100)
# grid the data.
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
# contour the gridded data, plotting dots at the randomly spaced data points.
pl.matplotlib.pyplot.autumn()

CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
#~ CK = pl.contourf(xi,yi,ziK,2)
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')

pl.clabel(CS, inline=1, fontsize=8)
430  pl.title("Heat Refraction across a clinal structure.")  pl.title("Heat Refraction across a clinal structure.")
431  pl.xlabel("Horizontal Displacement (m)")  pl.xlabel("Horizontal Displacement (m)")
432  pl.ylabel("Depth (m)")  pl.ylabel("Depth (m)")
#~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))

QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
433  pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")  pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
434  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))  pl.savefig(os.path.join(saved_path,"heatrefractionflux.png"))
435  \end{python}  \end{python}
436  The data has now been plotted.  to get the desired result, see \reffig{fig:hr001qumodel}.
\begin{figure}[ht]
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}}
\caption{Heat refraction model with gradient indicated by quivers.}
\label{fig:hr001qumodel}
\end{figure}
437
\newpage
\subsection{Fault and Overburden Model}
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.
438  \begin{figure}[ht]  \begin{figure}[ht]
439  \centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}}  \centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}}
440  \caption{Heat refraction model with three blocks and gradient quivers.}  \caption{Heat refraction model with three blocks and heat flux.}
441    \label{fig:hr002qumodel}
442  \end{figure}  \end{figure}
443
444    \section{Fault and Overburden Model}
445    \sslist{heatrefraction2_solver.py and cblib.py}
446    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.
447
448

Legend:
 Removed from v.2947 changed lines Added in v.2948

 ViewVC Help Powered by ViewVC 1.1.26