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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2597 - (hide annotations)
Thu Aug 6 03:30:09 2009 UTC (11 years, 5 months ago) by ahallam
Original Path: trunk/doc/cookbook/heatrefraction002.py
File MIME type: text/x-python
File size: 5196 byte(s)
Updates to: cookbook documentation; heat refraction problems; mencoder investigation for *.png animation
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     # 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, Poisson # This defines LinearPDE as LinearPDE
29     from esys.finley import Rectangle, ReadMesh, Domain # This imports the rectangle domain function from finley
30     import os #This package is necessary to handle saving our data.
31     from esys.escript.unitsSI import * # A useful unit handling package which will make sure all our units match up in the equations.
32     import numpy as np
33     import pylab as pl
34    
35     from esys.escript.pdetools import *
36    
37     ##ESTABLISHING VARIABLES
38 ahallam 2597 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
39 ahallam 2588 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
40    
41 ahallam 2597 # the folder to gett our outputs from, leave blank "" for script path -
42     # note these depen. are generated from heatrefraction_mesher001.py
43     saved_path = "data/heatrefrac002"
44    
45     ###### 2 BLOCK MODEL #########
46     ## DOMAIN
47     ## Anticline
48     mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh003.fly"))
49     tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
50     tpgx = tpg[:,0]
51     tpgy = tpg[:,1]
52     bpgl = np.loadtxt(os.path.join(saved_path,"botpgl"))
53     bpglx = bpgl[:,0]
54     bpgly = bpgl[:,1]
55     bpgr = np.loadtxt(os.path.join(saved_path,"botpgr"))
56     bpgrx = bpgr[:,0]
57     bpgry = bpgr[:,1]
58    
59 ahallam 2588 # set up kappa (thermal conductivity across domain) using tags
60     kappa=Scalar(0,Function(mymesh))
61 ahallam 2597 kappa.setTaggedValue("top",2.0)
62     kappa.setTaggedValue("bottomleft",18.0)
63     kappa.setTaggedValue("bottomright",6.0)
64 ahallam 2588
65     #... generate functionspace...
66     #... open PDE ...
67     mypde=LinearPDE(mymesh)
68 ahallam 2597 mypde.setSymmetryOn()
69 ahallam 2588 #define first coefficient
70     mypde.setValue(A=kappa*kronecker(mymesh))
71    
72     # ... set initial temperature ....
73     x=mymesh.getX()
74    
75 ahallam 2597 qH=Scalar(0,FunctionOnBoundary(mymesh))
76     qH.setTaggedValue("linebottom",qin)
77     mypde.setValue(q=whereZero(x[1]),r=Ti)
78     mypde.setValue(y=qH)
79 ahallam 2588
80     # get steady state solution and export to vtk.
81     T=mypde.getSolution()
82 ahallam 2597 #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
83 ahallam 2588
84     # rearrage mymesh to suit solution function space
85     oldspacecoords=mymesh.getX()
86     coords=Data(oldspacecoords, T.getFunctionSpace())
87    
88     # calculate gradient of solution for quiver plot
89     qu=-kappa*grad(T)
90    
91     # create quiver locations
92     quivshape = [20,20] #quivers x and quivers y
93     numquiv = quivshape[0]*quivshape[1] # total number of quivers
94     maxx = 5000. # maximum x displacement
95     dx = maxx/quivshape[0]+1. # quiver x spacing
96     maxy = -6000. # maximum y displacement
97     dy = maxy/quivshape[1]+1. # quiver y spacing
98     qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
99     # fill qulocs
100     for i in range(0,quivshape[0]-1):
101     for j in range(0,quivshape[1]-1):
102     qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
103     # retreive values for quivers direction and shape from qu
104     quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
105     qu = quL(qu) #array of dx,dy for quivers
106     qu = np.array(qu) #turn into a numpy array
107     qulocs = quL.getX() #returns actual locations from data
108     qulocs = np.array(qulocs) #turn into a numpy array
109    
110     kappaT = kappa.toListOfTuples(scalarastuple=False)
111     coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
112     coordsK = np.array(coordsK.toListOfTuples())
113     coordKX = coordsK[:,0]
114     coordKY = coordsK[:,1]
115    
116     coords = np.array(coords.toListOfTuples())
117     tempT = T.toListOfTuples(scalarastuple=False)
118    
119     coordX = coords[:,0]
120     coordY = coords[:,1]
121    
122     xi = np.linspace(0.0,5000.0,100)
123     yi = np.linspace(-6000.0,0.0,100)
124     # grid the data.
125     zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
126     ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
127     # contour the gridded data, plotting dots at the randomly spaced data points.
128    
129     pl.matplotlib.pyplot.autumn()
130 ahallam 2597 CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
131     #~ CK = pl.contourf(xi,yi,ziK,2)
132 ahallam 2588 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
133     pl.clabel(CS, inline=1, fontsize=8)
134 ahallam 2597 pl.title("Heat Refraction across an anisotropic structure.")
135     pl.xlabel("Horizontal Displacement (m)")
136     pl.ylabel("Depth (m)")
137     #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
138     pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
139 ahallam 2588
140 ahallam 2597 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
141     pl.title("Heat Refraction across an anisotropic structure \n with gradient quivers.")
142     pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))

  ViewVC Help
Powered by ViewVC 1.1.26