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

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

```
 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 # heatrefraction002.py 28 # Model steady state temperature distribution in three block model, 29 # mesh from heatrefraction_mesher002.py 30 31 #######################################################EXTERNAL MODULES 32 # To solve the problem it is necessary to import the modules we 33 # require. 34 # This imports everything from the escript library 35 from esys.escript import * 36 # This defines LinearPDE as LinearPDE 37 from esys.escript.linearPDEs import LinearPDE, Poisson 38 # This imports the rectangle domain function from finley 39 from esys.finley import Rectangle, ReadMesh, Domain 40 # This package is necessary to handle saving our data. 41 import os 42 # A useful unit handling package which will make sure all our units 43 # match up in the equations. 44 from esys.escript.unitsSI import * 45 # numpy for array handling 46 import numpy as np 47 import matplotlib 48 49 #For interactive use, you can comment out the next line 50 matplotlib.use('agg') #It's just here for automated testing 51 52 # pylab for matplotlib and plotting 53 import pylab as pl 54 # cblib functions 55 from cblib2 import toQuivLocs, toXYTuple 56 from cblib1 import needdirs 57 58 ########################################################MPI WORLD CHECK 59 if getMPISizeWorld() > 1: 60 import sys 61 print "This example will not run in an MPI world." 62 sys.exit(0) 63 64 #################################################ESTABLISHING VARIABLES 65 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero 66 Ti=290.15*K # Kelvin #the starting temperature of our iron bar 67 width=5000.0*m 68 depth=-6000.0*m 69 70 # the folder to gett our outputs from, leave blank "" for script path - 71 # note these depen. are generated from heatrefraction_mesher001.py 72 saved_path = save_path= os.path.join("data","heatrefrac002") 73 needdirs([saved_path]) 74 75 ################################################ESTABLISHING PARAMETERS 76 ## DOMAIN 77 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly")) 78 tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 79 tpgx = tpg[:,0] 80 tpgy = tpg[:,1] 81 bpgl = np.loadtxt(os.path.join(saved_path,"botpgl")) 82 bpglx = bpgl[:,0] 83 bpgly = bpgl[:,1] 84 bpgr = np.loadtxt(os.path.join(saved_path,"botpgr")) 85 bpgrx = bpgr[:,0] 86 bpgry = bpgr[:,1] 87 88 # set up kappa (thermal conductivity across domain) using tags 89 kappa=Scalar(0,Function(mymesh)) 90 kappa.setTaggedValue("top",2.0) 91 kappa.setTaggedValue("bottomleft",4.0) 92 kappa.setTaggedValue("bottomright",6.0) 93 94 #... generate functionspace... 95 #... open PDE ... 96 mypde=LinearPDE(mymesh) 97 mypde.setSymmetryOn() 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 116 # rearrage mymesh to suit solution function space 117 oldspacecoords=mymesh.getX() 118 coords=Data(oldspacecoords, T.getFunctionSpace()) 119 120 quivshape = [20,20] #quivers x and quivers y 121 # function to calculate quiver locations 122 qu,qulocs = toQuivLocs(quivshape,width,depth,qu) 123 124 kappaT = kappa.toListOfTuples(scalarastuple=False) 125 coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 126 coordKX, coordKY = toXYTuple(coordsK) 127 128 tempT = T.toListOfTuples(scalarastuple=False) 129 coordX, coordY = toXYTuple(coords) 130 131 xi = np.linspace(0.0,width,100) 132 yi = np.linspace(depth,0.0,100) 133 # grid the data. 134 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 135 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 136 # contour the gridded data, plotting dots at the randomly 137 # spaced data points. 138 139 pl.matplotlib.pyplot.autumn() 140 CKL = pl.fill(tpgx,tpgy,'brown',label='2 W/m/k',zorder=-1000) 141 CKM = pl.fill(bpglx,bpgly,'red',label='4 W/m/k',zorder=-1000) 142 CKN = pl.fill(bpgrx,bpgry,'orange',label='6 W/m/k',zorder=-1000) 143 144 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 145 pl.clabel(CS, inline=1, fontsize=8) 146 pl.title("Heat Refraction across an anisotropic structure.") 147 pl.xlabel("Horizontal Displacement (m)") 148 pl.ylabel("Depth (m)") 149 pl.legend() 150 if getMPIRankWorld() == 0: #check for MPI processing 151 pl.savefig(os.path.join(saved_path,"heatrefraction002_cont.png")) 152 153 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\ 154 angles='xy',color="white") 155 pl.title("Heat Refraction across an anisotropic structure\ 156 \n with gradient quivers.") 157 if getMPIRankWorld() == 0: #check for MPI processing 158 pl.savefig(os.path.join(saved_path,"heatrefraction002_contqu.png"))

 ViewVC Help Powered by ViewVC 1.1.26