# Contents of /trunk/doc/examples/cookbook/heatrefraction001.py

Revision 2681 - (show annotations)
Thu Sep 24 03:04:04 2009 UTC (11 years, 6 months ago) by ahallam
File MIME type: text/x-python
File size: 6661 byte(s)
```proof reading and review of cookbook, minor changes only + updated depth profiles to cb from last week
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2009 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 20 __url__= 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 toQuivLocs, toXYTuple, needdirs, toRegGrid, gradtoRegGrid 61 62 ########################################################MPI WORLD CHECK 63 if getMPISizeWorld() > 1: 64 import sys 65 print "This example will not run in an MPI world." 66 sys.exit(0) 67 68 #################################################ESTABLISHING VARIABLES 69 qin=70*Milli*W/(m*m) #our heat source temperature is now zero 70 Ti=290.15*K # Kelvin #the starting temperature of our iron bar 71 width=5000.0*m 72 depth=-6000.0*m 73 74 # the folder to gett our outputs from, leave blank "" for script path - 75 # note these depen. are generated from heatrefraction_mesher001.py 76 saved_path = save_path= os.path.join("data","heatrefrac001" ) 77 needdirs([saved_path]) 78 79 ################################################ESTABLISHING PARAMETERS 80 ## DOMAIN 81 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 82 tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 83 tpgx = tpg[:,0] 84 tpgy = tpg[:,1] 85 bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 86 bpgx = bpg[:,0] 87 bpgy = bpg[:,1] 88 89 # set up kappa (thermal conductivity across domain) using tags 90 kappa=Scalar(0,Function(mymesh)) 91 kappa.setTaggedValue("top",2.0) 92 kappa.setTaggedValue("bottom",4.0) 93 94 #... generate functionspace... 95 #... open PDE ... 96 mypde=LinearPDE(mymesh) 97 #define first coefficient 98 mypde.setValue(A=kappa*kronecker(mymesh)) 99 100 # ... set initial temperature .... 101 x=mymesh.getX() 102 103 qH=Scalar(0,FunctionOnBoundary(mymesh)) 104 qH.setTaggedValue("linebottom",qin) 105 mypde.setValue(q=whereZero(x[1]),r=Ti) 106 mypde.setValue(y=qH) 107 108 ###########################################################GET SOLUTION 109 T=mypde.getSolution() 110 111 ##########################################################VISUALISATION 112 # calculate gradient of solution for quiver plot 113 qu=-kappa*grad(T) 114 quT=qu.toListOfTuples() 115 116 #Projector is used to smooth the data. 117 proj=Projector(mymesh) 118 smthT=proj(T) 119 120 #move data to a regular grid for plotting 121 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth) 122 123 #Temperature Depth Profile along x[50] 124 cut=int(len(xi)/2) 125 pl.clf() 126 pl.plot(zi[:,cut],yi) 127 pl.title("Heat Refraction Temperature Depth Profile") 128 pl.xlabel("Temperature (K)") 129 pl.ylabel("Depth (m)") 130 if getMPIRankWorld() == 0: #check for MPI processing 131 pl.savefig(os.path.join(saved_path,"heatrefraction001_tdp.png")) 132 133 # Heat flow depth profile. 134 pl.clf() 135 # grid the data. 136 qu=proj(-kappa*grad(T)) 137 xiq,yiq,ziq = gradtoRegGrid(qu,mymesh,200,200,width,depth,1) 138 cut=int(len(xiq)/2) 139 pl.plot(ziq[:,cut]*1000.,yiq) 140 pl.title("Heat Flow Depth Profile") 141 pl.xlabel("Heat Flow (mW/m^2)") 142 pl.ylabel("Depth (m)") 143 if getMPIRankWorld() == 0: #check for MPI processing 144 pl.savefig(os.path.join(saved_path,"heatrefraction001_hf.png")) 145 146 # Temperature Gradient Depth Profile at x[50] 147 pl.clf() 148 zT=proj(-grad(T)) 149 xt,yt,zt=gradtoRegGrid(zT,mymesh,200,200,width,depth,1) 150 cut=int(len(xt)/2) 151 pl.plot(zt[:,cut]*1000.,yt) 152 pl.title("Heat Refraction Temperature Gradient \n Depth Profile") 153 pl.xlabel("Temperature (K/Km)") 154 pl.ylabel("Depth (m)") 155 if getMPIRankWorld() == 0: #check for MPI processing 156 pl.savefig(os.path.join(saved_path,"heatrefraction001_tgdp.png")) 157 158 # Thermal Conditions Depth Profile 159 pl.clf() 160 xk,yk,zk = toRegGrid(kappa,mymesh,200,200,width,depth) 161 cut=int(len(xk)/2) 162 pl.plot(zk[:,cut],yk) 163 pl.title("Heat Refraction Thermal Conductivity Depth Profile") 164 pl.xlabel("Conductivity (W/K/m)") 165 pl.ylabel("Depth (m)") 166 pl.axis([1,5,-6000,0]) 167 if getMPIRankWorld() == 0: #check for MPI processing 168 pl.savefig(os.path.join(saved_path,"heatrefraction001_tcdp.png")) 169 170 # contour the gridded data, 171 # plotting dots at the randomly spaced data points. 172 quivshape = [20,20] #quivers x and quivers y 173 # function to calculate quiver locations 174 qu,qulocs = toQuivLocs(quivshape,width,depth,qu) 175 176 # select colour 177 pl.matplotlib.pyplot.autumn() 178 pl.clf() 179 # plot polygons for boundaries 180 CKL = pl.fill(tpgx,tpgy,'brown',label='2 W/m/k',zorder=-1000) 181 CKM = pl.fill(bpgx,bpgy,'red',label='4 W/m/k',zorder=-1000) 182 # contour temperature 183 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 184 # labels and formatting 185 pl.clabel(CS, inline=1, fontsize=8) 186 pl.title("Heat Refraction across a clinal structure.") 187 pl.xlabel("Horizontal Displacement (m)") 188 pl.ylabel("Depth (m)") 189 pl.legend() 190 if getMPIRankWorld() == 0: #check for MPI processing 191 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png")) 192 193 #Quiver Plot qulocs -> tail location, qu -> quiver length/direction 194 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\ 195 angles='xy',color="white") 196 pl.title("Heat Refraction across a clinal structure\n with\ 197 gradient quivers.") 198 if getMPIRankWorld() == 0: #check for MPI processing 199 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png")) 200 201

 ViewVC Help Powered by ViewVC 1.1.26