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

revision 2534 by caltinay, Thu Jul 16 06:49:19 2009 UTC revision 2606 by ahallam, Thu Aug 13 04:32:23 2009 UTC
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  # twodheatdiff002.py
28  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  # Model temperature diffusion between a granite intrusion and sandstone
29  from esys.finley import Rectangle # This imports the rectangle domain function from finley  # country rock. This is a two dimensional problem with the granite as a
30    # heat source.
31
32    #######################################################EXTERNAL MODULES
33    #To solve the problem it is necessary to import the modules we require.
34    #This imports everything from the escript library
35    from esys.escript import *
36    # This defines the LinearPDE module as LinearPDE
37    from esys.escript.linearPDEs import LinearPDE
38    # This imports the rectangle domain function from finley.
39    from esys.finley import Rectangle
40    # A useful unit handling package which will make sure all our units
41    # match up in the equations under SI.
42    from esys.escript.unitsSI import *
43    import pylab as pl #Plotting package.
44    import numpy as np #Array package.
45  import os #This package is necessary to handle saving our data.  import os #This package is necessary to handle saving our data.
46
47    #################################################ESTABLISHING VARIABLES

##ESTABLISHING VARIABLES
48  #PDE related  #PDE related
49  mx = 600 # model lenght  mx = 600*m #meters - model length
50  my = 600 # model width  my = 600*m #meters - model width
51  ndx = 100 # steps in x direction  ndx = 100 #mesh steps in x direction
52  ndy = 100 # steps in y direction  ndy = 100 #mesh steps in y direction
53  r = 200 # radius of intrusion  r = 200*m #meters - radius of intrusion
54  ic = [300, 0] #centre of intrusion  ic = [300, 0] #centre of intrusion (meters)
55    q=0.*Celsius #our heat source temperature is now zero
q=0 #our heat source temperature is now zero
56
57  ## Intrusion Variables - Granite  ## Intrusion Variables - Granite
58  Ti=2273 # Kelvin #the starting temperature of our intrusion  Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
59  rhoi = 2750 #kg/m^{3} density  rhoi = 2750*kg/m**3 #kg/m^{3} density of granite
60  cpi = 790 #j/kg.k specific heat  cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity
61  rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT  rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT
63  kappai=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY  kappai=2.2*W/m/K #watts/m.K thermal conductivity
64    ## Country Rock Variables - Sandstone
65  Tc = 473 # Kelvin #the starting temperature of our country rock  Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
66  rhoc = 2000 #kg/m^{3} density  rhoc = 2000*kg/m**3 #kg/m^{3} density
67  cpc = 920 #j/kg.k specific heat  cpc = 920.*J/(kg*K) #j/kg.k specific heat
68  rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT  rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT
69  kappac = 1.9 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY  kappac = 1.9*W/m/K #watts/m.K thermal conductivity

70
71  #Script/Iteration Related  #Script/Iteration Related
72  t=0. #our start time, usually zero  t=0. #our start time, usually zero
# Line 63  tday=100*365. #the time we want to end t Line 74  tday=100*365. #the time we want to end t
74  tend=tday*24*60*60  tend=tday*24*60*60
75  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
76  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
77    #user warning
78  print "Expected Number of Output Files is: ", outputs  print "Expected Number of Output Files is: ", outputs
79  print "Step size is: ", h/(24.*60*60), "days"  print "Step size is: ", h/(24.*60*60), "days"

80  i=0 #loop counter  i=0 #loop counter
81  save_path = "data/twodheatdiff" #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
82    save_path="data/twodheatdiff"
83    ########## note this folder path must exist to work ###################
84
85  #... generate domain ...  ################################################ESTABLISHING PARAMETERS
86    #generate domain using rectangle
87  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
88  # extract finite points  #extract finite points - the solution points
89  x=model.getX()  x=model.getX()
90    #create the PDE
91  #... open PDE ...  mypde=LinearPDE(model) #assigns a domain to our PDE
92  mypde=LinearPDE(model)  mypde.setSymmetryOn() #set the fast solver on for symmetry
93  mypde.setSymmetryOn()  #establish location of boundary between two materials

94  bound = length(x-ic)-r #where the boundary will be located  bound = length(x-ic)-r #where the boundary will be located

95  A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)  A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)
96  D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)  D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)
97    #define our PDE coeffs
98  mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)  mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)
99    #set initial temperature
100    T= Ti*whereNegative(bound)+Tc*wherePositive(bound)
101
102  # ... set initial temperature ....  # rearrage mymesh to suit solution function space for contouring
103    oldspacecoords=model.getX()
104  T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures.  coords=Data(oldspacecoords, T.getFunctionSpace())
105  saveVTK(os.path.join(save_path,"dataedge.vtu"), sol=bound)  coords = np.array(coords.toListOfTuples())
106  saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)  coordX = coords[:,0]
107    coordY = coords[:,1]
108    # create regular grid
109    xi = np.linspace(0.0,600.0,100)
110    yi = np.linspace(0.0,600.0,100)
111
112  #... start iteration:  #... start iteration:
113  while t<=tend:  while t<=tend:
# Line 100  while t<=tend: Line 116  while t<=tend:
116        Y = T*D        Y = T*D
117        mypde.setValue(Y=Y)        mypde.setValue(Y=Y)
118        T=mypde.getSolution()        T=mypde.getSolution()
119        saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)        tempT = T.toListOfTuples(scalarastuple=False)
120          # grid the data.
121          zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
122          # contour the gridded data, plotting dots at the randomly spaced data points.
123          pl.matplotlib.pyplot.autumn()
124          pl.contourf(xi,yi,zi,10)
125          CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
126          pl.clabel(CS, inline=1, fontsize=8)
127          pl.axis([0,600,0,600])
128          pl.title("Heat diffusion from an intrusion.")
129          pl.xlabel("Horizontal Displacement (m)")
130          pl.ylabel("Depth (m)")
131          pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)
132          pl.clf()
133    #     saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
134
135    # compile the *.png files to create an *.avi video that shows T change
136    # with time. This opperation uses linux mencoder.
137    os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
138    w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
139    twodheatdiff001tempT.avi")

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