/[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 2648 - (hide annotations)
Mon Sep 7 00:06:15 2009 UTC (11 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 5246 byte(s)
ecording some changes related to Anthony's cookbook.
These tests still do not run.

test_heatref.py is the current problem.
Current problems include:
undefined symbols
a dependence on X somehow.

I've already addressed the matplotlib X thing somewhere else.
Find and apply.


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

  ViewVC Help
Powered by ViewVC 1.1.26