/[escript]/trunk/doc/examples/cookbook/example03a.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example03a.py

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

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 41  matplotlib.use('agg') #It's just here fo Line 41  matplotlib.use('agg') #It's just here fo
41  from esys.escript import *  from esys.escript import *
42  # This defines the LinearPDE module as LinearPDE  # This defines the LinearPDE module as LinearPDE
43  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
 # This imports the rectangle domain  
 from esys.finley import Rectangle  
44  # A useful unit handling package which will make sure all our units  # A useful unit handling package which will make sure all our units
45  # match up in the equations under SI.  # match up in the equations under SI.
46  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
# Line 50  import pylab as pl #Plotting package. Line 48  import pylab as pl #Plotting package.
48  import numpy as np #Array package.  import numpy as np #Array package.
49  import os #This package is necessary to handle saving our data.  import os #This package is necessary to handle saving our data.
50  from cblib import toXYTuple  from cblib import toXYTuple
51    try:
52        from esys.finley import Rectangle
53        HAVE_FINLEY = True
54    except ImportError:
55        print("Finley module not available")
56        HAVE_FINLEY = False
57  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
58  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
59          import sys          import sys
60          print("This example will not run in an MPI world.")          print("This example will not run in an MPI world.")
61          sys.exit(0)          sys.exit(0)
62    
63  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
64  #PDE related      #################################################ESTABLISHING VARIABLES
65  mx = 600*m #meters - model length      #PDE related
66  my = 600*m #meters - model width      mx = 600*m #meters - model length
67  ndx = 150 #mesh steps in x direction      my = 600*m #meters - model width
68  ndy = 150 #mesh steps in y direction      ndx = 150 #mesh steps in x direction
69  r = 200*m #meters - radius of intrusion      ndy = 150 #mesh steps in y direction
70  ic = [300*m, 0] #centre of intrusion (meters)      r = 200*m #meters - radius of intrusion
71  qH=0.*J/(sec*m**3) #our heat source temperature is now zero      ic = [300*m, 0] #centre of intrusion (meters)
72        qH=0.*J/(sec*m**3) #our heat source temperature is now zero
73  ## Intrusion Variables - Granite  
74  Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block      ## Intrusion Variables - Granite
75  rhoi = 2750*kg/m**3 #kg/m^{3} density of granite      Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block
76  cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity      rhoi = 2750*kg/m**3 #kg/m^{3} density of granite
77  rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT      cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity
78  kappai=2.2*W/m/K #watts/m.K thermal conductivity      rhocpi = rhoi*cpi   #DENSITY * SPECIFIC HEAT
79  ## Country Rock Variables - Sandstone      kappai=2.2*W/m/K #watts/m.K thermal conductivity
80  Tc = 473*Celsius # Kelvin #the starting temperature of our country rock      ## Country Rock Variables - Sandstone
81  rhoc = 2000*kg/m**3 #kg/m^{3} density      Tc = 473*Celsius # Kelvin #the starting temperature of our country rock
82  cpc = 920.*J/(kg*K) #j/kg.k specific heat      rhoc = 2000*kg/m**3 #kg/m^{3} density
83  rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT      cpc = 920.*J/(kg*K) #j/kg.k specific heat
84  kappac = 1.9*W/m/K #watts/m.K thermal conductivity      rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT
85        kappac = 1.9*W/m/K #watts/m.K thermal conductivity
86  #Script/Iteration Related  
87  t=0. #our start time, usually zero      #Script/Iteration Related
88  tend=200.* yr #the time we want to end the simulation      t=0. #our start time, usually zero
89  outputs = 200 # number of time steps required.      tend=200.* yr #the time we want to end the simulation
90  h=(tend-t)/outputs #size of time step      outputs = 200 # number of time steps required.
91  #user warning      h=(tend-t)/outputs #size of time step
92  print("Expected Number of Output Files is: ", outputs)      #user warning
93  print("Step size is: ", h/day, "days")      print("Expected Number of Output Files is: ", outputs)
94  i=0 #loop counter      print("Step size is: ", h/day, "days")
95  #the folder to put our outputs in, leave blank "" for script path      i=0 #loop counter
96  save_path= os.path.join("data","example03")      #the folder to put our outputs in, leave blank "" for script path
97  mkDir(save_path)      save_path= os.path.join("data","example03")
98  ########## note this folder path must exist to work ###################      mkDir(save_path)
99        ########## note this folder path must exist to work ###################
100  ################################################ESTABLISHING PARAMETERS  
101  #generate domain using rectangle      ################################################ESTABLISHING PARAMETERS
102  model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)      #generate domain using rectangle
103  #extract finite points - the solution points      model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
104  #create the PDE      #extract finite points - the solution points
105  mypde=LinearPDE(model) #assigns a domain to our PDE      #create the PDE
106  mypde.setSymmetryOn() #set the fast solver on for symmetry      mypde=LinearPDE(model) #assigns a domain to our PDE
107  #establish location of boundary between two materials      mypde.setSymmetryOn() #set the fast solver on for symmetry
108  x=Function(model).getX()      #establish location of boundary between two materials
109  bound = length(x-ic)-r #where the boundary will be located      x=Function(model).getX()
110  kappa = kappai*whereNegative(bound)+kappac*(1-whereNegative(bound))      bound = length(x-ic)-r #where the boundary will be located
111  rhocp = rhocpi*whereNegative(bound)+rhocpc*(1-whereNegative(bound))      kappa = kappai*whereNegative(bound)+kappac*(1-whereNegative(bound))
112  #define our PDE coeffs      rhocp = rhocpi*whereNegative(bound)+rhocpc*(1-whereNegative(bound))
113  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)      #define our PDE coeffs
114  #set initial temperature (make sure we use the right sample points)      mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
115  x=Solution(model).getX()      #set initial temperature (make sure we use the right sample points)
116  bound = length(x-ic)-r #where the boundary will be located      x=Solution(model).getX()
117  T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))      bound = length(x-ic)-r #where the boundary will be located
118        T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
119  # rearrage mymesh to suit solution function space for contouring        
120  coordX, coordY = toXYTuple(T.getFunctionSpace().getX())      # rearrage mymesh to suit solution function space for contouring      
121  # create regular grid      coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
122  xi = np.linspace(0.0,mx,75)      # create regular grid
123  yi = np.linspace(0.0,my, 75)      xi = np.linspace(0.0,mx,75)
124        yi = np.linspace(0.0,my, 75)
125  ########################################################START ITERATION  
126  while t<=tend:      ########################################################START ITERATION
127        i+=1 #counter      while t<=tend:
128        t+=h #current time            i+=1 #counter
129        mypde.setValue(Y=qH+T*rhocp/h)            t+=h #current time
130        T=mypde.getSolution()            mypde.setValue(Y=qH+T*rhocp/h)
131        tempT = T.toListOfTuples()            T=mypde.getSolution()
132        # grid the data.            tempT = T.toListOfTuples()
133        zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)            # grid the data.
134        # contour the gridded data, plotting dots at the            zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
135        # randomly spaced data points.            # contour the gridded data, plotting dots at the
136        pl.matplotlib.pyplot.autumn()            # randomly spaced data points.
137        pl.contourf(xi,yi,zi,10)            pl.matplotlib.pyplot.autumn()
138        CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')            pl.contourf(xi,yi,zi,10)
139        pl.clabel(CS, inline=1, fontsize=8)            CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
140        pl.axis([0,600,0,600])            pl.clabel(CS, inline=1, fontsize=8)
141        pl.title("Heat diffusion from an intrusion.")            pl.axis([0,600,0,600])
142        pl.xlabel("Horizontal Displacement (m)")            pl.title("Heat diffusion from an intrusion.")
143        pl.ylabel("Depth (m)")            pl.xlabel("Horizontal Displacement (m)")
144        pl.savefig(os.path.join(save_path,"temp%03d.png"%i))            pl.ylabel("Depth (m)")
145        pl.clf()                        pl.savefig(os.path.join(save_path,"temp%03d.png"%i))
146        print("time step %s at t=%e days completed."%(i,t/day))            pl.clf()            
147              print("time step %s at t=%e days completed."%(i,t/day))
148  #########################################################CREATE A MOVIE  
149  # compile the *.png files to create an *.avi video that shows T change      #########################################################CREATE A MOVIE
150  # with time. This opperation uses linux mencoder.      # compile the *.png files to create an *.avi video that shows T change
151  os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\      # with time. This opperation uses linux mencoder.
152        os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\
153  w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \  w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \
154  example03tempT.avi")  example03tempT.avi")

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26