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
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	qin=300.MilliW/(m*m) #our heat source temperature is now zero
39	Ti=290.15*K # Kelvin #the starting temperature of our iron bar
40
41	# 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	# set up kappa (thermal conductivity across domain) using tags
60	kappa=Scalar(0,Function(mymesh))
61	kappa.setTaggedValue("top",2.0)
62	kappa.setTaggedValue("bottomleft",18.0)
63	kappa.setTaggedValue("bottomright",6.0)
64
65	#... generate functionspace...
66	#... open PDE ...
67	mypde=LinearPDE(mymesh)
68	mypde.setSymmetryOn()
69	#define first coefficient
70	mypde.setValue(A=kappa*kronecker(mymesh))
71
72	# ... set initial temperature ....
73	x=mymesh.getX()
74
75	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
80	# get steady state solution and export to vtk.
81	T=mypde.getSolution()
82	#saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
83
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	CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
131	#~ CK = pl.contourf(xi,yi,ziK,2)
132	CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
133	pl.clabel(CS, inline=1, fontsize=8)
134	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
140	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"))