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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2706 - (hide annotations)
Fri Oct 2 02:15:04 2009 UTC (12 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 6680 byte(s)
Perhaps this will stop the cookbook from hating me.


1 ahallam 2588
2     ########################################################
3     #
4 jfenwick 2667 # Copyright (c) 2009 by University of Queensland
5 ahallam 2588 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14 jfenwick 2667 __copyright__="""Copyright (c) 2009 by University of Queensland
15 ahallam 2588 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="https://launchpad.net/escript-finley"
21    
22     """
23     Author: Antony Hallam antony.hallam@uqconnect.edu.au
24     """
25    
26 ahallam 2658 ############################################################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 ahallam 2634 # To solve the problem it is necessary to import the modules we
33     # require.
34 jfenwick 2667
35    
36    
37 ahallam 2634 # This imports everything from the escript library
38     from esys.escript import *
39     # This defines LinearPDE as LinearPDE
40 ahallam 2681 from esys.escript.linearPDEs import LinearPDE, Poisson
41     # smoothing operator
42 ahallam 2672 from esys.escript.pdetools import Projector
43 ahallam 2634 # This imports the rectangle domain function from finley
44     from esys.finley import Rectangle, ReadMesh, Domain
45     # This package is necessary to handle saving our data.
46     import os
47     # A useful unit handling package which will make sure all our units
48     # match up in the equations.
49     from esys.escript.unitsSI import *
50     # numpy for array handling
51 ahallam 2588 import numpy as np
52 jfenwick 2667
53 ahallam 2658 import matplotlib
54 jfenwick 2667 #For interactive use, you can comment out the next line
55 ahallam 2658 matplotlib.use('agg') #It's just here for automated testing
56 jfenwick 2667
57 ahallam 2634 # pylab for matplotlib and plotting
58 ahallam 2588 import pylab as pl
59 ahallam 2634 # cblib functions
60 jfenwick 2706 from cblib2 import toQuivLocs, toXYTuple, toRegGrid, gradtoRegGrid
61     from cblib1 import needdirs
62 ahallam 2588
63 ahallam 2658 ########################################################MPI WORLD CHECK
64     if getMPISizeWorld() > 1:
65     import sys
66     print "This example will not run in an MPI world."
67     sys.exit(0)
68    
69     #################################################ESTABLISHING VARIABLES
70 ahallam 2588 qin=70*Milli*W/(m*m) #our heat source temperature is now zero
71     Ti=290.15*K # Kelvin #the starting temperature of our iron bar
72 ahallam 2634 width=5000.0*m
73     depth=-6000.0*m
74 ahallam 2588
75 ahallam 2597 # the folder to gett our outputs from, leave blank "" for script path -
76     # note these depen. are generated from heatrefraction_mesher001.py
77 ahallam 2658 saved_path = save_path= os.path.join("data","heatrefrac001" )
78 jfenwick 2648 needdirs([saved_path])
79    
80 ahallam 2658 ################################################ESTABLISHING PARAMETERS
81 ahallam 2588 ## DOMAIN
82 ahallam 2597 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
83     tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
84 ahallam 2588 tpgx = tpg[:,0]
85     tpgy = tpg[:,1]
86 ahallam 2597 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
87 ahallam 2588 bpgx = bpg[:,0]
88     bpgy = bpg[:,1]
89    
90     # set up kappa (thermal conductivity across domain) using tags
91     kappa=Scalar(0,Function(mymesh))
92     kappa.setTaggedValue("top",2.0)
93     kappa.setTaggedValue("bottom",4.0)
94    
95     #... generate functionspace...
96     #... open PDE ...
97     mypde=LinearPDE(mymesh)
98     #define first coefficient
99     mypde.setValue(A=kappa*kronecker(mymesh))
100    
101     # ... set initial temperature ....
102     x=mymesh.getX()
103    
104 ahallam 2597 qH=Scalar(0,FunctionOnBoundary(mymesh))
105     qH.setTaggedValue("linebottom",qin)
106 ahallam 2588 mypde.setValue(q=whereZero(x[1]),r=Ti)
107 ahallam 2658 mypde.setValue(y=qH)
108 ahallam 2588
109 ahallam 2658 ###########################################################GET SOLUTION
110 ahallam 2588 T=mypde.getSolution()
111    
112 ahallam 2658 ##########################################################VISUALISATION
113     # calculate gradient of solution for quiver plot
114     qu=-kappa*grad(T)
115 ahallam 2672 quT=qu.toListOfTuples()
116 ahallam 2658
117 ahallam 2672 #Projector is used to smooth the data.
118     proj=Projector(mymesh)
119     smthT=proj(T)
120 ahallam 2588
121 ahallam 2672 #move data to a regular grid for plotting
122     xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
123 ahallam 2588
124 ahallam 2672 #Temperature Depth Profile along x[50]
125     cut=int(len(xi)/2)
126     pl.clf()
127     pl.plot(zi[:,cut],yi)
128     pl.title("Heat Refraction Temperature Depth Profile")
129     pl.xlabel("Temperature (K)")
130     pl.ylabel("Depth (m)")
131     if getMPIRankWorld() == 0: #check for MPI processing
132     pl.savefig(os.path.join(saved_path,"heatrefraction001_tdp.png"))
133    
134 ahallam 2681 # Heat flow depth profile.
135 ahallam 2672 pl.clf()
136     # grid the data.
137     qu=proj(-kappa*grad(T))
138     xiq,yiq,ziq = gradtoRegGrid(qu,mymesh,200,200,width,depth,1)
139     cut=int(len(xiq)/2)
140     pl.plot(ziq[:,cut]*1000.,yiq)
141     pl.title("Heat Flow Depth Profile")
142     pl.xlabel("Heat Flow (mW/m^2)")
143     pl.ylabel("Depth (m)")
144     if getMPIRankWorld() == 0: #check for MPI processing
145     pl.savefig(os.path.join(saved_path,"heatrefraction001_hf.png"))
146    
147 ahallam 2681 # Temperature Gradient Depth Profile at x[50]
148 ahallam 2672 pl.clf()
149     zT=proj(-grad(T))
150     xt,yt,zt=gradtoRegGrid(zT,mymesh,200,200,width,depth,1)
151     cut=int(len(xt)/2)
152     pl.plot(zt[:,cut]*1000.,yt)
153     pl.title("Heat Refraction Temperature Gradient \n Depth Profile")
154     pl.xlabel("Temperature (K/Km)")
155     pl.ylabel("Depth (m)")
156     if getMPIRankWorld() == 0: #check for MPI processing
157     pl.savefig(os.path.join(saved_path,"heatrefraction001_tgdp.png"))
158 ahallam 2681
159     # Thermal Conditions Depth Profile
160 ahallam 2672 pl.clf()
161     xk,yk,zk = toRegGrid(kappa,mymesh,200,200,width,depth)
162     cut=int(len(xk)/2)
163     pl.plot(zk[:,cut],yk)
164     pl.title("Heat Refraction Thermal Conductivity Depth Profile")
165     pl.xlabel("Conductivity (W/K/m)")
166     pl.ylabel("Depth (m)")
167     pl.axis([1,5,-6000,0])
168     if getMPIRankWorld() == 0: #check for MPI processing
169     pl.savefig(os.path.join(saved_path,"heatrefraction001_tcdp.png"))
170 ahallam 2681
171     # contour the gridded data,
172     # plotting dots at the randomly spaced data points.
173     quivshape = [20,20] #quivers x and quivers y
174     # function to calculate quiver locations
175     qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
176    
177     # select colour
178     pl.matplotlib.pyplot.autumn()
179     pl.clf()
180     # plot polygons for boundaries
181     CKL = pl.fill(tpgx,tpgy,'brown',label='2 W/m/k',zorder=-1000)
182     CKM = pl.fill(bpgx,bpgy,'red',label='4 W/m/k',zorder=-1000)
183     # contour temperature
184     CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
185     # labels and formatting
186     pl.clabel(CS, inline=1, fontsize=8)
187     pl.title("Heat Refraction across a clinal structure.")
188     pl.xlabel("Horizontal Displacement (m)")
189     pl.ylabel("Depth (m)")
190     pl.legend()
191     if getMPIRankWorld() == 0: #check for MPI processing
192     pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
193    
194     #Quiver Plot qulocs -> tail location, qu -> quiver length/direction
195     QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
196     angles='xy',color="white")
197     pl.title("Heat Refraction across a clinal structure\n with\
198     gradient quivers.")
199     if getMPIRankWorld() == 0: #check for MPI processing
200     pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
201    
202    

  ViewVC Help
Powered by ViewVC 1.1.26