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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 12106 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26