examples/cookbook/heatrefraction002.py


########################################################
#
# Copyright (c) 2003-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) 2003-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
"""

# To solve the problem it is necessary to import the modules we require.
from esys.escript import * # This imports everything from the escript library
from esys.escript.linearPDEs import LinearPDE, Poisson # This defines LinearPDE as LinearPDE
from esys.finley import Rectangle, ReadMesh, Domain # This imports the rectangle domain function from finley
import os #This package is necessary to handle saving our data.
from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
import numpy as np
import pylab as pl

from esys.escript.pdetools import *
from cblib import needdirs

##ESTABLISHING VARIABLES
qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
Ti=290.15*K # Kelvin #the starting temperature of our iron bar

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

###### 2 BLOCK MODEL #########
## DOMAIN
## Anticline
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))
tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
tpgx = tpg[:,0]
tpgy = tpg[:,1]
bpgl = np.loadtxt(os.path.join(saved_path,"botpgl"))
bpglx = bpgl[:,0]
bpgly = bpgl[:,1]
bpgr = np.loadtxt(os.path.join(saved_path,"botpgr"))
bpgrx = bpgr[:,0]
bpgry = bpgr[:,1]

# set up kappa (thermal conductivity across domain) using tags
kappa=Scalar(0,Function(mymesh))
kappa.setTaggedValue("top",2.0)
kappa.setTaggedValue("bottomleft",18.0)
kappa.setTaggedValue("bottomright",6.0)

#... generate functionspace...
#... open PDE ...
mypde=LinearPDE(mymesh)
mypde.setSymmetryOn()
#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 steady state solution and export to vtk.
T=mypde.getSolution()
#saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))

# rearrage mymesh to suit solution function space      
oldspacecoords=mymesh.getX()
coords=Data(oldspacecoords, T.getFunctionSpace())

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

# create quiver locations
quivshape = [20,20] #quivers x and quivers y
numquiv = quivshape[0]*quivshape[1] # total number of quivers
maxx = 5000. # maximum x displacement
dx = maxx/quivshape[0]+1. # quiver x spacing
maxy = -6000. # maximum y displacement
dy = maxy/quivshape[1]+1. # quiver y spacing
qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
# fill qulocs
for i in range(0,quivshape[0]-1):
    for j in range(0,quivshape[1]-1):
        qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
# retreive values for quivers direction and shape from qu
quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
qu = quL(qu) #array of dx,dy for quivers
qu = np.array(qu) #turn into a numpy array
qulocs = quL.getX() #returns actual locations from data
qulocs = np.array(qulocs) #turn into a numpy array

kappaT = kappa.toListOfTuples(scalarastuple=False)
coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
coordsK = np.array(coordsK.toListOfTuples())
coordKX = coordsK[:,0]
coordKY = coordsK[:,1]
      
coords = np.array(coords.toListOfTuples())
tempT = T.toListOfTuples(scalarastuple=False)

coordX = coords[:,0]
coordY = coords[:,1]

xi = np.linspace(0.0,5000.0,100)
yi = np.linspace(-6000.0,0.0,100)
# grid the data.
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
# contour the gridded data, plotting dots at the randomly spaced data points.

pl.matplotlib.pyplot.autumn()
CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
#~ CK = pl.contourf(xi,yi,ziK,2)
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
pl.clabel(CS, inline=1, fontsize=8)
pl.title("Heat Refraction across an anisotropic structure.")
pl.xlabel("Horizontal Displacement (m)")
pl.ylabel("Depth (m)")
#~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))

QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
pl.title("Heat Refraction across an anisotropic structure \n with gradient quivers.")
pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
1	ahallam	2588
2			########################################################
3			#
4			# Copyright (c) 2003-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) 2003-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			# To solve the problem it is necessary to import the modules we require.
27			from esys.escript import * # This imports everything from the escript library
28			from esys.escript.linearPDEs import LinearPDE, Poisson # This defines LinearPDE as LinearPDE
29			from esys.finley import Rectangle, ReadMesh, Domain # This imports the rectangle domain function from finley
30			import os #This package is necessary to handle saving our data.
31			from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
32			import numpy as np
33			import pylab as pl
34
35			from esys.escript.pdetools import *
36	jfenwick	2648	from cblib import needdirs
37	ahallam	2588
38			##ESTABLISHING VARIABLES
39	ahallam	2597	qin=300.MilliW/(m*m) #our heat source temperature is now zero
40	ahallam	2588	Ti=290.15*K # Kelvin #the starting temperature of our iron bar
41
42	ahallam	2597	# the folder to gett our outputs from, leave blank "" for script path -
43			# note these depen. are generated from heatrefraction_mesher001.py
44			saved_path = "data/heatrefrac002"
45	jfenwick	2648	needdirs([saved_path])
46	ahallam	2597
47			###### 2 BLOCK MODEL #########
48			## DOMAIN
49			## Anticline
50			mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))
51			tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
52			tpgx = tpg[:,0]
53			tpgy = tpg[:,1]
54			bpgl = np.loadtxt(os.path.join(saved_path,"botpgl"))
55			bpglx = bpgl[:,0]
56			bpgly = bpgl[:,1]
57			bpgr = np.loadtxt(os.path.join(saved_path,"botpgr"))
58			bpgrx = bpgr[:,0]
59			bpgry = bpgr[:,1]
60
61	ahallam	2588	# set up kappa (thermal conductivity across domain) using tags
62			kappa=Scalar(0,Function(mymesh))
63	ahallam	2597	kappa.setTaggedValue("top",2.0)
64			kappa.setTaggedValue("bottomleft",18.0)
65			kappa.setTaggedValue("bottomright",6.0)
66	ahallam	2588
67			#... generate functionspace...
68			#... open PDE ...
69			mypde=LinearPDE(mymesh)
70	ahallam	2597	mypde.setSymmetryOn()
71	ahallam	2588	#define first coefficient
72			mypde.setValue(A=kappa*kronecker(mymesh))
73
74			# ... set initial temperature ....
75			x=mymesh.getX()
76
77	ahallam	2597	qH=Scalar(0,FunctionOnBoundary(mymesh))
78			qH.setTaggedValue("linebottom",qin)
79			mypde.setValue(q=whereZero(x[1]),r=Ti)
80			mypde.setValue(y=qH)
81	ahallam	2588
82			# get steady state solution and export to vtk.
83			T=mypde.getSolution()
84	ahallam	2597	#saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
85	ahallam	2588
86			# rearrage mymesh to suit solution function space
87			oldspacecoords=mymesh.getX()
88			coords=Data(oldspacecoords, T.getFunctionSpace())
89
90			# calculate gradient of solution for quiver plot
91			qu=-kappa*grad(T)
92
93			# create quiver locations
94			quivshape = [20,20] #quivers x and quivers y
95			numquiv = quivshape[0]*quivshape[1] # total number of quivers
96			maxx = 5000. # maximum x displacement
97			dx = maxx/quivshape[0]+1. # quiver x spacing
98			maxy = -6000. # maximum y displacement
99			dy = maxy/quivshape[1]+1. # quiver y spacing
100			qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
101			# fill qulocs
102			for i in range(0,quivshape[0]-1):
103			for j in range(0,quivshape[1]-1):
104			qulocs[iquivshape[0]+j,:] = [dx+dxi,dy+dy*j]
105			# retreive values for quivers direction and shape from qu
106			quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
107			qu = quL(qu) #array of dx,dy for quivers
108			qu = np.array(qu) #turn into a numpy array
109			qulocs = quL.getX() #returns actual locations from data
110			qulocs = np.array(qulocs) #turn into a numpy array
111
112			kappaT = kappa.toListOfTuples(scalarastuple=False)
113			coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
114			coordsK = np.array(coordsK.toListOfTuples())
115			coordKX = coordsK[:,0]
116			coordKY = coordsK[:,1]
117
118			coords = np.array(coords.toListOfTuples())
119			tempT = T.toListOfTuples(scalarastuple=False)
120
121			coordX = coords[:,0]
122			coordY = coords[:,1]
123
124			xi = np.linspace(0.0,5000.0,100)
125			yi = np.linspace(-6000.0,0.0,100)
126			# grid the data.
127			zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
128			ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
129			# contour the gridded data, plotting dots at the randomly spaced data points.
130
131			pl.matplotlib.pyplot.autumn()
132	ahallam	2597	CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
133			#~ CK = pl.contourf(xi,yi,ziK,2)
134	ahallam	2588	CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
135			pl.clabel(CS, inline=1, fontsize=8)
136	ahallam	2597	pl.title("Heat Refraction across an anisotropic structure.")
137			pl.xlabel("Horizontal Displacement (m)")
138			pl.ylabel("Depth (m)")
139			#~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
140			pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
141	ahallam	2588
142	ahallam	2597	QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
143			pl.title("Heat Refraction across an anisotropic structure \n with gradient quivers.")
144			pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))