/[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 2657 by jfenwick, Mon Sep 7 00:06:15 2009 UTC revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC
# Line 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    
26    ############################################################FILE HEADER
27    # heatrefraction001.py
28    # Model steady state temperature distribution in two block model, mesh
29    # from heatrefraction_mesher001.py
30    
31    #######################################################EXTERNAL MODULES
32  # To solve the problem it is necessary to import the modules we  # To solve the problem it is necessary to import the modules we
33  # require.  # require.
34  # This imports everything from the escript library  # This imports everything from the escript library
# Line 38  import os Line 44  import os
44  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
45  # numpy for array handling  # numpy for array handling
46  import numpy as np  import numpy as np
47    import matplotlib
48    #For interactive use, you can comment out the next two lines
49    import matplotlib
50    matplotlib.use('agg') #It's just here for automated testing
51  # pylab for matplotlib and plotting  # pylab for matplotlib and plotting
52  import pylab as pl  import pylab as pl
53  # cblib functions  # cblib functions
54  from cblib import toQuivLocs, toXYTuple, needdirs  from cblib import toQuivLocs, toXYTuple, needdirs
55    
56  ##ESTABLISHING VARIABLES  ########################################################MPI WORLD CHECK
57    if getMPISizeWorld() > 1:
58        import sys
59        print "This example will not run in an MPI world."
60        sys.exit(0)
61    
62    #################################################ESTABLISHING VARIABLES
63  qin=70*Milli*W/(m*m) #our heat source temperature is now zero  qin=70*Milli*W/(m*m) #our heat source temperature is now zero
64  Ti=290.15*K # Kelvin #the starting temperature of our iron bar  Ti=290.15*K # Kelvin #the starting temperature of our iron bar
65  width=5000.0*m  width=5000.0*m
# Line 51  depth=-6000.0*m Line 67  depth=-6000.0*m
67    
68  # the folder to gett our outputs from, leave blank "" for script path -  # the folder to gett our outputs from, leave blank "" for script path -
69  # note these depen. are generated from heatrefraction_mesher001.py  # note these depen. are generated from heatrefraction_mesher001.py
70  saved_path = "data/heatrefrac001"  saved_path = save_path= os.path.join("data","heatrefrac001" )
   
71  needdirs([saved_path])  needdirs([saved_path])
72    
73  ###### 2 BLOCK MODEL #########  ################################################ESTABLISHING PARAMETERS
74  ## DOMAIN  ## DOMAIN
 ## Anticline  
75  mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))  mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
76  tpg = np.loadtxt(os.path.join(saved_path,"toppg"))  tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
77  tpgx = tpg[:,0]  tpgx = tpg[:,0]
# Line 65  tpgy = tpg[:,1] Line 79  tpgy = tpg[:,1]
79  bpg = np.loadtxt(os.path.join(saved_path,"botpg"))  bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
80  bpgx = bpg[:,0]  bpgx = bpg[:,0]
81  bpgy = bpg[:,1]  bpgy = bpg[:,1]
 ## Syncline  
 #~ mymesh = ReadMesh("heatrefraction_mesh002.fly")  
82    
83  # set up kappa (thermal conductivity across domain) using tags  # set up kappa (thermal conductivity across domain) using tags
84  kappa=Scalar(0,Function(mymesh))  kappa=Scalar(0,Function(mymesh))
85  kappa.setTaggedValue("top",2.0)  kappa.setTaggedValue("top",2.0)
86  kappa.setTaggedValue("bottom",4.0)  kappa.setTaggedValue("bottom",4.0)
87    
 ###### 3 BLOCK MODEL #########  
 # set up kappa (thermal conductivity across domain) using tags  
 #~ mymesh = ReadMesh("heatrefraction_mesh003.fly")  
 #~ kappa=Scalar(0,Function(mymesh))  
 #~ kappa.setTaggedValue("top",2.0)  
 #~ kappa.setTaggedValue("bottomleft",3.0)  
 #~ kappa.setTaggedValue("bottomright",4.0)  
   
88  #... generate functionspace...  #... generate functionspace...
89  #... open PDE ...  #... open PDE ...
90  mypde=LinearPDE(mymesh)  mypde=LinearPDE(mymesh)
# Line 93  x=mymesh.getX() Line 97  x=mymesh.getX()
97  qH=Scalar(0,FunctionOnBoundary(mymesh))  qH=Scalar(0,FunctionOnBoundary(mymesh))
98  qH.setTaggedValue("linebottom",qin)  qH.setTaggedValue("linebottom",qin)
99  mypde.setValue(q=whereZero(x[1]),r=Ti)  mypde.setValue(q=whereZero(x[1]),r=Ti)
100  mypde.setValue(y=qH)#,r=17*Celsius)  mypde.setValue(y=qH)
101    
102  # get steady state solution and export to vtk.  ###########################################################GET SOLUTION
103  T=mypde.getSolution()  T=mypde.getSolution()
104  #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))  
105    ##########################################################VISUALISATION
106    # calculate gradient of solution for quiver plot
107    qu=-kappa*grad(T)
108    
109  # rearrage mymesh to suit solution function space        # rearrage mymesh to suit solution function space      
110  oldspacecoords=mymesh.getX()  oldspacecoords=mymesh.getX()
111  coords=Data(oldspacecoords, T.getFunctionSpace())  coords=Data(oldspacecoords, T.getFunctionSpace())
112    
 # calculate gradient of solution for quiver plot  
 qu=-kappa*grad(T)  
113  quivshape = [20,20] #quivers x and quivers y  quivshape = [20,20] #quivers x and quivers y
114  # function to calculate quiver locations  # function to calculate quiver locations
115  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)  qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
# Line 121  yi = np.linspace(depth,0.0,100) Line 126  yi = np.linspace(depth,0.0,100)
126  # grid the data.  # grid the data.
127  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
128  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)  ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
129  # contour the gridded data, plotting dots at the randomly spaced data points.  # contour the gridded data,
130    # plotting dots at the randomly spaced data points.
131    
132    # select colour
133  pl.matplotlib.pyplot.autumn()  pl.matplotlib.pyplot.autumn()
134    # plot polygons for boundaries
135  CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)  CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
136  #~ CK = pl.contourf(xi,yi,ziK,2)  # contour temperature
137  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')  CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
138    # labels and formatting
139  pl.clabel(CS, inline=1, fontsize=8)  pl.clabel(CS, inline=1, fontsize=8)
140  pl.title("Heat Refraction across a clinal structure.")  pl.title("Heat Refraction across a clinal structure.")
141  pl.xlabel("Horizontal Displacement (m)")  pl.xlabel("Horizontal Displacement (m)")
142  pl.ylabel("Depth (m)")  pl.ylabel("Depth (m)")
143  #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')  if getMPIRankWorld() == 0: #check for MPI processing
144  pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))      pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
145    
146  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")  #Quiver Plot qulocs -> tail location, qu -> quiver length/direction
147  pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")  QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
148  pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))                  angles='xy',color="white")
149    pl.title("Heat Refraction across a clinal structure \n with\
150                                                        gradient quivers.")
151    if getMPIRankWorld() == 0: #check for MPI processing
152        pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))

Legend:
Removed from v.2657  
changed lines
  Added in v.2658

  ViewVC Help
Powered by ViewVC 1.1.26