/[escript]/branches/stage3.1/doc/examples/cookbook/twodheatdiff001.py
ViewVC logotype

Diff of /branches/stage3.1/doc/examples/cookbook/twodheatdiff001.py

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

revision 2944 by jfenwick, Thu Feb 4 01:42:47 2010 UTC revision 2945 by jfenwick, Wed Feb 24 00:17:46 2010 UTC
# Line 46  matplotlib.use('agg') #It's just here fo Line 46  matplotlib.use('agg') #It's just here fo
46  import pylab as pl #Plotting package.  import pylab as pl #Plotting package.
47  import numpy as np #Array package.  import numpy as np #Array package.
48  import os #This package is necessary to handle saving our data.  import os #This package is necessary to handle saving our data.
49  from cblib2 import toXYTuple  from cblib import toXYTuple
50    
51  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
52  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
# Line 58  if getMPISizeWorld() > 1: Line 58  if getMPISizeWorld() > 1:
58  #PDE related  #PDE related
59  mx = 600*m #meters - model length  mx = 600*m #meters - model length
60  my = 600*m #meters - model width  my = 600*m #meters - model width
61  ndx = 100 #mesh steps in x direction  ndx = 150 #mesh steps in x direction
62  ndy = 100 #mesh steps in y direction  ndy = 150 #mesh steps in y direction
63  r = 200*m #meters - radius of intrusion  r = 200*m #meters - radius of intrusion
64  ic = [300, 0] #centre of intrusion (meters)  ic = [300*m, 0] #centre of intrusion (meters)
65  q=0.*Celsius #our heat source temperature is now zero  qH=0.*J/(sec*m**3) #our heat source temperature is now zero
66    
67  ## Intrusion Variables - Granite  ## Intrusion Variables - Granite
68  Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block  Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
69  rhoi = 2750*kg/m**3 #kg/m^{3} density of granite  rhoi = 2750*kg/m**3 #kg/m^{3} density of granite
70  cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity  cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity
71  rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT  rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT
 eta=0.  # RADIATION CONDITION  
72  kappai=2.2*W/m/K #watts/m.K thermal conductivity  kappai=2.2*W/m/K #watts/m.K thermal conductivity
73  ## Country Rock Variables - Sandstone  ## Country Rock Variables - Sandstone
74  Tc = 473*Celsius # Kelvin #the starting temperature of our country rock  Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
# Line 80  kappac = 1.9*W/m/K #watts/m.K thermal co Line 79  kappac = 1.9*W/m/K #watts/m.K thermal co
79    
80  #Script/Iteration Related  #Script/Iteration Related
81  t=0. #our start time, usually zero  t=0. #our start time, usually zero
82  tday=100*365. #the time we want to end the simulation in days  tend=200.* yr #the time we want to end the simulation
 tend=tday*24*60*60  
83  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
84  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
85  #user warning  #user warning
86  print "Expected Number of Output Files is: ", outputs  print "Expected Number of Output Files is: ", outputs
87  print "Step size is: ", h/(24.*60*60), "days"  print "Step size is: ", h/day, "days"
88  i=0 #loop counter  i=0 #loop counter
89  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
90  save_path= os.path.join("data","twodheatdiff")  save_path= os.path.join("data","twodheatdiff")
# Line 97  mkDir(save_path) Line 95  mkDir(save_path)
95  #generate domain using rectangle  #generate domain using rectangle
96  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
97  #extract finite points - the solution points  #extract finite points - the solution points
 x=model.getX()  
98  #create the PDE  #create the PDE
99  mypde=LinearPDE(model) #assigns a domain to our PDE  mypde=LinearPDE(model) #assigns a domain to our PDE
100  mypde.setSymmetryOn() #set the fast solver on for symmetry  mypde.setSymmetryOn() #set the fast solver on for symmetry
101  #establish location of boundary between two materials  #establish location of boundary between two materials
102    x=Function(model).getX()
103  bound = length(x-ic)-r #where the boundary will be located  bound = length(x-ic)-r #where the boundary will be located
104  A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)  kappa = kappai*whereNegative(bound)+kappac*(1-whereNegative(bound))
105  D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)  rhocp = rhocpi*whereNegative(bound)+rhocpc*(1-whereNegative(bound))
106  #define our PDE coeffs  #define our PDE coeffs
107  mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
108  #set initial temperature  #set initial temperature (make sure we use the right sample points)
109  T= Ti*whereNegative(bound)+Tc*wherePositive(bound)  x=Solution(model).getX()
110    bound = length(x-ic)-r #where the boundary will be located
111    T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
112    
113  # rearrage mymesh to suit solution function space for contouring        # rearrage mymesh to suit solution function space for contouring      
114  oldspacecoords=model.getX()  coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
 coords=Data(oldspacecoords, T.getFunctionSpace())  
 #coords = np.array(coords.toListOfTuples())  
 coordX, coordY = toXYTuple(coords)  
115  # create regular grid  # create regular grid
116  xi = np.linspace(0.0,mx,100)  xi = np.linspace(0.0,mx,75)
117  yi = np.linspace(0.0,my,100)  yi = np.linspace(0.0,my, 75)
118    
119  ########################################################START ITERATION  ########################################################START ITERATION
120  while t<=tend:  while t<=tend:
121        i+=1 #counter        i+=1 #counter
122        t+=h #curretn time        t+=h #current time
123        Y = T*D #        mypde.setValue(Y=qH+T*rhocp/h)
       mypde.setValue(Y=Y)  
124        T=mypde.getSolution()        T=mypde.getSolution()
125        tempT = T.toListOfTuples(scalarastuple=False)        print T
126          tempT = T.toListOfTuples()
127        # grid the data.        # grid the data.
128        zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)        zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
129        # contour the gridded data, plotting dots at the        # contour the gridded data, plotting dots at the
# Line 139  while t<=tend: Line 136  while t<=tend:
136        pl.title("Heat diffusion from an intrusion.")        pl.title("Heat diffusion from an intrusion.")
137        pl.xlabel("Horizontal Displacement (m)")        pl.xlabel("Horizontal Displacement (m)")
138        pl.ylabel("Depth (m)")        pl.ylabel("Depth (m)")
139        if getMPIRankWorld() == 0:        pl.savefig(os.path.join(save_path,"heatrefraction%03d.png"%i))
           pl.savefig(os.path.join(save_path,\  
                                    "heatrefraction%03d.png"%i))  
140        pl.clf()                    pl.clf()            
141          print "time step %s at t=%e days completed."%(i,t/day)
142    
143  # compile the *.png files to create an *.avi video that shows T change  # compile the *.png files to create an *.avi video that shows T change
144  # with time. This opperation uses linux mencoder.  # with time. This opperation uses linux mencoder.

Legend:
Removed from v.2944  
changed lines
  Added in v.2945

  ViewVC Help
Powered by ViewVC 1.1.26