examples/cookbook/onedheatdiff002.py


########################################################
#
# Copyright (c) 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) 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
# onedheatdiff002.py
# Model temperature diffusion between two granite blocks. This is a one
# dimensional problem with no heat source and a single heat disparity.

#######################################################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 *
#For interactive use, you can comment out the next two lines
import matplotlib
matplotlib.use('agg') #It's just here for automated testing
import pylab as pl #Plotting package.
import numpy as np #Array package.
import os #This package is necessary to handle saving our data.
from cblib1 import needdirs

########################################################MPI WORLD CHECK
if getMPISizeWorld() > 1:
    import sys
    print "This example will not run in an MPI world."
    sys.exit(0)

#################################################ESTABLISHING VARIABLES
#PDE related
mx = 500*m #meters - model length
my = 100*m #meters - model width
ndx = 500 # mesh steps in x direction 
ndy = 1 # mesh steps in y direction
boundloc = mx/2 # location of boundary between two blocks
q=0.*Celsius #our heat source temperature is now zero
Tref=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
rho = 2750*kg/m**3 #kg/m^{3} density of granite
cp = 790.*J/(kg*K) #j/Kg.K thermal capacity
rhocp = rho*cp  #DENSITY * SPECIFIC HEAT
eta=0.  # RADIATION CONDITION
kappa=2.2*W/m/K #watts/m.K thermal conductivity

#Script/Iteration Related
t=0. #our start time, usually zero
tend=10*yr #the time we want to end the simulation in years
outputs = 40 # number of time steps required.
h=(tend-t)/outputs #size of time step
#user warning statement
print "Expected Number of Output Files is: ", outputs
print "Step size is: ", h/day, "days"
i=0 #loop counter 
#the folder to put our outputs in, leave blank "" for script path
save_path= os.path.join("data","onedheatdiff002")
needdirs([save_path])
########## 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
#define our PDE coeffs
mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
#establish location of boundary between two blocks
bound = x[0]-boundloc
#set initial temperature
T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)

#convert solution points for plotting
plx = x.toListOfTuples() 
plx = np.array(plx) #convert to tuple to numpy array
plx = plx[:,0] #extract x locations

########################################################START ITERATION
while t<=tend:
    i+=1 #increment the counter
    t+=h #increment the current time
        print "time step ",i," at time  ",t/day," days."
    mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
    T=mypde.getSolution() #find temperature solution
    #set up for plotting
    tempT = T.toListOfTuples(scalarastuple=False)
    pl.figure(1)
    pl.plot(plx,tempT)
    pl.axis([0,500,0,2500])
    pl.title("Temperature accross Interface")
    if getMPIRankWorld() == 0:
        pl.savefig(os.path.join(save_path,"intpyplot%03d.png" %i))
    pl.clf()
    
# compile the *.png files to create two *.avi videos that show T change
# with time. This opperation uses linux mencoder. For other operating 
# systems it may be possible to use your favourite video compiler to
# convert image files to videos. To enable this step uncomment the
# following lines.

#os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
#w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
#onedheatdiff002tempT.avi")
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/osl-3.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/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	# onedheatdiff002.py
28	# Model temperature diffusion between two granite blocks. This is a one
29	# dimensional problem with no heat source and a single heat disparity.
30
31	#######################################################EXTERNAL MODULES
32	#To solve the problem it is necessary to import the modules we require.
33	#This imports everything from the escript library
34	from esys.escript import *
35	# This defines the LinearPDE module as LinearPDE
36	from esys.escript.linearPDEs import LinearPDE
37	# This imports the rectangle domain function from finley.
38	from esys.finley import Rectangle
39	# A useful unit handling package which will make sure all our units
40	# match up in the equations under SI.
41	from esys.escript.unitsSI import *
42	#For interactive use, you can comment out the next two lines
43	import matplotlib
44	matplotlib.use('agg') #It's just here for automated testing
45	import pylab as pl #Plotting package.
46	import numpy as np #Array package.
47	import os #This package is necessary to handle saving our data.
48	from cblib1 import needdirs
49
50	########################################################MPI WORLD CHECK
51	if getMPISizeWorld() > 1:
52	import sys
53	print "This example will not run in an MPI world."
54	sys.exit(0)
55
56	#################################################ESTABLISHING VARIABLES
57	#PDE related
58	mx = 500*m #meters - model length
59	my = 100*m #meters - model width
60	ndx = 500 # mesh steps in x direction
61	ndy = 1 # mesh steps in y direction
62	boundloc = mx/2 # location of boundary between two blocks
63	q=0.*Celsius #our heat source temperature is now zero
64	Tref=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
65	rho = 2750kg/m*3 #kg/m^{3} density of granite
66	cp = 790.J/(kgK) #j/Kg.K thermal capacity
67	rhocp = rhocp #DENSITY SPECIFIC HEAT
68	eta=0. # RADIATION CONDITION
69	kappa=2.2*W/m/K #watts/m.K thermal conductivity
70
71	#Script/Iteration Related
72	t=0. #our start time, usually zero
73	tend=10*yr #the time we want to end the simulation in years
74	outputs = 40 # number of time steps required.
75	h=(tend-t)/outputs #size of time step
76	#user warning statement
77	print "Expected Number of Output Files is: ", outputs
78	print "Step size is: ", h/day, "days"
79	i=0 #loop counter
80	#the folder to put our outputs in, leave blank "" for script path
81	save_path= os.path.join("data","onedheatdiff002")
82	needdirs([save_path])
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	#define our PDE coeffs
94	mypde.setValue(A=kappakronecker(model),D=rhocp/h,d=eta,y=etaTref)
95	#establish location of boundary between two blocks
96	bound = x[0]-boundloc
97	#set initial temperature
98	T= 0TrefwhereNegative(bound)+Tref*wherePositive(bound)
99
100	#convert solution points for plotting
101	plx = x.toListOfTuples()
102	plx = np.array(plx) #convert to tuple to numpy array
103	plx = plx[:,0] #extract x locations
104
105	########################################################START ITERATION
106	while t<=tend:
107	i+=1 #increment the counter
108	t+=h #increment the current time
109	print "time step ",i," at time ",t/day," days."
110	mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
111	T=mypde.getSolution() #find temperature solution
112	#set up for plotting
113	tempT = T.toListOfTuples(scalarastuple=False)
114	pl.figure(1)
115	pl.plot(plx,tempT)
116	pl.axis([0,500,0,2500])
117	pl.title("Temperature accross Interface")
118	if getMPIRankWorld() == 0:
119	pl.savefig(os.path.join(save_path,"intpyplot%03d.png" %i))
120	pl.clf()
121
122	# compile the .png files to create two .avi videos that show T change
123	# with time. This opperation uses linux mencoder. For other operating
124	# systems it may be possible to use your favourite video compiler to
125	# convert image files to videos. To enable this step uncomment the
126	# following lines.
127
128	#os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
129	#w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
130	#onedheatdiff002tempT.avi")