/[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 842 by gross, Thu Sep 7 05:40:22 2006 UTC revision 850 by gross, Sun Sep 17 23:27:00 2006 UTC
# Line 21  __date__="$Date$" Line 21  __date__="$Date$"
21  from esys.escript import *  from esys.escript import *
22  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
23  from esys.finley import Brick  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  resolution=4000.  # number of elements per m
30  l=80000.          # width and length m  l=100000.          # width and length m
31  h=30000.          # height in m  h=30000.          # height in m
32    o=1               # element order
33    
34  rho_bedrock=8e3  rho_bedrock=8e3
35  mu_bedrock=1.7e11  mu_bedrock=1.7e11
# Line 56  event=numarray.array([1.,0.,0.])*1.e6 Line 61  event=numarray.array([1.,0.,0.])*1.e6
61  # time and length of the event  # time and length of the event
62  tc_length=2.  tc_length=2.
63  tc=sqrt(5.*tc_length)  tc=sqrt(5.*tc_length)
64  print "event location = ",xc  if output:
65  print "radius of event = ",src_radius     print "event location = ",xc
66  print "time of event = ",tc     print "radius of event = ",src_radius
67  print "length of event = ",tc_length     print "time of event = ",tc
68  print "direction = ",event     print "length of event = ",tc_length
69       print "direction = ",event
70    
71  t_end=22.  t_end=20.
72    
73  def getDomain():  def getDomain():
74      """      """
# Line 70  def getDomain(): Line 76  def getDomain():
76    
77                
78      """      """
79      ne_l=int(l/resolution+0.5)      global netotal
80      ne_h=int(h/resolution+0.5)      v_p_bedrock=sqrt((2*mu_bedrock+lambda_bedrock)/rho_bedrock)
81      print "grid : %s x %s x %s"%(ne_l,ne_l,ne_h)      v_p_sand=sqrt((2*mu_sand+lambda_sand)/rho_sand)
82      dom=Brick(ne_l,ne_l,ne_h,l0=l,l1=l,l2=h,order=1)      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)      fs=Function(dom)
134      x=Function(dom).getX()      x=Function(dom).getX()
# Line 110  def wavePropagation(dom,rho,lame_mu,lame Line 165  def wavePropagation(dom,rho,lame_mu,lame
165     mypde.setValue(D=k*rho)     mypde.setValue(D=k*rho)
166    
167     v_p=sqrt((2*lame_mu+lame_lambda)/rho)     v_p=sqrt((2*lame_mu+lame_lambda)/rho)
168     print "v_p=",v_p     if output: print "v_p=",v_p
169     v_s=sqrt(lame_mu/rho)     v_s=sqrt(lame_mu/rho)
170     print "v_s=",v_s     if output: print "v_s=",v_s
171     dt=(1./5.)*inf(dom.getSize()/v_p)     dt=(1./5.)*inf(dom.getSize()/v_p)
172     print "time step size = ",dt     if output: print "time step size = ",dt
173     # ... set initial values ....     # ... set initial values ....
174     n=0     n=0
175     t=0     t=0
# Line 123  def wavePropagation(dom,rho,lame_mu,lame Line 178  def wavePropagation(dom,rho,lame_mu,lame
178     u     =Vector(0.,Solution(dom))     u     =Vector(0.,Solution(dom))
179     u_last=Vector(0.,Solution(dom))     u_last=Vector(0.,Solution(dom))
180    
181     while t<t_end:     starttime = time.clock()
182       print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),     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 ....       # ... get current stress ....
185       eps=symmetric(grad(u))       eps=symmetric(grad(u))
186       stress=lame_lambda*trace(eps)*k+2*lame_mu*eps       stress=lame_lambda*trace(eps)*k+2*lame_mu*eps
187       # ... force due to event:       # ... force due to event:
188       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
189       print Lsup(F)       if output: print Lsup(F)
190       # ... get new acceleration ....       # ... get new acceleration ....
191       mypde.setValue(X=-stress,Y=F)       mypde.setValue(X=-stress,Y=F)
192       a=mypde.getSolution()       a=mypde.getSolution()
# Line 139  def wavePropagation(dom,rho,lame_mu,lame Line 195  def wavePropagation(dom,rho,lame_mu,lame
195       # ... shift displacements ....       # ... shift displacements ....
196       u_last,u=u,u_new       u_last,u=u,u_new
197       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
198       if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))       if output:
199              if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))
200    
201       t+=dt       t+=dt
202       n+=1       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__":  if __name__ =="__main__":
209     dom=getDomain()     dom=getDomain()
210     rho,lame_mu,lame_lambda=getMaterialProperties(dom)     rho,lame_mu,lame_lambda=getMaterialProperties(dom)

Legend:
Removed from v.842  
changed lines
  Added in v.850

  ViewVC Help
Powered by ViewVC 1.1.26