/[escript]/branches/symbolic_from_3470/dudley/test/python/seismic_wave.py
ViewVC logotype

Annotation of /branches/symbolic_from_3470/dudley/test/python/seismic_wave.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (hide annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 11935 byte(s)
Fast forward to latest trunk revision 3788.

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

  ViewVC Help
Powered by ViewVC 1.1.26