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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (10 years, 2 months ago) by jfenwick
Original Path: trunk/doc/examples/cookbook/onedheatdiff001.py
File MIME type: text/x-python
File size: 5570 byte(s)
Renamed the main cookbook tex file to match our convention.
Replaced doc/cookbook/figures/heatrefraction002contqu.pdf with
a version which is actually pdf. However it needs to be regnerated since
it it sideways.

The examples have had their copyright notices fixed (dates were too early).
sb2.py has been removed since it uses pyvisi.

scons will now build the cookbook as parts of a docs build.
Also in reposnse to :
scons cookbook_pdf


1
2 ########################################################
3 #
4 # Copyright (c) 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) 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 ############################################################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
30 #######################################################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 #For interactive use, you can comment out the next two lines
42 import matplotlib
43 matplotlib.use('agg') #It's just here for automated testing
44 import pylab as pl #Plotting package.
45 import numpy as np #Array package.
46 import os #This package is necessary to handle saving our data.
47 import cblib
48
49 ########################################################MPI WORLD CHECK
50 if getMPISizeWorld() > 1:
51 import sys
52 print "This example will not run in an MPI world."
53 sys.exit(0)
54
55 #################################################ESTABLISHING VARIABLES
56 #Domain related.
57 mx = 1*m #meters - model length
58 my = .1*m #meters - model width
59 ndx = 100 # mesh steps in x direction
60 ndy = 1 # mesh steps in y direction - one dimension means one element
61
62 #PDE related
63 q=200. * Celsius #Kelvin - our heat source temperature
64 Tref = 0. * Celsius #Kelvin - starting temp of iron bar
65 rho = 7874. *kg/m**3 #kg/m^{3} density of iron
66 cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
67 rhocp = rho*cp
68 kappa = 80.*W/m/K #watts/m.Kthermal conductivity
69
70 #Script/Iteration Related
71 t=0 #our start time, usually zero
72 tend=5.*minute #seconds - time to end simulation
73 outputs = 200 # number of time steps required.
74 h=(tend-t)/outputs #size of time step
75 #user warning statement
76 print "Expected Number of time outputs is: ", (tend-t)/h
77 i=0 #loop counter
78 #the folder to put our outputs in, leave blank "" for script path
79 save_path= os.path.join("data","onedheatdiff001")
80
81 #ensure the dir exists
82 cblib.needdirs([save_path, os.path.join(save_path,"tempT"),\
83 os.path.join(save_path, "totT")])
84
85 ################################################ESTABLISHING PARAMETERS
86 #generate domain using rectangle
87 rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
88 #extract finite points - the solution points
89 x=rod.getX()
90 #create the PDE
91 mypde=LinearPDE(rod) #assigns a domain to our PDE
92 mypde.setSymmetryOn() #set the fast solver on for symmetry
93 mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) #define our PDE coeffs
94 qH=q*whereZero(x[0]) #set heat source
95 T=Tref # set initial temperature
96
97 #convert solution points for plotting
98 plx = x.toListOfTuples()
99 plx = np.array(plx) #convert to tuple to numpy array
100 plx = plx[:,0] #extract x locations
101
102 ########################################################START ITERATION
103 while t<=tend:
104 i+=1 #increment the counter
105 t+=h #increment the current time
106 mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients
107 T=mypde.getSolution() #get the PDE solution
108 totT = rhocp*T #get the total heat solution in the system
109
110 #establish figure 1 for temperature vs x plots
111 tempT = T.toListOfTuples(scalarastuple=False)
112 pl.figure(1) #current figure
113 pl.plot(plx,tempT) #plot solution
114 #define axis extents and title
115 pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
116 pl.title("Temperature accross Rod")
117 #save figure to file
118 if getMPIRankWorld() == 0:
119 pl.savefig(os.path.join(save_path,"tempT",\
120 "rodpyplot%03d.png"%i))
121 pl.clf() #clear figure
122
123 #establish figure 2 for total temperature vs x plots and repeat
124 tottempT = totT.toListOfTuples(scalarastuple=False)
125 pl.figure(2)
126 pl.plot(plx,tottempT)
127 pl.axis([0,1.0,9.657E08,12000+9.657E08])
128 pl.title("Total temperature accross Rod")
129 if getMPIRankWorld() == 0:
130 pl.savefig(os.path.join(save_path,"totT",\
131 "ttrodpyplot%03d.png"%i))
132 pl.clf()
133
134 # compile the *.png files to create two *.avi videos that show T change
135 # with time. This opperation uses linux mencoder. For other operating
136 # systems it may be possible to use your favourite video compiler to
137 # convert image files to videos. To enable this step uncomment the
138 # following lines.
139
140 #os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
141 #w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
142 #onedheatdiff001tempT.avi")
143
144 #os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
145 #w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
146 #onedheatdiff001totT.avi")

  ViewVC Help
Powered by ViewVC 1.1.26