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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 864 by gross, Tue Oct 3 06:01:30 2006 UTC revision 865 by gross, Fri Oct 6 10:02:18 2006 UTC
# Line 26  import time Line 26  import time
26  output=True  output=True
27  n_end=10000  n_end=10000
28    
29  resolution=5000.  # number of elements per m  resolution=1000.  # number of elements per m in the finest region
30  l=100000.          # width and length m  o=1              # element order
 h=30000.          # height in m  
 o=1               # element order  
   
 rho_bedrock=8e3  
 mu_bedrock=1.7e11  
 lambda_bedrock=1.7e11  
   
 l_x_water=10000    # length of water in x  
 l_y_water=10000    # length of water in y  
 h_water=max(2000,resolution)       # depth of water region  
   
 water_tag=2  
 rho_water=1e3  
 mu_water=0.  
 lambda_water=1.e9  
   
 d0_sand=10000      # thickness of sand region on surface  
 d_sand=max(2000,resolution)        # thickness of sand layer under the water  
   
 sand_tag=3  
 rho_sand=5e3  
 mu_sand=1.5e10  
 lambda_sand=1.5e10  
31    
32    l=100000.           # width and length m (without obsorber)
33    h=30000.            # height in m        (without obsorber)
34    d_absorber=l*0.10   # thickness of absorbing layer
35    
36  # location and geometrical size of event:  l_sand=5000.          # thickness of sand region on surface
37  xc=[l*0.5,l*0.5,h*0.3]  h_sand=2000.           # thickness of sand layer under the water
38    
39    l_x_water=10000.       # length of water in x
40    l_y_water=10000.       # length of water in y
41    h_water=2000.          # depth of water region
42    
43    x_sand=l/2-l_x_water/2-l_sand # x coordinate of location of sand region (without obsorber)
44    y_sand=l/2-l_y_water/2-l_sand # y coordinate of location of sand region (without obsorber)
45    
46    
47    # origin
48    origin={"x": -d_absorber, "y" : -d_absorber , "z" : -h-d_absorber }
49    # location and geometrical size of event reltive to origin:
50    xc=[l*0.2,l*0.3,-h*0.7]
51  src_radius  = 2*resolution  src_radius  = 2*resolution
52  # direction of event:  # direction of event:
53  event=numarray.array([0.,0.,1.])*1.e6  event=numarray.array([0.,0.,1.])*1.e6
54  # time and length of the event  # time and length of the event
55  tc_length=2.  tc=2.
56  tc=sqrt(5.*tc_length)  tc_length=0.5
57    
58    # material properties:
59    bedrock=0
60    absorber=1
61    water=2
62    sand=3
63    
64    rho_tab={}
65    rho_tab[bedrock]=8e3
66    rho_tab[absorber]=rho_tab[bedrock]
67    rho_tab[water]=1e3
68    rho_tab[sand]=5e3
69    
70    mu_tab={}
71    mu_tab[bedrock]=1.7e11
72    mu_tab[absorber]=mu_tab[bedrock]
73    mu_tab[water]=0.
74    mu_tab[sand]=1.5e10
75    
76    lmbd_tab={}
77    lmbd_tab[bedrock]=1.7e11
78    lmbd_tab[absorber]=lmbd_tab[bedrock]
79    lmbd_tab[water]=1.e9
80    lmbd_tab[sand]=1.5e10
81    
82    eta_tab={}
83    eta_tab[absorber]=-log(0.05)*sqrt(rho_tab[absorber]*(lmbd_tab[absorber]+2*mu_tab[absorber]))/d_absorber
84    eta_tab[sand]=eta_tab[absorber]/40.
85    eta_tab[water]=eta_tab[absorber]/40.
86    eta_tab[bedrock]=eta_tab[absorber]/40.
87    
88  if output:  if output:
89     print "event location = ",xc     print "event location = ",xc
90     print "radius of event = ",src_radius     print "radius of event = ",src_radius
# Line 69  if output: Line 93  if output:
93     print "direction = ",event     print "direction = ",event
94    
95  t_end=30.  t_end=30.
96    dt_write=0.1
97    
98    
99  def getDomain():  def getDomain():
100      """      """
# Line 77  def getDomain(): Line 103  def getDomain():
103                
104      """      """
105      global netotal      global netotal
     v_p_bedrock=sqrt((2*mu_bedrock+lambda_bedrock)/rho_bedrock)  
     v_p_sand=sqrt((2*mu_sand+lambda_sand)/rho_sand)  
     v_p_water=sqrt((2*mu_water+lambda_water)/rho_water)  
   
     print v_p_bedrock,v_p_sand,v_p_water  
106            
107      ne_l_x_bed=int((l-d0_sand-l_x_water)/resolution+0.5)      v_p={}
108      ne_l_y_bed=int((l-d0_sand-l_y_water)/resolution+0.5)      for tag in rho_tab.keys():
109      ne_h_bed=int((h-d_sand-h_water)/resolution+0.5)         v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])
110        v_p_ref=min(v_p.values())
111      ne_l_sand=int(d0_sand*v_p_bedrock/v_p_sand/resolution+0.5)  
112      ne_h_sand=int(d_sand*v_p_bedrock/v_p_sand/resolution+0.5)      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)
113    
114      ne_l_x_water=int(l_x_water*v_p_bedrock/v_p_water/resolution+0.5)      sections={}
115      ne_l_y_water=int(l_y_water*v_p_bedrock/v_p_water/resolution+0.5)      sections["x"]=[d_absorber, x_sand, l_sand, l_x_water, l_sand, l-x_sand-2*l_sand-l_x_water, d_absorber]
116      ne_h_water=int(h_water*v_p_bedrock/v_p_water/resolution+0.5)      sections["y"]=[d_absorber, y_sand, l_sand, l_y_water, l_sand, l-y_sand-2*l_sand-l_y_water, d_absorber]
117        sections["z"]=[d_absorber,h-h_water-h_sand,h_sand,h_water]
118      ne_l_x=ne_l_x_bed+ne_l_sand+ne_l_x_water      if output:
119      ne_l_y=ne_l_y_bed+ne_l_sand+ne_l_y_water        print "sections x = ",sections["x"]
120      ne_h=ne_h_bed+ne_h_sand+ne_h_water        print "sections y = ",sections["y"]
121          print "sections z = ",sections["z"]
122      print ne_l_x,ne_l_x_bed,ne_l_sand,ne_l_x_water  
123      print ne_l_y,ne_l_y_bed,ne_l_sand,ne_l_y_water      mats= [
124      print ne_h,ne_h_bed,ne_h_sand,ne_h_water              [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
125                  [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
126      netotal=ne_l_x*ne_l_y*ne_h                [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
127      if output: print "grid : %s x %s x %s (%s elements)"%(ne_l_x,ne_l_y,ne_h,netotal)                [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
128      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)                [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
129      x=dom.getX()                [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
130                  [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
131    
132                [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
133                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
134                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
135                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
136                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
137                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
138                  [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
139    
140                [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
141                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
142                  [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
143                  [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
144                  [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
145                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
146                  [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ],
147    
148                [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
149                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
150                  [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
151                  [absorber, bedrock , sand    , water   , sand    , bedrock , absorber],
152                  [absorber, bedrock , sand    , sand    , sand    , bedrock , absorber],
153                  [absorber, bedrock , bedrock , bedrock , bedrock , bedrock , absorber],
154                  [absorber, absorber, absorber, absorber, absorber, absorber, absorber] ] ]
155            
156      x0=x[0]      num_elem={}
157      x0_new = (l-d0_sand-l_x_water)/(o*ne_l_x_bed)*x0      for d in sections:
158      msk=whereNonPositive(x0-o*ne_l_x_bed)         num_elem[d]=[]
159      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)         for i in range(len(sections[d])):
160      msk=whereNonPositive(x0-o*(ne_l_x_bed+ne_l_sand))             if d=="x":
161      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)                v_p_min=v_p[mats[0][0][i]]
162                  for q in range(len(sections["y"])):
163      x1=x[1]                   for r in range(len(sections["z"])):
164      x1_new = (l-d0_sand-l_y_water)/(o*ne_l_y_bed)*x1                      v_p_min=min(v_p[mats[r][q][i]],v_p_min)
165      msk=whereNonPositive(x1-o*ne_l_y_bed)             elif d=="y":
166      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)                v_p_min=v_p[mats[0][i][0]]
167      msk=whereNonPositive(x1-o*(ne_l_y_bed+ne_l_sand))                for q in range(len(sections["x"])):
168      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)                   for r in range(len(sections["z"])):
169                        v_p_min=min(v_p[mats[r][i][q]],v_p_min)
170      x2=x[2]             elif d=="z":
171      x2_new = (h-d_sand-h_water)/(o*ne_h_bed)*x2                v_p_min=v_p[mats[i][0][0]]
172      msk=whereNonPositive(x2-o*ne_h_bed)                for q in range(len(sections["x"])):
173      x2_new = x2_new * msk + (d_sand/(o*ne_h_sand)*(x2-o*ne_h_bed)+h-d_sand-h_water)*(1.-msk)                   for r in range(len(sections["y"])):
174      msk=whereNonPositive(x2-o*(ne_h_bed+ne_h_sand))                      v_p_min=min(v_p[mats[i][r][q]],v_p_min)
175      x2_new = x2_new * msk + (h_water/(o*ne_h_water)*(x2-o*(ne_h_bed+ne_h_sand))+h-h_water)*(1.-msk)             num_elem[d].append(max(1,int(sections[d][i] * v_p_ref/v_p_min /resolution+0.5)))
176          
177      dom.setX(x0_new*[1,0,0]+x1_new*[0,1,0]+x2_new*[0,0,1])      ne_x=sum(num_elem["x"])
178        ne_y=sum(num_elem["y"])
179        ne_z=sum(num_elem["z"])
180        netotal=ne_x*ne_y*ne_z
181        if output: print "grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal)
182        dom=Brick(ne_x,ne_y,ne_z,l0=o*ne_x,l1=o*ne_y,l2=o*ne_z,order=o)
183        x_old=dom.getX()
184        x_new=0
185    
186        for d in sections:
187            if d=="x":
188                i=0
189                f=[1,0,0]
190            if d=="y":
191                i=1
192                f=[0,1,0]
193            if d=="z":
194                i=2
195                f=[0,0,1]
196            x=x_old[i]
197    
198            p=origin[d]
199            ne=0
200            s=0.
201    
202            for i in range(len(sections[d])-1):
203                msk=whereNonPositive(x-o*ne+0.5)
204                s=s*msk + (sections[d][i]/(o*num_elem[d][i])*(x-o*ne)+p)*(1.-msk)
205                ne+=num_elem[d][i]
206                p+=sections[d][i]
207            x_new=x_new + s * f
208        dom.setX(x_new)
209    
210      fs=Function(dom)      fs=Function(dom)
211      x=Function(dom).getX()      x=Function(dom).getX()
212      fs.setTags(sand_tag,wherePositive(x[0]-(l-l_x_water-d0_sand)) \      x0=x[0]
213                         *wherePositive(x[1]-(l-l_y_water-d0_sand)) \      x1=x[1]
214                         *wherePositive(x[2]-(h-h_water-d_sand)))      x2=x[2]
215      fs.setTags(water_tag,wherePositive(x[0]-(l-l_x_water)) \      p_z=origin["z"]
216                          *wherePositive(x[1]-(l-l_y_water)) \      for i in range(len(mats)):
217                          *wherePositive(x[2]-(h-h_water)))         f_z=wherePositive(x2-p_z)*wherePositive(x2-p_z+sections["z"][i])
218           p_y=origin["y"]
219           for j in range(len(mats[i])):
220             f_y=wherePositive(x1-p_y)*wherePositive(x1-p_z+sections["y"][j])
221             p_x=origin["x"]
222             for k in range(len(mats[i][j])):
223                 f_x=wherePositive(x0-p_x)*wherePositive(x0-p_x+sections["x"][k])
224                 fs.setTags(mats[i][j][k],f_x*f_y*f_z)
225                 p_x+=sections["x"][k]
226             p_y+=sections["y"][j]
227           p_z+=sections["z"][i]
228      return dom      return dom
229    
230  def getMaterialProperties(dom):  def getMaterialProperties(dom):
231     rho        =Scalar(rho_bedrock,Function(dom))     rho =Scalar(rho_tab[bedrock],Function(dom))
232     rho.setTaggedValue(sand_tag,rho_sand)     eta =Scalar(eta_tab[bedrock],Function(dom))
233     rho.setTaggedValue(water_tag,rho_water)     mu  =Scalar(mu_tab[bedrock],Function(dom))
234       lmbd=Scalar(lmbd_tab[bedrock],Function(dom))
235     lame_mu    =Scalar(mu_bedrock,Function(dom))     tags=Scalar(bedrock,Function(dom))
236     lame_mu.setTaggedValue(sand_tag,mu_sand)    
237     lame_mu.setTaggedValue(water_tag,mu_water)     for tag in rho_tab.keys():
238          rho.setTaggedValue(tag,rho_tab[tag])
239     lame_lambda=Scalar(lambda_bedrock,Function(dom))        eta.setTaggedValue(tag,eta_tab[tag])
240     lame_lambda.setTaggedValue(sand_tag,lambda_sand)        mu.setTaggedValue(tag,mu_tab[tag])
241     lame_lambda.setTaggedValue(water_tag,lambda_water)        lmbd.setTaggedValue(tag,lmbd_tab[tag])
242          tags.setTaggedValue(tag,tag)
243     return rho,lame_mu,lame_lambda     return rho,mu,lmbd,eta
244    
245    
246  def wavePropagation(dom,rho,lame_mu,lame_lambda):  def wavePropagation(dom,rho,mu,lmbd,eta):
247     x=Function(dom).getX()     x=Function(dom).getX()
248     # ... open new PDE ...     # ... open new PDE ...
249     mypde=LinearPDE(dom)     mypde=LinearPDE(dom)
# Line 164  def wavePropagation(dom,rho,lame_mu,lame Line 251  def wavePropagation(dom,rho,lame_mu,lame
251     k=kronecker(Function(dom))     k=kronecker(Function(dom))
252     mypde.setValue(D=k*rho)     mypde.setValue(D=k*rho)
253    
254     v_p=sqrt((2*lame_mu+lame_lambda)/rho)     dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))
    if output: print "v_p=",v_p  
    v_s=sqrt(lame_mu/rho)  
    if output: print "v_s=",v_s  
    dt=(1./5.)*inf(dom.getSize()/v_p)  
255     if output: print "time step size = ",dt     if output: print "time step size = ",dt
256     # ... set initial values ....     # ... set initial values ....
257     n=0     n=0
258     t=0     t=0
259     t_write=0     t_write=0.
260       n_write=0
261     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
262     # for first two time steps     # for first two time steps
263     u     =Vector(0.,Solution(dom))     u     =Vector(0.,Solution(dom))
264     u_last=Vector(0.,Solution(dom))     u_last=Vector(0.,Solution(dom))
265       v=Vector(0.,Solution(dom))
266    
267     starttime = time.clock()     starttime = time.clock()
268     while t<t_end and n<n_end:     while t<t_end and n<n_end:
269       if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),       if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
270       # ... get current stress ....       # ... get current stress ....
271       eps=symmetric(grad(u))       eps=symmetric(grad(u))
272       stress=lame_lambda*trace(eps)*k+2*lame_mu*eps       stress=lmbd*trace(eps)*k+2*mu*eps
273       # ... force due to event:       # ... force due to event:
274       F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event       F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
275       if output: print Lsup(F)       if output: print Lsup(F)
276       # ... get new acceleration ....       # ... get new acceleration ....
277       mypde.setValue(X=-stress,Y=F)       mypde.setValue(X=-stress,Y=F-eta*v)
278       a=mypde.getSolution()       a=mypde.getSolution()
279       # ... get new displacement ...       # ... get new displacement ...
280       u_new=2*u-u_last+dt**2*a       u_new=2*u-u_last+dt**2*a
281       # ... shift displacements ....       # ... shift displacements ....
282         v=(u_new-u)/dt
283       u_last,u=u,u_new       u_last,u=u,u_new
284       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
285       if output:       if output:
286            if t>=t_write:            if t>=t_write:
287               saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))               saveVTK("disp.%i.vtu"%n,displacement=u, amplitude=length(u))
288               t_write+=0.5               t_write+=dt_write
289                 n_write+=1
290       t+=dt       t+=dt
291       n+=1       n+=1
292    
# Line 209  def wavePropagation(dom,rho,lame_mu,lame Line 296  def wavePropagation(dom,rho,lame_mu,lame
296     print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)     print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
297  if __name__ =="__main__":  if __name__ =="__main__":
298     dom=getDomain()     dom=getDomain()
299     rho,lame_mu,lame_lambda=getMaterialProperties(dom)     rho,mu,lmbd,eta=getMaterialProperties(dom)
300     wavePropagation(dom,rho,lame_mu,lame_lambda)     wavePropagation(dom,rho,mu,lmbd,eta)

Legend:
Removed from v.864  
changed lines
  Added in v.865

  ViewVC Help
Powered by ViewVC 1.1.26