/[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 865 by gross, Fri Oct 6 10:02:18 2006 UTC revision 866 by gross, Mon Oct 9 03:50:20 2006 UTC
# Line 33  l=100000.           # width and length m Line 33  l=100000.           # width and length m
33  h=30000.            # height in m        (without obsorber)  h=30000.            # height in m        (without obsorber)
34  d_absorber=l*0.10   # thickness of absorbing layer  d_absorber=l*0.10   # thickness of absorbing layer
35    
36  l_sand=5000.          # thickness of sand region on surface  l_sand=20000.          # thickness of sand region on surface
37  h_sand=2000.           # thickness of sand layer under the water  h_sand=5000.           # thickness of sand layer under the water
38    
39  l_x_water=10000.       # length of water in x  l_x_water=10000.       # length of water in x
40  l_y_water=10000.       # length of water in y  l_y_water=10000.       # length of water in y
# Line 85  eta_tab[sand]=eta_tab[absorber]/40. Line 85  eta_tab[sand]=eta_tab[absorber]/40.
85  eta_tab[water]=eta_tab[absorber]/40.  eta_tab[water]=eta_tab[absorber]/40.
86  eta_tab[bedrock]=eta_tab[absorber]/40.  eta_tab[bedrock]=eta_tab[absorber]/40.
87    
88    
89    # material properties:
90    bedrock=0
91    absorber=1
92    water=2
93    sand=3
94    
95    rho={}
96    rho[bedrock]=8e3
97    rho[absorber]=rho[bedrock]
98    rho[water]=1e3
99    rho[sand]=5e3
100    
101    mu={}
102    mu[bedrock]=1.7e11
103    mu[absorber]=mu[bedrock]
104    mu[water]=0.
105    mu[sand]=1.5e10
106    
107    lmbd={}
108    lmbd[bedrock]=1.7e11
109    lmbd_absorber=lmbd[bedrock]
110    lmbd[water]=1.e9
111    lmbd[sand]=1.5e10
112    
113    eta={}
114    eta[absorber]=-log(0.05)*sqrt(rho[absorber]*(lmbd_absorber+2*mu[absorber]))/d_absorber
115    eta[sand]=eta[absorber]/40.
116    eta[water]=eta[absorber]/40.
117    eta[bedrock]=eta[absorber]/40.
118    
119  if output:  if output:
120     print "event location = ",xc     print "event location = ",xc
121     print "radius of event = ",src_radius     print "radius of event = ",src_radius
# Line 108  def getDomain(): Line 139  def getDomain():
139      for tag in rho_tab.keys():      for tag in rho_tab.keys():
140         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])
141      v_p_ref=min(v_p.values())      v_p_ref=min(v_p.values())
   
142      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)
143    
144      sections={}      sections={}
# Line 242  def getMaterialProperties(dom): Line 272  def getMaterialProperties(dom):
272        tags.setTaggedValue(tag,tag)        tags.setTaggedValue(tag,tag)
273     return rho,mu,lmbd,eta     return rho,mu,lmbd,eta
274    
   
275  def wavePropagation(dom,rho,mu,lmbd,eta):  def wavePropagation(dom,rho,mu,lmbd,eta):
276     x=Function(dom).getX()     x=Function(dom).getX()
277     # ... open new PDE ...     # ... open new PDE ...
# Line 260  def wavePropagation(dom,rho,mu,lmbd,eta) Line 289  def wavePropagation(dom,rho,mu,lmbd,eta)
289     n_write=0     n_write=0
290     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
291     # for first two time steps     # for first two time steps
292     u     =Vector(0.,Solution(dom))     u=Vector(0.,Solution(dom))
293     u_last=Vector(0.,Solution(dom))     v=Vector(0.,Solution(dom))
294       a=Vector(0.,Solution(dom))
295       a2=Vector(0.,Solution(dom))
296     v=Vector(0.,Solution(dom))     v=Vector(0.,Solution(dom))
297    
298     starttime = time.clock()     starttime = time.clock()
299     while t<t_end and n<n_end:     while t<t_end and n<n_end:
300       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),
301         # prediction:
302         u_pr=u+dt*v+(dt**2/2)*a+(dt**3/6)*a2
303         v_pr=v+dt*a+(dt**2/2)*a2
304         a_pr=a+dt*a2
305       # ... get current stress ....       # ... get current stress ....
306       eps=symmetric(grad(u))       eps=symmetric(grad(u_pr))
307       stress=lmbd*trace(eps)*k+2*mu*eps       stress=lmbd*trace(eps)*k+2*mu*eps
308       # ... force due to event:       # ... force due to event:
309       F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event       if abs(t-tc)<5*tc_length:
310       if output: print Lsup(F)          F=exp(-((t-tc)/tc_length)**2)*exp(-(length(x-xc)/src_radius)**2)*event
311            if output: print Lsup(F)
312         else:
313            if output: print 0.
314       # ... get new acceleration ....       # ... get new acceleration ....
315       mypde.setValue(X=-stress,Y=F-eta*v)       mypde.setValue(X=-stress,Y=F-eta*v_pr)
316       a=mypde.getSolution()       a=mypde.getSolution()
317       # ... get new displacement ...       # ... get new displacement ...
318       u_new=2*u-u_last+dt**2*a       da=a-a_pr
319       # ... shift displacements ....       u=u_pr+(dt**2/12.)*da
320       v=(u_new-u)/dt       v=v_pr+(5*dt/12.)*da
321       u_last,u=u,u_new       a2+=da/dt
322       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
323       if output:       if output:
324            if t>=t_write:            if t>=t_write:
Legend:
Removed from v.865  
changed lines
  Added in v.866
  ViewVC Help
Powered by ViewVC 1.1.26