/[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 843 by gross, Thu Sep 7 05:55:12 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=False
27    n_end=100
28    
29  resolution=4000.  # number of elements per m  resolution=4000.  # number of elements per m
30  l=80000.          # width and length m  l=80000.          # width and length m
# Line 56  event=numarray.array([1.,0.,0.])*1.e6 Line 60  event=numarray.array([1.,0.,0.])*1.e6
60  # time and length of the event  # time and length of the event
61  tc_length=2.  tc_length=2.
62  tc=sqrt(5.*tc_length)  tc=sqrt(5.*tc_length)
63  print "event location = ",xc  if output:
64  print "radius of event = ",src_radius     print "event location = ",xc
65  print "time of event = ",tc     print "radius of event = ",src_radius
66  print "length of event = ",tc_length     print "time of event = ",tc
67  print "direction = ",event     print "length of event = ",tc_length
68       print "direction = ",event
69    
70  t_end=22.  t_end=20.
71    
72  def getDomain():  def getDomain():
73      """      """
# Line 70  def getDomain(): Line 75  def getDomain():
75    
76                
77      """      """
78        global netotal
79      ne_l=int(l/resolution+0.5)      ne_l=int(l/resolution+0.5)
80      ne_h=int(h/resolution+0.5)      ne_h=int(h/resolution+0.5)
81      print "grid : %s x %s x %s"%(ne_l,ne_l,ne_h)      netotal=ne_l*ne_l*ne_h
82        if output: print "grid : %s x %s x %s"%(ne_l,ne_l,ne_h)
83      dom=Brick(ne_l,ne_l,ne_h,l0=l,l1=l,l2=h,order=1)      dom=Brick(ne_l,ne_l,ne_h,l0=l,l1=l,l2=h,order=1)
84    
85      fs=Function(dom)      fs=Function(dom)
# Line 110  def wavePropagation(dom,rho,lame_mu,lame Line 117  def wavePropagation(dom,rho,lame_mu,lame
117     mypde.setValue(D=k*rho)     mypde.setValue(D=k*rho)
118    
119     v_p=sqrt((2*lame_mu+lame_lambda)/rho)     v_p=sqrt((2*lame_mu+lame_lambda)/rho)
120     print "v_p=",v_p     if output: print "v_p=",v_p
121     v_s=sqrt(lame_mu/rho)     v_s=sqrt(lame_mu/rho)
122     print "v_s=",v_s     if output: print "v_s=",v_s
123     dt=(1./5.)*inf(dom.getSize()/v_p)     dt=(1./5.)*inf(dom.getSize()/v_p)
124     print "time step size = ",dt     if output: print "time step size = ",dt
125     # ... set initial values ....     # ... set initial values ....
126     n=0     n=0
127     t=0     t=0
# Line 123  def wavePropagation(dom,rho,lame_mu,lame Line 130  def wavePropagation(dom,rho,lame_mu,lame
130     u     =Vector(0.,Solution(dom))     u     =Vector(0.,Solution(dom))
131     u_last=Vector(0.,Solution(dom))     u_last=Vector(0.,Solution(dom))
132    
133     while t<t_end:     starttime = time.clock()
134       print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),     while t<t_end and n<n_end:
135         if output: print n+1,"-th time step t ",t+dt," max u and F: ",Lsup(u),
136       # ... get current stress ....       # ... get current stress ....
137       eps=symmetric(grad(u))       eps=symmetric(grad(u))
138       stress=lame_lambda*trace(eps)*k+2*lame_mu*eps       stress=lame_lambda*trace(eps)*k+2*lame_mu*eps
139       # ... force due to event:       # ... force due to event:
140       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
141       print Lsup(F)       if output: print Lsup(F)
142       # ... get new acceleration ....       # ... get new acceleration ....
143       mypde.setValue(X=-stress,Y=F)       mypde.setValue(X=-stress,Y=F)
144       a=mypde.getSolution()       a=mypde.getSolution()
# Line 139  def wavePropagation(dom,rho,lame_mu,lame Line 147  def wavePropagation(dom,rho,lame_mu,lame
147       # ... shift displacements ....       # ... shift displacements ....
148       u_last,u=u,u_new       u_last,u=u,u_new
149       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
150       if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))       if output:
151              if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))
152    
153       t+=dt       t+=dt
154       n+=1       n+=1
155    
156       endtime = time.clock()
157       totaltime = endtime-starttime
158       global netotal
159       print ">>number of elements: %s, total time: %s, per time step: %s <<"%(netotal,totaltime,totaltime/n)
160  if __name__ =="__main__":  if __name__ =="__main__":
161     dom=getDomain()     dom=getDomain()
162     rho,lame_mu,lame_lambda=getMaterialProperties(dom)     rho,lame_mu,lame_lambda=getMaterialProperties(dom)
Legend:
Removed from v.842  
changed lines
  Added in v.843
  ViewVC Help
Powered by ViewVC 1.1.26