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

Contents of /trunk/doc/cookbook/twodheatdiff001.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2401 - (show annotations)
Wed Apr 29 04:23:07 2009 UTC (11 years, 7 months ago) by ahallam
File MIME type: text/x-python
File size: 3371 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
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 = 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
45 ## Intrusion Variables - Granite
46 Ti=2273 # Kelvin #the starting temperature of our intrusion
47 rhoi = 2750 #kg/m^{3} density
48 cpi = 790 #j/kg.k specific heat
49 rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT
50 eta=0. # RADIATION CONDITION
51 kappai=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY
52
53 Tc = 473 # Kelvin #the starting temperature of our country rock
54 rhoc = 2000 #kg/m^{3} density
55 cpc = 920 #j/kg.k specific heat
56 rhocpc = rhoc*cpc #DENSITY * SPECIFIC HEAT
57 kappac = 1.9 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY
58
59
60 #Script/Iteration Related
61 t=0. #our start time, usually zero
62 tday=100*365. #the time we want to end the simulation in days
63 tend=tday*24*60*60
64 outputs = 200 # number of time steps required.
65 h=(tend-t)/outputs #size of time step
66
67 print "Expected Number of Output Files is: ", outputs
68 print "Step size is: ", h/(24.*60*60), "days"
69
70
71 i=0 #loop counter
72 save_path = "data/twodheatdiff" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work
73
74 #... generate domain ...
75 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
76 # extract finite points
77 x=model.getX()
78
79 #... open PDE ...
80 mypde=LinearPDE(model)
81 mypde.setSymmetryOn()
82
83 bound = length(x-ic)-r #where the boundary will be located
84
85 A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)
86 D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)
87
88 mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)
89
90 # ... set initial temperature ....
91
92 T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures.
93 saveVTK(os.path.join(save_path,"dataedge.xml"), sol=bound)
94 saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
95
96 #... start iteration:
97 while t<=tend:
98 i+=1
99 t+=h
100 Y = T*D
101 mypde.setValue(Y=Y)
102 T=mypde.getSolution()
103 saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)

  ViewVC Help
Powered by ViewVC 1.1.26