examples/cookbook/onedheatdiff001.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
# onedheatdiff001.py
# Model temperature diffusion in an Iron bar. This is a one dimensional
# problem with a single heat source at the LHS

#######################################################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
#Domain related.
mx = 1*m #meters - model length
my = .1*m #meters - model width
ndx = 100 # mesh steps in x direction 
ndy = 1 # mesh steps in y direction - one dimension means one element

#PDE related
q=200. * Celsius #Kelvin - our heat source temperature
Tref = 0. * Celsius #Kelvin - starting temp of iron bar
rho = 7874. *kg/m**3 #kg/m^{3} density of iron
cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
rhocp = rho*cp
kappa = 80.*W/m/K #watts/m.Kthermal conductivity

#Script/Iteration Related
t=0 #our start time, usually zero
tend=5.*minute #seconds - time to end simulation
outputs = 200 # number of time steps required.
h=(tend-t)/outputs #size of time step
#user warning statement
print "Expected Number of time outputs is: ", (tend-t)/h
i=0 #loop counter
#the folder to put our outputs in, leave blank "" for script path 
save_path="data/onedheatdiff001"
########## note this folder path must exist to work ###################

################################################ESTABLISHING PARAMETERS
#generate domain using rectangle
rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
#extract finite points - the solution points
x=rod.getX()
#create the PDE
mypde=LinearPDE(rod) #assigns a domain to our PDE
mypde.setSymmetryOn() #set the fast solver on for symmetry
mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) #define our PDE coeffs
qH=q*whereZero(x[0]) #set heat source
T=Tref # set initial temperature

#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
    mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients
    T=mypde.getSolution() #get the PDE solution
    totT = rhocp*T #get the total heat solution in the system   
    
    #establish figure 1 for temperature vs x plots
    tempT = T.toListOfTuples(scalarastuple=False)
    pl.figure(1) #current figure
    pl.plot(plx,tempT) #plot solution
        #define axis extents and title
    pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
    pl.title("Temperature accross Rod")
        #save figure to file
    pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i)
    pl.clf() #clear figure
    
    #establish figure 2 for total temperature vs x plots and repeat
    tottempT = totT.toListOfTuples(scalarastuple=False)
    pl.figure(2)
    pl.plot(plx,tottempT)
    pl.axis([0,1.0,9.657E08,12000+9.657E08])
    pl.title("Total temperature accross Rod")
    pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%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 is possible to use your favourite video compiler to
# convert image files to videos.

os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
onedheatdiff001tempT.avi")

os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
onedheatdiff001totT.avi")
1	ahallam	2401
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	ahallam	2606	############################################################FILE HEADER
26			# onedheatdiff001.py
27			# Model temperature diffusion in an Iron bar. This is a one dimensional
28			# problem with a single heat source at the LHS
29	ahallam	2401
30	ahallam	2606	#######################################################EXTERNAL MODULES
31			#To solve the problem it is necessary to import the modules we require.
32			#This imports everything from the escript library
33			from esys.escript import *
34			# This defines the LinearPDE module as LinearPDE
35			from esys.escript.linearPDEs import LinearPDE
36			# This imports the rectangle domain function from finley.
37			from esys.finley import Rectangle
38			# A useful unit handling package which will make sure all our units
39			# match up in the equations under SI.
40			from esys.escript.unitsSI import *
41			import pylab as pl #Plotting package.
42			import numpy as np #Array package.
43	ahallam	2401	import os #This package is necessary to handle saving our data.
44
45	ahallam	2606	#################################################ESTABLISHING VARIABLES
46	ahallam	2401	#Domain related.
47	ahallam	2606	mx = 1*m #meters - model length
48	ahallam	2494	my = .1*m #meters - model width
49	ahallam	2606	ndx = 100 # mesh steps in x direction
50			ndy = 1 # mesh steps in y direction - one dimension means one element
51	ahallam	2401
52			#PDE related
53	ahallam	2494	q=200. * Celsius #Kelvin - our heat source temperature
54	ahallam	2606	Tref = 0. * Celsius #Kelvin - starting temp of iron bar
55	ahallam	2494	rho = 7874. kg/m*3 #kg/m^{3} density of iron
56			cp = 449.J/(kgK) #j/Kg.K thermal capacity
57	ahallam	2401	rhocp = rho*cp
58	ahallam	2494	kappa = 80.*W/m/K #watts/m.Kthermal conductivity
59	ahallam	2606
60	ahallam	2401	#Script/Iteration Related
61			t=0 #our start time, usually zero
62	ahallam	2494	tend=5.*minute #seconds - time to end simulation
63	ahallam	2401	outputs = 200 # number of time steps required.
64			h=(tend-t)/outputs #size of time step
65	ahallam	2606	#user warning statement
66	ahallam	2494	print "Expected Number of time outputs is: ", (tend-t)/h
67	ahallam	2401	i=0 #loop counter
68			#the folder to put our outputs in, leave blank "" for script path
69	ahallam	2494	save_path="data/onedheatdiff001"
70	ahallam	2606	########## note this folder path must exist to work ###################
71	ahallam	2401
72	ahallam	2606	################################################ESTABLISHING PARAMETERS
73			#generate domain using rectangle
74	ahallam	2401	rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
75	ahallam	2606	#extract finite points - the solution points
76	ahallam	2401	x=rod.getX()
77	ahallam	2606	#create the PDE
78			mypde=LinearPDE(rod) #assigns a domain to our PDE
79			mypde.setSymmetryOn() #set the fast solver on for symmetry
80			mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) #define our PDE coeffs
81			qH=q*whereZero(x[0]) #set heat source
82			T=Tref # set initial temperature
83	ahallam	2401
84	ahallam	2606	#convert solution points for plotting
85			plx = x.toListOfTuples()
86			plx = np.array(plx) #convert to tuple to numpy array
87			plx = plx[:,0] #extract x locations
88	ahallam	2401
89	ahallam	2606	########################################################START ITERATION
90	ahallam	2401	while t<=tend:
91	ahallam	2606	i+=1 #increment the counter
92			t+=h #increment the current time
93			mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients
94			T=mypde.getSolution() #get the PDE solution
95			totT = rhocp*T #get the total heat solution in the system
96
97			#establish figure 1 for temperature vs x plots
98	ahallam	2589	tempT = T.toListOfTuples(scalarastuple=False)
99	ahallam	2606	pl.figure(1) #current figure
100			pl.plot(plx,tempT) #plot solution
101			#define axis extents and title
102			pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
103			pl.title("Temperature accross Rod")
104			#save figure to file
105			pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i)
106			pl.clf() #clear figure
107
108			#establish figure 2 for total temperature vs x plots and repeat
109			tottempT = totT.toListOfTuples(scalarastuple=False)
110			pl.figure(2)
111			pl.plot(plx,tottempT)
112			pl.axis([0,1.0,9.657E08,12000+9.657E08])
113			pl.title("Total temperature accross Rod")
114			pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i)
115			pl.clf()
116	ahallam	2401
117	ahallam	2606	# compile the .png files to create two .avi videos that show T change
118			# with time. This opperation uses linux mencoder. For other operating
119			# systems it is possible to use your favourite video compiler to
120			# convert image files to videos.
121
122			os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
123			w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
124			onedheatdiff001tempT.avi")
125
126			os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
127			w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
128			onedheatdiff001totT.avi")