/[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 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 4990 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
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 cblib 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 = 400 # 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/(24.*60*60), "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 mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
110 T=mypde.getSolution() #find temperature solution
111 #set up for plotting
112 tempT = T.toListOfTuples(scalarastuple=False)
113 pl.figure(1)
114 pl.plot(plx,tempT)
115 pl.axis([0,500,0,2500])
116 pl.title("Temperature accross Interface")
117 if getMPIRankWorld() == 0:
118 pl.savefig(os.path.join(save_path,"intpyplot%03d.png" %i))
119 pl.clf()
120
121 # compile the *.png files to create two *.avi videos that show T change
122 # with time. This opperation uses linux mencoder. For other operating
123 # systems it may be possible to use your favourite video compiler to
124 # convert image files to videos. To enable this step uncomment the
125 # following lines.
126
127 #os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
128 #w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
129 #onedheatdiff002tempT.avi")

  ViewVC Help
Powered by ViewVC 1.1.26