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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2632 - (show annotations)
Wed Aug 26 22:18:19 2009 UTC (9 years, 7 months ago) by ahallam
File MIME type: text/x-python
File size: 5475 byte(s)
Regigger of cookbook directory structure. Examlples->examples/cookbook TEXT->doc/cookbook Figures-> doc/cookbook/figures
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=70*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/heatrefrac001"
44
45 ###### 2 BLOCK MODEL #########
46 ## DOMAIN
47 ## Anticline
48 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
49 tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
50 tpgx = tpg[:,0]
51 tpgy = tpg[:,1]
52 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
53 bpgx = bpg[:,0]
54 bpgy = bpg[:,1]
55 ## Syncline
56 #~ mymesh = ReadMesh("heatrefraction_mesh002.fly")
57
58 # set up kappa (thermal conductivity across domain) using tags
59 kappa=Scalar(0,Function(mymesh))
60 kappa.setTaggedValue("top",2.0)
61 kappa.setTaggedValue("bottom",4.0)
62
63 ###### 3 BLOCK MODEL #########
64 # set up kappa (thermal conductivity across domain) using tags
65 #~ mymesh = ReadMesh("heatrefraction_mesh003.fly")
66 #~ kappa=Scalar(0,Function(mymesh))
67 #~ kappa.setTaggedValue("top",2.0)
68 #~ kappa.setTaggedValue("bottomleft",3.0)
69 #~ kappa.setTaggedValue("bottomright",4.0)
70
71 #... generate functionspace...
72 #... open PDE ...
73 mypde=LinearPDE(mymesh)
74 #define first coefficient
75 mypde.setValue(A=kappa*kronecker(mymesh))
76
77 # ... set initial temperature ....
78 x=mymesh.getX()
79
80 qH=Scalar(0,FunctionOnBoundary(mymesh))
81 qH.setTaggedValue("linebottom",qin)
82 mypde.setValue(q=whereZero(x[1]),r=Ti)
83 mypde.setValue(y=qH)#,r=17*Celsius)
84
85 # get steady state solution and export to vtk.
86 T=mypde.getSolution()
87 #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
88
89 # rearrage mymesh to suit solution function space
90 oldspacecoords=mymesh.getX()
91 coords=Data(oldspacecoords, T.getFunctionSpace())
92
93 # calculate gradient of solution for quiver plot
94 qu=-kappa*grad(T)
95
96 # create quiver locations
97 quivshape = [20,20] #quivers x and quivers y
98 numquiv = quivshape[0]*quivshape[1] # total number of quivers
99 maxx = 5000. # maximum x displacement
100 dx = maxx/quivshape[0]+1. # quiver x spacing
101 maxy = -6000. # maximum y displacement
102 dy = maxy/quivshape[1]+1. # quiver y spacing
103 qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
104 # fill qulocs
105 for i in range(0,quivshape[0]-1):
106 for j in range(0,quivshape[1]-1):
107 qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
108 # retreive values for quivers direction and shape from qu
109 quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
110 qu = quL(qu) #array of dx,dy for quivers
111 qu = np.array(qu) #turn into a numpy array
112 qulocs = quL.getX() #returns actual locations from data
113 qulocs = np.array(qulocs) #turn into a numpy array
114
115 kappaT = kappa.toListOfTuples(scalarastuple=False)
116 coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
117 coordsK = np.array(coordsK.toListOfTuples())
118 coordKX = coordsK[:,0]
119 coordKY = coordsK[:,1]
120
121 coords = np.array(coords.toListOfTuples())
122 tempT = T.toListOfTuples(scalarastuple=False)
123
124 coordX = coords[:,0]
125 coordY = coords[:,1]
126
127 xi = np.linspace(0.0,5000.0,100)
128 yi = np.linspace(-6000.0,0.0,100)
129 # grid the data.
130 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
131 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
132 # contour the gridded data, plotting dots at the randomly spaced data points.
133
134 pl.matplotlib.pyplot.autumn()
135 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
136 #~ CK = pl.contourf(xi,yi,ziK,2)
137 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
138 pl.clabel(CS, inline=1, fontsize=8)
139 pl.title("Heat Refraction across a clinal structure.")
140 pl.xlabel("Horizontal Displacement (m)")
141 pl.ylabel("Depth (m)")
142 #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
143 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
144
145 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
146 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
147 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))
148
149 pl.clf()
150 CKL = pl.fill(tpgx,tpgy,'white',bpgx,bpgy,'white',zorder=-1000)
151 pl.savefig('model.png')

  ViewVC Help
Powered by ViewVC 1.1.26