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

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

Parent Directory Parent Directory | Revision Log Revision Log


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 # 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 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 # 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