/[escript]/trunk/downunder/py_src/seismic.py
ViewVC logotype

Diff of /trunk/downunder/py_src/seismic.py

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

revision 4555 by gross, Fri Nov 29 05:24:40 2013 UTC revision 4556 by sshaw, Tue Dec 3 05:40:53 2013 UTC
# Line 35  from esys.escript.linearPDEs import Line Line 35  from esys.escript.linearPDEs import Line
35    
36                            
37  class Wavelet(object):  class Wavelet(object):
38      """          """
39      place holder for source wavelet          place holder for source wavelet
40      """          """
41      pass          pass
42                
43  class Ricker(Wavelet):  class Ricker(Wavelet):
44      """          """
45      The Ricker Wavelet w=f(t)          The Ricker Wavelet w=f(t)
46      """          """
47      def __init__(self, f_dom=40, t_dom=None):          def __init__(self, f_dom=40, t_dom=None):
48              """                  """
49              set up Ricker wavelet wih dominant frequence f_dom and center at time t_dom. If t_dom                  set up Ricker wavelet wih dominant frequence f_dom and center at time t_dom. If t_dom
50              is not given an estimate for suitable t_dom is calculated so f(0)~0.                  is not given an estimate for suitable t_dom is calculated so f(0)~0.
51              """                  """
52          drop=18                  drop=18
53          self.__f=f_dom                  self.__f=f_dom
54          self.__f_max=sqrt(7)*f_dom                  self.__f_max=sqrt(7)*f_dom
55          self.__s=pi*self.__f                  self.__s=pi*self.__f
56          if t_dom == None:                  if t_dom == None:
57              t_dom=sqrt(drop)/self.__s                          t_dom=sqrt(drop)/self.__s
58          self.__t0=t_dom                  self.__t0=t_dom
59                            
60      def getCenter(self):          def getCenter(self):
61             """                 """
62             return value of wavelet center                 return value of wavelet center
63             """                 """
64                 return self.__t0                 return self.__t0
65                    
66          def getTimeScale(self):          def getTimeScale(self):
67                  """                  """
68                  returns the time scale which is the inverse of the largest freqence with a significant                  returns the time scale which is the inverse of the largest freqence with a significant
69                  spectral component.                  spectral component.
70                  """                  """
71          return 1/self.__f_max                  return 1/self.__f_max
72    
73      def getValue(self, t):          def getValue(self, t):
74              """                  """
75              get value of wavelet at time t                  get value of wavelet at time t
76              """                  """
77          e2=(self.__s*(t-self.__t0))**2                  e2=(self.__s*(t-self.__t0))**2
78          return (1-2*e2)*exp(-e2)                  return (1-2*e2)*exp(-e2)
79    
80      def getAcceleration(self, t):          def getAcceleration(self, t):
81              """                  """
82              get the acceleration f''(t) at time t                  get the acceleration f''(t) at time t
83              """                  """
84          e2=(self.__s*(t-self.__t0))**2                  e2=(self.__s*(t-self.__t0))**2
85                  return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)                  return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)
86                                    
87    
88  class SimpleSEGYWriter(object):  class SimpleSEGYWriter(object):
89      """          """
90      as simple writer for 2D and 3D seimic lines in particular for synthetic data          as simple writer for 2D and 3D seimic lines in particular for synthetic data
91    
92      Typical usage:          Typical usage:
93    
94             from esys.escript import unitsSI as U             from esys.escript import unitsSI as U
95         sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)             sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)
96         while n < 10:             while n < 10:
97              sw.addRecord([i*2., i*0.67, i**2, -i*7])                  sw.addRecord([i*2., i*0.67, i**2, -i*7])
98         sw.write('example.segy')             sw.write('example.segy')
99    
100      :note: the writer uses `obsy`          :note: the writer uses `obsy`
101      """          """
102      COORDINATE_SCALE = 1000.          COORDINATE_SCALE = 1000.
103      def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):          def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):
104          """                  """
105          initalize writer                  initalize writer
106    
107          :param receiver_group: list of receiver coordinates (in meters). For the 2D case a list of floats is given, for the 3D case a list of coordinate tuples                  :param receiver_group: list of receiver coordinates (in meters). For the 2D case a list of floats is given, for the 3D case a list of coordinate tuples
108                                 are given                                         are given
109          :param source: coordinates of the source (in meters).  For the 2D case a single floats is given, for the 3D case a coordinate tuples                  :param source: coordinates of the source (in meters).  For the 2D case a single floats is given, for the 3D case a coordinate tuples
110          :param sampling_interval: sample rate in seconds                  :param sampling_interval: sample rate in seconds
111          :param text: a text for the header file (e.g a description)                  :param text: a text for the header file (e.g a description)
112          """                  """
113          if isinstance(source, float) or  isinstance(source, int) :                  if isinstance(source, float) or  isinstance(source, int) :
114              DIM=1                          DIM=1
115              source = (source, 0.)                          source = (source, 0.)
116          elif hasattr(source, '__len__'):                  elif hasattr(source, '__len__'):
117              DIM=len(source)                          DIM=len(source)
118              if DIM == 1:                          if DIM == 1:
119                  source = (source[0], 0.)                              source = (source[0], 0.)
120                  elif DIM==2:                          elif DIM==2:
121                  source = (source[0], source[1] )                              source = (source[0], source[1] )
122              else:                          else:
123                      raise ValueError("Only 1D or 2D arrays accepted.")                              raise ValueError("Only 1D or 2D arrays accepted.")
124              else:                  else:
125              raise TypeError("illegal source type.")                      raise TypeError("illegal source type.")
126              if receiver_group== None:                  if receiver_group== None:
127              if DIM == 1:                          if DIM == 1:
128                  receiver_group=[0.]                                  receiver_group=[0.]
129              else:                          else:
130                  receiver_group=[(0.,0.)]                                  receiver_group=[(0.,0.)]
131    
132          if isinstance(receiver_group, float) or  isinstance(receiver_group, int) :                  if isinstance(receiver_group, float) or  isinstance(receiver_group, int) :
133              if DIM == 1:                          if DIM == 1:
134                  rg = [ (receiver_group, 0. ) ]                              rg = [ (receiver_group, 0. ) ]
135              else:                          else:
136                      raise TypeError("illegal receiver_group type.")                              raise TypeError("illegal receiver_group type.")
137          elif hasattr(receiver_group, '__len__'):                  elif hasattr(receiver_group, '__len__'):
138              if DIM == 1:                          if DIM == 1:
139                  rg = [(c,0)  for c in receiver_group]                              rg = [(c,0)  for c in receiver_group]
140              else:                          else:
141                  rg = [(c[0],c[1])  for c in receiver_group]                              rg = [(c[0],c[1])  for c in receiver_group]
142              else:                  else:
143              raise TypeError("illegal receiver_group type.")                      raise TypeError("illegal receiver_group type.")
144                                                            
145          self.__source=source                  self.__source=source
146          self.__receiver_group=rg                  self.__receiver_group=rg
147          self.__text=text                  self.__text=text
148                  self.__trace = [ [] for i in self.__receiver_group ]                  self.__trace = [ [] for i in self.__receiver_group ]
149                  self.__sampling_interval = sampling_interval                    self.__sampling_interval = sampling_interval  
150    
151      def addRecord(self, record):          def addRecord(self, record):
152             """                 """
153             adds a record to the traces. a time difference of sample_interval between two records is assumed.                 adds a record to the traces. a time difference of sample_interval between two records is assumed.
154             The record mast be a list of as many values as given receivers or a float if a single receiver is used.                 The record mast be a list of as many values as given receivers or a float if a single receiver is used.
155    
156                 :param record: list of tracks to be added to the record.                 :param record: list of tracks to be added to the record.
157             """                 """
158             if not len(self.__receiver_group) == len(record):                 if not len(self.__receiver_group) == len(record):
159                       raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))                           raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))
160             if len(self.__receiver_group) == 1:                 if len(self.__receiver_group) == 1:
161                 if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):                         if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
162                       self.__trace[0].append(record)                                   self.__trace[0].append(record)
163                 else:                         else:
164                       self.__trace[0].append(record[0])                                   self.__trace[0].append(record[0])
165             else:                 else:
166                 for i in xrange(len(record)):                         for i in xrange(len(record)):
167                     self.__trace[i].append(record[i])                                 self.__trace[i].append(record[i])
168    
169      def getSamplingInterval(self):          def getSamplingInterval(self):
170             """                 """
171             returns the sampling interval in seconds.                 returns the sampling interval in seconds.
172             """                 """
173                 return self.__sampling_interval                 return self.__sampling_interval
174    
175      def write(self, filename):          def write(self, filename):
176              """              """
177              writes to segy file              writes to segy file
178                            
# Line 180  class SimpleSEGYWriter(object): Line 180  class SimpleSEGYWriter(object):
180              :note: the function uses the `obspy` module.              :note: the function uses the `obspy` module.
181              """              """
182              from obspy import Trace, Stream, UTCDateTime              from obspy import Trace, Stream, UTCDateTime
183          from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader              from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
184              from obspy.core import AttribDict              from obspy.core import AttribDict
185    
186          stream=Stream()              stream=Stream()
187                        
188          for i in xrange(len(self.__receiver_group)):              for i in xrange(len(self.__receiver_group)):
189              trace = Trace(data=np.array(self.__trace[i], dtype='float32'))                      trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
190                      # Attributes in trace.stats will overwrite everything in                      # Attributes in trace.stats will overwrite everything in
191                      # trace.stats.segy.trace_header (in Hz)                      # trace.stats.segy.trace_header (in Hz)
192                      trace.stats.sampling_rate = 1./self.getSamplingInterval()                      trace.stats.sampling_rate = 1./self.getSamplingInterval()
# Line 245  class WaveBase(object): Line 245  class WaveBase(object):
245           self.__dt=dt           self.__dt=dt
246                    
247        def getTimeStepSize(self):        def getTimeStepSize(self):
248             return self.__dt             return self.__dt
249                        
250        def _getAcceleration(self, t, u):        def _getAcceleration(self, t, u):
251             """             """
252             returns the acceleraton for time t and solution u at time t             returns the acceleraton for time t and solution u at time t
# Line 266  class WaveBase(object): Line 266  class WaveBase(object):
266               # apply Verlet scheme until               # apply Verlet scheme until
267               while self.t < t:               while self.t < t:
268                    a=self._getAcceleration(self.t, self.u)                    a=self._getAcceleration(self.t, self.u)
269            self.v += dt * a                    self.v += dt * a
270            self.u += dt * self.v                    self.u += dt * self.v
271                    self.t += dt                    self.t += dt
272                            
273               # now we work backwards               # now we work backwards
# Line 295  def createAbsorbtionLayerFunction(x, abs Line 295  def createAbsorbtionLayerFunction(x, abs
295      for i in xrange(DIM):      for i in xrange(DIM):
296          x_i=x[i]          x_i=x[i]
297          x_l=x_i-(bb[i][0]+absorption_zone)          x_l=x_i-(bb[i][0]+absorption_zone)
298      m_l=whereNegative(x_l)          m_l=whereNegative(x_l)
299      f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )          f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
300      if not DIM-1 == i:          if not DIM-1 == i:
301          x_r=(bb[i][1]-absorption_zone)-x_i              x_r=(bb[i][1]-absorption_zone)-x_i
302          m_r=whereNegative(x_r)              m_r=whereNegative(x_r)
303          f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )              f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
304      return f      return f
305            
306  class SonicWave(WaveBase):  class SonicWave(WaveBase):
307      """          """
308      Solving the sonic wave equation          Solving the sonic wave equation
309                
310      p_tt = (v_p**2 * p_i)_i  + f(t) * delta_s   where (p-) velocity v_p.          p_tt = (v_p**2 * p_i)_i  + f(t) * delta_s   where (p-) velocity v_p.
311                
312      f(t) is wavelet acting at a point source term at positon s          f(t) is wavelet acting at a point source term at positon s
313      """          """
314      def __init__(self, domain, v_p, wavelet, source_tag, dt=None, p0=None, p0_t=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):          def __init__(self, domain, v_p, wavelet, source_tag, dt=None, p0=None, p0_t=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
315             """             """
316             initialize the sonic wave solver             initialize the sonic wave solver
317                        
# Line 333  class SonicWave(WaveBase): Line 333  class SonicWave(WaveBase):
333             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
334             v_p=v_p*f             v_p=v_p*f
335    
336         if p0 == None:             if p0 == None:
337            p0=Scalar(0.,Solution(domain))                p0=Scalar(0.,Solution(domain))
338         else:             else:
339            p0=interpolate(p0, Solution(domain ))                p0=interpolate(p0, Solution(domain ))
340                            
341         if p0_t == None:             if p0_t == None:
342            p0_t=Scalar(0.,Solution(domain))                p0_t=Scalar(0.,Solution(domain))
343         else:             else:
344            p0_t=interpolate(p0_t, Solution(domain ))                p0_t=interpolate(p0_t, Solution(domain ))
345                    
346         if dt == None:             if dt == None:
347                    dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())                    dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
348                            
349             super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)             super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
# Line 353  class SonicWave(WaveBase): Line 353  class SonicWave(WaveBase):
353             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
354             self.__mypde.setSymmetryOn()             self.__mypde.setSymmetryOn()
355             self.__mypde.setValue(D=1.)             self.__mypde.setValue(D=1.)
356         self.__source_tag=source_tag             self.__source_tag=source_tag
357             self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))             self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
358         self.__vp2=v_p**2             self.__vp2=v_p**2
359    
360    
361          def  _getAcceleration(self, t, u):          def  _getAcceleration(self, t, u):

Legend:
Removed from v.4555  
changed lines
  Added in v.4556

  ViewVC Help
Powered by ViewVC 1.1.26