doc/cookbook/twodheatdiff001.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
"""

############################################################FILE HEADER
# twodheatdiff002.py
# Model temperature diffusion between a granite intrusion and sandstone 
# country rock. This is a two dimensional problem with the granite as a
# heat source.

#######################################################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 the LinearPDE module as LinearPDE
from esys.escript.linearPDEs import LinearPDE 
# This imports the rectangle domain function from finley.
from esys.finley import Rectangle 
# A useful unit handling package which will make sure all our units
# match up in the equations under SI.
from esys.escript.unitsSI import * 
import pylab as pl #Plotting package.
import numpy as np #Array package.
import os #This package is necessary to handle saving our data.

#################################################ESTABLISHING VARIABLES
#PDE related
mx = 600*m #meters - model length
my = 600*m #meters - model width
ndx = 100 #mesh steps in x direction 
ndy = 100 #mesh steps in y direction
r = 200*m #meters - radius of intrusion
ic = [300, 0] #centre of intrusion (meters)
q=0.*Celsius #our heat source temperature is now zero

## Intrusion Variables - Granite
Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
rhoi = 2750*kg/m**3 #kg/m^{3} density of granite
cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity
rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT
eta=0.  # RADIATION CONDITION
kappai=2.2*W/m/K #watts/m.K thermal conductivity
## Country Rock Variables - Sandstone
Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
rhoc = 2000*kg/m**3 #kg/m^{3} density
cpc = 920.*J/(kg*K) #j/kg.k specific heat
rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT
kappac = 1.9*W/m/K #watts/m.K thermal conductivity

#Script/Iteration Related
t=0. #our start time, usually zero
tday=100*365. #the time we want to end the simulation in days
tend=tday*24*60*60
outputs = 200 # number of time steps required.
h=(tend-t)/outputs #size of time step
#user warning
print "Expected Number of Output Files is: ", outputs
print "Step size is: ", h/(24.*60*60), "days"
i=0 #loop counter 
#the folder to put our outputs in, leave blank "" for script path 
save_path="data/twodheatdiff"
########## note this folder path must exist to work ###################

################################################ESTABLISHING PARAMETERS
#generate domain using rectangle
model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
#extract finite points - the solution points
x=model.getX()
#create the PDE
mypde=LinearPDE(model) #assigns a domain to our PDE
mypde.setSymmetryOn() #set the fast solver on for symmetry
#establish location of boundary between two materials
bound = length(x-ic)-r #where the boundary will be located
A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)
D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)
#define our PDE coeffs
mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)
#set initial temperature
T= Ti*whereNegative(bound)+Tc*wherePositive(bound)

# rearrage mymesh to suit solution function space for contouring      
oldspacecoords=model.getX()
coords=Data(oldspacecoords, T.getFunctionSpace())
coords = np.array(coords.toListOfTuples())
coordX = coords[:,0]
coordY = coords[:,1]
# create regular grid
xi = np.linspace(0.0,600.0,100)
yi = np.linspace(0.0,600.0,100)

#... start iteration:
while t<=tend:
      i+=1
      t+=h
      Y = T*D
      mypde.setValue(Y=Y)
      T=mypde.getSolution()
      tempT = T.toListOfTuples(scalarastuple=False)
      # grid the data.
      zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
      # contour the gridded data, plotting dots at the randomly spaced data points.
      pl.matplotlib.pyplot.autumn()
      pl.contourf(xi,yi,zi,10)
      CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
      pl.clabel(CS, inline=1, fontsize=8)
      pl.axis([0,600,0,600])
      pl.title("Heat diffusion from an intrusion.")
      pl.xlabel("Horizontal Displacement (m)")
      pl.ylabel("Depth (m)")
      pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)
      pl.clf()            
#     saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)

# compile the *.png files to create an *.avi video that shows T change
# with time. This opperation uses linux mencoder.
os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
twodheatdiff001tempT.avi")
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	############################################################FILE HEADER
27	# twodheatdiff002.py
28	# Model temperature diffusion between a granite intrusion and sandstone
29	# country rock. This is a two dimensional problem with the granite as a
30	# heat source.
31
32	#######################################################EXTERNAL MODULES
33	#To solve the problem it is necessary to import the modules we require.
34	#This imports everything from the escript library
35	from esys.escript import *
36	# This defines the LinearPDE module as LinearPDE
37	from esys.escript.linearPDEs import LinearPDE
38	# This imports the rectangle domain function from finley.
39	from esys.finley import Rectangle
40	# A useful unit handling package which will make sure all our units
41	# match up in the equations under SI.
42	from esys.escript.unitsSI import *
43	import pylab as pl #Plotting package.
44	import numpy as np #Array package.
45	import os #This package is necessary to handle saving our data.
46
47	#################################################ESTABLISHING VARIABLES
48	#PDE related
49	mx = 600*m #meters - model length
50	my = 600*m #meters - model width
51	ndx = 100 #mesh steps in x direction
52	ndy = 100 #mesh steps in y direction
53	r = 200*m #meters - radius of intrusion
54	ic = [300, 0] #centre of intrusion (meters)
55	q=0.*Celsius #our heat source temperature is now zero
56
57	## Intrusion Variables - Granite
58	Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
59	rhoi = 2750kg/m*3 #kg/m^{3} density of granite
60	cpi = 790.J/(kgK) #j/Kg.K thermal capacity
61	rhocpi = rhoicpi #DENSITY SPECIFIC HEAT
62	eta=0. # RADIATION CONDITION
63	kappai=2.2*W/m/K #watts/m.K thermal conductivity
64	## Country Rock Variables - Sandstone
65	Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
66	rhoc = 2000kg/m*3 #kg/m^{3} density
67	cpc = 920.J/(kgK) #j/kg.k specific heat
68	rhocpc = rhoccpc #DENSITY SPECIFIC HEAT
69	kappac = 1.9*W/m/K #watts/m.K thermal conductivity
70
71	#Script/Iteration Related
72	t=0. #our start time, usually zero
73	tday=100*365. #the time we want to end the simulation in days
74	tend=tday2460*60
75	outputs = 200 # number of time steps required.
76	h=(tend-t)/outputs #size of time step
77	#user warning
78	print "Expected Number of Output Files is: ", outputs
79	print "Step size is: ", h/(24.6060), "days"
80	i=0 #loop counter
81	#the folder to put our outputs in, leave blank "" for script path
82	save_path="data/twodheatdiff"
83	########## note this folder path must exist to work ###################
84
85	################################################ESTABLISHING PARAMETERS
86	#generate domain using rectangle
87	model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
88	#extract finite points - the solution points
89	x=model.getX()
90	#create the PDE
91	mypde=LinearPDE(model) #assigns a domain to our PDE
92	mypde.setSymmetryOn() #set the fast solver on for symmetry
93	#establish location of boundary between two materials
94	bound = length(x-ic)-r #where the boundary will be located
95	A = (kappai)whereNegative(bound)+(kappac)wherePositive(bound)
96	D = (rhocpi/h)whereNegative(bound)+(rhocpc/h)wherePositive(bound)
97	#define our PDE coeffs
98	mypde.setValue(A=Akronecker(model),D=D,d=eta,y=etaTc)
99	#set initial temperature
100	T= TiwhereNegative(bound)+TcwherePositive(bound)
101
102	# rearrage mymesh to suit solution function space for contouring
103	oldspacecoords=model.getX()
104	coords=Data(oldspacecoords, T.getFunctionSpace())
105	coords = np.array(coords.toListOfTuples())
106	coordX = coords[:,0]
107	coordY = coords[:,1]
108	# create regular grid
109	xi = np.linspace(0.0,600.0,100)
110	yi = np.linspace(0.0,600.0,100)
111
112	#... start iteration:
113	while t<=tend:
114	i+=1
115	t+=h
116	Y = T*D
117	mypde.setValue(Y=Y)
118	T=mypde.getSolution()
119	tempT = T.toListOfTuples(scalarastuple=False)
120	# grid the data.
121	zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
122	# contour the gridded data, plotting dots at the randomly spaced data points.
123	pl.matplotlib.pyplot.autumn()
124	pl.contourf(xi,yi,zi,10)
125	CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
126	pl.clabel(CS, inline=1, fontsize=8)
127	pl.axis([0,600,0,600])
128	pl.title("Heat diffusion from an intrusion.")
129	pl.xlabel("Horizontal Displacement (m)")
130	pl.ylabel("Depth (m)")
131	pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)
132	pl.clf()
133	# saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
134
135	# compile the .png files to create an .avi video that shows T change
136	# with time. This opperation uses linux mencoder.
137	os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
138	w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
139	twodheatdiff001tempT.avi")