/[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 4575 by gross, Mon Dec 9 03:42:50 2013 UTC revision 4576 by sshaw, Mon Dec 9 23:35:30 2013 UTC
# Line 165  class SimpleSEGYWriter(object): Line 165  class SimpleSEGYWriter(object):
165                         else:                         else:
166                                   self.__trace[0].append(record[0])                                   self.__trace[0].append(record[0])
167                 else:                 else:
168                         for i in xrange(len(record)):                         for i in range(len(record)):
169                                 self.__trace[i].append(record[i])                                 self.__trace[i].append(record[i])
170    
171          def getSamplingInterval(self):          def getSamplingInterval(self):
# Line 187  class SimpleSEGYWriter(object): Line 187  class SimpleSEGYWriter(object):
187    
188              stream=Stream()              stream=Stream()
189                            
190              for i in xrange(len(self.__receiver_group)):              for i in range(len(self.__receiver_group)):
191                      trace = Trace(data=np.array(self.__trace[i], dtype='float32'))                      trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
192                      # Attributes in trace.stats will overwrite everything in                      # Attributes in trace.stats will overwrite everything in
193                      # trace.stats.segy.trace_header (in Hz)                      # trace.stats.segy.trace_header (in Hz)
# Line 294  def createAbsorbtionLayerFunction(x, abs Line 294  def createAbsorbtionLayerFunction(x, abs
294      DIM=dom.getDim()      DIM=dom.getDim()
295      decay=-log(absorption_cut)/absorption_zone**2      decay=-log(absorption_cut)/absorption_zone**2
296      f=1      f=1
297      for i in xrange(DIM):      for i in range(DIM):
298          x_i=x[i]          x_i=x[i]
299          x_l=x_i-(bb[i][0]+absorption_zone)          x_l=x_i-(bb[i][0]+absorption_zone)
300          m_l=whereNegative(x_l)          m_l=whereNegative(x_l)
# Line 436  class VTIWave(WaveBase): Line 436  class VTIWave(WaveBase):
436                        
437    
438             if DIM ==2 :             if DIM ==2 :
439            source_vector= [source_vector[0],source_vector[2]]                source_vector= [source_vector[0],source_vector[2]]
440    
441                    
442         self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))             self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
443             self.__r.setTaggedValue(self.__source_tag, source_vector)             self.__r.setTaggedValue(self.__source_tag, source_vector)
444    
445                        
446             self.c33=v_p**2 * rho             self.c33=v_p**2 * rho
447         self.c44=v_s**2 * rho             self.c44=v_s**2 * rho
448         self.c11=(1+2*eps) * self.c33             self.c11=(1+2*eps) * self.c33
449         self.c66=(1+2*gamma) * self.c44             self.c66=(1+2*gamma) * self.c44
450         self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44             self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
451         self.c12=self.c11-2*self.c66             self.c12=self.c11-2*self.c66
452                    
453          def  _getAcceleration(self, t, u):          def  _getAcceleration(self, t, u):
454               """               """
455               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 458  class VTIWave(WaveBase): Line 458  class VTIWave(WaveBase):
458               sigma=self.__mypde.getCoefficient('X')               sigma=self.__mypde.getCoefficient('X')
459                            
460               if self.__mypde.getDim() == 3:               if self.__mypde.getDim() == 3:
461          e11=du[0,0]                  e11=du[0,0]
462          e22=du[1,1]                  e22=du[1,1]
463          e33=du[2,2]                  e33=du[2,2]
464                            
465          sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33                  sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
466          sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33                  sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
467          sigma[2,2]=self.c13*(e11+e22)+self.c33*e33                  sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
468    
469          s=self.c44*(du[2,1]+du[1,2])                  s=self.c44*(du[2,1]+du[1,2])
470          sigma[1,2]=s                  sigma[1,2]=s
471          sigma[2,1]=s                              sigma[2,1]=s            
472    
473          s=self.c44*(du[2,0]+du[0,2])                  s=self.c44*(du[2,0]+du[0,2])
474          sigma[0,2]=s                  sigma[0,2]=s
475          sigma[2,0]=s                  sigma[2,0]=s
476                            
477          s=self.c66*(du[0,1]+du[1,0])                  s=self.c66*(du[0,1]+du[1,0])
478          sigma[0,1]=s                  sigma[0,1]=s
479          sigma[1,0]=s                  sigma[1,0]=s
480                            
481    
482               else:               else:
483          e11=du[0,0]                  e11=du[0,0]
484          e22=du[1,1]                  e22=du[1,1]
485                            
486          sigma[0,0]=self.c11*e11+self.c13*e22                  sigma[0,0]=self.c11*e11+self.c13*e22
487          sigma[1,1]=self.c13*e11+self.c33*e22                  sigma[1,1]=self.c13*e11+self.c33*e22
488                            
489          s=self.c44*(du[1,0]+du[0,1])                  s=self.c44*(du[1,0]+du[0,1])
490          sigma[0,1]=s                  sigma[0,1]=s
491          sigma[1,0]=s                  sigma[1,0]=s
492                            
493               self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))               self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
494               return self.__mypde.getSolution()               return self.__mypde.getSolution()
# Line 557  class HTIWave(WaveBase): Line 557  class HTIWave(WaveBase):
557                        
558    
559             if DIM ==2 :             if DIM ==2 :
560            source_vector= [source_vector[0],source_vector[1]]                source_vector= [source_vector[0],source_vector[1]]
561    
562                    
563         self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))             self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
564             self.__r.setTaggedValue(self.__source_tag, source_vector)             self.__r.setTaggedValue(self.__source_tag, source_vector)
565    
566                        
567             self.c33=v_p**2 * rho             self.c33=v_p**2 * rho
568         self.c44=v_s**2 * rho             self.c44=v_s**2 * rho
569         self.c11=(1+2*eps) * self.c33             self.c11=(1+2*eps) * self.c33
570         self.c66=(1+2*gamma) * self.c44             self.c66=(1+2*gamma) * self.c44
571         self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44             self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
572         self.c23=self.c33-2*self.c66             self.c23=self.c33-2*self.c66
573                    
574          def  _getAcceleration(self, t, u):          def  _getAcceleration(self, t, u):
575               """               """
576               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 579  class HTIWave(WaveBase): Line 579  class HTIWave(WaveBase):
579               sigma=self.__mypde.getCoefficient('X')               sigma=self.__mypde.getCoefficient('X')
580                            
581               if self.__mypde.getDim() == 3:               if self.__mypde.getDim() == 3:
582          e11=du[0,0]                  e11=du[0,0]
583          e22=du[1,1]                  e22=du[1,1]
584          e33=du[2,2]                  e33=du[2,2]
585                            
586          sigma[0,0]=self.c11*e11+self.c13*(e22+e33)                  sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
587          sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33                  sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
588          sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33                  sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
589    
590          s=self.c44*(du[2,1]+du[1,2])                  s=self.c44*(du[2,1]+du[1,2])
591          sigma[1,2]=s                  sigma[1,2]=s
592          sigma[2,1]=s                              sigma[2,1]=s            
593    
594          s=self.c66*(du[2,0]+du[0,2])                  s=self.c66*(du[2,0]+du[0,2])
595          sigma[0,2]=s                  sigma[0,2]=s
596          sigma[2,0]=s                  sigma[2,0]=s
597                            
598          s=self.c66*(du[0,1]+du[1,0])                  s=self.c66*(du[0,1]+du[1,0])
599          sigma[0,1]=s                  sigma[0,1]=s
600          sigma[1,0]=s                  sigma[1,0]=s
601    
602               else:               else:
603          e11=du[0,0]                  e11=du[0,0]
604          e22=du[1,1]                  e22=du[1,1]
605                            
606          sigma[0,0]=self.c11*e11+self.c13*e22                  sigma[0,0]=self.c11*e11+self.c13*e22
607          sigma[1,1]=self.c13*e11+self.c33*e22                  sigma[1,1]=self.c13*e11+self.c33*e22
608                            
609          s=self.c66*(du[0,1]+du[1,0])                  s=self.c66*(du[0,1]+du[1,0])
610          sigma[0,1]=s                  sigma[0,1]=s
611          sigma[1,0]=s                  sigma[1,0]=s
612                            
613               self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))               self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
614               return self.__mypde.getSolution()                           return self.__mypde.getSolution()            

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

  ViewVC Help
Powered by ViewVC 1.1.26