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

trunk/doc/examples/cookbook/onedheatdiff.py revision 2697 by jfenwick, Wed Sep 30 01:29:22 2009 UTC trunk/doc/examples/cookbook/onedheatdiffbase.py revision 2904 by gross, Tue Feb 2 04:21:52 2010 UTC
# Line 25  Author: Antony Hallam antony.hallam@uqco Line 25  Author: Antony Hallam antony.hallam@uqco
25
26  # To solve the problem it is necessary to import the modules we require.  # To solve the problem it is necessary to import the modules we require.
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.unitsSI import *
29  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
30  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
import os #This package is necessary to handle saving our data.
#from cblib import needdirs
from cblib1 import needdirs
31
32  ##ESTABLISHING VARIABLES  ##ESTABLISHING VARIABLES
33    #Domain related.
34    mx = 500*m #meters - model length
35    my = 100*m #meters - model width
36    ndx = 100 # mesh steps in x direction
37    ndy = 1 # mesh steps in y direction - one dimension means one element
38    boundloc = mx/2 # location of boundary between the two blocks
39  #PDE related  #PDE related
40  q=50.e6 #our heat source temperature  rho = 7874. *kg/m**3 #kg/m^{3} density of iron
41  Tref=0. #the starting temperature of our iron bar  cp = 449.*J/(kg*K) # J/Kg.K thermal capacity
42  rho=2.6e6  rhocp = rho*cp
43  eta=0#75.  kappa = 80.*W/m/K   # watts/m.Kthermal conductivity
44  kappa=240.  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source
45  #Script/Iteration Related  T1=20 * Celsius # initial temperature at Block 1
46  t=0 #our start time, usually zero  T2=2273. * Celsius # initial temperature at Block 2
47  tend=5.#the time we want to end the simulation
48  h=0.05 #size of time step  t=0 * day  # our start time, usually zero
49    tend=50 * yr # - time to end simulation
50  print "Expected Number of Output Files is: ", (tend-t)/h  outputs = 200 # number of time steps required.
51    h=(tend-t)/outputs #size of time step
52  i=0 #loop counter  #user warning statement
53  save_path = "data/onedheatdiff" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work  print "Expected Number of time outputs is: ", (tend-t)/h
54  needdirs([save_path])  i=0 #loop counter
55    #the folder to put our outputs in, leave blank "" for script path
56    save_path= os.path.join("data","onedheatdiff001")
57    #ensure the dir exists
58    mkDir(save_path, os.path.join(save_path,"tempT"))
59
60  #... generate domain ...  #... generate domain ...
61  rod = Rectangle(l0=0.05,l1=.01,n0=500, n1=1)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
62  # extract finite points  #... open PDE and set coefficients ...
63  x=rod.getX()  mypde=LinearPDE(blocks)
#... open PDE ...
mypde=LinearPDE(rod)
64  mypde.setSymmetryOn()  mypde.setSymmetryOn()
65  mypde.setValue(A=kappa*kronecker(rod),D=rho/h,d=eta,y=eta*Tref)  A=zeros((2,2))
66  # ... set heat source: ....  A[0,0]=kappa
67    mypde.setValue(A=A,D=rhocp/h)
qH=q*whereZero(x[0])
68  # ... set initial temperature ....  # ... set initial temperature ....
69  T=Tref  x=Solution(blocks).getX()
70    T= T1*whereNegative(x[0]-boundloc)+T2*(1-whereNegative(x[0]-boundloc))
71
72  # ... start iteration:  # ... start iteration:
73  while t<=tend:  while t<tend:
74        i+=1        i+=1
75        t+=h        t+=h
76        mypde.setValue(Y=qH+rho/h*T)        mypde.setValue(Y=qH+rhocp/h*T)
77        T=mypde.getSolution()        T=mypde.getSolution()
78        print T        totE=integrate(rhocp*T)
79        saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)        print "time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE)

#~ command = ('mencoder',
#~ 'mf://*.png',
#~ '-mf',
#~ 'type=png:w=800:h=600:fps=25',
#~ '-ovc',
#~ 'lavc',
#~ '-lavcopts',
#~ 'vcodec=mpeg4',
#~ '-oac',
#~ 'copy',
#~ '-o',
#~ 'output.avi')

Legend:
 Removed from v.2697 changed lines Added in v.2904