/[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

trunk/doc/cookbook/onedheatdiff001.py revision 2534 by caltinay, Thu Jul 16 06:49:19 2009 UTC trunk/doc/examples/cookbook/onedheatdiff001.py revision 2648 by jfenwick, Mon Sep 7 00:06:15 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  from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.  #######################################################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  #plotting tools  import cblib
 #import matplotlib as ptool  
45    
46    
47  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
48  #Domain related.  #Domain related.
49  mx = 1*m #meters - model lenght  mx = 1*m #meters - model length
50  my = .1*m #meters - model width  my = .1*m #meters - model width
51  ndx = 100 # steps in x direction  ndx = 100 # mesh steps in x direction
52  ndy = 1 # steps in y direction  ndy = 1 # mesh steps in y direction - one dimension means one element
53    
54  #PDE related  #PDE related
55  q=200. * Celsius #Kelvin - our heat source temperature  q=200. * Celsius #Kelvin - our heat source temperature
56  Tref = 0. * Celsius # Kelvin - starting temp of iron bar  Tref = 0. * Celsius #Kelvin - starting temp of iron bar
57  rho = 7874. *kg/m**3 #kg/m^{3} density of iron  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
58  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
59  rhocp = rho*cp  rhocp = rho*cp
60  kappa = 80.*W/m/K #watts/m.Kthermal conductivity  kappa = 80.*W/m/K #watts/m.Kthermal conductivity
61    
62  #Script/Iteration Related  #Script/Iteration Related
63  t=0 #our start time, usually zero  t=0 #our start time, usually zero
64  tend=5.*minute #seconds - time to end simulation  tend=5.*minute #seconds - time to end simulation
65  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
66  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
67    #user warning statement
68  print "Expected Number of time outputs is: ", (tend-t)/h  print "Expected Number of time outputs is: ", (tend-t)/h
69  i=0 #loop counter  i=0 #loop counter
70  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
71  save_path="data/onedheatdiff001"  save_path="data/onedheatdiff001"
 #note this folder path must exist to work  
72    
73    #ensure the dir exists
74    cblib.needdirs([save_path, os.path.join(save_path,"tempT"), os.path.join(save_path, "totT")])
75    
76  #... generate domain ...  ################################################ESTABLISHING PARAMETERS
77    #generate domain using rectangle
78  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
79  # extract finite points  #extract finite points - the solution points
80  x=rod.getX()  x=rod.getX()
81  #... open PDE ...  #create the PDE
82  mypde=LinearPDE(rod)  mypde=LinearPDE(rod) #assigns a domain to our PDE
83  mypde.setSymmetryOn()  mypde.setSymmetryOn() #set the fast solver on for symmetry
84  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) #define our PDE coeffs
85    qH=q*whereZero(x[0]) #set heat source
86  # ... set heat source: ....  T=Tref # set initial temperature
87  qH=q*whereZero(x[0])  
88  # ... set initial temperature ....  #convert solution points for plotting
89  T=Tref  plx = x.toListOfTuples()
90    plx = np.array(plx) #convert to tuple to numpy array
91    plx = plx[:,0] #extract x locations
92    
93  # ... start iteration:  ########################################################START ITERATION
94  while t<=tend:  while t<=tend:
95        i+=1      i+=1 #increment the counter
96        t+=h      t+=h #increment the current time
97        mypde.setValue(Y=qH+rhocp/h*T)      mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients
98        T=mypde.getSolution()      T=mypde.getSolution() #get the PDE solution
99        #ptool.plot(T)      totT = rhocp*T #get the total heat solution in the system  
100        saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)      
101              #establish figure 1 for temperature vs x plots
102        tempT = T.toListOfTuples(scalarastuple=False)
103        pl.figure(1) #current figure
104        pl.plot(plx,tempT) #plot solution
105            #define axis extents and title
106        pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
107        pl.title("Temperature accross Rod")
108            #save figure to file
109        pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i)
110        pl.clf() #clear figure
111        
112        #establish figure 2 for total temperature vs x plots and repeat
113        tottempT = totT.toListOfTuples(scalarastuple=False)
114        pl.figure(2)
115        pl.plot(plx,tottempT)
116        pl.axis([0,1.0,9.657E08,12000+9.657E08])
117        pl.title("Total temperature accross Rod")
118        pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%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 is possible to use your favourite video compiler to
124    # convert image files to videos.
125    
126    os.system("mencoder mf://"+save_path+"/tempT"+"/*.png -mf type=png:\
127    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
128    onedheatdiff001tempT.avi")
129    
130    os.system("mencoder mf://"+save_path+"/totT"+"/*.png -mf type=png:\
131    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
132    onedheatdiff001totT.avi")

Legend:
Removed from v.2534  
changed lines
  Added in v.2648

  ViewVC Help
Powered by ViewVC 1.1.26