examples/cookbook/heatrefraction001.py


########################################################
#
# Copyright (c) 2009 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2009 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

"""
Author: Antony Hallam antony.hallam@uqconnect.edu.au
"""

############################################################FILE HEADER
# heatrefraction001.py
# Model steady state temperature distribution in two block model, mesh
# from heatrefraction_mesher001.py 

#######################################################EXTERNAL MODULES
# To solve the problem it is necessary to import the modules we 
# require.


# This imports everything from the escript library
from esys.escript import * 
# This defines LinearPDE as LinearPDE
from esys.escript.linearPDEs import LinearPDE, Poisson
# smoothing operator 
from esys.escript.pdetools import Projector
# This imports the rectangle domain function from finley
from esys.finley import Rectangle, ReadMesh, Domain 
# This package is necessary to handle saving our data.
import os
# A useful unit handling package which will make sure all our units
# match up in the equations.
from esys.escript.unitsSI import * 
# numpy for array handling
import numpy as np

import matplotlib
#For interactive use, you can comment out the next line
matplotlib.use('agg') #It's just here for automated testing

# pylab for matplotlib and plotting
import pylab as pl
# cblib functions
from cblib2 import toQuivLocs, toXYTuple, toRegGrid, gradtoRegGrid
from cblib1 import needdirs

########################################################MPI WORLD CHECK
if getMPISizeWorld() > 1:
    import sys
    print "This example will not run in an MPI world."
    sys.exit(0)

#################################################ESTABLISHING VARIABLES
qin=70*Milli*W/(m*m) #our heat source temperature is now zero
Ti=290.15*K # Kelvin #the starting temperature of our iron bar
width=5000.0*m
depth=-6000.0*m

# the folder to gett our outputs from, leave blank "" for script path - 
# note these depen. are generated from heatrefraction_mesher001.py
saved_path = save_path= os.path.join("data","heatrefrac001" )
needdirs([saved_path])

################################################ESTABLISHING PARAMETERS
## DOMAIN
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
tpgx = tpg[:,0]
tpgy = tpg[:,1]
bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
bpgx = bpg[:,0]
bpgy = bpg[:,1]

# set up kappa (thermal conductivity across domain) using tags
kappa=Scalar(0,Function(mymesh))
kappa.setTaggedValue("top",2.0)
kappa.setTaggedValue("bottom",4.0)

#... generate functionspace...
#... open PDE ...
mypde=LinearPDE(mymesh)
#define first coefficient
mypde.setValue(A=kappa*kronecker(mymesh))

# ... set initial temperature ....
x=mymesh.getX()

qH=Scalar(0,FunctionOnBoundary(mymesh))
qH.setTaggedValue("linebottom",qin)
mypde.setValue(q=whereZero(x[1]),r=Ti)
mypde.setValue(y=qH)

###########################################################GET SOLUTION
T=mypde.getSolution()

##########################################################VISUALISATION
# calculate gradient of solution for quiver plot
qu=-kappa*grad(T)
quT=qu.toListOfTuples()

#Projector is used to smooth the data.
proj=Projector(mymesh)
smthT=proj(T)

#move data to a regular grid for plotting
xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)

#Temperature Depth Profile along x[50]
cut=int(len(xi)/2)
pl.clf()
pl.plot(zi[:,cut],yi)
pl.title("Heat Refraction Temperature Depth Profile")
pl.xlabel("Temperature (K)")
pl.ylabel("Depth (m)")
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_tdp.png"))
    
# Heat flow depth profile.
pl.clf()
# grid the data.
qu=proj(-kappa*grad(T))
xiq,yiq,ziq = gradtoRegGrid(qu,mymesh,200,200,width,depth,1)
cut=int(len(xiq)/2)
pl.plot(ziq[:,cut]*1000.,yiq)
pl.title("Heat Flow Depth Profile")
pl.xlabel("Heat Flow (mW/m^2)")
pl.ylabel("Depth (m)")
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_hf.png"))

# Temperature Gradient Depth Profile at x[50]
pl.clf()
zT=proj(-grad(T))
xt,yt,zt=gradtoRegGrid(zT,mymesh,200,200,width,depth,1)
cut=int(len(xt)/2)
pl.plot(zt[:,cut]*1000.,yt)
pl.title("Heat Refraction Temperature Gradient \n Depth Profile")
pl.xlabel("Temperature (K/Km)")
pl.ylabel("Depth (m)")
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_tgdp.png"))

# Thermal Conditions Depth Profile    
pl.clf()
xk,yk,zk = toRegGrid(kappa,mymesh,200,200,width,depth)
cut=int(len(xk)/2)
pl.plot(zk[:,cut],yk)
pl.title("Heat Refraction Thermal Conductivity Depth Profile")
pl.xlabel("Conductivity (W/K/m)")
pl.ylabel("Depth (m)")
pl.axis([1,5,-6000,0])
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_tcdp.png"))

# contour the gridded data, 
# plotting dots at the randomly spaced data points.
quivshape = [20,20] #quivers x and quivers y
# function to calculate quiver locations
qu,qulocs = toQuivLocs(quivshape,width,depth,qu)

# select colour
pl.matplotlib.pyplot.autumn()
pl.clf()
# plot polygons for boundaries
CKL = pl.fill(tpgx,tpgy,'brown',label='2 W/m/k',zorder=-1000)
CKM = pl.fill(bpgx,bpgy,'red',label='4 W/m/k',zorder=-1000)
# contour temperature
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
# labels and formatting
pl.clabel(CS, inline=1, fontsize=8)
pl.title("Heat Refraction across a clinal structure.")
pl.xlabel("Horizontal Displacement (m)")
pl.ylabel("Depth (m)")
pl.legend()
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))

#Quiver Plot qulocs -> tail location, qu -> quiver length/direction
QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
                angles='xy',color="white")
pl.title("Heat Refraction across a clinal structure\n with\
gradient quivers.")
if getMPIRankWorld() == 0: #check for MPI processing
    pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))


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