/[escript]/trunk/doc/examples/cookbook/onedheatdiff001.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/onedheatdiff001.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2589 - (hide annotations)
Wed Aug 5 23:06:35 2009 UTC (11 years, 5 months ago) by ahallam
Original Path: trunk/doc/cookbook/onedheatdiff001.py
File MIME type: text/x-python
File size: 3067 byte(s)


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 ahallam 2494 from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
31 ahallam 2589 import pylab as pl
32     import numpy as np
33 ahallam 2401 import os #This package is necessary to handle saving our data.
34 ahallam 2494 #plotting tools
35     #import matplotlib as ptool
36 ahallam 2401
37    
38     ##ESTABLISHING VARIABLES
39     #Domain related.
40 ahallam 2494 mx = 1*m #meters - model lenght
41     my = .1*m #meters - model width
42 ahallam 2401 ndx = 100 # steps in x direction
43     ndy = 1 # steps in y direction
44    
45     #PDE related
46 ahallam 2494 q=200. * Celsius #Kelvin - our heat source temperature
47     Tref = 0. * Celsius # Kelvin - starting temp of iron bar
48     rho = 7874. *kg/m**3 #kg/m^{3} density of iron
49     cp = 449.*J/(kg*K) #j/Kg.K thermal capacity
50 ahallam 2401 rhocp = rho*cp
51 ahallam 2494 kappa = 80.*W/m/K #watts/m.Kthermal conductivity
52 ahallam 2401 #Script/Iteration Related
53     t=0 #our start time, usually zero
54 ahallam 2494 tend=5.*minute #seconds - time to end simulation
55 ahallam 2401 outputs = 200 # number of time steps required.
56     h=(tend-t)/outputs #size of time step
57 ahallam 2494 print "Expected Number of time outputs is: ", (tend-t)/h
58 ahallam 2401 i=0 #loop counter
59     #the folder to put our outputs in, leave blank "" for script path
60 ahallam 2494 save_path="data/onedheatdiff001"
61 ahallam 2401 #note this folder path must exist to work
62    
63 ahallam 2494
64 ahallam 2401 #... generate domain ...
65     rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
66     # extract finite points
67     x=rod.getX()
68     #... open PDE ...
69     mypde=LinearPDE(rod)
70     mypde.setSymmetryOn()
71 ahallam 2494 mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h)
72 ahallam 2401
73     # ... set heat source: ....
74     qH=q*whereZero(x[0])
75     # ... set initial temperature ....
76     T=Tref
77    
78     # ... start iteration:
79     while t<=tend:
80 ahallam 2589 i+=1
81     t+=h
82     mypde.setValue(Y=qH+rhocp/h*T)
83     T=mypde.getSolution()
84     #ptool.plot(T)
85     saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)
86     plx = x.toListOfTuples()
87     plx = np.array(plx)
88     plx = plx[:,0]
89     tempT = T.toListOfTuples(scalarastuple=False)
90     pl.plot(plx,tempT)
91     pl.title("Temperatuer accross Rod")
92     #~ pl.axis([0.0,1.0,0.0,0.004*2.73146e2])
93     pl.savefig(os.path.join(save_path,"rodpyplot%03d.png") %i)
94     pl.clf()
95 ahallam 2401

  ViewVC Help
Powered by ViewVC 1.1.26