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

  ViewVC Help
Powered by ViewVC 1.1.26