examples/cookbook/twodheatdiff.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
import os #This package is necessary to handle saving our data.

from cblib import needdirs


##ESTABLISHING VARIABLES
#PDE related
mx = 600 # model lenght
my = 600 # model width
ndx = 100 # steps in x direction 
ndy = 100 # steps in y direction
r = 200 # radius of intrusion
ic = [300, 0] #centre of intrusion

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

Tc = 200
rhoc = 2200
cpc = 400
rhocpc = rhoc*cpc
kappac = 0.1


#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

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


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

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

#... open PDE ...
mypde=LinearPDE(model)
mypde.setSymmetryOn()

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)

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

# ... set initial temperature ....

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

#... start iteration:
while t<=tend:
      i+=1
      t+=h
      Y = T*D
      mypde.setValue(Y=Y)
      T=mypde.getSolution()
      saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
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	# 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 = rhoicpi #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=tday2460*60
63	outputs = 200 # number of time steps required.
64	h=(tend-t)/outputs #size of time step
65
66	print "Expected Number of Output Files is: ", outputs
67	print "Step size is: ", h/(24.6060), "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(x-ic)-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=Akronecker(model),D=D,d=eta,y=etaTc)
89
90	# ... set initial temperature ....
91
92	T= TiwhereNegative(bound)+TcwherePositive(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)