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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2672 - (show annotations)
Fri Sep 18 01:46:47 2009 UTC (11 years, 5 months ago) by ahallam
File MIME type: text/x-python
File size: 10911 byte(s)
Heat refraction updates : alex's downhole profiles
1 ########################################################
2 #
3 # Copyright (c) 2009 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2009 by University of Queensland
14 Earth Systems Science Computational Center (ESSCC)
15 http://www.uq.edu.au/esscc
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="https://launchpad.net/escript-finley"
20
21 """
22 Author: Antony Hallam antony.hallam@uqconnect.edu.au
23 """
24
25 from esys.pycad import CurveLoop
26
27 # numpy for array handling
28 import numpy as np
29 import pylab as pl
30 # tools for dealing with PDEs - contains locator
31 from esys.escript.pdetools import Locator, Projector
32
33 from esys.escript import *
34 from esys.escript.linearPDEs import LinearPDE
35
36 # routine to find consecutive coordinates of a loop in pycad
37 def getLoopCoords(loop):
38 # return all construction points of input
39 temp = loop.getConstructionPoints()
40 #create a numpy array for xyz components or construction points
41 coords = np.zeros([len(temp),3],float)
42 #place construction points in array
43 for i in range(0,len(temp)):
44 coords[i,:]=temp[i].getCoordinates()
45 #return a numpy array
46 return coords
47
48 # Calculate the location of quivers for a matplotlib plot
49 # quivshape :: [x,y] :: number of quivers in x and y direction
50 # lenxax :: length of model along x
51 # lenyax :: length of model along y
52 # qu :: gradient of escript solution ie grad(T)
53 def toQuivLocs(quivshape,lenxax,lenyax,qu):
54 numquiv = quivshape[0]*quivshape[1] # total number of quivers
55 dx = lenxax/quivshape[0]+1. # quiver x spacing
56 dy = lenyax/quivshape[1]+1. # quiver y spacing
57 qulocs=np.zeros([numquiv,2],float) # memory for quiver locations
58 # fill qulocs
59 for i in range(0,quivshape[0]-1):
60 for j in range(0,quivshape[1]-1):
61 qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j]
62 # retreive values for quivers direction and shape from qu
63 quL = Locator(qu.getFunctionSpace(),qulocs.tolist())
64 qu = quL(qu) #array of dx,dy for quivers
65 qu = np.array(qu) #turn into a numpy array
66 qulocs = quL.getX() #returns actual locations from data
67 qulocs = np.array(qulocs) #turn into a numpy array
68 return qu,qulocs
69
70
71
72
73 # Extract the X and Y coordinates of an array
74 # coords :: escript coordiantes from .getX function
75 def toXYTuple(coords):
76 coords = np.array(coords.toListOfTuples()) #convert to Tuple
77 coordX = coords[:,0]; coordY = coords[:,1] #X and Y components.
78 return coordX,coordY
79
80 def toRegGrid(grid,domain,newx,newy,width,depth):
81 oldspacecoords=domain.getX()
82 gridT = grid.toListOfTuples(scalarastuple=False)
83 cxy = Data(oldspacecoords, grid.getFunctionSpace())
84 cX, cY = toXYTuple(cxy)
85 xi = np.linspace(0.0,width,newx)
86 yi = np.linspace(depth,0.0,newy)
87 # grid the data.
88 zi = pl.matplotlib.mlab.griddata(cX,cY,gridT,xi,yi)
89 return xi,yi,zi
90
91 def gradtoRegGrid(grid,domain,newx,newy,width,depth,dim):
92 cxy = grid.getFunctionSpace().getX()
93 gridT = grid.toListOfTuples()#(scalarastuple=False)
94 #cxy = Data(oldspacecoords, grid.getFunctionSpace())
95 cX, cY = toXYTuple(cxy)
96 xi = np.linspace(0.0,width,newx)
97 yi = np.linspace(depth,0.0,newy)
98
99 gridT = np.array(gridT)
100 gridT = gridT[:,dim]
101
102 # grid the data.
103
104 zi = pl.matplotlib.mlab.griddata(cX,cY,gridT,xi,yi)
105 return xi,yi,zi
106 ########################################################
107 # subroutine: cbphones
108 # Allows us to record the values of a PDE at various
109 # specified locations in the model.
110 # Arguments:
111 # domain : domain of model
112 # U : Current time state displacement solution.
113 # phones : Geophone Locations
114 # dim : model dimesions
115 # savepath: where to output the data files local is default
116 ########################################################
117 def cbphones(domain,U,phones,dim,savepath=""):
118 #find the number of geophones
119 nphones = len(phones)
120 u_pot = np.zeros([nphones,dim],float)
121
122 for i in range(0,nphones):
123 # define the location of the phone source
124 L=Locator(domain,numpy.array(phones[i]))
125 # find potential at point source.
126 temp = L.getValue(U)
127 for j in range(0,dim):
128 u_pot[i,j]=temp[j]
129
130 # open file to save displacement at point source
131 return u_pot
132
133 ########################################################
134 # subroutine: wavesolver2d
135 # Can solve a generic 2D wave propagation problem with a
136 # point source in a homogeneous medium.
137 # Arguments:
138 # domain : domain to solve over
139 # h : time step
140 # tend : end time
141 # lam, mu : lames constants for domain
142 # rho : density of domain
143 # U0 : magnitude of source
144 # xc : source location in domain (Vector)
145 # savepath: where to output the data files
146 ########################################################
147 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
148 from esys.escript.linearPDEs import LinearPDE
149 x=domain.getX()
150 # ... open new PDE ...
151 mypde=LinearPDE(domain)
152 #mypde.setSolverMethod(LinearPDE.LUMPING)
153 mypde.setSymmetryOn()
154 kmat = kronecker(domain)
155 mypde.setValue(D=kmat*rho)
156
157 # define small radius around point xc
158 # Lsup(x) returns the maximum value of the argument x
159 src_radius = 50#2*Lsup(domain.getSize())
160 print "src_radius = ",src_radius
161
162 dunit=numpy.array([0.,1.]) # defines direction of point source
163
164
165 # ... set initial values ....
166 n=0
167 # initial value of displacement at point source is constant (U0=0.01)
168 # for first two time steps
169 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
170 u_m1=u
171 t=0
172
173 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
174 u_pc_x1 = u_pot[0,0]
175 u_pc_y1 = u_pot[0,1]
176 u_pc_x2 = u_pot[1,0]
177 u_pc_y2 = u_pot[1,1]
178 u_pc_x3 = u_pot[2,0]
179 u_pc_y3 = u_pot[2,1]
180
181 # open file to save displacement at point source
182 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
183 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
184
185 while t<tend:
186 # ... get current stress ....
187
188 ##OLD WAY
189 g=grad(u)
190 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
191 ### ... get new acceleration ....
192 #mypde.setValue(X=-stress)
193 #a=mypde.getSolution()
194 ### ... get new displacement ...
195 #u_p1=2*u-u_m1+h*h*a
196 ###NEW WAY
197 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
198 u_p1 = mypde.getSolution()
199 # ... shift displacements ....
200 u_m1=u
201 u=u_p1
202 #stress =
203 t+=h
204 n+=1
205 print n,"-th time step t ",t
206 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
207
208 # print "u at point charge=",u_pc
209 u_pc_x1 = u_pot[0,0]
210 u_pc_y1 = u_pot[0,1]
211 u_pc_x2 = u_pot[1,0]
212 u_pc_y2 = u_pot[1,1]
213 u_pc_x3 = u_pot[2,0]
214 u_pc_y3 = u_pot[2,1]
215
216 # save displacements at point source to file for t > 0
217 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
218
219 # ... save current acceleration in units of gravity and displacements
220 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
221 #displacement = length(u), tensor = stress, Ux = u[0] )
222 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
223
224 u_pc_data.close()
225
226
227 ########################################################
228 # subroutine: wavesolver2d
229 # Can solve a generic 2D wave propagation problem with a
230 # point source in a homogeneous medium with friction.
231 # Arguments:
232 # domain : domain to solve over
233 # h : time step
234 # tend : end time
235 # lam, mu : lames constants for domain
236 # rho : density of domain
237 # U0 : magnitude of source
238 # xc : source location in domain (Vector)
239 # savepath: where to output the data files
240 ########################################################
241 def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
242 x=domain.getX()
243 # ... open new PDE ...
244 mypde=LinearPDE(domain)
245 #mypde.setSolverMethod(LinearPDE.LUMPING)
246 mypde.setSymmetryOn()
247 kmat = kronecker(domain)
248 mypde.setValue(D=kmat)
249 b=0.9
250
251 # define small radius around point xc
252 # Lsup(x) returns the maximum value of the argument x
253 src_radius = 50#2*Lsup(domain.getSize())
254 print "src_radius = ",src_radius
255
256 dunit=numpy.array([0.,1.]) # defines direction of point source
257
258
259 # ... set initial values ....
260 n=0
261 # initial value of displacement at point source is constant (U0=0.01)
262 # for first two time steps
263 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
264 u_m1=u
265 t=0
266
267 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
268 u_pc_x1 = u_pot[0,0]
269 u_pc_y1 = u_pot[0,1]
270 u_pc_x2 = u_pot[1,0]
271 u_pc_y2 = u_pot[1,1]
272 u_pc_x3 = u_pot[2,0]
273 u_pc_y3 = u_pot[2,1]
274
275 # open file to save displacement at point source
276 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
277 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
278
279 while t<tend:
280 # ... get current stress ....
281
282 ##OLD WAY
283 g=grad(u)
284 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
285 ### ... get new acceleration ....
286 #mypde.setValue(X=-stress)
287 #a=mypde.getSolution()
288 ### ... get new displacement ...
289 #u_p1=2*u-u_m1+h*h*a
290 ###NEW WAY
291 y = ((rho/(-rho-b*h))*(u_m1-2*u))+(((b*h)/(-rho-(b*h)))*-u)
292 mypde.setValue(X=-stress*((h*h)/(-rho-h*b)),Y=y)
293 u_p1 = mypde.getSolution()
294 # ... shift displacements ....
295 u_m1=u
296 u=u_p1
297 #stress =
298 t+=h
299 n+=1
300 print n,"-th time step t ",t
301 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
302
303 # print "u at point charge=",u_pc
304 u_pc_x1 = u_pot[0,0]
305 u_pc_y1 = u_pot[0,1]
306 u_pc_x2 = u_pot[1,0]
307 u_pc_y2 = u_pot[1,1]
308 u_pc_x3 = u_pot[2,0]
309 u_pc_y3 = u_pot[2,1]
310
311 # save displacements at point source to file for t > 0
312 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
313
314 # ... save current acceleration in units of gravity and displacements
315 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
316 #displacement = length(u), tensor = stress, Ux = u[0] )
317 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
318
319 u_pc_data.close()
320
321 # joel wrote this to create directories for you
322 import os
323 def needdirs(dirlist):
324 for name in dirlist:
325 if name == '':
326 continue
327 if not os.path.exists(name):
328 try:
329 os.makedirs(name)
330 except OSError:
331 if not os.path.exists(save_path):
332 raise

  ViewVC Help
Powered by ViewVC 1.1.26