/[escript]/trunk/doc/examples/cookbook/heatrefraction001.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/heatrefraction001.py

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

revision 2671 by jfenwick, Thu Sep 17 01:49:11 2009 UTC revision 2672 by ahallam, Fri Sep 18 01:46:47 2009 UTC
# Line 38  Author: Antony Hallam antony.hallam@uqco Line 38  Author: Antony Hallam antony.hallam@uqco
38  from esys.escript import *  from esys.escript import *
39  # This defines LinearPDE as LinearPDE  # This defines LinearPDE as LinearPDE
40  from esys.escript.linearPDEs import LinearPDE, Poisson  from esys.escript.linearPDEs import LinearPDE, Poisson
41    from esys.escript.pdetools import Projector
42  # This imports the rectangle domain function from finley  # This imports the rectangle domain function from finley
43  from esys.finley import Rectangle, ReadMesh, Domain  from esys.finley import Rectangle, ReadMesh, Domain
44  # This package is necessary to handle saving our data.  # This package is necessary to handle saving our data.
# Line 55  matplotlib.use('agg') #It's just here fo Line 56  matplotlib.use('agg') #It's just here fo
56  # pylab for matplotlib and plotting  # pylab for matplotlib and plotting
57  import pylab as pl  import pylab as pl
58  # cblib functions  # cblib functions
59  from cblib import toQuivLocs, toXYTuple, needdirs  from cblib import toQuivLocs, toXYTuple, needdirs, toRegGrid, gradtoRegGrid
60    
61  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
62  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
# Line 109  T=mypde.getSolution() Line 110  T=mypde.getSolution()
110  ##########################################################VISUALISATION  ##########################################################VISUALISATION
111  # calculate gradient of solution for quiver plot  # calculate gradient of solution for quiver plot
112  qu=-kappa*grad(T)  qu=-kappa*grad(T)
113    quT=qu.toListOfTuples()
114    
115  # rearrage mymesh to suit solution function space        #Projector is used to smooth the data.
116  oldspacecoords=mymesh.getX()  proj=Projector(mymesh)
117  coords=Data(oldspacecoords, T.getFunctionSpace())  smthT=proj(T)
118    
119  quivshape = [20,20] #quivers x and quivers y  quivshape = [20,20] #quivers x and quivers y
120  # function to calculate quiver locations  # function to calculate quiver locations
121  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
122    
123  kappaT = kappa.toListOfTuples(scalarastuple=False)  #move data to a regular grid for plotting
124  coordsK = Data(oldspacecoords, kappa.getFunctionSpace())  xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
 coordKX, coordKY = toXYTuple(coordsK)  
         
 tempT = T.toListOfTuples(scalarastuple=False)  
 coordX, coordY = toXYTuple(coords)  
125    
 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)  
126  # contour the gridded data,  # contour the gridded data,
127  # plotting dots at the randomly spaced data points.  # plotting dots at the randomly spaced data points.
128    
129  # select colour  # select colour
130  pl.matplotlib.pyplot.autumn()  pl.matplotlib.pyplot.autumn()
131  # plot polygons for boundaries  # plot polygons for boundaries
132  CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)  CKL = pl.fill(tpgx,tpgy,'brown',label='2 W/m/k',zorder=-1000)
133    CKM = pl.fill(bpgx,bpgy,'red',label='4 W/m/k',zorder=-1000)
134  # contour temperature  # contour temperature
135  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
136  # labels and formatting  # labels and formatting
# Line 144  pl.clabel(CS, inline=1, fontsize=8) Line 138  pl.clabel(CS, inline=1, fontsize=8)
138  pl.title("Heat Refraction across a clinal structure.")  pl.title("Heat Refraction across a clinal structure.")
139  pl.xlabel("Horizontal Displacement (m)")  pl.xlabel("Horizontal Displacement (m)")
140  pl.ylabel("Depth (m)")  pl.ylabel("Depth (m)")
141    pl.legend()
142  if getMPIRankWorld() == 0: #check for MPI processing  if getMPIRankWorld() == 0: #check for MPI processing
143      pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))      pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
144    
145  #Quiver Plot qulocs -> tail location, qu -> quiver length/direction  #Quiver Plot qulocs -> tail location, qu -> quiver length/direction
146  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
147                  angles='xy',color="white")                  angles='xy',color="white")
148  pl.title("Heat Refraction across a clinal structure \n with\  pl.title("Heat Refraction across a clinal structure\n with\
149                                                      gradient quivers.")  gradient quivers.")
150  if getMPIRankWorld() == 0: #check for MPI processing  if getMPIRankWorld() == 0: #check for MPI processing
151      pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))      pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
152    
153    #Temperature Depth Profile along x[50]
154    cut=int(len(xi)/2)
155    pl.clf()
156    pl.plot(zi[:,cut],yi)
157    pl.title("Heat Refraction Temperature Depth Profile")
158    pl.xlabel("Temperature (K)")
159    pl.ylabel("Depth (m)")
160    if getMPIRankWorld() == 0: #check for MPI processing
161        pl.savefig(os.path.join(saved_path,"heatrefraction001_tdp.png"))
162        
163    #Temperature Gradient Profile along x[50]
164    
165    
166    
167    pl.clf()
168    # grid the data.
169    qu=proj(-kappa*grad(T))
170    xiq,yiq,ziq = gradtoRegGrid(qu,mymesh,200,200,width,depth,1)
171    cut=int(len(xiq)/2)
172    pl.plot(ziq[:,cut]*1000.,yiq)
173    pl.title("Heat Flow Depth Profile")
174    pl.xlabel("Heat Flow (mW/m^2)")
175    pl.ylabel("Depth (m)")
176    if getMPIRankWorld() == 0: #check for MPI processing
177        pl.savefig(os.path.join(saved_path,"heatrefraction001_hf.png"))
178    
179    pl.clf()
180    zT=proj(-grad(T))
181    
182    xt,yt,zt=gradtoRegGrid(zT,mymesh,200,200,width,depth,1)
183    cut=int(len(xt)/2)
184    pl.plot(zt[:,cut]*1000.,yt)
185    pl.title("Heat Refraction Temperature Gradient \n Depth Profile")
186    pl.xlabel("Temperature (K/Km)")
187    pl.ylabel("Depth (m)")
188    if getMPIRankWorld() == 0: #check for MPI processing
189        pl.savefig(os.path.join(saved_path,"heatrefraction001_tgdp.png"))
190        
191    pl.clf()
192    xk,yk,zk = toRegGrid(kappa,mymesh,200,200,width,depth)
193    cut=int(len(xk)/2)
194    pl.plot(zk[:,cut],yk)
195    pl.title("Heat Refraction Thermal Conductivity Depth Profile")
196    pl.xlabel("Conductivity (W/K/m)")
197    pl.ylabel("Depth (m)")
198    pl.axis([1,5,-6000,0])
199    if getMPIRankWorld() == 0: #check for MPI processing
200        pl.savefig(os.path.join(saved_path,"heatrefraction001_tcdp.png"))
201        

Legend:
Removed from v.2671  
changed lines
  Added in v.2672

  ViewVC Help
Powered by ViewVC 1.1.26