/[escript]/trunk/doc/cookbook/onedheatdiff002.py
ViewVC logotype

Contents of /trunk/doc/cookbook/onedheatdiff002.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2606 - (show annotations)
Thu Aug 13 04:32:23 2009 UTC (9 years, 10 months ago) by ahallam
File MIME type: text/x-python
File size: 4431 byte(s)
Updates to first chapter in cookbook -> approaching completion.
New figures added for model descriptions.
Scripts updated to use pylab/matplotlib and movie making scripts added for *.avi creation.
Still to finalise in chapter 1
- twodheatdiff model of intrusion
	-contouring

- heat refraction scripts and begin cookbook sections
	- pycad discussion
	- steady state pde discussion
	- quiver plot discussion
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 ############################################################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 import pylab as pl #Plotting package.
43 import numpy as np #Array package.
44 import os #This package is necessary to handle saving our data.
45
46 #################################################ESTABLISHING VARIABLES
47 #PDE related
48 mx = 500*m #meters - model length
49 my = 100*m #meters - model width
50 ndx = 500 # mesh steps in x direction
51 ndy = 1 # mesh steps in y direction
52 boundloc = mx/2 # location of boundary between two blocks
53 q=0.*Celsius #our heat source temperature is now zero
54 Tref=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
55 rho = 2750*kg/m**3 #kg/m^{3} density of granite
56 cp = 790.*J/(kg*K) #j/Kg.K thermal capacity
57 rhocp = rho*cp #DENSITY * SPECIFIC HEAT
58 eta=0. # RADIATION CONDITION
59 kappa=2.2*W/m/K #watts/m.K thermal conductivity
60
61 #Script/Iteration Related
62 t=0. #our start time, usually zero
63 tday=10*365. #the time we want to end the simulation in days
64 tend=tday*24*60*60
65 outputs = 400 # number of time steps required.
66 h=(tend-t)/outputs #size of time step
67 #user warning statement
68 print "Expected Number of Output Files is: ", outputs
69 print "Step size is: ", h/(24.*60*60), "days"
70 i=0 #loop counter
71 #the folder to put our outputs in, leave blank "" for script path
72 save_path="data/onedheatdiff002"
73 ########## note this folder path must exist to work ###################
74
75 ################################################ESTABLISHING PARAMETERS
76 #generate domain using rectangle
77 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
78 #extract finite points - the solution points
79 x=model.getX()
80 #create the PDE
81 mypde=LinearPDE(model) #assigns a domain to our PDE
82 mypde.setSymmetryOn() #set the fast solver on for symmetry
83 #define our PDE coeffs
84 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
85 #establish location of boundary between two blocks
86 bound = x[0]-boundloc
87 #set initial temperature
88 T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
89
90 #convert solution points for plotting
91 plx = x.toListOfTuples()
92 plx = np.array(plx) #convert to tuple to numpy array
93 plx = plx[:,0] #extract x locations
94
95 ########################################################START ITERATION
96 while t<=tend:
97 i+=1 #increment the counter
98 t+=h #increment the current time
99 mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
100 T=mypde.getSolution() #find temperature solution
101 #set up for plotting
102 tempT = T.toListOfTuples(scalarastuple=False)
103 pl.figure(1)
104 pl.plot(plx,tempT)
105 pl.axis([0,500,0,2500])
106 pl.title("Temperature accross Interface")
107 pl.savefig(os.path.join(save_path,"intpyplot%03d.png") %i)
108 pl.clf()
109
110 # compile the *.png files to create an *.avi video that shows T change
111 # with time. This opperation uses linux mencoder.
112 os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
113 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
114 onedheatdiff002tempT.avi")

  ViewVC Help
Powered by ViewVC 1.1.26