/[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 850 by gross, Sun Sep 17 23:27:00 2006 UTC revision 862 by gross, Tue Oct 3 06:01:30 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=4000.  # number of elements per m  resolution=5000.  # number of elements per m
30  l=100000.          # 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  o=1               # element order
# Line 54  lambda_sand=1.5e10 Line 54  lambda_sand=1.5e10
54    
55    
56  # location and geometrical size of event:  # location and geometrical size of event:
57  xc=[30000.,40000.,10000.]  xc=[l*0.5,l*0.5,h*0.3]
58  src_radius = 0.1*h  src_radius  = 2*resolution
59  # direction of event:  # direction of event:
60  event=numarray.array([1.,0.,0.])*1.e6  event=numarray.array([0.,0.,1.])*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)
# Line 68  if output: Line 68  if output:
68     print "length of event = ",tc_length     print "length of event = ",tc_length
69     print "direction = ",event     print "direction = ",event
70    
71  t_end=20.  t_end=30.
72    
73  def getDomain():  def getDomain():
74      """      """
# Line 116  def getDomain(): Line 116  def getDomain():
116    
117      x1=x[1]      x1=x[1]
118      x1_new = (l-d0_sand-l_y_water)/(o*ne_l_y_bed)*x1      x1_new = (l-d0_sand-l_y_water)/(o*ne_l_y_bed)*x1
119      msk=whereNonPositive(x1-(o*ne_l_y_bed))      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)      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))      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)      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]      x2=x[2]
125      x2_new = (h-d_sand-h_water)/(o*ne_h_bed)*x2      x2_new = (h-d_sand-h_water)/(o*ne_h_bed)*x2
126      msk=whereNonPositive(x2-(o*ne_h_bed+1))      msk=whereNonPositive(x2-o*ne_h_bed)
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)      x2_new = x2_new * msk + (d_sand/(o*ne_h_sand)*(x2-o*ne_h_bed)+h-d_sand-h_water)*(1.-msk)
128      msk=whereNonPositive(x2-(o*(ne_h_bed+ne_h_sand)+1))      msk=whereNonPositive(x2-o*(ne_h_bed+ne_h_sand))
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)      x2_new = x2_new * msk + (h_water/(o*ne_h_water)*(x2-o*(ne_h_bed+ne_h_sand))+h-h_water)*(1.-msk)
130    
131      dom.setX(x0_new*[1,0,0]+x1_new*[0,1,0]+x2_new*[0,0,1])      dom.setX(x0_new*[1,0,0]+x1_new*[0,1,0]+x2_new*[0,0,1])
132    
# Line 173  def wavePropagation(dom,rho,lame_mu,lame Line 173  def wavePropagation(dom,rho,lame_mu,lame
173     # ... set initial values ....     # ... set initial values ....
174     n=0     n=0
175     t=0     t=0
176       t_write=0
177     # initial value of displacement at point source is constant (U0=0.01)     # initial value of displacement at point source is constant (U0=0.01)
178     # for first two time steps     # for first two time steps
179     u     =Vector(0.,Solution(dom))     u     =Vector(0.,Solution(dom))
# Line 196  def wavePropagation(dom,rho,lame_mu,lame Line 197  def wavePropagation(dom,rho,lame_mu,lame
197       u_last,u=u,u_new       u_last,u=u,u_new
198       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
199       if output:       if output:
200            if n%10==0: saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))            if t>=t_write:
201                 saveVTK("disp.%i.vtu"%(n/10),displacement=u, amplitude=length(u))
202                 t_write+=0.5
203       t+=dt       t+=dt
204       n+=1       n+=1
205    
Legend:
Removed from v.850  
changed lines
  Added in v.862
  ViewVC Help
Powered by ViewVC 1.1.26