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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2633 by ahallam, Wed Aug 26 22:18:19 2009 UTC revision 2634 by ahallam, Thu Aug 27 04:03:32 2009 UTC
# Line 14  Line 14 
14  \section{Steady-state heat refraction}  \section{Steady-state heat refraction}
15  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 develop finite element grids 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 escript to define properties like material constants and source locations.  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 develop finite element grids 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 escript to define properties like material constants and source locations.
16    
17    \subsection{Creating the model with \pycad}
18    
19  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 infact warped by the structure of model and is heavily dependant on the material properties and shape.  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 infact warped by the structure of model and is heavily dependant on the material properties and shape.
20    
21  \begin{figure}[h!]  \begin{figure}[h!]
# Line 22  Our first heat refraction model will be Line 24  Our first heat refraction model will be
24  \label{fig:anticlinehrmodel}  \label{fig:anticlinehrmodel}
25  \end{figure}  \end{figure}
26    
27  We will start by defining our domain and material boundaries using \pycad. This example is located in the file _________  and a labeled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have:  We will start by defining our domain and material boundaries using \pycad. This example is located in the \esc directory under \exf and a labeled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have:
28  \begin{verbatim}  \begin{verbatim}
29  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
30  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
# Line 42  Line() #Creates a line from a number of Line 44  Line() #Creates a line from a number of
44  CurveLoop() #Creates a closed loop from a number of lines.  CurveLoop() #Creates a closed loop from a number of lines.
45  Surface() #Creates a surface based on a CurveLoop  Surface() #Creates a surface based on a CurveLoop
46  \end{verbatim}  \end{verbatim}
47  We start by inputting 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 and 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.  We start by inputting the variables we need to construct the model.
48    \begin{verbatim}
49    #Model Parameters
50    width=5000.0   #width of model
51    depth=-6000.0  #depth of model
52    sspl=51 #number of discrete points in spline
53    dsp=width/(sspl-1) #dx of spline steps for width
54    dep_sp=2500.0 #avg depth of spline
55    h_sp=1500.0 #heigh of spline
56    orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
57    \end{verbatim}
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.
59  \begin{verbatim}  \begin{verbatim}
60  # Overall Domain  # Overall Domain
61  p0=Point(0.0,        0.0, 0.0)  p0=Point(0.0,      0.0, 0.0)
62  p1=Point(0.0,    -6000.0, 0.0)  p1=Point(0.0,    depth, 0.0)
63  p2=Point(5000.0, -6000.0, 0.0)  p2=Point(width, depth, 0.0)
64  p3=Point(5000.0,     0.0, 0.0)  p3=Point(width,   0.0, 0.0)
65  \end{verbatim}  \end{verbatim}
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 in an \textbf{anti-clockwise} direction using our points. This forms a rectangle around our domain.
67  \begin{verbatim}  \begin{verbatim}
# Line 63  c=CurveLoop(l01, l12, l23, l30) Line 76  c=CurveLoop(l01, l12, l23, l30)
76  \end{verbatim}  \end{verbatim}
77  The curve defining our clinal structure is located in approximately the middle of the domain and has a cosinusoidal 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.  The curve defining our clinal structure is located in approximately the middle of the domain and has a cosinusoidal 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.
78  \begin{verbatim}  \begin{verbatim}
79  x=[ Point(i*100.0,-2500+modal*1500.*cos(pi*i*100.0/2500.0+pi)) for i in range(0,51) ]  # Material Boundary
80    x=[ Point(i*dsp\
81        ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
82         for i in range(0,sspl)\
83        ]
84  mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple  mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple
85  \end{verbatim}  \end{verbatim}
86  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.
# Line 112  domain.write(os.path.join(save_path,"hea Line 129  domain.write(os.path.join(save_path,"hea
129  \end{verbatim}  \end{verbatim}
130  The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.  The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem.
131    
132    \subsection{Steady-state PDE solution}
133    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.
134    \begin{equation}
135    \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
136    \label{eqn:hd2}
137    \end{equation}
138    In the steady state PDE however, there is no time dependance. This means that \refEq{eqn:hd2} reduces to
139    \begin{equation}
140    - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
141    \label{eqn:sshd}
142    \end{equation}
143    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
144    \begin{equation}
145    u=r \textit{    where    } q>0
146    \label{eqn:bdr}
147    \end{equation}
148    using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true.
149    The structure of our solution script is similar to those we have used already. Firstly, the dependencies must be imported. Then the variables must be defined. 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.
150    \begin{verbatim}
151    ##ESTABLISHING VARIABLES
152    qin=70*Milli*W/(m*m) #our heat source temperature is now zero
153    Ti=290.15*K # Kelvin #the starting temperature of our iron bar
154    width=5000.0*m
155    depth=-6000.0*m
156    \end{verbatim}
157    The mesh is now imported via the \verb ReadMesh()  function and loaded into a domain variable \verb mymesh  . The structural polygons are next. 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.
158    \begin{verbatim}
159    mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
160    tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
161    tpgx = tpg[:,0]
162    tpgy = tpg[:,1]
163    bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
164    bpgx = bpg[:,0]
165    bpgy = bpg[:,1]
166    
167    # set up kappa (thermal conductivity across domain) using tags
168    kappa=Scalar(0,Function(mymesh))
169    kappa.setTaggedValue("top",2.0)
170    kappa.setTaggedValue("bottom",4.0)
171    \end{verbatim}
172    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()  .
173    \begin{verbatim}
174    #... generate functionspace...
175    #... open PDE ...
176    mypde=LinearPDE(mymesh)
177    #define first coefficient
178    mypde.setValue(A=kappa*kronecker(mymesh))
179    
180    # ... set initial temperature ....
181    x=mymesh.getX()
182    
183    qH=Scalar(0,FunctionOnBoundary(mymesh))
184    qH.setTaggedValue("linebottom",qin)
185    mypde.setValue(q=whereZero(x[1]),r=Ti)
186    mypde.setValue(y=qH)#,r=17*Celsius)
187    
188    # get steady state solution and export to vtk.
189    T=mypde.getSolution()
190    \end{verbatim}
191    
192    \subsection{Quiver plots in \mpl}
193    \begin{verbatim}
194    # rearrage mymesh to suit solution function space      
195    oldspacecoords=mymesh.getX()
196    coords=Data(oldspacecoords, T.getFunctionSpace())
197    
198    # calculate gradient of solution for quiver plot
199    qu=-kappa*grad(T)
200    
201    # create quiver locations
202    quivshape = [20,20] #quivers x and quivers y
203    numquiv = quivshape[0]*quivshape[1] # total number of quivers
204    maxx = 5000. # maximum x displacement
205    dx = maxx/quivshape[0]+1. # quiver x spacing
206    maxy = -6000. # maximum y displacement
207    dy = maxy/quivshape[1]+1. # quiver y spacing
208    qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
209    # fill qulocs
210    for i in range(0,quivshape[0]-1):
211        for j in range(0,quivshape[1]-1):
212            qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
213    # retreive values for quivers direction and shape from qu
214    quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
215    qu = quL(qu) #array of dx,dy for quivers
216    qu = np.array(qu) #turn into a numpy array
217    qulocs = quL.getX() #returns actual locations from data
218    qulocs = np.array(qulocs) #turn into a numpy array
219    
220    kappaT = kappa.toListOfTuples(scalarastuple=False)
221    coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
222    coordsK = np.array(coordsK.toListOfTuples())
223    coordKX = coordsK[:,0]
224    coordKY = coordsK[:,1]
225          
226    coords = np.array(coords.toListOfTuples())
227    tempT = T.toListOfTuples(scalarastuple=False)
228    
229    coordX = coords[:,0]
230    coordY = coords[:,1]
231    
232    xi = np.linspace(0.0,5000.0,100)
233    yi = np.linspace(-6000.0,0.0,100)
234    # grid the data.
235    zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
236    ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
237    # contour the gridded data, plotting dots at the randomly spaced data points.
238    
239    pl.matplotlib.pyplot.autumn()
240    CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
241    #~ CK = pl.contourf(xi,yi,ziK,2)
242    CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
243    pl.clabel(CS, inline=1, fontsize=8)
244    pl.title("Heat Refraction across a clinal structure.")
245    pl.xlabel("Horizontal Displacement (m)")
246    pl.ylabel("Depth (m)")
247    #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
248    pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
249    
250    QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
251    pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
252    pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
253    \end{verbatim}
254    
255    
256    \begin{figure}
257    \centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}}
258    \caption{Heat refraction model with gradient indicated by quivers.}
259    \label{fig:hr001qumodel}
260    \end{figure}
261    

Legend:
Removed from v.2633  
changed lines
  Added in v.2634

  ViewVC Help
Powered by ViewVC 1.1.26