/[escript]/trunk/doc/cookbook/onedheatdiff002.py
ViewVC logotype

Diff of /trunk/doc/cookbook/onedheatdiff002.py

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

revision 2605 by caltinay, Thu Jul 16 06:49:19 2009 UTC revision 2606 by ahallam, Thu Aug 13 04:32:23 2009 UTC
# Line 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    
26  # To solve the problem it is necessary to import the modules we require.  ############################################################FILE HEADER
27  from esys.escript import * # This imports everything from the escript library  # onedheatdiff002.py
28  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  # Model temperature diffusion between two granite blocks. This is a one
29  from esys.finley import Rectangle # This imports the rectangle domain function from finley  # 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.  import os #This package is necessary to handle saving our data.
45    
46    #################################################ESTABLISHING VARIABLES
   
 ##ESTABLISHING VARIABLES  
47  #PDE related  #PDE related
48  mx = 500 # model lenght  mx = 500*m #meters - model length
49  my = 100 # model width  my = 100*m #meters - model width
50  ndx = 500 # steps in x direction  ndx = 500 # mesh steps in x direction
51  ndy = 1 # steps in y direction  ndy = 1 # mesh steps in y direction
52    boundloc = mx/2 # location of boundary between two blocks
53  q=0 #our heat source temperature is now zero  q=0.*Celsius #our heat source temperature is now zero
54  Tref=2273 # Kelvin #the starting temperature of our intrusion  Tref=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
55  rho = 2750 #kg/m^{3} density  rho = 2750*kg/m**3 #kg/m^{3} density of granite
56  cp = 790 #j/kg specific heat  cp = 790.*J/(kg*K) #j/Kg.K thermal capacity
57  rhocp = rho*cp  #DENSITY * SPECIFIC HEAT  rhocp = rho*cp  #DENSITY * SPECIFIC HEAT
58  eta=0.  # RADIATION CONDITION  eta=0.  # RADIATION CONDITION
59  kappa=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY  kappa=2.2*W/m/K #watts/m.K thermal conductivity
60    
61  #Script/Iteration Related  #Script/Iteration Related
62  t=0. #our start time, usually zero  t=0. #our start time, usually zero
63  tday=10*365. #the time we want to end the simulation in days  tday=10*365. #the time we want to end the simulation in days
64  tend=tday*24*60*60  tend=tday*24*60*60
65  outputs = 400 # number of time steps required.  outputs = 400 # 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 Output Files is: ", outputs  print "Expected Number of Output Files is: ", outputs
69  print "Step size is: ", h/(24.*60*60), "days"  print "Step size is: ", h/(24.*60*60), "days"
   
   
70  i=0 #loop counter  i=0 #loop counter
71  save_path = "data/onedheatdiff_var001" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work  #the folder to put our outputs in, leave blank "" for script path
72    save_path="data/onedheatdiff002"
73    ########## note this folder path must exist to work ###################
74    
75  #... generate domain ...  ################################################ESTABLISHING PARAMETERS
76    #generate domain using rectangle
77  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
78  # extract finite points  #extract finite points - the solution points
79  x=model.getX()  x=model.getX()
80  #... open PDE ...  #create the PDE
81  mypde=LinearPDE(model)  mypde=LinearPDE(model) #assigns a domain to our PDE
82  mypde.setSymmetryOn()  mypde.setSymmetryOn() #set the fast solver on for symmetry
83    #define our PDE coeffs
84  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
85    #establish location of boundary between two blocks
86  # ... set initial temperature ....  bound = x[0]-boundloc
87  bound = x[0]-mx/(ndx/250.)  #set initial temperature
88  T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)  T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
89    
90  saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)  #convert solution points for plotting
91    plx = x.toListOfTuples()
92    plx = np.array(plx) #convert to tuple to numpy array
93    plx = plx[:,0] #extract x locations
94    
95  # ... start iteration:  ########################################################START ITERATION
96  while t<=tend:  while t<=tend:
97        i+=1      i+=1 #increment the counter
98        t+=h      t+=h #increment the current time
99        mypde.setValue(Y=rhocp/h*T)      mypde.setValue(Y=rhocp/h*T) #reset variable PDE coefficients
100        T=mypde.getSolution()      T=mypde.getSolution() #find temperature solution
101        saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)      #set up for plotting
102              tempT = T.toListOfTuples(scalarastuple=False)
103  # iteration var 2      pl.figure(1)
104  #Tl = 0      pl.plot(plx,tempT)
105  #Tr = Tref      pl.axis([0,500,0,2500])
106  #while Tl < Tr*0.8:      pl.title("Temperature accross Interface")
107        #mypde.setValue(Y=rhocp/h*T)      pl.savefig(os.path.join(save_path,"intpyplot%03d.png") %i)
108        #T=mypde.getSolution()      pl.clf()
109        #i+=1      
110        #x=rod.getX()  # compile the *.png files to create an *.avi video that shows T change
111        #Tl= x[0]  # with time. This opperation uses linux mencoder.
112        #Tr= x[  os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
113    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
114  #print "Finish temp balance in:", i*h/(24.*60*60), " days"  onedheatdiff002tempT.avi")

Legend:
Removed from v.2605  
changed lines
  Added in v.2606

  ViewVC Help
Powered by ViewVC 1.1.26