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

Contents of /trunk/doc/examples/cookbook/twodheatdiff001.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: 5896 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 # twodheatdiff002.py
28 # Model temperature diffusion between a granite intrusion and sandstone
29 # country rock. This is a two dimensional problem with the granite as a
30 # heat source.
31
32 #######################################################EXTERNAL MODULES
33 #To solve the problem it is necessary to import the modules we require.
34 #This imports everything from the escript library
35 from esys.escript import *
36 # This defines the LinearPDE module as LinearPDE
37 from esys.escript.linearPDEs import LinearPDE
38 # This imports the rectangle domain function from finley.
39 from esys.finley import Rectangle
40 # A useful unit handling package which will make sure all our units
41 # match up in the equations under SI.
42 from esys.escript.unitsSI import *
43 #For interactive use, you can comment out the next two lines
44 import matplotlib
45 matplotlib.use('agg') #It's just here for automated testing
46 import pylab as pl #Plotting package.
47 import numpy as np #Array package.
48 import os #This package is necessary to handle saving our data.
49 from cblib import toXYTuple, needdirs
50
51 ########################################################MPI WORLD CHECK
52 if getMPISizeWorld() > 1:
53 import sys
54 print "This example will not run in an MPI world."
55 sys.exit(0)
56
57 #################################################ESTABLISHING VARIABLES
58 #PDE related
59 mx = 600*m #meters - model length
60 my = 600*m #meters - model width
61 ndx = 100 #mesh steps in x direction
62 ndy = 100 #mesh steps in y direction
63 r = 200*m #meters - radius of intrusion
64 ic = [300, 0] #centre of intrusion (meters)
65 q=0.*Celsius #our heat source temperature is now zero
66
67 ## Intrusion Variables - Granite
68 Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
69 rhoi = 2750*kg/m**3 #kg/m^{3} density of granite
70 cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity
71 rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT
72 eta=0. # RADIATION CONDITION
73 kappai=2.2*W/m/K #watts/m.K thermal conductivity
74 ## Country Rock Variables - Sandstone
75 Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
76 rhoc = 2000*kg/m**3 #kg/m^{3} density
77 cpc = 920.*J/(kg*K) #j/kg.k specific heat
78 rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT
79 kappac = 1.9*W/m/K #watts/m.K thermal conductivity
80
81 #Script/Iteration Related
82 t=0. #our start time, usually zero
83 tday=100*365. #the time we want to end the simulation in days
84 tend=tday*24*60*60
85 outputs = 200 # number of time steps required.
86 h=(tend-t)/outputs #size of time step
87 #user warning
88 print "Expected Number of Output Files is: ", outputs
89 print "Step size is: ", h/(24.*60*60), "days"
90 i=0 #loop counter
91 #the folder to put our outputs in, leave blank "" for script path
92 save_path= os.path.join("data","twodheatdiff")
93 needdirs([save_path])
94 ########## note this folder path must exist to work ###################
95
96 ################################################ESTABLISHING PARAMETERS
97 #generate domain using rectangle
98 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
99 #extract finite points - the solution points
100 x=model.getX()
101 #create the PDE
102 mypde=LinearPDE(model) #assigns a domain to our PDE
103 mypde.setSymmetryOn() #set the fast solver on for symmetry
104 #establish location of boundary between two materials
105 bound = length(x-ic)-r #where the boundary will be located
106 A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)
107 D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)
108 #define our PDE coeffs
109 mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)
110 #set initial temperature
111 T= Ti*whereNegative(bound)+Tc*wherePositive(bound)
112
113 # rearrage mymesh to suit solution function space for contouring
114 oldspacecoords=model.getX()
115 coords=Data(oldspacecoords, T.getFunctionSpace())
116 #coords = np.array(coords.toListOfTuples())
117 coordX, coordY = toXYTuple(coords)
118 # create regular grid
119 xi = np.linspace(0.0,mx,100)
120 yi = np.linspace(0.0,my,100)
121
122 ########################################################START ITERATION
123 while t<=tend:
124 i+=1 #counter
125 t+=h #curretn time
126 Y = T*D #
127 mypde.setValue(Y=Y)
128 T=mypde.getSolution()
129 tempT = T.toListOfTuples(scalarastuple=False)
130 # grid the data.
131 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
132 # contour the gridded data, plotting dots at the
133 # randomly spaced data points.
134 pl.matplotlib.pyplot.autumn()
135 pl.contourf(xi,yi,zi,10)
136 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
137 pl.clabel(CS, inline=1, fontsize=8)
138 pl.axis([0,600,0,600])
139 pl.title("Heat diffusion from an intrusion.")
140 pl.xlabel("Horizontal Displacement (m)")
141 pl.ylabel("Depth (m)")
142 if getMPIRankWorld() == 0:
143 pl.savefig(os.path.join(save_path,\
144 "heatrefraction%03d.png"%i))
145 pl.clf()
146
147 # compile the *.png files to create an *.avi video that shows T change
148 # with time. This opperation uses linux mencoder.
149 os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
150 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
151 twodheatdiff001tempT.avi")

  ViewVC Help
Powered by ViewVC 1.1.26