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 *

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

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