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

Contents of /trunk/doc/examples/cookbook/twodheatdiff.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (9 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 3192 byte(s)
Renamed the main cookbook tex file to match our convention.
Replaced doc/cookbook/figures/heatrefraction002contqu.pdf with
a version which is actually pdf. However it needs to be regnerated since
it it sideways.

The examples have had their copyright notices fixed (dates were too early).
sb2.py has been removed since it uses pyvisi.

scons will now build the cookbook as parts of a docs build.
Also in reposnse to :
scons cookbook_pdf


1
2 ########################################################
3 #
4 # Copyright (c) 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) 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 from cblib import needdirs
33
34
35 ##ESTABLISHING VARIABLES
36 #PDE related
37 mx = 600 # model lenght
38 my = 600 # model width
39 ndx = 100 # steps in x direction
40 ndy = 100 # steps in y direction
41 r = 200 # radius of intrusion
42 ic = [300, 0] #centre of intrusion
43
44 q=0 #our heat source temperature is now zero
45 Ti=2273 # Kelvin #the starting temperature of our iron bar
46 rhoi = 2750 #kg/m^{3} density
47 cpi = 790 #j/kg specific heat
48 rhocpi = rhoi*cpi #DENSITY * SPECIFIC HEAT
49 eta=0. # RADIATION CONDITION
50 kappai=2.2 # Watts/(meter*Kelvin) DIFFUSION CONSTANT, HEAT PERMEABILITY
51
52 Tc = 200
53 rhoc = 2200
54 cpc = 400
55 rhocpc = rhoc*cpc
56 kappac = 0.1
57
58
59 #Script/Iteration Related
60 t=0. #our start time, usually zero
61 tday=100*365. #the time we want to end the simulation in days
62 tend=tday*24*60*60
63 outputs = 200 # number of time steps required.
64 h=(tend-t)/outputs #size of time step
65
66 print "Expected Number of Output Files is: ", outputs
67 print "Step size is: ", h/(24.*60*60), "days"
68
69
70 i=0 #loop counter
71 save_path = "data/twodheatdiff" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work
72 needdirs([save_path])
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.vtu"), sol=bound)
94 saveVTK(os.path.join(save_path,"data%03d.vtu") %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.vtu") %i,sol=T)

  ViewVC Help
Powered by ViewVC 1.1.26