1 

2 
######################################################## 
3 
# 
4 
# Copyright (c) 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/osl3.0.php 
11 
# 
12 
######################################################## 
13 

14 
__copyright__="""Copyright (c) 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/osl3.0.php""" 
20 
__url__="https://launchpad.net/escriptfinley" 
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 # This defines LinearPDE as LinearPDE 
29 
from esys.finley import Rectangle # This imports the rectangle domain function from finley 
30 
import os #This package is necessary to handle saving our data. 
31 

32 
from cblib import needdirs 
33 

34 

35 
##ESTABLISHING VARIABLES 
36 
#PDE related 
37 
mx = 600 # model lenght 
38 
my = 600 # model width 
39 
ndx = 100 # steps in x direction 
40 
ndy = 100 # steps in y direction 
41 
r = 200 # radius of intrusion 
42 
ic = [300, 0] #centre of intrusion 
43 

44 
q=0 #our heat source temperature is now zero 
45 
Ti=2273 # Kelvin #the starting temperature of our iron bar 
46 
rhoi = 2750 #kg/m^{3} density 
47 
cpi = 790 #j/kg specific heat 
48 
rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT 
49 
eta=0. # RADIATION CONDITION 
50 
kappai=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY 
51 

52 
Tc = 200 
53 
rhoc = 2200 
54 
cpc = 400 
55 
rhocpc = rhoc*cpc 
56 
kappac = 0.1 
57 

58 

59 
#Script/Iteration Related 
60 
t=0. #our start time, usually zero 
61 
tday=100*365. #the time we want to end the simulation in days 
62 
tend=tday*24*60*60 
63 
outputs = 200 # number of time steps required. 
64 
h=(tendt)/outputs #size of time step 
65 

66 
print "Expected Number of Output Files is: ", outputs 
67 
print "Step size is: ", h/(24.*60*60), "days" 
68 

69 

70 
i=0 #loop counter 
71 
save_path = "data/twodheatdiff" #the folder to put our outputs in, leave blank "" for script path  note this folder path must exist to work 
72 
needdirs([save_path]) 
73 

74 
#... generate domain ... 
75 
model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
76 
# extract finite points 
77 
x=model.getX() 
78 

79 
#... open PDE ... 
80 
mypde=LinearPDE(model) 
81 
mypde.setSymmetryOn() 
82 

83 
bound = length(xic)r #where the boundary will be located 
84 

85 
A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound) 
86 
D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound) 
87 

88 
mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc) 
89 

90 
# ... set initial temperature .... 
91 

92 
T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures. 
93 
saveVTK(os.path.join(save_path,"dataedge.vtu"), sol=bound) 
94 
saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T) 
95 

96 
#... start iteration: 
97 
while t<=tend: 
98 
i+=1 
99 
t+=h 
100 
Y = T*D 
101 
mypde.setValue(Y=Y) 
102 
T=mypde.getSolution() 
103 
saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T) 