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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2658 - (hide annotations)
Thu Sep 10 02:58:44 2009 UTC (9 years, 11 months ago) by ahallam
File MIME type: text/x-python
File size: 5394 byte(s)
Updates to all files scripts to support MPI testing proceedure. Updates to cookbook, new section on functino spaces/domains (needs work). Finalising first 3 chapters for editing.
1 ahallam 2588
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 ahallam 2658 ############################################################FILE HEADER
27     # heatrefraction001.py
28     # Model steady state temperature distribution in two block model, mesh
29     # from heatrefraction_mesher001.py
30    
31     #######################################################EXTERNAL MODULES
32 ahallam 2634 # To solve the problem it is necessary to import the modules we
33     # require.
34     # This imports everything from the escript library
35     from esys.escript import *
36     # This defines LinearPDE as LinearPDE
37     from esys.escript.linearPDEs import LinearPDE, Poisson
38     # This imports the rectangle domain function from finley
39     from esys.finley import Rectangle, ReadMesh, Domain
40     # This package is necessary to handle saving our data.
41     import os
42     # A useful unit handling package which will make sure all our units
43     # match up in the equations.
44     from esys.escript.unitsSI import *
45     # numpy for array handling
46 ahallam 2588 import numpy as np
47 ahallam 2658 import matplotlib
48     #For interactive use, you can comment out the next two lines
49     import matplotlib
50     matplotlib.use('agg') #It's just here for automated testing
51 ahallam 2634 # pylab for matplotlib and plotting
52 ahallam 2588 import pylab as pl
53 ahallam 2634 # cblib functions
54 jfenwick 2648 from cblib import toQuivLocs, toXYTuple, needdirs
55 ahallam 2588
56 ahallam 2658 ########################################################MPI WORLD CHECK
57     if getMPISizeWorld() > 1:
58     import sys
59     print "This example will not run in an MPI world."
60     sys.exit(0)
61    
62     #################################################ESTABLISHING VARIABLES
63 ahallam 2588 qin=70*Milli*W/(m*m) #our heat source temperature is now zero
64     Ti=290.15*K # Kelvin #the starting temperature of our iron bar
65 ahallam 2634 width=5000.0*m
66     depth=-6000.0*m
67 ahallam 2588
68 ahallam 2597 # the folder to gett our outputs from, leave blank "" for script path -
69     # note these depen. are generated from heatrefraction_mesher001.py
70 ahallam 2658 saved_path = save_path= os.path.join("data","heatrefrac001" )
71 jfenwick 2648 needdirs([saved_path])
72    
73 ahallam 2658 ################################################ESTABLISHING PARAMETERS
74 ahallam 2588 ## DOMAIN
75 ahallam 2597 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
76     tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
77 ahallam 2588 tpgx = tpg[:,0]
78     tpgy = tpg[:,1]
79 ahallam 2597 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
80 ahallam 2588 bpgx = bpg[:,0]
81     bpgy = bpg[:,1]
82    
83     # set up kappa (thermal conductivity across domain) using tags
84     kappa=Scalar(0,Function(mymesh))
85     kappa.setTaggedValue("top",2.0)
86     kappa.setTaggedValue("bottom",4.0)
87    
88     #... generate functionspace...
89     #... open PDE ...
90     mypde=LinearPDE(mymesh)
91     #define first coefficient
92     mypde.setValue(A=kappa*kronecker(mymesh))
93    
94     # ... set initial temperature ....
95     x=mymesh.getX()
96    
97 ahallam 2597 qH=Scalar(0,FunctionOnBoundary(mymesh))
98     qH.setTaggedValue("linebottom",qin)
99 ahallam 2588 mypde.setValue(q=whereZero(x[1]),r=Ti)
100 ahallam 2658 mypde.setValue(y=qH)
101 ahallam 2588
102 ahallam 2658 ###########################################################GET SOLUTION
103 ahallam 2588 T=mypde.getSolution()
104    
105 ahallam 2658 ##########################################################VISUALISATION
106     # calculate gradient of solution for quiver plot
107     qu=-kappa*grad(T)
108    
109 ahallam 2588 # rearrage mymesh to suit solution function space
110     oldspacecoords=mymesh.getX()
111     coords=Data(oldspacecoords, T.getFunctionSpace())
112    
113     quivshape = [20,20] #quivers x and quivers y
114 ahallam 2634 # function to calculate quiver locations
115     qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
116 ahallam 2588
117     kappaT = kappa.toListOfTuples(scalarastuple=False)
118     coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
119 ahallam 2634 coordKX, coordKY = toXYTuple(coordsK)
120 ahallam 2588
121     tempT = T.toListOfTuples(scalarastuple=False)
122 ahallam 2634 coordX, coordY = toXYTuple(coords)
123 ahallam 2588
124 ahallam 2634 xi = np.linspace(0.0,width,100)
125     yi = np.linspace(depth,0.0,100)
126 ahallam 2588 # grid the data.
127     zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
128     ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
129 ahallam 2658 # contour the gridded data,
130     # plotting dots at the randomly spaced data points.
131 ahallam 2588
132 ahallam 2658 # select colour
133 ahallam 2588 pl.matplotlib.pyplot.autumn()
134 ahallam 2658 # plot polygons for boundaries
135 ahallam 2588 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
136 ahallam 2658 # contour temperature
137 ahallam 2588 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
138 ahallam 2658 # labels and formatting
139 ahallam 2588 pl.clabel(CS, inline=1, fontsize=8)
140 ahallam 2597 pl.title("Heat Refraction across a clinal structure.")
141     pl.xlabel("Horizontal Displacement (m)")
142     pl.ylabel("Depth (m)")
143 ahallam 2658 if getMPIRankWorld() == 0: #check for MPI processing
144     pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
145 ahallam 2588
146 ahallam 2658 #Quiver Plot qulocs -> tail location, qu -> quiver length/direction
147     QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],\
148     angles='xy',color="white")
149     pl.title("Heat Refraction across a clinal structure \n with\
150     gradient quivers.")
151     if getMPIRankWorld() == 0: #check for MPI processing
152     pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))

  ViewVC Help
Powered by ViewVC 1.1.26