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

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

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

revision 3788 by caltinay, Tue Mar 15 04:23:54 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 133  eta[water]=eta[absorber]/40. Line 133  eta[water]=eta[absorber]/40.
133  eta[bedrock]=eta[absorber]/40.  eta[bedrock]=eta[absorber]/40.
134    
135  if output:  if output:
136     print "event location = ",xc     print("event location = ",xc)
137     print "radius of event = ",src_radius     print("radius of event = ",src_radius)
138     print "time of event = ",tc     print("time of event = ",tc)
139     print "length of event = ",tc_length     print("length of event = ",tc_length)
140     print "direction = ",event     print("direction = ",event)
141    
142  t_end=30.  t_end=30.
143  dt_write=0.1  dt_write=0.1
# Line 152  def getDomain(): Line 152  def getDomain():
152      global netotal      global netotal
153            
154      v_p={}      v_p={}
155      for tag in rho_tab.keys():      for tag in list(rho_tab.keys()):
156         v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])         v_p[tag]=sqrt((2*mu_tab[tag]+lmbd_tab[tag])/rho_tab[tag])
157      v_p_ref=min(v_p.values())      v_p_ref=min(v_p.values())
158      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)      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    
160      sections={}      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]      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]      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]      sections["z"]=[d_absorber,h-h_water-h_sand,h_sand,h_water]
164      if output:      if output:
165        print "sections x = ",sections["x"]        print("sections x = ",sections["x"])
166        print "sections y = ",sections["y"]        print("sections y = ",sections["y"])
167        print "sections z = ",sections["z"]        print("sections z = ",sections["z"])
168    
169      mats= [      mats= [
170              [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],              [ [absorber, absorber, absorber, absorber, absorber, absorber, absorber],
# Line 224  def getDomain(): Line 224  def getDomain():
224      ne_y=sum(num_elem["y"])      ne_y=sum(num_elem["y"])
225      ne_z=sum(num_elem["z"])      ne_z=sum(num_elem["z"])
226      netotal=ne_x*ne_y*ne_z      netotal=ne_x*ne_y*ne_z
227      if output: print "grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal)      if output: print("grid : %s x %s x %s (%s elements)"%(ne_x,ne_y,ne_z,netotal))
228      dom=Brick(ne_x,ne_y,ne_z,l0=o*ne_x,l1=o*ne_y,l2=o*ne_z,order=o)      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()      x_old=dom.getX()
230      x_new=0      x_new=0
# Line 280  def getMaterialProperties(dom): Line 280  def getMaterialProperties(dom):
280     lmbd=Scalar(lmbd_tab[bedrock],Function(dom))     lmbd=Scalar(lmbd_tab[bedrock],Function(dom))
281     tags=Scalar(bedrock,Function(dom))     tags=Scalar(bedrock,Function(dom))
282        
283     for tag in rho_tab.keys():     for tag in list(rho_tab.keys()):
284        rho.setTaggedValue(tag,rho_tab[tag])        rho.setTaggedValue(tag,rho_tab[tag])
285        eta.setTaggedValue(tag,eta_tab[tag])        eta.setTaggedValue(tag,eta_tab[tag])
286        mu.setTaggedValue(tag,mu_tab[tag])        mu.setTaggedValue(tag,mu_tab[tag])
# Line 297  def wavePropagation(dom,rho,mu,lmbd,eta) Line 297  def wavePropagation(dom,rho,mu,lmbd,eta)
297     mypde.setValue(D=k*rho)     mypde.setValue(D=k*rho)
298    
299     dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))     dt=(1./5.)*inf(dom.getSize()/sqrt((2*mu+lmbd)/rho))
300     if output: print "time step size = ",dt     if output: print("time step size = ",dt)
301     # ... set initial values ....     # ... set initial values ....
302     n=0     n=0
303     t=0     t=0
# Line 315  def wavePropagation(dom,rho,mu,lmbd,eta) Line 315  def wavePropagation(dom,rho,mu,lmbd,eta)
315    
316     starttime = time.clock()     starttime = time.clock()
317     while t<t_end and n<n_end:     while t<t_end and n<n_end:
318       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), end=' ')
319       # prediction:       # prediction:
320       u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2       u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2
321       v_pr=v+dt*a+(dt**2/2)*a2       v_pr=v+dt*a+(dt**2/2)*a2
# Line 326  def wavePropagation(dom,rho,mu,lmbd,eta) Line 326  def wavePropagation(dom,rho,mu,lmbd,eta)
326       # ... force due to event:       # ... force due to event:
327       if abs(t-tc)<5*tc_length:       if abs(t-tc)<5*tc_length:
328          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
329          if output: print Lsup(F)          if output: print(Lsup(F))
330       else:       else:
331          if output: print 0.          if output: print(0.)
332       # ... get new acceleration ....       # ... get new acceleration ....
333       mypde.setValue(X=-stress,Y=F-eta*v_pr)       mypde.setValue(X=-stress,Y=F-eta*v_pr)
334       a=mypde.getSolution()       a=mypde.getSolution()
# Line 349  def wavePropagation(dom,rho,mu,lmbd,eta) Line 349  def wavePropagation(dom,rho,mu,lmbd,eta)
349     endtime = time.clock()     endtime = time.clock()
350     totaltime = endtime-starttime     totaltime = endtime-starttime
351     global netotal     global netotal
352     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))
353  if __name__ =="__main__":  if __name__ =="__main__":
354     dom=getDomain()     dom=getDomain()
355     rho,mu,lmbd,eta=getMaterialProperties(dom)     rho,mu,lmbd,eta=getMaterialProperties(dom)

Legend:
Removed from v.3788  
changed lines
  Added in v.3789

  ViewVC Help
Powered by ViewVC 1.1.26