/[escript]/trunk/finley/test/python/seismic_wave.py
ViewVC logotype

Contents of /trunk/finley/test/python/seismic_wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 4 months ago) by ksteube
File MIME type: text/x-python
File size: 12256 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1 #
2 # $Id$
3 #
4 #######################################################
5 #
6 # Copyright 2003-2007 by ACceSS MNRF
7 # Copyright 2007 by University of Queensland
8 #
9 # http://esscc.uq.edu.au
10 # Primary Business: Queensland, Australia
11 # Licensed under the Open Software License version 3.0
12 # http://www.opensource.org/licenses/osl-3.0.php
13 #
14 #######################################################
15 #
16
17 """
18 seismic wave propagation
19
20 @var __author__: name of author
21 @var __licence__: licence agreement
22 @var __url__: url entry point on documentation
23 @var __version__: version
24 @var __date__: date of the version
25 """
26
27 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
28 http://www.access.edu.au
29 Primary Business: Queensland, Australia"""
30 __license__="""Licensed under the Open Software License version 3.0
31 http://www.opensource.org/licenses/osl-3.0.php"""
32 __author__="Lutz Gross, l.gross@uq.edu.au"
33 __url__="http://www.iservo.edu.au/esys/escript"
34 __version__="$Revision$"
35 __date__="$Date$"
36
37 from esys.escript import *
38 from esys.escript.linearPDEs import LinearPDE
39 from esys.finley import Brick
40 import time
41
42 WORKDIR="/raid2/lutz/waves/"
43 output=True
44 n_end=10000
45
46 resolution=1000. # number of elements per m in the finest region
47 resolution=400. # number of elements per m in the finest region
48 o=1 # element order
49
50 l=100000. # width and length m (without obsorber)
51 h=30000. # height in m (without obsorber)
52 d_absorber=l*0.10 # thickness of absorbing layer
53
54 l_sand=20000. # thickness of sand region on surface
55 h_sand=5000. # thickness of sand layer under the water
56
57 l_x_water=10000. # length of water in x
58 l_y_water=10000. # length of water in y
59 h_water=2000. # depth of water region
60
61 x_sand=l/2-l_x_water/2-l_sand # x coordinate of location of sand region (without obsorber)
62 y_sand=l/2-l_y_water/2-l_sand # y coordinate of location of sand region (without obsorber)
63
64
65 # origin
66 origin={"x": -d_absorber, "y" : -d_absorber , "z" : -h-d_absorber }
67 # location and geometrical size of event reltive to origin:
68 xc=[l*0.2,l*0.3,-h*0.7]
69 src_radius = 2*resolution
70 # direction of event:
71 event=numarray.array([0.,0.,1.])*1.e6
72 # time and length of the event
73 tc=2.
74 tc_length=0.5
75
76 # material properties:
77 bedrock=0
78 absorber=1
79 water=2
80 sand=3
81
82 rho_tab={}
83 rho_tab[bedrock]=8e3
84 rho_tab[absorber]=rho_tab[bedrock]
85 rho_tab[water]=1e3
86 rho_tab[sand]=5e3
87
88 mu_tab={}
89 mu_tab[bedrock]=1.7e11
90 mu_tab[absorber]=mu_tab[bedrock]
91 mu_tab[water]=0.
92 mu_tab[sand]=1.5e10
93
94 lmbd_tab={}
95 lmbd_tab[bedrock]=1.7e11
96 lmbd_tab[absorber]=lmbd_tab[bedrock]
97 lmbd_tab[water]=1.e9
98 lmbd_tab[sand]=1.5e10
99
100 eta_tab={}
101 eta_tab[absorber]=-log(0.05)*sqrt(rho_tab[absorber]*(lmbd_tab[absorber]+2*mu_tab[absorber]))/d_absorber
102 eta_tab[sand]=eta_tab[absorber]/40.
103 eta_tab[water]=eta_tab[absorber]/40.
104 eta_tab[bedrock]=eta_tab[absorber]/40.
105
106
107 # material properties:
108 bedrock=0
109 absorber=1
110 water=2
111 sand=3
112
113 rho={}
114 rho[bedrock]=8e3
115 rho[absorber]=rho[bedrock]
116 rho[water]=1e3
117 rho[sand]=5e3
118
119 mu={}
120 mu[bedrock]=1.7e11
121 mu[absorber]=mu[bedrock]
122 mu[water]=0.
123 mu[sand]=1.5e10
124
125 lmbd={}
126 lmbd[bedrock]=1.7e11
127 lmbd_absorber=lmbd[bedrock]
128 lmbd[water]=1.e9
129 lmbd[sand]=1.5e10
130
131 eta={}
132 eta[absorber]=-log(0.05)*sqrt(rho[absorber]*(lmbd_absorber+2*mu[absorber]))/d_absorber
133 eta[sand]=eta[absorber]/40.
134 eta[water]=eta[absorber]/40.
135 eta[bedrock]=eta[absorber]/40.
136
137 if output:
138 print "event location = ",xc
139 print "radius of event = ",src_radius
140 print "time of event = ",tc
141 print "length of event = ",tc_length
142 print "direction = ",event
143
144 t_end=30.
145 dt_write=0.1
146
147
148 def getDomain():
149 """
150 this defines a dom as a brick of length and width l and hight h
151
152
153 """
154 global netotal
155
156 v_p={}
157 for tag in rho_tab.keys():
158 v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])
159 v_p_ref=min(v_p.values())
160 print "velocities: bedrock = %s, sand = %s, water =%s, absorber =%s, reference =%s"%(v_p[bedrock],v_p[sand],v_p[water],v_p[absorber],v_p_ref)
161
162 sections={}
163 sections["x"]=[d_absorber, x_sand, l_sand, l_x_water, l_sand, l-x_sand-2*l_sand-l_x_water, d_absorber]
164 sections["y"]=[d_absorber, y_sand, l_sand, l_y_water, l_sand, l-y_sand-2*l_sand-l_y_water, d_absorber]
165 sections["z"]=[d_absorber,h-h_water-h_sand,h_sand,h_water]
166 if output:
167 print "sections x = ",sections["x"]
168 print "sections y = ",sections["y"]
169 print "sections z = ",sections["z"]
170
171 mats= [
172 [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
173 [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
174 [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
175 [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
176 [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
177 [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
178 [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
179
180 [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
181 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
182 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
183 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
184 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
185 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
186 [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
187
188 [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
189 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
190 [absorber, bedrock , sand , sand , sand , bedrock , absorber],
191 [absorber, bedrock , sand , sand , sand , bedrock , absorber],
192 [absorber, bedrock , sand , sand , sand , bedrock , absorber],
193 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
194 [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
195
196 [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
197 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
198 [absorber, bedrock , sand , sand , sand , bedrock , absorber],
199 [absorber, bedrock , sand , water , sand , bedrock , absorber],
200 [absorber, bedrock , sand , sand , sand , bedrock , absorber],
201 [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
202 [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ] ]
203
204 num_elem={}
205 for d in sections:
206 num_elem[d]=[]
207 for i in range(len(sections[d])):
208 if d=="x":
209 v_p_min=v_p[mats[0][0][i]]
210 for q in range(len(sections["y"])):
211 for r in range(len(sections["z"])):
212 v_p_min=min(v_p[mats[r][q][i]],v_p_min)
213 elif d=="y":
214 v_p_min=v_p[mats[0][i][0]]
215 for q in range(len(sections["x"])):
216 for r in range(len(sections["z"])):
217 v_p_min=min(v_p[mats[r][i][q]],v_p_min)
218 elif d=="z":
219 v_p_min=v_p[mats[i][0][0]]
220 for q in range(len(sections["x"])):
221 for r in range(len(sections["y"])):
222 v_p_min=min(v_p[mats[i][r][q]],v_p_min)
223 num_elem[d].append(max(1,int(sections[d][i] * v_p_ref/v_p_min /resolution+0.5)))
224
225 ne_x=sum(num_elem["x"])
226 ne_y=sum(num_elem["y"])
227 ne_z=sum(num_elem["z"])
228 netotal=ne_x*ne_y*ne_z
229 if output: print "grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal)
230 dom=Brick(ne_x,ne_y,ne_z,l0=o*ne_x,l1=o*ne_y,l2=o*ne_z,order=o)
231 x_old=dom.getX()
232 x_new=0
233
234 for d in sections:
235 if d=="x":
236 i=0
237 f=[1,0,0]
238 if d=="y":
239 i=1
240 f=[0,1,0]
241 if d=="z":
242 i=2
243 f=[0,0,1]
244 x=x_old[i]
245
246 p=origin[d]
247 ne=0
248 s=0.
249
250 for i in range(len(sections[d])-1):
251 msk=whereNonPositive(x-o*ne+0.5)
252 s=s*msk + (sections[d][i]/(o*num_elem[d][i])*(x-o*ne)+p)*(1.-msk)
253 ne+=num_elem[d][i]
254 p+=sections[d][i]
255 x_new=x_new + s * f
256 dom.setX(x_new)
257
258 fs=Function(dom)
259 x=Function(dom).getX()
260 x0=x[0]
261 x1=x[1]
262 x2=x[2]
263 p_z=origin["z"]
264 for i in range(len(mats)):
265 f_z=wherePositive(x2-p_z)*wherePositive(x2-p_z+sections["z"][i])
266 p_y=origin["y"]
267 for j in range(len(mats[i])):
268 f_y=wherePositive(x1-p_y)*wherePositive(x1-p_z+sections["y"][j])
269 p_x=origin["x"]
270 for k in range(len(mats[i][j])):
271 f_x=wherePositive(x0-p_x)*wherePositive(x0-p_x+sections["x"][k])
272 fs.setTags(mats[i][j][k],f_x*f_y*f_z)
273 p_x+=sections["x"][k]
274 p_y+=sections["y"][j]
275 p_z+=sections["z"][i]
276 return dom
277
278 def getMaterialProperties(dom):
279 rho =Scalar(rho_tab[bedrock],Function(dom))
280 eta =Scalar(eta_tab[bedrock],Function(dom))
281 mu =Scalar(mu_tab[bedrock],Function(dom))
282 lmbd=Scalar(lmbd_tab[bedrock],Function(dom))
283 tags=Scalar(bedrock,Function(dom))
284
285 for tag in rho_tab.keys():
286 rho.setTaggedValue(tag,rho_tab[tag])
287 eta.setTaggedValue(tag,eta_tab[tag])
288 mu.setTaggedValue(tag,mu_tab[tag])
289 lmbd.setTaggedValue(tag,lmbd_tab[tag])
290 tags.setTaggedValue(tag,tag)
291 return rho,mu,lmbd,eta
292
293 def wavePropagation(dom,rho,mu,lmbd,eta):
294 x=Function(dom).getX()
295 # ... open new PDE ...
296 mypde=LinearPDE(dom)
297 mypde.setSolverMethod(LinearPDE.LUMPING)
298 k=kronecker(Function(dom))
299 mypde.setValue(D=k*rho)
300
301 dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))
302 if output: print "time step size = ",dt
303 # ... set initial values ....
304 n=0
305 t=0
306 t_write=0.
307 n_write=0
308 # initial value of displacement at point source is constant (U0=0.01)
309 # for first two time steps
310 u=Vector(0.,Solution(dom))
311 v=Vector(0.,Solution(dom))
312 a=Vector(0.,Solution(dom))
313 a2=Vector(0.,Solution(dom))
314 v=Vector(0.,Solution(dom))
315
316 starttime = time.clock()
317 while t<t_end and n<n_end:
318 if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
319 # prediction:
320 u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2
321 v_pr=v+dt*a+(dt**2/2)*a2
322 a_pr=a+dt*a2
323 # ... get current stress ....
324 eps=symmetric(grad(u_pr))
325 stress=lmbd*trace(eps)*k+2*mu*eps
326 # ... force due to event:
327 if abs(t-tc)<5*tc_length:
328 F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
329 if output: print Lsup(F)
330 else:
331 if output: print 0.
332 # ... get new acceleration ....
333 mypde.setValue(X=-stress,Y=F-eta*v_pr)
334 a=mypde.getSolution()
335 # ... get new displacement ...
336 da=a-a_pr
337 u=u_pr+(dt**2/12.)*da
338 v=v_pr+(5*dt/12.)*da
339 a2+=da/dt
340 # ... save current acceleration in units of gravity and displacements
341 if output:
342 if t>=t_write:
343 saveVTK(WORKDIR+"disp.%i.vtu"%n_write,displacement=u, amplitude=length(u))
344 t_write+=dt_write
345 n_write+=1
346 t+=dt
347 n+=1
348
349 endtime = time.clock()
350 totaltime = endtime-starttime
351 global netotal
352 print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
353 if __name__ =="__main__":
354 dom=getDomain()
355 rho,mu,lmbd,eta=getMaterialProperties(dom)
356 wavePropagation(dom,rho,mu,lmbd,eta)

  ViewVC Help
Powered by ViewVC 1.1.26