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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 5366 byte(s)
Renamed the main cookbook tex file to match our convention.
Replaced doc/cookbook/figures/heatrefraction002contqu.pdf with
a version which is actually pdf. However it needs to be regnerated since
it it sideways.

The examples have had their copyright notices fixed (dates were too early).
sb2.py has been removed since it uses pyvisi.

scons will now build the cookbook as parts of a docs build.
Also in reposnse to :
scons cookbook_pdf


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 # 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 # This imports the rectangle domain function from finley
42 from esys.finley import Rectangle, ReadMesh, Domain
43 # This package is necessary to handle saving our data.
44 import os
45 # A useful unit handling package which will make sure all our units
46 # match up in the equations.
47 from esys.escript.unitsSI import *
48 # numpy for array handling
49 import numpy as np
50
51 import matplotlib
52 #For interactive use, you can comment out the next line
53 matplotlib.use('agg') #It's just here for automated testing
54
55 # pylab for matplotlib and plotting
56 import pylab as pl
57 # cblib functions
58 from cblib import toQuivLocs, toXYTuple, needdirs
59
60 ########################################################MPI WORLD CHECK
61 if getMPISizeWorld() > 1:
62 import sys
63 print "This example will not run in an MPI world."
64 sys.exit(0)
65
66 #################################################ESTABLISHING VARIABLES
67 qin=70*Milli*W/(m*m) #our heat source temperature is now zero
68 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
69 width=5000.0*m
70 depth=-6000.0*m
71
72 # the folder to gett our outputs from, leave blank "" for script path -
73 # note these depen. are generated from heatrefraction_mesher001.py
74 saved_path = save_path= os.path.join("data","heatrefrac001" )
75 needdirs([saved_path])
76
77 ################################################ESTABLISHING PARAMETERS
78 ## DOMAIN
79 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
80 tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
81 tpgx = tpg[:,0]
82 tpgy = tpg[:,1]
83 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
84 bpgx = bpg[:,0]
85 bpgy = bpg[:,1]
86
87 # set up kappa (thermal conductivity across domain) using tags
88 kappa=Scalar(0,Function(mymesh))
89 kappa.setTaggedValue("top",2.0)
90 kappa.setTaggedValue("bottom",4.0)
91
92 #... generate functionspace...
93 #... open PDE ...
94 mypde=LinearPDE(mymesh)
95 #define first coefficient
96 mypde.setValue(A=kappa*kronecker(mymesh))
97
98 # ... set initial temperature ....
99 x=mymesh.getX()
100
101 qH=Scalar(0,FunctionOnBoundary(mymesh))
102 qH.setTaggedValue("linebottom",qin)
103 mypde.setValue(q=whereZero(x[1]),r=Ti)
104 mypde.setValue(y=qH)
105
106 ###########################################################GET SOLUTION
107 T=mypde.getSolution()
108
109 ##########################################################VISUALISATION
110 # calculate gradient of solution for quiver plot
111 qu=-kappa*grad(T)
112
113 # rearrage mymesh to suit solution function space
114 oldspacecoords=mymesh.getX()
115 coords=Data(oldspacecoords, T.getFunctionSpace())
116
117 quivshape = [20,20] #quivers x and quivers y
118 # function to calculate quiver locations
119 qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
120
121 kappaT = kappa.toListOfTuples(scalarastuple=False)
122 coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
123 coordKX, coordKY = toXYTuple(coordsK)
124
125 tempT = T.toListOfTuples(scalarastuple=False)
126 coordX, coordY = toXYTuple(coords)
127
128 xi = np.linspace(0.0,width,100)
129 yi = np.linspace(depth,0.0,100)
130 # grid the data.
131 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
132 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
133 # contour the gridded data,
134 # plotting dots at the randomly spaced data points.
135
136 # select colour
137 pl.matplotlib.pyplot.autumn()
138 # plot polygons for boundaries
139 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
140 # contour temperature
141 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
142 # labels and formatting
143 pl.clabel(CS, inline=1, fontsize=8)
144 pl.title("Heat Refraction across a clinal structure.")
145 pl.xlabel("Horizontal Displacement (m)")
146 pl.ylabel("Depth (m)")
147 if getMPIRankWorld() == 0: #check for MPI processing
148 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
149
150 #Quiver Plot qulocs -> tail location, qu -> quiver length/direction
151 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
152 angles='xy',color="white")
153 pl.title("Heat Refraction across a clinal structure \n with\
154 gradient quivers.")
155 if getMPIRankWorld() == 0: #check for MPI processing
156 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))

  ViewVC Help
Powered by ViewVC 1.1.26