/[escript]/trunk/doc/examples/inversion/synthetic_HTI.py
ViewVC logotype

Diff of /trunk/doc/examples/inversion/synthetic_HTI.py

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

revision 4575 by gross, Mon Dec 9 04:29:44 2013 UTC revision 4576 by sshaw, Mon Dec 9 23:35:30 2013 UTC
# Line 79  width_y=width_x Line 79  width_y=width_x
79  #  #
80  # create array  # create array
81  #  #
82  receiver_line=[2*absorption_zone + i * (rangeRcv/(numRcvPerLine-1)) for i in xrange(numRcvPerLine) ]  receiver_line=[2*absorption_zone + i * (rangeRcv/(numRcvPerLine-1)) for i in range(numRcvPerLine) ]
83  #  #
84  #   set source location with tag "source""  #   set source location with tag "source""
85  #  #
# Line 104  rcv_locations=[] Line 104  rcv_locations=[]
104  rg=[]  rg=[]
105  mid_point=receiver_line[len(receiver_line)/2]  mid_point=receiver_line[len(receiver_line)/2]
106    
107  for ix in xrange(len(receiver_line)):  for ix in range(len(receiver_line)):
108          if DIM == 2:          if DIM == 2:
109              rcv_locations.append((receiver_line[ix], mid_point))              rcv_locations.append((receiver_line[ix], mid_point))
110              rg.append( ( receiver_line[ix], mid_point) )              rg.append( ( receiver_line[ix], mid_point) )
# Line 112  for ix in xrange(len(receiver_line)): Line 112  for ix in xrange(len(receiver_line)):
112             rcv_locations.append((receiver_line[ix], mid_point, depth))             rcv_locations.append((receiver_line[ix], mid_point, depth))
113             rg.append( ( receiver_line[ix], mid_point) )             rg.append( ( receiver_line[ix], mid_point) )
114  # North-south line of receiver  # North-south line of receiver
115  for iy in xrange(len(receiver_line)):  for iy in range(len(receiver_line)):
116          if DIM == 2:          if DIM == 2:
117              rcv_locations.append((mid_point, receiver_line[iy]))              rcv_locations.append((mid_point, receiver_line[iy]))
118              rg.append( ( mid_point, receiver_line[iy]) )              rg.append( ( mid_point, receiver_line[iy]) )
# Line 124  for iy in xrange(len(receiver_line)): Line 124  for iy in xrange(len(receiver_line)):
124  #  #
125  if DIM == 2:  if DIM == 2:
126     domain=Rectangle(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),l0=width_x,l1=width_y,     domain=Rectangle(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),l0=width_x,l1=width_y,
127          diracPoints=src_locations, diracTags=src_tags)          diracPoints=src_locations, diracTags=src_tags)
128  else:  else:
129     domain=Brick(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,     domain=Brick(ceil(ne_z*width_x/depth),ceil(ne_z*width_y/depth),ne_z,l0=width_x,l1=width_y,l2=depth,
130          diracPoints=src_locations, diracTags=src_tags)          diracPoints=src_locations, diracTags=src_tags)
131  wl=Ricker(frq)  wl=Ricker(frq)
132    
133  #======================================================================  #======================================================================
# Line 139  if DIM == 3: Line 139  if DIM == 3:
139     vareps=0     vareps=0
140     gamma=0     gamma=0
141     rho=0     rho=0
142     for l in xrange(len(layers)):     for l in range(len(layers)):
143         m=wherePositive(z-z_bottom)*whereNonNegative(z-(z_bottom+layers[l]))         m=wherePositive(z-z_bottom)*whereNonNegative(z-(z_bottom+layers[l]))
144         v_p=v_p*(1-m)+v_Ps[l]*m         v_p=v_p*(1-m)+v_Ps[l]*m
145         v_s=v_s*(1-m)+v_Ss[l]*m         v_s=v_s*(1-m)+v_Ss[l]*m
# Line 163  mkDir('tmp') Line 163  mkDir('tmp')
163  n=0  n=0
164  k=0  k=0
165  while t < t_end:  while t < t_end:
166      t,u = sw.update(t+sampling_interval)      t,u = sw.update(t+sampling_interval)
167      tracer_x.addRecord(loc(u[0]))      tracer_x.addRecord(loc(u[0]))
168      tracer_y.addRecord(loc(u[1]))      tracer_y.addRecord(loc(u[1]))
169      if DIM==3:      if DIM==3:
170             tracer_z.addRecord(loc(u[2]))          tracer_z.addRecord(loc(u[2]))
171      print t, loc(u[0])[len(rg)/2-4:len(rg)/2+1], wl.getValue(t)      print t, loc(u[0])[len(rg)/2-4:len(rg)/2+1], wl.getValue(t)
172      #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)      #if n%5 == 0 : saveSilo("tmp/u_%d.silo"%(n/5,), u=u)
173      if t>0.3 and t< 0.5:      if t>0.3 and t< 0.5:
174         saveSilo("tmp/u_%d.silo"%(k,), u=u)          saveSilo("tmp/u_%d.silo"%(k,), u=u)
175         k+=1          k+=1
176          n+=1          n+=1
177  tracer_x.write('line_x.sgy')  tracer_x.write('line_x.sgy')
178  tracer_y.write('line_y.sgy')  tracer_y.write('line_y.sgy')

Legend:
Removed from v.4575  
changed lines
  Added in v.4576

  ViewVC Help
Powered by ViewVC 1.1.26