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

Revision 2697 - (show annotations)
Wed Sep 30 01:29:22 2009 UTC (11 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 3063 byte(s)
```Trying to get cookbook scripts passing.
Split cblib in two to avoid importing pylab in scripts which
don't otheriwse need it.

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 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) 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 from cblib1 import needdirs 33 34 35 ##ESTABLISHING VARIABLES 36 #PDE related 37 mx = 500 # model lenght 38 my = 100 # model width 39 ndx = 500 # steps in x direction 40 ndy = 1 # steps in y direction 41 42 q=0 #our heat source temperature is now zero 43 Tref=2273 # Kelvin #the starting temperature of our intrusion 44 rho = 2750 #kg/m^{3} density 45 cp = 790 #j/kg specific heat 46 rhocp = rho*cp #DENSITY * SPECIFIC HEAT 47 eta=0. # RADIATION CONDITION 48 kappa=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY 49 #Script/Iteration Related 50 t=0. #our start time, usually zero 51 tday=10*365. #the time we want to end the simulation in days 52 tend=tday*24*60*60 53 outputs = 400 # number of time steps required. 54 h=(tend-t)/outputs #size of time step 55 56 print "Expected Number of Output Files is: ", outputs 57 print "Step size is: ", h/(24.*60*60), "days" 58 59 60 i=0 #loop counter 61 save_path = "data/onedheatdiff_var001" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work 62 needdirs([save_path]) 63 #... generate domain ... 64 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 65 # extract finite points 66 x=model.getX() 67 #... open PDE ... 68 mypde=LinearPDE(model) 69 mypde.setSymmetryOn() 70 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref) 71 72 # ... set initial temperature .... 73 bound = x[0]-mx/(ndx/250.) 74 T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound) 75 76 saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T) 77 78 # ... start iteration: 79 while t<=tend: 80 i+=1 81 t+=h 82 mypde.setValue(Y=rhocp/h*T) 83 T=mypde.getSolution() 84 saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T) 85 86 # iteration var 2 87 #Tl = 0 88 #Tr = Tref 89 #while Tl < Tr*0.8: 90 #mypde.setValue(Y=rhocp/h*T) 91 #T=mypde.getSolution() 92 #i+=1 93 #x=rod.getX() 94 #Tl= x[0] 95 #Tr= x[ 96 97 #print "Finish temp balance in:", i*h/(24.*60*60), " days"

 ViewVC Help Powered by ViewVC 1.1.26