# Contents of /trunk/doc/cookbook/twodheatdiff.py

Revision 2534 - (show annotations)
Thu Jul 16 06:49:19 2009 UTC (9 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 3153 byte(s)
```Changed examples, tests and tutorials to save VTK files as .vtu instead .xml.
Visit doesn't know what to do with xml's and vtu is the proper extension
anyway.

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2009 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 22 """ 23 Author: Antony Hallam antony.hallam@uqconnect.edu.au 24 """ 25 26 # 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 28 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE 29 from esys.finley import Rectangle # This imports the rectangle domain function from finley 30 import os #This package is necessary to handle saving our data. 31 32 33 34 ##ESTABLISHING VARIABLES 35 #PDE related 36 mx = 600 # model lenght 37 my = 600 # model width 38 ndx = 100 # steps in x direction 39 ndy = 100 # steps in y direction 40 r = 200 # radius of intrusion 41 ic = [300, 0] #centre of intrusion 42 43 q=0 #our heat source temperature is now zero 44 Ti=2273 # Kelvin #the starting temperature of our iron bar 45 rhoi = 2750 #kg/m^{3} density 46 cpi = 790 #j/kg specific heat 47 rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT 48 eta=0. # RADIATION CONDITION 49 kappai=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY 50 51 Tc = 200 52 rhoc = 2200 53 cpc = 400 54 rhocpc = rhoc*cpc 55 kappac = 0.1 56 57 58 #Script/Iteration Related 59 t=0. #our start time, usually zero 60 tday=100*365. #the time we want to end the simulation in days 61 tend=tday*24*60*60 62 outputs = 200 # number of time steps required. 63 h=(tend-t)/outputs #size of time step 64 65 print "Expected Number of Output Files is: ", outputs 66 print "Step size is: ", h/(24.*60*60), "days" 67 68 69 i=0 #loop counter 70 save_path = "data/twodheatdiff" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work 71 72 #... generate domain ... 73 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 74 # extract finite points 75 x=model.getX() 76 77 #... open PDE ... 78 mypde=LinearPDE(model) 79 mypde.setSymmetryOn() 80 81 bound = length(x-ic)-r #where the boundary will be located 82 83 A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound) 84 D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound) 85 86 mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc) 87 88 # ... set initial temperature .... 89 90 T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures. 91 saveVTK(os.path.join(save_path,"dataedge.vtu"), sol=bound) 92 saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T) 93 94 #... start iteration: 95 while t<=tend: 96 i+=1 97 t+=h 98 Y = T*D 99 mypde.setValue(Y=Y) 100 T=mypde.getSolution() 101 saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)