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

Diff of /trunk/doc/examples/cookbook/onedheatdiff001.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2401 by ahallam, Wed Apr 29 04:23:07 2009 UTC revision 2606 by ahallam, Thu Aug 13 04:32:23 2009 UTC
# Line 22  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    ############################################################FILE HEADER
26  # To solve the problem it is necessary to import the modules we require.  # onedheatdiff001.py
27  from esys.escript import * # This imports everything from the escript library  # Model temperature diffusion in an Iron bar. This is a one dimensional
28  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  # problem with a single heat source at the LHS
29  from esys.finley import Rectangle # This imports the rectangle domain function from finley  
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    import pylab as pl #Plotting package.
42    import numpy as np #Array package.
43  import os #This package is necessary to handle saving our data.  import os #This package is necessary to handle saving our data.
44    
45    #################################################ESTABLISHING VARIABLES
   
 ##ESTABLISHING VARIABLES  
46  #Domain related.  #Domain related.
47  mx = 1 #meters - model lenght  mx = 1*m #meters - model length
48  my = .1 #meters - model width  my = .1*m #meters - model width
49  ndx = 100 # steps in x direction  ndx = 100 # mesh steps in x direction
50  ndy = 1 # steps in y direction  ndy = 1 # mesh steps in y direction - one dimension means one element
51    
52  #PDE related  #PDE related
53  q=473. #Kelvin - our heat source temperature  q=200. * Celsius #Kelvin - our heat source temperature
54  Tref = 273. # Kelvin - starting temp of iron bar  Tref = 0. * Celsius #Kelvin - starting temp of iron bar
55  rho = 7874. #kg/m^{3} density of iron  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
56  cp = 449. #j/Kg.K  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
57  rhocp = rho*cp  rhocp = rho*cp
58  eta = 0 #radiation condition  kappa = 80.*W/m/K #watts/m.Kthermal conductivity
59  kappa = 68. #temperature diffusion constant  
60  #Script/Iteration Related  #Script/Iteration Related
61  t=0 #our start time, usually zero  t=0 #our start time, usually zero
62  tend=5.*60. #seconds - time to end simulation  tend=5.*minute #seconds - time to end simulation
63  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
64  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
65  print "Expected Number of Output Files is: ", (tend-t)/h  #user warning statement
66    print "Expected Number of time outputs is: ", (tend-t)/h
67  i=0 #loop counter  i=0 #loop counter
68  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
69  #note this folder path must exist to work  save_path="data/onedheatdiff001"
70  save_path = "data/onedheatdiff001"  ########## note this folder path must exist to work ###################
71    
72  #... generate domain ...  ################################################ESTABLISHING PARAMETERS
73    #generate domain using rectangle
74  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
75  # extract finite points  #extract finite points - the solution points
76  x=rod.getX()  x=rod.getX()
77  #... open PDE ...  #create the PDE
78  mypde=LinearPDE(rod)  mypde=LinearPDE(rod) #assigns a domain to our PDE
79  mypde.setSymmetryOn()  mypde.setSymmetryOn() #set the fast solver on for symmetry
80  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h,d=eta,y=eta*Tref)  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) #define our PDE coeffs
81    qH=q*whereZero(x[0]) #set heat source
82  # ... set heat source: ....  T=Tref # set initial temperature
83  qH=q*whereZero(x[0])  
84  # ... set initial temperature ....  #convert solution points for plotting
85  T=Tref  plx = x.toListOfTuples()
86    plx = np.array(plx) #convert to tuple to numpy array
87  #saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)  plx = plx[:,0] #extract x locations
88    
89  # ... start iteration:  ########################################################START ITERATION
90  while t<=tend:  while t<=tend:
91        i+=1      i+=1 #increment the counter
92        t+=h      t+=h #increment the current time
93        mypde.setValue(Y=qH+rhocp/h*T)      mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients
94        T=mypde.getSolution()      T=mypde.getSolution() #get the PDE solution
95        saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)      totT = rhocp*T #get the total heat solution in the system  
96              
97        #establish figure 1 for temperature vs x plots
98        tempT = T.toListOfTuples(scalarastuple=False)
99        pl.figure(1) #current figure
100        pl.plot(plx,tempT) #plot solution
101            #define axis extents and title
102        pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
103        pl.title("Temperature accross Rod")
104            #save figure to file
105        pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i)
106        pl.clf() #clear figure
107        
108        #establish figure 2 for total temperature vs x plots and repeat
109        tottempT = totT.toListOfTuples(scalarastuple=False)
110        pl.figure(2)
111        pl.plot(plx,tottempT)
112        pl.axis([0,1.0,9.657E08,12000+9.657E08])
113        pl.title("Total temperature accross Rod")
114        pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i)
115        pl.clf()
116    
117    # compile the *.png files to create two *.avi videos that show T change
118    # with time. This opperation uses linux mencoder. For other operating
119    # systems it is possible to use your favourite video compiler to
120    # convert image files to videos.
121    
122    os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
123    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
124    onedheatdiff001tempT.avi")
125    
126    os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
127    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
128    onedheatdiff001totT.avi")
Legend:
Removed from v.2401  
changed lines
  Added in v.2606
  ViewVC Help
Powered by ViewVC 1.1.26