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

Contents of /trunk/doc/cookbook/heatrefraction002.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2597 - (show annotations)
Thu Aug 6 03:30:09 2009 UTC (10 years, 1 month ago) by ahallam
File MIME type: text/x-python
File size: 5196 byte(s)
Updates to: cookbook documentation; heat refraction problems; mencoder investigation for *.png animation
1
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 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
39 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
40
41 # 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 # set up kappa (thermal conductivity across domain) using tags
60 kappa=Scalar(0,Function(mymesh))
61 kappa.setTaggedValue("top",2.0)
62 kappa.setTaggedValue("bottomleft",18.0)
63 kappa.setTaggedValue("bottomright",6.0)
64
65 #... generate functionspace...
66 #... open PDE ...
67 mypde=LinearPDE(mymesh)
68 mypde.setSymmetryOn()
69 #define first coefficient
70 mypde.setValue(A=kappa*kronecker(mymesh))
71
72 # ... set initial temperature ....
73 x=mymesh.getX()
74
75 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
80 # get steady state solution and export to vtk.
81 T=mypde.getSolution()
82 #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
83
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 CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
131 #~ CK = pl.contourf(xi,yi,ziK,2)
132 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
133 pl.clabel(CS, inline=1, fontsize=8)
134 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
140 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