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

revision 4572 by gross, Fri Dec 6 02:22:37 2013 UTC revision 4573 by gross, Mon Dec 9 03:42:50 2013 UTC
22
23  __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave' ]  __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave', 'HTIWave' ]
24
25
26  from math import pi  from math import pi
# Line 48  class Ricker(Wavelet): Line 48  class Ricker(Wavelet):
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                    :note: maximum frequence is about 2 x the dominant frequence.
53                  """                  """
54                  drop=18                  drop=18
55                  self.__f=f_dom                  self.__f=f_dom
# Line 274  class WaveBase(object): Line 276  class WaveBase(object):
276               self.t_last=t               self.t_last=t
277               return t, self.u + self.v * (t-self.t)               return t, self.u + self.v * (t-self.t)
278
279  def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2):  def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
280      """      """
281      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
282      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'
# Line 297  def createAbsorbtionLayerFunction(x, abs Line 299  def createAbsorbtionLayerFunction(x, abs
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)
301          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 )
302          if not DIM-1 == i:          if  top_absorbation or not DIM-1 == i:
303              x_r=(bb[i][1]-absorption_zone)-x_i              x_r=(bb[i][1]-absorption_zone)-x_i
304              m_r=whereNegative(x_r)              m_r=whereNegative(x_r)
305              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 )
# Line 371  class SonicWave(WaveBase): Line 373  class SonicWave(WaveBase):
373  class VTIWave(WaveBase):  class VTIWave(WaveBase):
374          """          """
375          Solving the VTI wave equation          Solving the VTI wave equation
376
377            :note: In case of a two dimensional domain the second spatial dimenion is depth.
378          """          """
379          def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,              def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
380                       dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):                       dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
# Line 379  class VTIWave(WaveBase): Line 383  class VTIWave(WaveBase):
383
384             :param domain: domain of the problem             :param domain: domain of the problem
385             :type domain: `Doamin`                     :type domain: `Doamin`
386             :param v_p: p-velocity field                 :param v_p: vertical p-velocity field
387             :type v_p: `Scalar`                   :type v_p: `Scalar`
388               :param v_s: vertical s-velocity field
389               :type v_s: `Scalar`
390             :param wavelet: wavelet to describe the time evolution of source term             :param wavelet: wavelet to describe the time evolution of source term
391             :type wavelet: `Wavelet`                       :type wavelet: `Wavelet`
392             :param source_tag: tag of the source location             :param source_tag: tag of the source location
393             :type source_tag: 'str' or 'int'                       :type source_tag: 'str' or 'int'
394               :param source_vector: source orientation vector
395               :param eps: first Thompsen parameter
396               :param delta: second Thompsen parameter
397               :param gamma: third Thompsen parameter
398               :param rho: density
399             :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.
400             :param p0: initial solution. If not present zero is used.                       :param u0: initial solution. If not present zero is used.
401             :param p0_t: initial solution change rate. If not present zero is used.                       :param v0: initial solution change rate. If not present zero is used.
402             :param absorption_zone: thickness of absorption zone                       :param absorption_zone: thickness of absorption zone
403             :param absorption_cut: boundary value of absorption decay factor             :param absorption_cut: boundary value of absorption decay factor
404             :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.
# Line 455  class VTIWave(WaveBase): Line 466  class VTIWave(WaveBase):
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
s=self.c66*(du[0,1]+du[1,0])
sigma[0,1]=s
sigma[1,0]=s

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])
478            sigma[0,1]=s
479            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]
# Line 480  class VTIWave(WaveBase): Line 493  class VTIWave(WaveBase):
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()
495
496    class HTIWave(WaveBase):
497            """
498            Solving the HTI wave equation (along the x_0 axis)
499
500            :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
501            """
502            def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
503                         dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
504               """
505               initialize the VTI wave solver
506
507               :param domain: domain of the problem
508               :type domain: `Doamin`
509               :param v_p: vertical p-velocity field
510               :type v_p: `Scalar`
511               :param v_s: vertical s-velocity field
512               :type v_s: `Scalar`
513               :param wavelet: wavelet to describe the time evolution of source term
514               :type wavelet: `Wavelet`
515               :param source_tag: tag of the source location
516               :type source_tag: 'str' or 'int'
517               :param source_vector: source orientation vector
518               :param eps: first Thompsen parameter
519               :param delta: second Thompsen parameter
520               :param gamma: third Thompsen parameter
521               :param rho: density
522               :param dt: time step size. If not present a suitable time step size is calculated.
523               :param u0: initial solution. If not present zero is used.
524               :param v0: initial solution change rate. If not present zero is used.
525               :param absorption_zone: thickness of absorption zone
526               :param absorption_cut: boundary value of absorption decay factor
527               :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
528               """
529               DIM=domain.getDim()
530               f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut,  top_absorbation=True)
531
532               v_p=v_p*f
533               v_s=v_s*f
534
535               if u0 == None:
536                  u0=Vector(0.,Solution(domain))
537               else:
538                  u0=interpolate(p0, Solution(domain ))
539
540               if v0 == None:
541                  v0=Vector(0.,Solution(domain))
542               else:
543                  v0=interpolate(v0, Solution(domain ))
544
545               if dt == None:
546                      dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
547
548               super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
549
550               self.__wavelet=wavelet
551
552               self.__mypde=LinearPDESystem(domain)
553               if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
554               self.__mypde.setSymmetryOn()
555               self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
556               self.__source_tag=source_tag
557
558
559               if DIM ==2 :
560              source_vector= [source_vector[0],source_vector[1]]
561
562
563           self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
564               self.__r.setTaggedValue(self.__source_tag, source_vector)
565
566
567               self.c33=v_p**2 * rho
568           self.c44=v_s**2 * rho
569           self.c11=(1+2*eps) * self.c33
570           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
572           self.c23=self.c33-2*self.c66
573
574            def  _getAcceleration(self, t, u):
575                 """
576                 returns the acceleraton for time t and solution u at time t
577                 """
579                 sigma=self.__mypde.getCoefficient('X')
580
581                 if self.__mypde.getDim() == 3:
582            e11=du[0,0]
583            e22=du[1,1]
584            e33=du[2,2]
585
586            sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
587            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
589
590            s=self.c44*(du[2,1]+du[1,2])
591            sigma[1,2]=s
592            sigma[2,1]=s
593
594            s=self.c66*(du[2,0]+du[0,2])
595            sigma[0,2]=s
596            sigma[2,0]=s
597
598            s=self.c66*(du[0,1]+du[1,0])
599            sigma[0,1]=s
600            sigma[1,0]=s
601
602                 else:
603            e11=du[0,0]
604            e22=du[1,1]
605
606            sigma[0,0]=self.c11*e11+self.c13*e22
607            sigma[1,1]=self.c13*e11+self.c33*e22
608
609            s=self.c66*(du[0,1]+du[1,0])
610            sigma[0,1]=s
611            sigma[1,0]=s
612
613                 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
614                 return self.__mypde.getSolution()
615

Legend:
 Removed from v.4572 changed lines Added in v.4573