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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2648 - (show annotations)
Mon Sep 7 00:06:15 2009 UTC (11 years, 2 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
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 from cblib import needdirs
37
38 ##ESTABLISHING VARIABLES
39 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
40 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
41
42 # 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 needdirs([saved_path])
46
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 # set up kappa (thermal conductivity across domain) using tags
62 kappa=Scalar(0,Function(mymesh))
63 kappa.setTaggedValue("top",2.0)
64 kappa.setTaggedValue("bottomleft",18.0)
65 kappa.setTaggedValue("bottomright",6.0)
66
67 #... generate functionspace...
68 #... open PDE ...
69 mypde=LinearPDE(mymesh)
70 mypde.setSymmetryOn()
71 #define first coefficient
72 mypde.setValue(A=kappa*kronecker(mymesh))
73
74 # ... set initial temperature ....
75 x=mymesh.getX()
76
77 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
82 # get steady state solution and export to vtk.
83 T=mypde.getSolution()
84 #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
85
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 CKL = pl.fill(tpgx,tpgy,'brown',bpglx,bpgly,'red',bpgrx,bpgry,'orange',zorder=-1000)
133 #~ CK = pl.contourf(xi,yi,ziK,2)
134 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
135 pl.clabel(CS, inline=1, fontsize=8)
136 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
142 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