/[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 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: 4766 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
27 # require.
28 # This imports everything from the escript library
29 from esys.escript import *
30 # This defines LinearPDE as LinearPDE
31 from esys.escript.linearPDEs import LinearPDE, Poisson
32 # This imports the rectangle domain function from finley
33 from esys.finley import Rectangle, ReadMesh, Domain
34 # This package is necessary to handle saving our data.
35 import os
36 # A useful unit handling package which will make sure all our units
37 # match up in the equations.
38 from esys.escript.unitsSI import *
39 # numpy for array handling
40 import numpy as np
41 # pylab for matplotlib and plotting
42 import pylab as pl
43 # cblib functions
44 from cblib import toQuivLocs, toXYTuple, needdirs
45
46 ##ESTABLISHING VARIABLES
47 qin=70*Milli*W/(m*m) #our heat source temperature is now zero
48 Ti=290.15*K # Kelvin #the starting temperature of our iron bar
49 width=5000.0*m
50 depth=-6000.0*m
51
52 # the folder to gett our outputs from, leave blank "" for script path -
53 # note these depen. are generated from heatrefraction_mesher001.py
54 saved_path = "data/heatrefrac001"
55
56 needdirs([saved_path])
57
58 ###### 2 BLOCK MODEL #########
59 ## DOMAIN
60 ## Anticline
61 mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly"))
62 tpg = np.loadtxt(os.path.join(saved_path,"toppg"))
63 tpgx = tpg[:,0]
64 tpgy = tpg[:,1]
65 bpg = np.loadtxt(os.path.join(saved_path,"botpg"))
66 bpgx = bpg[:,0]
67 bpgy = bpg[:,1]
68 ## Syncline
69 #~ mymesh = ReadMesh("heatrefraction_mesh002.fly")
70
71 # set up kappa (thermal conductivity across domain) using tags
72 kappa=Scalar(0,Function(mymesh))
73 kappa.setTaggedValue("top",2.0)
74 kappa.setTaggedValue("bottom",4.0)
75
76 ###### 3 BLOCK MODEL #########
77 # set up kappa (thermal conductivity across domain) using tags
78 #~ mymesh = ReadMesh("heatrefraction_mesh003.fly")
79 #~ kappa=Scalar(0,Function(mymesh))
80 #~ kappa.setTaggedValue("top",2.0)
81 #~ kappa.setTaggedValue("bottomleft",3.0)
82 #~ kappa.setTaggedValue("bottomright",4.0)
83
84 #... generate functionspace...
85 #... open PDE ...
86 mypde=LinearPDE(mymesh)
87 #define first coefficient
88 mypde.setValue(A=kappa*kronecker(mymesh))
89
90 # ... set initial temperature ....
91 x=mymesh.getX()
92
93 qH=Scalar(0,FunctionOnBoundary(mymesh))
94 qH.setTaggedValue("linebottom",qin)
95 mypde.setValue(q=whereZero(x[1]),r=Ti)
96 mypde.setValue(y=qH)#,r=17*Celsius)
97
98 # get steady state solution and export to vtk.
99 T=mypde.getSolution()
100 #saveVTK("tempheatrefract.xml",sol=T, q=-kappa*grad(T))
101
102 # rearrage mymesh to suit solution function space
103 oldspacecoords=mymesh.getX()
104 coords=Data(oldspacecoords, T.getFunctionSpace())
105
106 # calculate gradient of solution for quiver plot
107 qu=-kappa*grad(T)
108 quivshape = [20,20] #quivers x and quivers y
109 # function to calculate quiver locations
110 qu,qulocs = toQuivLocs(quivshape,width,depth,qu)
111
112 kappaT = kappa.toListOfTuples(scalarastuple=False)
113 coordsK = Data(oldspacecoords, kappa.getFunctionSpace())
114 coordKX, coordKY = toXYTuple(coordsK)
115
116 tempT = T.toListOfTuples(scalarastuple=False)
117 coordX, coordY = toXYTuple(coords)
118
119 xi = np.linspace(0.0,width,100)
120 yi = np.linspace(depth,0.0,100)
121 # grid the data.
122 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
123 ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi)
124 # contour the gridded data, plotting dots at the randomly spaced data points.
125
126 pl.matplotlib.pyplot.autumn()
127
128 CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=-1000)
129 #~ CK = pl.contourf(xi,yi,ziK,2)
130 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
131
132 pl.clabel(CS, inline=1, fontsize=8)
133 pl.title("Heat Refraction across a clinal structure.")
134 pl.xlabel("Horizontal Displacement (m)")
135 pl.ylabel("Depth (m)")
136 #~ CB = pl.colorbar(CS, shrink=0.8, extend='both')
137 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
138
139 QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white")
140 pl.title("Heat Refraction across a clinal structure \n with gradient quivers.")
141 pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png"))

  ViewVC Help
Powered by ViewVC 1.1.26