/[escript]/trunk/doc/cookbook/onedheatdiff002.py
ViewVC logotype

Annotation of /trunk/doc/cookbook/onedheatdiff002.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2401 - (hide annotations)
Wed Apr 29 04:23:07 2009 UTC (12 years, 1 month ago) by ahallam
File MIME type: text/x-python
File size: 3024 byte(s)
End Week 5:
Renamed *.tex files to more logical, changed to chapter format as HEAT_DIFFUSION.tex
Editted and Reviewed text portion and some python scripts for heat diffusion.
Looking to push for finalisation on heat diffusion in the next couple of week to begin next modelling problem/Chapter.
1 ahallam 2401
2     ########################################################
3     #
4     # Copyright (c) 2003-2009 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
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     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="https://launchpad.net/escript-finley"
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 = 500 # model lenght
37     my = 100 # model width
38     ndx = 500 # steps in x direction
39     ndy = 1 # steps in y direction
40    
41     q=0 #our heat source temperature is now zero
42     Tref=2273 # Kelvin #the starting temperature of our intrusion
43     rho = 2750 #kg/m^{3} density
44     cp = 790 #j/kg specific heat
45     rhocp = rho*cp #DENSITY * SPECIFIC HEAT
46     eta=0. # RADIATION CONDITION
47     kappa=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY
48     #Script/Iteration Related
49     t=0. #our start time, usually zero
50     tday=10*365. #the time we want to end the simulation in days
51     tend=tday*24*60*60
52     outputs = 400 # number of time steps required.
53     h=(tend-t)/outputs #size of time step
54    
55     print "Expected Number of Output Files is: ", outputs
56     print "Step size is: ", h/(24.*60*60), "days"
57    
58    
59     i=0 #loop counter
60     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
61    
62     #... generate domain ...
63     model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
64     # extract finite points
65     x=model.getX()
66     #... open PDE ...
67     mypde=LinearPDE(model)
68     mypde.setSymmetryOn()
69     mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
70    
71     # ... set initial temperature ....
72     bound = x[0]-mx/(ndx/250.)
73     T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
74    
75     saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
76    
77     # ... start iteration:
78     while t<=tend:
79     i+=1
80     t+=h
81     mypde.setValue(Y=rhocp/h*T)
82     T=mypde.getSolution()
83     saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
84    
85     # iteration var 2
86     #Tl = 0
87     #Tr = Tref
88     #while Tl < Tr*0.8:
89     #mypde.setValue(Y=rhocp/h*T)
90     #T=mypde.getSolution()
91     #i+=1
92     #x=rod.getX()
93     #Tl= x[0]
94     #Tr= x[
95    
96     #print "Finish temp balance in:", i*h/(24.*60*60), " days"

  ViewVC Help
Powered by ViewVC 1.1.26