/[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 2648 - (show annotations)
Mon Sep 7 00:06:15 2009 UTC (11 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 4461 byte(s)
ecording some changes related to Anthony's cookbook.
These tests still do not run.

test_heatref.py is the current problem.
Current problems include:
undefined symbols
a dependence on X somehow.

I've already addressed the matplotlib X thing somewhere else.
Find and apply.


1
2 ########################################################
3 #
4 # Copyright (c) 2003-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) 2003-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 import pylab as pl #Plotting package.
43 import numpy as np #Array package.
44 import os #This package is necessary to handle saving our data.
45 from cblib import needdirs
46
47
48 #################################################ESTABLISHING VARIABLES
49 #PDE related
50 mx = 500*m #meters - model length
51 my = 100*m #meters - model width
52 ndx = 500 # mesh steps in x direction
53 ndy = 1 # mesh steps in y direction
54 boundloc = mx/2 # location of boundary between two blocks
55 q=0.*Celsius #our heat source temperature is now zero
56 Tref=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
57 rho = 2750*kg/m**3 #kg/m^{3} density of granite
58 cp = 790.*J/(kg*K) #j/Kg.K thermal capacity
59 rhocp = rho*cp #DENSITY * SPECIFIC HEAT
60 eta=0. # RADIATION CONDITION
61 kappa=2.2*W/m/K #watts/m.K thermal conductivity
62
63 #Script/Iteration Related
64 t=0. #our start time, usually zero
65 tend=10*yr #the time we want to end the simulation in years
66 outputs = 400 # number of time steps required.
67 h=(tend-t)/outputs #size of time step
68 #user warning statement
69 print "Expected Number of Output Files is: ", outputs
70 print "Step size is: ", h/(24.*60*60), "days"
71 i=0 #loop counter
72 #the folder to put our outputs in, leave blank "" for script path
73 save_path="data/onedheatdiff002"
74 needdirs([save_path])
75 ########## note this folder path must exist to work ###################
76
77 ################################################ESTABLISHING PARAMETERS
78 #generate domain using rectangle
79 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
80 #extract finite points - the solution points
81 x=model.getX()
82 #create the PDE
83 mypde=LinearPDE(model) #assigns a domain to our PDE
84 mypde.setSymmetryOn() #set the fast solver on for symmetry
85 #define our PDE coeffs
86 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
87 #establish location of boundary between two blocks
88 bound = x[0]-boundloc
89 #set initial temperature
90 T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
91
92 #convert solution points for plotting
93 plx = x.toListOfTuples()
94 plx = np.array(plx) #convert to tuple to numpy array
95 plx = plx[:,0] #extract x locations
96
97 ########################################################START ITERATION
98 while t<=tend:
99 i+=1 #increment the counter
100 t+=h #increment the current time
101 mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
102 T=mypde.getSolution() #find temperature solution
103 #set up for plotting
104 tempT = T.toListOfTuples(scalarastuple=False)
105 pl.figure(1)
106 pl.plot(plx,tempT)
107 pl.axis([0,500,0,2500])
108 pl.title("Temperature accross Interface")
109 pl.savefig(os.path.join(save_path,"intpyplot%03d.png") %i)
110 pl.clf()
111
112 # compile the *.png files to create an *.avi video that shows T change
113 # with time. This opperation uses linux mencoder.
114 os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
115 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
116 onedheatdiff002tempT.avi")

  ViewVC Help
Powered by ViewVC 1.1.26