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
"""

# 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 # This defines LinearPDE as LinearPDE
from esys.finley import Rectangle # This imports the rectangle domain function from finley
from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
import os #This package is necessary to handle saving our data.
#plotting tools
#import matplotlib as ptool


##ESTABLISHING VARIABLES
#Domain related.
mx = 1*m #meters - model lenght
my = .1*m #meters - model width
ndx = 100 # steps in x direction 
ndy = 1 # steps in y direction

#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
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 


#... generate domain ...
rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
# extract finite points
x=rod.getX()
#... open PDE ...
mypde=LinearPDE(rod)
mypde.setSymmetryOn()
mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)

# ... set heat source: ....
qH=q*whereZero(x[0])
# ... set initial temperature ....
T=Tref

# ... start iteration:
while t<=tend:
      i+=1
      t+=h
      mypde.setValue(Y=qH+rhocp/h*T)
      T=mypde.getSolution()
      #ptool.plot(T)
      saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
      

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
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	ahallam	2494	from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
31	ahallam	2401	import os #This package is necessary to handle saving our data.
32	ahallam	2494	#plotting tools
33			#import matplotlib as ptool
34	ahallam	2401
35
36			##ESTABLISHING VARIABLES
37			#Domain related.
38	ahallam	2494	mx = 1*m #meters - model lenght
39			my = .1*m #meters - model width
40	ahallam	2401	ndx = 100 # steps in x direction
41			ndy = 1 # steps in y direction
42
43			#PDE related
44	ahallam	2494	q=200. * Celsius #Kelvin - our heat source temperature
45			Tref = 0. * Celsius # Kelvin - starting temp of iron bar
46			rho = 7874. kg/m*3 #kg/m^{3} density of iron
47			cp = 449.J/(kgK) #j/Kg.K thermal capacity
48	ahallam	2401	rhocp = rho*cp
49	ahallam	2494	kappa = 80.*W/m/K #watts/m.Kthermal conductivity
50	ahallam	2401	#Script/Iteration Related
51			t=0 #our start time, usually zero
52	ahallam	2494	tend=5.*minute #seconds - time to end simulation
53	ahallam	2401	outputs = 200 # number of time steps required.
54			h=(tend-t)/outputs #size of time step
55	ahallam	2494	print "Expected Number of time outputs is: ", (tend-t)/h
56	ahallam	2401	i=0 #loop counter
57			#the folder to put our outputs in, leave blank "" for script path
58	ahallam	2494	save_path="data/onedheatdiff001"
59	ahallam	2401	#note this folder path must exist to work
60
61	ahallam	2494
62	ahallam	2401	#... generate domain ...
63			rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
64			# extract finite points
65			x=rod.getX()
66			#... open PDE ...
67			mypde=LinearPDE(rod)
68			mypde.setSymmetryOn()
69	ahallam	2494	mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)
70	ahallam	2401
71			# ... set heat source: ....
72			qH=q*whereZero(x[0])
73			# ... set initial temperature ....
74			T=Tref
75
76			# ... start iteration:
77			while t<=tend:
78			i+=1
79			t+=h
80			mypde.setValue(Y=qH+rhocp/h*T)
81			T=mypde.getSolution()
82	ahallam	2494	#ptool.plot(T)
83	caltinay	2534	saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
84	ahallam	2401
85