/[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 2494 by ahallam, Thu Jun 25 05:00:57 2009 UTC
# Line 27  Author: Antony Hallam antony.hallam@uqco Line 27  Author: Antony Hallam antony.hallam@uqco
27  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
28  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
29  from esys.finley import Rectangle # This imports the rectangle domain function from finley  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.
31  import os #This package is necessary to handle saving our data.  import os #This package is necessary to handle saving our data.
32    #plotting tools
33    #import matplotlib as ptool
34    
35    
36  ##ESTABLISHING VARIABLES  ##ESTABLISHING VARIABLES
37  #Domain related.  #Domain related.
38  mx = 1 #meters - model lenght  mx = 1*m #meters - model lenght
39  my = .1 #meters - model width  my = .1*m #meters - model width
40  ndx = 100 # steps in x direction  ndx = 100 # steps in x direction
41  ndy = 1 # steps in y direction  ndy = 1 # steps in y direction
42    
43  #PDE related  #PDE related
44  q=473. #Kelvin - our heat source temperature  q=200. * Celsius #Kelvin - our heat source temperature
45  Tref = 273. # Kelvin - starting temp of iron bar  Tref = 0. * Celsius # Kelvin - starting temp of iron bar
46  rho = 7874. #kg/m^{3} density of iron  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
47  cp = 449. #j/Kg.K  cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
48  rhocp = rho*cp  rhocp = rho*cp
49  eta = 0 #radiation condition  kappa = 80.*W/m/K #watts/m.Kthermal conductivity
 kappa = 68. #temperature diffusion constant  
50  #Script/Iteration Related  #Script/Iteration Related
51  t=0 #our start time, usually zero  t=0 #our start time, usually zero
52  tend=5.*60. #seconds - time to end simulation  tend=5.*minute #seconds - time to end simulation
53  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
54  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
55  print "Expected Number of Output Files is: ", (tend-t)/h  print "Expected Number of time outputs is: ", (tend-t)/h
56  i=0 #loop counter  i=0 #loop counter
57  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
58    save_path="data/onedheatdiff001"
59  #note this folder path must exist to work  #note this folder path must exist to work
60  save_path = "data/onedheatdiff001"  
61    
62  #... generate domain ...  #... generate domain ...
63  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
# Line 64  x=rod.getX() Line 66  x=rod.getX()
66  #... open PDE ...  #... open PDE ...
67  mypde=LinearPDE(rod)  mypde=LinearPDE(rod)
68  mypde.setSymmetryOn()  mypde.setSymmetryOn()
69  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h,d=eta,y=eta*Tref)  mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)
70    
71  # ... set heat source: ....  # ... set heat source: ....
72  qH=q*whereZero(x[0])  qH=q*whereZero(x[0])
73  # ... set initial temperature ....  # ... set initial temperature ....
74  T=Tref  T=Tref
75    
 #saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)  
   
76  # ... start iteration:  # ... start iteration:
77  while t<=tend:  while t<=tend:
78        i+=1        i+=1
79        t+=h        t+=h
80        mypde.setValue(Y=qH+rhocp/h*T)        mypde.setValue(Y=qH+rhocp/h*T)
81        T=mypde.getSolution()        T=mypde.getSolution()
82          #ptool.plot(T)
83        saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)        saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
84                
85    
Legend:
Removed from v.2401  
changed lines
  Added in v.2494
  ViewVC Help
Powered by ViewVC 1.1.26