/[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 4621 by gross, Thu Jan 9 05:30:27 2014 UTC revision 4622 by sshaw, Fri Jan 17 04:55:41 2014 UTC
# Line 20  __license__="""Licensed under the Open S Line 20  __license__="""Licensed under the Open S
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
21  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
22    
23  __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave', 'HTIWave', 'createAbsorbtionLayerFunction' ]  __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave',
24                'HTIWave', 'createAbsorbtionLayerFunction' ]
25    
26    
27  from math import pi  from math import pi
# Line 29  import sys Line 30  import sys
30  import time  import time
31  from esys.escript import *  from esys.escript import *
32  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
33  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem, VTIWavePDE
34    
   
   
               
35  class Wavelet(object):  class Wavelet(object):
36          """          """
37          place holder for source wavelet          place holder for source wavelet
38          """          """
39          pass          pass
40            
41  class Ricker(Wavelet):  class Ricker(Wavelet):
42          """          """
43          The Ricker Wavelet w=f(t)          The Ricker Wavelet w=f(t)
# Line 48  class Ricker(Wavelet): Line 46  class Ricker(Wavelet):
46                  """                  """
47                  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
48                  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.
49                    
50                  :note: maximum frequence is about 2 x the dominant frequence.                  :note: maximum frequence is about 2 x the dominant frequence.
51                  """                  """
52                  drop=18                  drop=18
# Line 58  class Ricker(Wavelet): Line 56  class Ricker(Wavelet):
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    
# Line 85  class Ricker(Wavelet): Line 83  class Ricker(Wavelet):
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  class SimpleSEGYWriter(object):  class SimpleSEGYWriter(object):
88          """          """
# Line 125  class SimpleSEGYWriter(object): Line 122  class SimpleSEGYWriter(object):
122                              raise ValueError("Only 1D or 2D arrays accepted.")                              raise ValueError("Only 1D or 2D arrays accepted.")
123                  else:                  else:
124                      raise TypeError("illegal source type.")                      raise TypeError("illegal source type.")
125                  if receiver_group== None:                  if receiver_group== None:
126                          if DIM == 1:                          if DIM == 1:
127                                  receiver_group=[0.]                                  receiver_group=[0.]
128                          else:                          else:
# Line 143  class SimpleSEGYWriter(object): Line 140  class SimpleSEGYWriter(object):
140                              rg = [(c[0],c[1])  for c in receiver_group]                              rg = [(c[0],c[1])  for c in receiver_group]
141                  else:                  else:
142                      raise TypeError("illegal receiver_group type.")                      raise TypeError("illegal receiver_group type.")
143                                    
144                  self.__source=source                  self.__source=source
145                  self.__receiver_group=rg                  self.__receiver_group=rg
146                  self.__text=text                  self.__text=text
147                  self.__trace = [ [] for i in self.__receiver_group ]                  self.__trace = [ [] for i in self.__receiver_group ]
148                  self.__sampling_interval = sampling_interval                    self.__sampling_interval = sampling_interval
149    
150          def addRecord(self, record):          def addRecord(self, record):
151                 """                 """
152                 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.
153                 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.
154    
155                 :param record: list of tracks to be added to the record.                 :param record: list of tracks to be added to the record.
156                 """                 """
157                 if not len(self.__receiver_group) == len(record):                 if not len(self.__receiver_group) == len(record):
158                           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)))
159                 if len(self.__receiver_group) == 1:                 if len(self.__receiver_group) == 1:
160                         if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):                         if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
# Line 170  class SimpleSEGYWriter(object): Line 167  class SimpleSEGYWriter(object):
167    
168          def getSamplingInterval(self):          def getSamplingInterval(self):
169                 """                 """
170                 returns the sampling interval in seconds.                 returns the sampling interval in seconds.
171                 """                 """
172                 return self.__sampling_interval                 return self.__sampling_interval
173    
174          def write(self, filename):          def write(self, filename):
175              """              """
176              writes to segy file              writes to segy file
177                
178              :param filename: file name              :param filename: file name
179              :note: the function uses the `obspy` module.              :note: the function uses the `obspy` module.
180              """              """
181              try:              try:
182                  from obspy import Trace, Stream, UTCDateTime                  from obspy import Trace, Stream, UTCDateTime
# Line 191  class SimpleSEGYWriter(object): Line 188  class SimpleSEGYWriter(object):
188                          "https://github.com/obspy/obspy for install guide")                          "https://github.com/obspy/obspy for install guide")
189    
190              stream=Stream()              stream=Stream()
191                
192              for i in range(len(self.__receiver_group)):              for i in range(len(self.__receiver_group)):
193                      trace = Trace(data=np.array(self.__trace[i], dtype='float32'))                      trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
194                      # Attributes in trace.stats will overwrite everything in                      # Attributes in trace.stats will overwrite everything in
# Line 224  class SimpleSEGYWriter(object): Line 221  class SimpleSEGYWriter(object):
221    
222  class WaveBase(object):  class WaveBase(object):
223        """        """
224        Base for wave propagation using the Verlet scheme.        Base for wave propagation using the Verlet scheme.
225                
226              u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0              u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
227    
228        with a given  acceleration force as function of time.        with a given  acceleration force as function of time.
229          
230        a_n=A(t_{n-1})        a_n=A(t_{n-1})
231        v_n=v_{n-1} +  dt *  a_n        v_n=v_{n-1} +  dt *  a_n
232        u_n=u_{n-1} +  dt *  v_n        u_n=u_{n-1} +  dt *  v_n
233          
234          
235        """        """
236        def __init__(self, dt, u0, v0, t0=0.):        def __init__(self, dt, u0, v0, t0=0.):
237           """           """
238           set up the wave base           set up the wave base
239            
240           :param dt: time step size (need to be sufficiently small)           :param dt: time step size (need to be sufficiently small)
241           :param u0: initial value           :param u0: initial value
242           :param v0: initial velocity           :param v0: initial velocity
243           :param t0: initial time           :param t0: initial time
244           """           """
245           self.u=u0           self.u=u0
246           self.v=v0           self.v=v0
247           self.t=t0           self.t=t0
248           self.t_last=t0           self.t_last=t0
249           self.__dt=dt           self.__dt=dt
250            
251        def getTimeStepSize(self):        def getTimeStepSize(self):
252             return self.__dt             return self.__dt
253              
254        def _getAcceleration(self, t, u):        def _getAcceleration(self, t, u):
255             """             """
256             returns the acceleraton for time t and solution u at time t             returns the acceleraton for time t and solution u at time t
257             : note: this function needs to be overwritten by a particular wave problem             : note: this function needs to be overwritten by a particular wave problem
258             """             """
259             pass             pass
260              
261        def update(self, t):        def update(self, t):
262               """               """
263               returns the solution for the next time marker t which needs to greater than the time marker from the               returns the solution for the next time marker t which needs to greater than the time marker from the
264               previous call.               previous call.
265               """               """
266               if not self.t_last <= t:               if not self.t_last <= t:
267                    raise ValueError("You can not move backward in time.")                    raise ValueError("You can not move backward in time.")
268                      
269               dt = self.getTimeStepSize()               dt = self.getTimeStepSize()
270               # apply Verlet scheme until               # apply Verlet scheme until
271               while self.t < t:               while self.t < t:
272                    a=self._getAcceleration(self.t, self.u)                    a=self._getAcceleration(self.t, self.u)
273                    self.v += dt * a                    self.v += dt * a
274                    self.u += dt * self.v                    self.u += dt * self.v
275                    self.t += dt                    self.t += dt
276                
277               # now we work backwards               # now we work backwards
278               self.t_last=t               self.t_last=t
279               return t, self.u + self.v * (t-self.t)               return t, self.u + self.v * (t-self.t)
280                
281  def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):  def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
282      """      """
283      creating a distribution which is one in the interior of the domain of x      creating a distribution which is one in the interior of the domain of x
284      and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'      and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
285      toward each boundary except the top of the domain.      toward each boundary except the top of the domain.
286        
287      :param x: location of points in the domain      :param x: location of points in the domain
288      :type x: `Data`      :type x: `Data`
289      :param absorption_zone: thickness of the aborption zone      :param absorption_zone: thickness of the aborption zone
290      :param absorption_cut: value of decay function on domain boundary      :param absorption_cut: value of decay function on domain boundary
291      :return: function on 'x' which is one in the iterior and decays to almost zero over a margin      :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
292               toward the boundary.                           toward the boundary.
293      """          """
294      dom=x.getDomain()        dom=x.getDomain()
295      bb=boundingBox(dom)      bb=boundingBox(dom)
296      DIM=dom.getDim()      DIM=dom.getDim()
297      decay=-log(absorption_cut)/absorption_zone**2      decay=-log(absorption_cut)/absorption_zone**2
# Line 309  def createAbsorbtionLayerFunction(x, abs Line 306  def createAbsorbtionLayerFunction(x, abs
306              m_r=whereNegative(x_r)              m_r=whereNegative(x_r)
307              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 )
308      return f      return f
309        
310  class SonicWave(WaveBase):  class SonicWave(WaveBase):
311          """          """
312          Solving the sonic wave equation          Solving the sonic wave equation
313            
314          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.
315            
316          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
317          """          """
318          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):
319             """             """
320             initialize the sonic wave solver             initialize the sonic wave solver
321              
322             :param domain: domain of the problem             :param domain: domain of the problem
323             :type domain: `Doamin`                     :type domain: `Doamin`
324             :param v_p: p-velocity field                 :param v_p: p-velocity field
325             :type v_p: `Scalar`                   :type v_p: `Scalar`
326             :param wavelet: wavelet to describe the time evolution of source term             :param wavelet: wavelet to describe the time evolution of source term
327             :type wavelet: `Wavelet`                       :type wavelet: `Wavelet`
328             :param source_tag: tag of the source location             :param source_tag: tag of the source location
329             :type source_tag: 'str' or 'int'                       :type source_tag: 'str' or 'int'
330             :param dt: time step size. If not present a suitable time step size is calculated.                       :param dt: time step size. If not present a suitable time step size is calculated.
331             :param p0: initial solution. If not present zero is used.                       :param p0: initial solution. If not present zero is used.
332             :param p0_t: initial solution change rate. If not present zero is used.                       :param p0_t: initial solution change rate. If not present zero is used.
333             :param absorption_zone: thickness of absorption zone                       :param absorption_zone: thickness of absorption zone
334             :param absorption_cut: boundary value of absorption decay factor             :param absorption_cut: boundary value of absorption decay factor
335             :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.             :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
336             """             """
337             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
338             v_p=v_p*f             v_p=v_p*f
# Line 344  class SonicWave(WaveBase): Line 341  class SonicWave(WaveBase):
341                p0=Scalar(0.,Solution(domain))                p0=Scalar(0.,Solution(domain))
342             else:             else:
343                p0=interpolate(p0, Solution(domain ))                p0=interpolate(p0, Solution(domain ))
344                  
345             if p0_t == None:             if p0_t == None:
346                p0_t=Scalar(0.,Solution(domain))                p0_t=Scalar(0.,Solution(domain))
347             else:             else:
348                p0_t=interpolate(p0_t, Solution(domain ))                p0_t=interpolate(p0_t, Solution(domain ))
349              
350             if dt == None:             if dt == None:
351                    dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())                    dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
352                
353             super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)             super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
354              
355             self.__wavelet=wavelet             self.__wavelet=wavelet
356             self.__mypde=LinearSinglePDE(domain)             self.__mypde=LinearSinglePDE(domain)
357             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
# Line 372  class SonicWave(WaveBase): Line 369  class SonicWave(WaveBase):
369               self.__mypde.setValue(X=-grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)               self.__mypde.setValue(X=-grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
370               return self.__mypde.getSolution()               return self.__mypde.getSolution()
371    
               
372    
373  class VTIWave(WaveBase):  class VTIWave(WaveBase):
374          """      """
375          Solving the VTI wave equation      Solving the VTI wave equation
           
         :note: In case of a two dimensional domain the second spatial dimenion is depth.  
         """  
         def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,      
                      dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):  
            """  
            initialize the VTI wave solver  
             
            :param domain: domain of the problem  
            :type domain: `Doamin`          
            :param v_p: vertical p-velocity field      
            :type v_p: `Scalar`  
            :param v_s: vertical s-velocity field      
            :type v_s: `Scalar`            
            :param wavelet: wavelet to describe the time evolution of source term  
            :type wavelet: `Wavelet`            
            :param source_tag: tag of the source location  
            :type source_tag: 'str' or 'int'  
            :param source_vector: source orientation vector  
            :param eps: first Thompsen parameter  
            :param delta: second Thompsen parameter  
            :param gamma: third Thompsen parameter  
            :param rho: density            
            :param dt: time step size. If not present a suitable time step size is calculated.            
            :param u0: initial solution. If not present zero is used.            
            :param v0: initial solution change rate. If not present zero is used.            
            :param absorption_zone: thickness of absorption zone            
            :param absorption_cut: boundary value of absorption decay factor  
            :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.  
            """  
            DIM=domain.getDim()  
            f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)  
             
            v_p=v_p*f  
            v_s=v_s*f  
   
             
             
            if u0 == None:  
               u0=Vector(0.,Solution(domain))  
            else:  
               u0=interpolate(p0, Solution(domain ))  
                 
            if v0 == None:  
               v0=Vector(0.,Solution(domain))  
            else:  
               v0=interpolate(v0, Solution(domain ))  
             
            if dt == None:  
                   dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())  
               
            super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)  
             
            self.__wavelet=wavelet  
             
            self.__mypde=LinearPDESystem(domain)  
            if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)  
            self.__mypde.setSymmetryOn()  
            self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))  
            self.__source_tag=source_tag  
             
   
            if DIM ==2 :  
               source_vector= [source_vector[0],source_vector[2]]  
   
             
            self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))  
            self.__r.setTaggedValue(self.__source_tag, source_vector)  
   
             
            self.c33=v_p**2 * rho  
            self.c44=v_s**2 * rho  
            self.c11=(1+2*eps) * self.c33  
            self.c66=(1+2*gamma) * self.c44  
            self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44  
            self.c12=self.c11-2*self.c66  
             
         def  _getAcceleration(self, t, u):  
              """  
              returns the acceleraton for time t and solution u at time t  
              """  
              du = grad(u)  
              sigma=self.__mypde.getCoefficient('X')  
               
              if self.__mypde.getDim() == 3:  
                 e11=du[0,0]  
                 e22=du[1,1]  
                 e33=du[2,2]  
               
                 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33  
                 sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33  
                 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33  
   
                 s=self.c44*(du[2,1]+du[1,2])  
                 sigma[1,2]=s  
                 sigma[2,1]=s              
376    
377                  s=self.c44*(du[2,0]+du[0,2])      :note: In case of a two dimensional domain the second spatial dimenion is depth.
378                  sigma[0,2]=s      """
379                  sigma[2,0]=s      def __init__(self, domain, v_p, v_s, wavelet, source_tag,
380                                source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
381                  s=self.c66*(du[0,1]+du[1,0])              dt=None, u0=None, v0=None, absorption_zone=300*U.m,
382                  sigma[0,1]=s              absorption_cut=1e-2, lumping=True):
383                  sigma[1,0]=s          """
384                            initialize the VTI wave solver
385    
386            :param domain: domain of the problem
387            :type domain: `Domain`
388            :param v_p: vertical p-velocity field
389            :type v_p: `Scalar`
390            :param v_s: vertical s-velocity field
391            :type v_s: `Scalar`
392            :param wavelet: wavelet to describe the time evolution of source term
393            :type wavelet: `Wavelet`
394            :param source_tag: tag of the source location
395            :type source_tag: 'str' or 'int'
396            :param source_vector: source orientation vector
397            :param eps: first Thompsen parameter
398            :param delta: second Thompsen parameter
399            :param gamma: third Thompsen parameter
400            :param rho: density
401            :param dt: time step size. If not present a suitable time step size is calculated.
402            :param u0: initial solution. If not present zero is used.
403            :param v0: initial solution change rate. If not present zero is used.
404            :param absorption_zone: thickness of absorption zone
405            :param absorption_cut: boundary value of absorption decay factor
406            :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
407            """
408            DIM=domain.getDim()
409            f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
410    
411            v_p=v_p*f
412            v_s=v_s*f
413    
414    
415    
416            if u0 == None:
417              u0=Vector(0.,Solution(domain))
418            else:
419              u0=interpolate(p0, Solution(domain ))
420    
421            if v0 == None:
422              v0=Vector(0.,Solution(domain))
423            else:
424              v0=interpolate(v0, Solution(domain ))
425    
426            if dt == None:
427                  dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
428    
429            super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
430    
431            self.__wavelet=wavelet
432    
433            c33=v_p**2 * rho
434            c44=v_s**2 * rho
435            c11=(1+2*eps) * c33
436            c66=(1+2*gamma) * c44
437            c13=sqrt(2*c33*(c33-c44) * delta + (c33-c44)**2)-c44
438            c12=c11-2*c66
439            self.__mypde=VTIWavePDE(domain, [c11, c12, c13, c33, c44, c66])
440    
441            if lumping:
442                self.__mypde.getSolverOptions().setSolverMethod(
443                        self.__mypde.getSolverOptions().HRZ_LUMPING)
444            self.__mypde.setSymmetryOn()
445            self.__mypde.setValue(D=rho*kronecker(DIM))
446            self.__source_tag=source_tag
447    
448            if DIM ==2 :
449              source_vector= [source_vector[0],source_vector[2]]
450    
451            self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
452            self.__r.setTaggedValue(self.__source_tag, source_vector)
453    
454        def _getAcceleration(self, t, u):
455            """
456            returns the acceleraton for time t and solution u at time t
457            """
458            du = grad(u)
459            self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
460            return self.__mypde.getSolution()
461    
              else:  
                 e11=du[0,0]  
                 e22=du[1,1]  
               
                 sigma[0,0]=self.c11*e11+self.c13*e22  
                 sigma[1,1]=self.c13*e11+self.c33*e22  
               
                 s=self.c44*(du[1,0]+du[0,1])  
                 sigma[0,1]=s  
                 sigma[1,0]=s  
               
              self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))  
              return self.__mypde.getSolution()  
462    
463  class HTIWave(WaveBase):  class HTIWave(WaveBase):
464          """          """
465          Solving the HTI wave equation (along the x_0 axis)          Solving the HTI wave equation (along the x_0 axis)
466            
467          :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.          :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
468          """          """
469          def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,              
470                       dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):          def __init__(self, domain, v_p, v_s,   wavelet, source_tag,
471                    source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
472                    dt=None, u0=None, v0=None, absorption_zone=300*U.m,
473                    absorption_cut=1e-2, lumping=True):
474             """             """
475             initialize the VTI wave solver             initialize the VTI wave solver
476              
477             :param domain: domain of the problem             :param domain: domain of the problem
478             :type domain: `Doamin`                     :type domain: `Doamin`
479             :param v_p: vertical p-velocity field                 :param v_p: vertical p-velocity field
480             :type v_p: `Scalar`             :type v_p: `Scalar`
481             :param v_s: vertical s-velocity field                 :param v_s: vertical s-velocity field
482             :type v_s: `Scalar`                       :type v_s: `Scalar`
483             :param wavelet: wavelet to describe the time evolution of source term             :param wavelet: wavelet to describe the time evolution of source term
484             :type wavelet: `Wavelet`                       :type wavelet: `Wavelet`
485             :param source_tag: tag of the source location             :param source_tag: tag of the source location
486             :type source_tag: 'str' or 'int'             :type source_tag: 'str' or 'int'
487             :param source_vector: source orientation vector             :param source_vector: source orientation vector
488             :param eps: first Thompsen parameter             :param eps: first Thompsen parameter
489             :param delta: second Thompsen parameter             :param delta: second Thompsen parameter
490             :param gamma: third Thompsen parameter             :param gamma: third Thompsen parameter
491             :param rho: density                       :param rho: density
492             :param dt: time step size. If not present a suitable time step size is calculated.                       :param dt: time step size. If not present a suitable time step size is calculated.
493             :param u0: initial solution. If not present zero is used.                       :param u0: initial solution. If not present zero is used.
494             :param v0: initial solution change rate. If not present zero is used.                       :param v0: initial solution change rate. If not present zero is used.
495             :param absorption_zone: thickness of absorption zone                       :param absorption_zone: thickness of absorption zone
496             :param absorption_cut: boundary value of absorption decay factor             :param absorption_cut: boundary value of absorption decay factor
497             :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.             :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
498             """             """
499             DIM=domain.getDim()             DIM=domain.getDim()
500             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)             f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
501              
502             v_p=v_p*f             v_p=v_p*f
503             v_s=v_s*f             v_s=v_s*f
504              
505             if u0 == None:             if u0 == None:
506                u0=Vector(0.,Solution(domain))                u0=Vector(0.,Solution(domain))
507             else:             else:
508                u0=interpolate(p0, Solution(domain ))                u0=interpolate(p0, Solution(domain ))
509                  
510             if v0 == None:             if v0 == None:
511                v0=Vector(0.,Solution(domain))                v0=Vector(0.,Solution(domain))
512             else:             else:
513                v0=interpolate(v0, Solution(domain ))                v0=interpolate(v0, Solution(domain ))
514              
515             if dt == None:             if dt == None:
516                    dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())                    dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
517                
518             super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)             super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
519              
520             self.__wavelet=wavelet             self.__wavelet=wavelet
521              
522             self.__mypde=LinearPDESystem(domain)             self.__mypde=LinearPDESystem(domain)
523             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)             if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
524             self.__mypde.setSymmetryOn()             self.__mypde.setSymmetryOn()
525             self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))             self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
526             self.__source_tag=source_tag             self.__source_tag=source_tag
             
527    
528             if DIM ==2 :             if DIM ==2 :
529                source_vector= [source_vector[0],source_vector[2]]                source_vector= [source_vector[0],source_vector[2]]
530              
531             self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))             self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
532             self.__r.setTaggedValue(self.__source_tag, source_vector)             self.__r.setTaggedValue(self.__source_tag, source_vector)
533    
             
534             self.c33=v_p**2 * rho             self.c33=v_p**2 * rho
535             self.c44=v_s**2 * rho             self.c44=v_s**2 * rho
536             self.c11=(1+2*eps) * self.c33             self.c11=(1+2*eps) * self.c33
537             self.c66=(1+2*gamma) * self.c44             self.c66=(1+2*gamma) * self.c44
538             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
539             self.c23=self.c33-2*self.c66             self.c23=self.c33-2*self.c66
540              
541          def  _getAcceleration(self, t, u):          def  _getAcceleration(self, t, u):
542               """               """
543               returns the acceleraton for time t and solution u at time t               returns the acceleraton for time t and solution u at time t
544               """               """
545               du = grad(u)               du = grad(u)
546               sigma=self.__mypde.getCoefficient('X')               sigma=self.__mypde.getCoefficient('X')
547                
548               if self.__mypde.getDim() == 3:               if self.__mypde.getDim() == 3:
549                  e11=du[0,0]                  e11=du[0,0]
550                  e22=du[1,1]                  e22=du[1,1]
551                  e33=du[2,2]                  e33=du[2,2]
552                
553                  sigma[0,0]=self.c11*e11+self.c13*(e22+e33)                  sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
554                  sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33                  sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
555                  sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33                  sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
556    
557                  s=self.c44*(du[2,1]+du[1,2])                  s=self.c44*(du[2,1]+du[1,2])
558                  sigma[1,2]=s                  sigma[1,2]=s
559                  sigma[2,1]=s                              sigma[2,1]=s
560    
561                  s=self.c66*(du[2,0]+du[0,2])                  s=self.c66*(du[2,0]+du[0,2])
562                  sigma[0,2]=s                  sigma[0,2]=s
563                  sigma[2,0]=s                  sigma[2,0]=s
564                    
565                  s=self.c66*(du[0,1]+du[1,0])                  s=self.c66*(du[0,1]+du[1,0])
566                  sigma[0,1]=s                  sigma[0,1]=s
567                  sigma[1,0]=s                  sigma[1,0]=s
# Line 611  class HTIWave(WaveBase): Line 575  class HTIWave(WaveBase):
575                  s=self.c66*(du[1,0]+du[0,1])                  s=self.c66*(du[1,0]+du[0,1])
576                  sigma[0,1]=s                  sigma[0,1]=s
577                  sigma[1,0]=s                  sigma[1,0]=s
578                
579               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))
580               return self.__mypde.getSolution()                           return self.__mypde.getSolution()
581    

Legend:
Removed from v.4621  
changed lines
  Added in v.4622

  ViewVC Help
Powered by ViewVC 1.1.26