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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 10001 byte(s)
Don't panic.
Updating copyright stamps

1 ########################################################
2 #
3 # Copyright (c) 2009-2010 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-2010 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 A collection of routines to use in cookbook examples.
23 Author: Antony Hallam antony.hallam@uqconnect.edu.au
24 """
25
26 from esys.pycad import CurveLoop
27
28 # numpy for array handling
29 import numpy as np
30
31 # tools for dealing with PDEs - contains locator
32 from esys.escript.pdetools import Locator, Projector
33
34 from esys.escript import *
35 from esys.escript.linearPDEs import LinearPDE
36
37 # routine to find consecutive coordinates of a loop in pycad
38 def getLoopCoords(loop):
39 # return all construction points of input
40 temp = loop.getConstructionPoints()
41 #create a numpy array for xyz components or construction points
42 coords = np.zeros([len(temp),3],float)
43 #place construction points in array
44 for i in range(0,len(temp)):
45 coords[i,:]=temp[i].getCoordinates()
46 #return a numpy array
47 return coords
48
49
50 ########################################################
51 # subroutine: cbphones
52 # Allows us to record the values of a PDE at various
53 # specified locations in the model.
54 # Arguments:
55 # domain : domain of model
56 # U : Current time state displacement solution.
57 # phones : Geophone Locations
58 # dim : model dimesions
59 # savepath: where to output the data files local is default
60 ########################################################
61 def cbphones(domain,U,phones,dim,savepath=""):
62 #find the number of geophones
63 nphones = len(phones)
64 u_pot = np.zeros([nphones,dim],float)
65
66 for i in range(0,nphones):
67 # define the location of the phone source
68 L=Locator(domain,numpy.array(phones[i]))
69 # find potential at point source.
70 temp = L.getValue(U)
71 for j in range(0,dim):
72 u_pot[i,j]=temp[j]
73
74 # open file to save displacement at point source
75 return u_pot
76
77 ########################################################
78 # subroutine: wavesolver2d
79 # Can solve a generic 2D wave propagation problem with a
80 # point source in a homogeneous medium.
81 # Arguments:
82 # domain : domain to solve over
83 # h : time step
84 # tend : end time
85 # lam, mu : lames constants for domain
86 # rho : density of domain
87 # U0 : magnitude of source
88 # xc : source location in domain (Vector)
89 # savepath: where to output the data files
90 ########################################################
91 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath,output="vtk"):
92 from esys.escript.linearPDEs import LinearPDE
93 x=domain.getX()
94 # ... open new PDE ...
95 mypde=LinearPDE(domain)
96 #mypde.setSolverMethod(LinearPDE.LUMPING)
97 mypde.setSymmetryOn()
98 kmat = kronecker(domain)
99 mypde.setValue(D=kmat*rho)
100
101 # define small radius around point xc
102 # Lsup(x) returns the maximum value of the argument x
103 src_radius = 50#2*Lsup(domain.getSize())
104 print "src_radius = ",src_radius
105
106 dunit=numpy.array([0.,1.]) # defines direction of point source
107
108
109 # ... set initial values ....
110 n=0
111 # initial value of displacement at point source is constant (U0=0.01)
112 # for first two time steps
113 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
114 u_m1=u
115 t=0
116
117 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
118 u_pc_x1 = u_pot[0,0]
119 u_pc_y1 = u_pot[0,1]
120 u_pc_x2 = u_pot[1,0]
121 u_pc_y2 = u_pot[1,1]
122 u_pc_x3 = u_pot[2,0]
123 u_pc_y3 = u_pot[2,1]
124
125 # open file to save displacement at point source
126 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
127 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))
128
129 # while t<tend:
130 while t<1.:
131
132 # ... get current stress ....
133 t=1.
134 ##OLD WAY
135 break
136 g=grad(u)
137 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
138 ### ... get new acceleration ....
139 #mypde.setValue(X=-stress)
140 #a=mypde.getSolution()
141 ### ... get new displacement ...
142 #u_p1=2*u-u_m1+h*h*a
143 ###NEW WAY
144 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
145 u_p1 = mypde.getSolution()
146 # ... shift displacements ....
147 u_m1=u
148 u=u_p1
149 #stress =
150 t+=h
151 n+=1
152 print n,"-th time step t ",t
153 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
154
155 # print "u at point charge=",u_pc
156 u_pc_x1 = u_pot[0,0]
157 u_pc_y1 = u_pot[0,1]
158 u_pc_x2 = u_pot[1,0]
159 u_pc_y2 = u_pot[1,1]
160 u_pc_x3 = u_pot[2,0]
161 u_pc_y3 = u_pot[2,1]
162
163 # save displacements at point source to file for t > 0
164 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))
165
166 # ... save current acceleration in units of gravity and displacements
167 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
168 #displacement = length(u), tensor = stress, Ux = u[0] )
169 if output == "vtk":
170 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
171 else:
172 quT=qu.toListOfTuples()
173 #Projector is used to smooth the data.
174 proj=Projector(mymesh)
175 smthT=proj(T)
176
177 #move data to a regular grid for plotting
178 xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth)
179
180 # contour the gridded data,
181 # select colour
182 pl.matplotlib.pyplot.autumn()
183 pl.clf()
184 # contour temperature
185 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
186 # labels and formatting
187 pl.clabel(CS, inline=1, fontsize=8)
188 pl.title("Heat Refraction across a clinal structure.")
189 pl.xlabel("Horizontal Displacement (m)")
190 pl.ylabel("Depth (m)")
191 pl.legend()
192 if getMPIRankWorld() == 0: #check for MPI processing
193 pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png"))
194
195 u_pc_data.close()
196
197
198 ########################################################
199 # subroutine: wavesolver2d
200 # Can solve a generic 2D wave propagation problem with a
201 # point source in a homogeneous medium with friction.
202 # Arguments:
203 # domain : domain to solve over
204 # h : time step
205 # tend : end time
206 # lam, mu : lames constants for domain
207 # rho : density of domain
208 # U0 : magnitude of source
209 # xc : source location in domain (Vector)
210 # savepath: where to output the data files
211 ########################################################
212 def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
213 x=domain.getX()
214 # ... open new PDE ...
215 mypde=LinearPDE(domain)
216 #mypde.setSolverMethod(LinearPDE.LUMPING)
217 mypde.setSymmetryOn()
218 kmat = kronecker(domain)
219 mypde.setValue(D=kmat)
220 b=0.9
221
222 # define small radius around point xc
223 # Lsup(x) returns the maximum value of the argument x
224 src_radius = 50#2*Lsup(domain.getSize())
225 print "src_radius = ",src_radius
226
227 dunit=numpy.array([0.,1.]) # defines direction of point source
228
229
230 # ... set initial values ....
231 n=0
232 # initial value of displacement at point source is constant (U0=0.01)
233 # for first two time steps
234 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
235 u_m1=u
236 t=0
237
238 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
239 u_pc_x1 = u_pot[0,0]
240 u_pc_y1 = u_pot[0,1]
241 u_pc_x2 = u_pot[1,0]
242 u_pc_y2 = u_pot[1,1]
243 u_pc_x3 = u_pot[2,0]
244 u_pc_y3 = u_pot[2,1]
245
246 # open file to save displacement at point source
247 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
248 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))
249
250 while t<tend:
251 # ... get current stress ....
252
253 ##OLD WAY
254 g=grad(u)
255 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
256 ### ... get new acceleration ....
257 #mypde.setValue(X=-stress)
258 #a=mypde.getSolution()
259 ### ... get new displacement ...
260 #u_p1=2*u-u_m1+h*h*a
261 ###NEW WAY
262 y = ((rho/(-rho-b*h))*(u_m1-2*u))+(((b*h)/(-rho-(b*h)))*-u)
263 mypde.setValue(X=-stress*((h*h)/(-rho-h*b)),Y=y)
264 u_p1 = mypde.getSolution()
265 # ... shift displacements ....
266 u_m1=u
267 u=u_p1
268 #stress =
269 t+=h
270 n+=1
271 print n,"-th time step t ",t
272 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
273
274 # print "u at point charge=",u_pc
275 u_pc_x1 = u_pot[0,0]
276 u_pc_y1 = u_pot[0,1]
277 u_pc_x2 = u_pot[1,0]
278 u_pc_y2 = u_pot[1,1]
279 u_pc_x3 = u_pot[2,0]
280 u_pc_y3 = u_pot[2,1]
281
282 # save displacements at point source to file for t > 0
283 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))
284
285 # ... save current acceleration in units of gravity and displacements
286 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
287 #displacement = length(u), tensor = stress, Ux = u[0] )
288 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
289
290 u_pc_data.close()
291
292 # joel wrote this to create directories for you
293 # It only needs to be this complicated to accomodate our
294 # automated tests
295 import os
296 def needdirs(dirlist):
297 for name in dirlist:
298 if name == '':
299 continue
300 if not os.path.exists(name):
301 try:
302 os.makedirs(name)
303 except OSError:
304 if not os.path.exists(save_path):
305 raise

  ViewVC Help
Powered by ViewVC 1.1.26