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

revision 2589 by ahallam, Wed Aug 5 23:06:35 2009 UTC revision 2606 by ahallam, Thu Aug 13 04:32:23 2009 UTC
22  """  """
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
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  import pylab as pl  #To solve the problem it is necessary to import the modules we require.
32  import numpy as np  #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.
#plotting tools
#import matplotlib as ptool

44
45  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
46  #Domain related.  #Domain related.
47  mx = 1*m #meters - model lenght  mx = 1*m #meters - model length
48  my = .1*m #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=200. * Celsius #Kelvin - our heat source temperature  q=200. * Celsius #Kelvin - our heat source temperature
54  Tref = 0. * Celsius # Kelvin - starting temp of iron bar  Tref = 0. * Celsius #Kelvin - starting temp of iron bar
55  rho = 7874. *kg/m**3 #kg/m^{3} density of iron  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
56  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
57  rhocp = rho*cp  rhocp = rho*cp
58  kappa = 80.*W/m/K #watts/m.Kthermal conductivity  kappa = 80.*W/m/K #watts/m.Kthermal conductivity
59
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.*minute #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    #user warning statement
66  print "Expected Number of time outputs is: ", (tend-t)/h  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  save_path="data/onedheatdiff001"  save_path="data/onedheatdiff001"
70  #note this folder path must exist to work  ########## note this folder path must exist to work ###################
71
72    ################################################ESTABLISHING PARAMETERS
73  #... generate domain ...  #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)  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    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      #ptool.plot(T)      totT = rhocp*T #get the total heat solution in the system
96      saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
97      plx = x.toListOfTuples()      #establish figure 1 for temperature vs x plots
plx = np.array(plx)
plx = plx[:,0]
98      tempT = T.toListOfTuples(scalarastuple=False)      tempT = T.toListOfTuples(scalarastuple=False)
99      pl.plot(plx,tempT)      pl.figure(1) #current figure
100      pl.title("Temperatuer accross Rod")      pl.plot(plx,tempT) #plot solution
101      #~ pl.axis([0.0,1.0,0.0,0.004*2.73146e2])          #define axis extents and title
102      pl.savefig(os.path.join(save_path,"rodpyplot%03d.png") %i)      pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499])
103      pl.clf()            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.2589 changed lines Added in v.2606