/[escript]/branches/stage3.1/doc/examples/cookbook/heatrefraction001.py
ViewVC logotype

Contents of /branches/stage3.1/doc/examples/cookbook/heatrefraction001.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2945 - (show annotations)
Wed Feb 24 00:17:46 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 6682 byte(s)
Bringing release stage up to trunk version 2944

1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 by University of Queensland
5 # 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 __copyright__="""Copyright (c) 2009-2010 by University of Queensland
15 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 ############################################################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
33 # require.
34
35
36
37 # This imports everything from the escript library
38 from esys.escript import *
39 # This defines LinearPDE as LinearPDE
40 from esys.escript.linearPDEs import LinearPDE, Poisson
41 # smoothing operator
42 from esys.escript.pdetools import Projector
43 # 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 import numpy as np
52
53 import matplotlib
54 #For interactive use, you can comment out the next line
55 matplotlib.use('agg') #It's just here for automated testing
56
57 # pylab for matplotlib and plotting
58 import pylab as pl
59 # cblib functions
60 from cblib import toXYTuple
61 from cblib2 import toQuivLocs, toRegGrid, gradtoRegGrid
62
63 ########################################################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 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 width=5000.0*m
73 depth=-6000.0*m
74
75 # the folder to gett our outputs from, leave blank "" for script path -
76 # note these depen. are generated from heatrefraction_mesher001.py
77 saved_path = save_path= os.path.join("data","heatrefrac001" )
78 mkDir(saved_path)
79
80 ################################################ESTABLISHING PARAMETERS
81 ## DOMAIN
82 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
83 tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
84 tpgx = tpg[:,0]
85 tpgy = tpg[:,1]
86 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
87 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 qH=Scalar(0,FunctionOnBoundary(mymesh))
105 qH.setTaggedValue("linebottom",qin)
106 mypde.setValue(q=whereZero(x[1]),r=Ti)
107 mypde.setValue(y=qH)
108
109 ###########################################################GET SOLUTION
110 T=mypde.getSolution()
111
112 ##########################################################VISUALISATION
113 # calculate gradient of solution for quiver plot
114 qu=-kappa*grad(T)
115 quT=qu.toListOfTuples()
116
117 #Projector is used to smooth the data.
118 proj=Projector(mymesh)
119 smthT=proj(T)
120
121 #move data to a regular grid for plotting
122 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
123
124 #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 # Heat flow depth profile.
135 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 # Temperature Gradient Depth Profile at x[50]
148 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
159 # Thermal Conditions Depth Profile
160 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
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 pl.clf()

  ViewVC Help
Powered by ViewVC 1.1.26