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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (hide 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 ahallam 2401
2     ########################################################
3     #
4 jfenwick 2667 # Copyright (c) 2009 by University of Queensland
5 ahallam 2401 # 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 jfenwick 2667 __copyright__="""Copyright (c) 2009 by University of Queensland
15 ahallam 2401 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 ahallam 2606 ############################################################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 ahallam 2401
30 ahallam 2606 #######################################################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 ahallam 2658 #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 ahallam 2606 import pylab as pl #Plotting package.
45     import numpy as np #Array package.
46 ahallam 2401 import os #This package is necessary to handle saving our data.
47 jfenwick 2648 import cblib
48 ahallam 2401
49 ahallam 2658 ########################################################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 jfenwick 2648
55 ahallam 2606 #################################################ESTABLISHING VARIABLES
56 ahallam 2401 #Domain related.
57 ahallam 2606 mx = 1*m #meters - model length
58 ahallam 2494 my = .1*m #meters - model width
59 ahallam 2606 ndx = 100 # mesh steps in x direction
60     ndy = 1 # mesh steps in y direction - one dimension means one element
61 ahallam 2401
62     #PDE related
63 ahallam 2494 q=200. * Celsius #Kelvin - our heat source temperature
64 ahallam 2606 Tref = 0. * Celsius #Kelvin - starting temp of iron bar
65 ahallam 2494 rho = 7874. *kg/m**3 #kg/m^{3} density of iron
66     cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
67 ahallam 2401 rhocp = rho*cp
68 ahallam 2494 kappa = 80.*W/m/K #watts/m.Kthermal conductivity
69 ahallam 2606
70 ahallam 2401 #Script/Iteration Related
71     t=0 #our start time, usually zero
72 ahallam 2494 tend=5.*minute #seconds - time to end simulation
73 ahallam 2401 outputs = 200 # number of time steps required.
74     h=(tend-t)/outputs #size of time step
75 ahallam 2606 #user warning statement
76 ahallam 2494 print "Expected Number of time outputs is: ", (tend-t)/h
77 ahallam 2401 i=0 #loop counter
78     #the folder to put our outputs in, leave blank "" for script path
79 ahallam 2658 save_path= os.path.join("data","onedheatdiff001")
80 ahallam 2401
81 jfenwick 2648 #ensure the dir exists
82 ahallam 2658 cblib.needdirs([save_path, os.path.join(save_path,"tempT"),\
83     os.path.join(save_path, "totT")])
84 jfenwick 2648
85 ahallam 2606 ################################################ESTABLISHING PARAMETERS
86     #generate domain using rectangle
87 ahallam 2401 rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
88 ahallam 2606 #extract finite points - the solution points
89 ahallam 2401 x=rod.getX()
90 ahallam 2606 #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 ahallam 2401
97 ahallam 2606 #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 ahallam 2401
102 ahallam 2606 ########################################################START ITERATION
103 ahallam 2401 while t<=tend:
104 ahallam 2606 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 ahallam 2589 tempT = T.toListOfTuples(scalarastuple=False)
112 ahallam 2606 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 ahallam 2658 if getMPIRankWorld() == 0:
119     pl.savefig(os.path.join(save_path,"tempT",\
120     "rodpyplot%03d.png"%i))
121 ahallam 2606 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 ahallam 2658 if getMPIRankWorld() == 0:
130     pl.savefig(os.path.join(save_path,"totT",\
131     "ttrodpyplot%03d.png"%i))
132 ahallam 2606 pl.clf()
133 ahallam 2401
134 ahallam 2606 # 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 jfenwick 2667 # systems it may be possible to use your favourite video compiler to
137 ahallam 2658 # convert image files to videos. To enable this step uncomment the
138     # following lines.
139 ahallam 2606
140 ahallam 2658 #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 ahallam 2606
144 ahallam 2658 #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