# Contents of /trunk/doc/examples/cookbook/example03a.py

Revision 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (4 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 5957 byte(s)
```python3ified things, replaced mixed whitespace and xrange calls
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2009-2013 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 __copyright__="""Copyright (c) 2009-2013 by University of Queensland 17 http://www.uq.edu.au 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 """ 24 Author: Antony Hallam antony.hallam@uqconnect.edu.au 25 """ 26 27 ############################################################FILE HEADER 28 # example03a.py 29 # Model temperature diffusion between a granite intrusion and sandstone 30 # country rock. This is a two dimensional problem with the granite as a 31 # heat source. 32 33 #######################################################EXTERNAL MODULES 34 #To solve the problem it is necessary to import the modules we require. 35 #For interactive use, you can comment out the next two lines 36 import matplotlib 37 matplotlib.use('agg') #It's just here for automated testing 38 #This imports everything from the escript library 39 from esys.escript import * 40 # This defines the LinearPDE module as LinearPDE 41 from esys.escript.linearPDEs import LinearPDE 42 # This imports the rectangle domain 43 from esys.finley import Rectangle 44 # A useful unit handling package which will make sure all our units 45 # match up in the equations under SI. 46 from esys.escript.unitsSI import * 47 import pylab as pl #Plotting package. 48 import numpy as np #Array package. 49 import os #This package is necessary to handle saving our data. 50 from cblib import toXYTuple 51 52 ########################################################MPI WORLD CHECK 53 if getMPISizeWorld() > 1: 54 import sys 55 print("This example will not run in an MPI world.") 56 sys.exit(0) 57 58 #################################################ESTABLISHING VARIABLES 59 #PDE related 60 mx = 600*m #meters - model length 61 my = 600*m #meters - model width 62 ndx = 150 #mesh steps in x direction 63 ndy = 150 #mesh steps in y direction 64 r = 200*m #meters - radius of intrusion 65 ic = [300*m, 0] #centre of intrusion (meters) 66 qH=0.*J/(sec*m**3) #our heat source temperature is now zero 67 68 ## Intrusion Variables - Granite 69 Ti=2273.*Celsius # Kelvin -the starting temperature of our RHS Block 70 rhoi = 2750*kg/m**3 #kg/m^{3} density of granite 71 cpi = 790.*J/(kg*K) #j/Kg.K thermal capacity 72 rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT 73 kappai=2.2*W/m/K #watts/m.K thermal conductivity 74 ## Country Rock Variables - Sandstone 75 Tc = 473*Celsius # Kelvin #the starting temperature of our country rock 76 rhoc = 2000*kg/m**3 #kg/m^{3} density 77 cpc = 920.*J/(kg*K) #j/kg.k specific heat 78 rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT 79 kappac = 1.9*W/m/K #watts/m.K thermal conductivity 80 81 #Script/Iteration Related 82 t=0. #our start time, usually zero 83 tend=200.* yr #the time we want to end the simulation 84 outputs = 200 # number of time steps required. 85 h=(tend-t)/outputs #size of time step 86 #user warning 87 print("Expected Number of Output Files is: ", outputs) 88 print("Step size is: ", h/day, "days") 89 i=0 #loop counter 90 #the folder to put our outputs in, leave blank "" for script path 91 save_path= os.path.join("data","example03") 92 mkDir(save_path) 93 ########## note this folder path must exist to work ################### 94 95 ################################################ESTABLISHING PARAMETERS 96 #generate domain using rectangle 97 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 98 #extract finite points - the solution points 99 #create the PDE 100 mypde=LinearPDE(model) #assigns a domain to our PDE 101 mypde.setSymmetryOn() #set the fast solver on for symmetry 102 #establish location of boundary between two materials 103 x=Function(model).getX() 104 bound = length(x-ic)-r #where the boundary will be located 105 kappa = kappai*whereNegative(bound)+kappac*(1-whereNegative(bound)) 106 rhocp = rhocpi*whereNegative(bound)+rhocpc*(1-whereNegative(bound)) 107 #define our PDE coeffs 108 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 109 #set initial temperature (make sure we use the right sample points) 110 x=Solution(model).getX() 111 bound = length(x-ic)-r #where the boundary will be located 112 T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound)) 113 114 # rearrage mymesh to suit solution function space for contouring 115 coordX, coordY = toXYTuple(T.getFunctionSpace().getX()) 116 # create regular grid 117 xi = np.linspace(0.0,mx,75) 118 yi = np.linspace(0.0,my, 75) 119 120 ########################################################START ITERATION 121 while t<=tend: 122 i+=1 #counter 123 t+=h #current time 124 mypde.setValue(Y=qH+T*rhocp/h) 125 T=mypde.getSolution() 126 tempT = T.toListOfTuples() 127 # grid the data. 128 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 129 # contour the gridded data, plotting dots at the 130 # randomly spaced data points. 131 pl.matplotlib.pyplot.autumn() 132 pl.contourf(xi,yi,zi,10) 133 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 134 pl.clabel(CS, inline=1, fontsize=8) 135 pl.axis([0,600,0,600]) 136 pl.title("Heat diffusion from an intrusion.") 137 pl.xlabel("Horizontal Displacement (m)") 138 pl.ylabel("Depth (m)") 139 pl.savefig(os.path.join(save_path,"temp%03d.png"%i)) 140 pl.clf() 141 print("time step %s at t=%e days completed."%(i,t/day)) 142 143 #########################################################CREATE A MOVIE 144 # compile the *.png files to create an *.avi video that shows T change 145 # with time. This opperation uses linux mencoder. 146 os.system("mencoder mf://"+save_path+"/*.png -mf type=png:\ 147 w=800:h=600:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o \ 148 example03tempT.avi")

 ViewVC Help Powered by ViewVC 1.1.26