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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2397 - (hide annotations)
Fri Apr 24 04:24:28 2009 UTC (10 years, 6 months ago) by ahallam
File MIME type: text/x-python
File size: 3024 byte(s)
End Week 4, Finalised 1D examples, need work on text but script ok.
Begining of 2D example builds on 1D. Script semi-finalised, text to complete.
1 ahallam 2393
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 ahallam 2397 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 ahallam 2393 q=0 #our heat source temperature is now zero
42 ahallam 2397 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 ahallam 2393 #Script/Iteration Related
49 ahallam 2397 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 ahallam 2393
55 ahallam 2397 print "Expected Number of Output Files is: ", outputs
56     print "Step size is: ", h/(24.*60*60), "days"
57 ahallam 2393
58 ahallam 2397
59 ahallam 2393 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 ahallam 2397 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
64 ahallam 2393 # extract finite points
65 ahallam 2397 x=model.getX()
66 ahallam 2393 #... open PDE ...
67 ahallam 2397 mypde=LinearPDE(model)
68 ahallam 2393 mypde.setSymmetryOn()
69 ahallam 2397 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h,d=eta,y=eta*Tref)
70 ahallam 2393
71     # ... set initial temperature ....
72 ahallam 2397 bound = x[0]-mx/(ndx/250.)
73     T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound)
74 ahallam 2393
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 ahallam 2397 mypde.setValue(Y=rhocp/h*T)
82 ahallam 2393 T=mypde.getSolution()
83     saveVTK(os.path.join(save_path,"data%03d.xml") %i,sol=T)
84    
85 ahallam 2397 # 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 ahallam 2393
96 ahallam 2397 #print "Finish temp balance in:", i*h/(24.*60*60), " days"

  ViewVC Help
Powered by ViewVC 1.1.26