/[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 850 - (show annotations)
Sun Sep 17 23:27:00 2006 UTC (13 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 7425 byte(s)
some problems with vtk writer fixed
1 """
2 seismic wave propagation
3
4 @var __author__: name of author
5 @var __licence__: licence agreement
6 @var __url__: url entry point on documentation
7 @var __version__: version
8 @var __date__: date of the version
9 """
10
11 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12 http://www.access.edu.au
13 Primary Business: Queensland, Australia"""
14 __license__="""Licensed under the Open Software License version 3.0
15 http://www.opensource.org/licenses/osl-3.0.php"""
16 __author__="Lutz Gross, l.gross@uq.edu.au"
17 __url__="http://www.iservo.edu.au/esys/escript"
18 __version__="$Revision$"
19 __date__="$Date$"
20
21 from esys.escript import *
22 from esys.escript.linearPDEs import LinearPDE
23 from esys.finley import Brick
24 import time
25
26 output=True
27 n_end=10000
28
29 resolution=4000. # number of elements per m
30 l=100000. # width and length m
31 h=30000. # height in m
32 o=1 # element order
33
34 rho_bedrock=8e3
35 mu_bedrock=1.7e11
36 lambda_bedrock=1.7e11
37
38 l_x_water=10000 # length of water in x
39 l_y_water=10000 # length of water in y
40 h_water=max(2000,resolution) # depth of water region
41
42 water_tag=2
43 rho_water=1e3
44 mu_water=0.
45 lambda_water=1.e9
46
47 d0_sand=10000 # thickness of sand region on surface
48 d_sand=max(2000,resolution) # thickness of sand layer under the water
49
50 sand_tag=3
51 rho_sand=5e3
52 mu_sand=1.5e10
53 lambda_sand=1.5e10
54
55
56 # location and geometrical size of event:
57 xc=[30000.,40000.,10000.]
58 src_radius = 0.1*h
59 # direction of event:
60 event=numarray.array([1.,0.,0.])*1.e6
61 # time and length of the event
62 tc_length=2.
63 tc=sqrt(5.*tc_length)
64 if output:
65 print "event location = ",xc
66 print "radius of event = ",src_radius
67 print "time of event = ",tc
68 print "length of event = ",tc_length
69 print "direction = ",event
70
71 t_end=20.
72
73 def getDomain():
74 """
75 this defines a dom as a brick of length and width l and hight h
76
77
78 """
79 global netotal
80 v_p_bedrock=sqrt((2*mu_bedrock+lambda_bedrock)/rho_bedrock)
81 v_p_sand=sqrt((2*mu_sand+lambda_sand)/rho_sand)
82 v_p_water=sqrt((2*mu_water+lambda_water)/rho_water)
83
84 print v_p_bedrock,v_p_sand,v_p_water
85
86 ne_l_x_bed=int((l-d0_sand-l_x_water)/resolution+0.5)
87 ne_l_y_bed=int((l-d0_sand-l_y_water)/resolution+0.5)
88 ne_h_bed=int((h-d_sand-h_water)/resolution+0.5)
89
90 ne_l_sand=int(d0_sand*v_p_bedrock/v_p_sand/resolution+0.5)
91 ne_h_sand=int(d_sand*v_p_bedrock/v_p_sand/resolution+0.5)
92
93 ne_l_x_water=int(l_x_water*v_p_bedrock/v_p_water/resolution+0.5)
94 ne_l_y_water=int(l_y_water*v_p_bedrock/v_p_water/resolution+0.5)
95 ne_h_water=int(h_water*v_p_bedrock/v_p_water/resolution+0.5)
96
97 ne_l_x=ne_l_x_bed+ne_l_sand+ne_l_x_water
98 ne_l_y=ne_l_y_bed+ne_l_sand+ne_l_y_water
99 ne_h=ne_h_bed+ne_h_sand+ne_h_water
100
101 print ne_l_x,ne_l_x_bed,ne_l_sand,ne_l_x_water
102 print ne_l_y,ne_l_y_bed,ne_l_sand,ne_l_y_water
103 print ne_h,ne_h_bed,ne_h_sand,ne_h_water
104
105 netotal=ne_l_x*ne_l_y*ne_h
106 if output: print "grid : %s x %s x %s (%s elements)"%(ne_l_x,ne_l_y,ne_h,netotal)
107 dom=Brick(ne_l_x,ne_l_y,ne_h,l0=o*ne_l_x,l1=o*ne_l_y,l2=o*ne_h,order=o)
108 x=dom.getX()
109
110 x0=x[0]
111 x0_new = (l-d0_sand-l_x_water)/(o*ne_l_x_bed)*x0
112 msk=whereNonPositive(x0-o*ne_l_x_bed)
113 x0_new = x0_new * msk + (d0_sand/(o*ne_l_sand)*(x0-o*ne_l_x_bed)+l-d0_sand-l_x_water)*(1.-msk)
114 msk=whereNonPositive(x0-o*(ne_l_x_bed+ne_l_sand))
115 x0_new = x0_new * msk + (l_x_water/(o*ne_l_x_water)*(x0-o*(ne_l_x_bed+ne_l_sand))+l-l_x_water)*(1.-msk)
116
117 x1=x[1]
118 x1_new = (l-d0_sand-l_y_water)/(o*ne_l_y_bed)*x1
119 msk=whereNonPositive(x1-(o*ne_l_y_bed))
120 x1_new = x1_new * msk + (d0_sand/(o*ne_l_sand)*(x1-o*ne_l_y_bed)+l-d0_sand-l_y_water)*(1.-msk)
121 msk=whereNonPositive(x1-o*(ne_l_y_bed+ne_l_sand))
122 x1_new = x1_new * msk + (l_y_water/(o*ne_l_y_water)*(x1-o*(ne_l_y_bed+ne_l_sand))+l-l_y_water)*(1.-msk)
123
124 x2=x[2]
125 x2_new = (h-d_sand-h_water)/(o*ne_h_bed)*x2
126 msk=whereNonPositive(x2-(o*ne_h_bed+1))
127 x2_new = x2_new * msk + (d_sand/(o*ne_h_sand)*(x2-(o*ne_h_bed+1))+h-d_sand-h_water)*(1.-msk)
128 msk=whereNonPositive(x2-(o*(ne_h_bed+ne_h_sand)+1))
129 x2_new = x2_new * msk + (h_water/(o*ne_h_water)*(x2-(o*(ne_h_bed+ne_h_sand)+1))+h-h_water)*(1.-msk)
130
131 dom.setX(x0_new*[1,0,0]+x1_new*[0,1,0]+x2_new*[0,0,1])
132
133 fs=Function(dom)
134 x=Function(dom).getX()
135 fs.setTags(sand_tag,wherePositive(x[0]-(l-l_x_water-d0_sand)) \
136 *wherePositive(x[1]-(l-l_y_water-d0_sand)) \
137 *wherePositive(x[2]-(h-h_water-d_sand)))
138 fs.setTags(water_tag,wherePositive(x[0]-(l-l_x_water)) \
139 *wherePositive(x[1]-(l-l_y_water)) \
140 *wherePositive(x[2]-(h-h_water)))
141 return dom
142
143 def getMaterialProperties(dom):
144 rho =Scalar(rho_bedrock,Function(dom))
145 rho.setTaggedValue(sand_tag,rho_sand)
146 rho.setTaggedValue(water_tag,rho_water)
147
148 lame_mu =Scalar(mu_bedrock,Function(dom))
149 lame_mu.setTaggedValue(sand_tag,mu_sand)
150 lame_mu.setTaggedValue(water_tag,mu_water)
151
152 lame_lambda=Scalar(lambda_bedrock,Function(dom))
153 lame_lambda.setTaggedValue(sand_tag,lambda_sand)
154 lame_lambda.setTaggedValue(water_tag,lambda_water)
155
156 return rho,lame_mu,lame_lambda
157
158
159 def wavePropagation(dom,rho,lame_mu,lame_lambda):
160 x=Function(dom).getX()
161 # ... open new PDE ...
162 mypde=LinearPDE(dom)
163 mypde.setSolverMethod(LinearPDE.LUMPING)
164 k=kronecker(Function(dom))
165 mypde.setValue(D=k*rho)
166
167 v_p=sqrt((2*lame_mu+lame_lambda)/rho)
168 if output: print "v_p=",v_p
169 v_s=sqrt(lame_mu/rho)
170 if output: print "v_s=",v_s
171 dt=(1./5.)*inf(dom.getSize()/v_p)
172 if output: print "time step size = ",dt
173 # ... set initial values ....
174 n=0
175 t=0
176 # initial value of displacement at point source is constant (U0=0.01)
177 # for first two time steps
178 u =Vector(0.,Solution(dom))
179 u_last=Vector(0.,Solution(dom))
180
181 starttime = time.clock()
182 while t<t_end and n<n_end:
183 if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
184 # ... get current stress ....
185 eps=symmetric(grad(u))
186 stress=lame_lambda*trace(eps)*k+2*lame_mu*eps
187 # ... force due to event:
188 F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
189 if output: print Lsup(F)
190 # ... get new acceleration ....
191 mypde.setValue(X=-stress,Y=F)
192 a=mypde.getSolution()
193 # ... get new displacement ...
194 u_new=2*u-u_last+dt**2*a
195 # ... shift displacements ....
196 u_last,u=u,u_new
197 # ... save current acceleration in units of gravity and displacements
198 if output:
199 if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))
200
201 t+=dt
202 n+=1
203
204 endtime = time.clock()
205 totaltime = endtime-starttime
206 global netotal
207 print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
208 if __name__ =="__main__":
209 dom=getDomain()
210 rho,lame_mu,lame_lambda=getMaterialProperties(dom)
211 wavePropagation(dom,rho,lame_mu,lame_lambda)

  ViewVC Help
Powered by ViewVC 1.1.26