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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 5049 byte(s)
Don't panic.
Updating copyright stamps

1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 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-2010 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 = 2750*kg/m**3 #kg/m^{3} density of granite
66 cp = 790.*J/(kg*K) #j/Kg.K thermal capacity
67 rhocp = rho*cp #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=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
95 #establish location of boundary between two blocks
96 bound = x[0]-boundloc
97 #set initial temperature
98 T= 0*Tref*whereNegative(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")

  ViewVC Help
Powered by ViewVC 1.1.26