/[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 4561 by sshaw, Tue Dec 3 05:40:53 2013 UTC revision 4562 by gross, Wed Dec 4 05:41:10 2013 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' ]  __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave' ]
24    
25    
26  from math import pi  from math import pi
# Line 29  import sys Line 29  import sys
29  import time  import time
30  from esys.escript import *  from esys.escript import *
31  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
32  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
33    
34    
35    
# Line 349  class SonicWave(WaveBase): Line 349  class SonicWave(WaveBase):
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.)
350                        
351             self.__wavelet=wavelet             self.__wavelet=wavelet
352             self.__mypde=LinearPDE(domain)             self.__mypde=LinearSinglePDE(domain)
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.)
# Line 367  class SonicWave(WaveBase): Line 367  class SonicWave(WaveBase):
367               return self.__mypde.getSolution()               return self.__mypde.getSolution()
368    
369                            
370    
371    class VTIWave(WaveBase):
372            """
373            Solving the VTI wave equation
374            """
375            def __init__(self, domain, v_p, v_s,   wavelet, source_tag, source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,    
376                         dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
377               """
378               initialize the VTI wave solver
379              
380               :param domain: domain of the problem
381               :type domain: `Doamin`        
382               :param v_p: p-velocity field    
383               :type v_p: `Scalar`      
384               :param wavelet: wavelet to describe the time evolution of source term
385               :type wavelet: `Wavelet`          
386               :param source_tag: tag of the source location
387               :type source_tag: 'str' or 'int'          
388               :param dt: time step size. If not present a suitable time step size is calculated.          
389               :param p0: initial solution. If not present zero is used.          
390               :param p0_t: initial solution change rate. If not present zero is used.          
391               :param absorption_zone: thickness of absorption zone          
392               :param absorption_cut: boundary value of absorption decay factor
393               :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
394               """
395               DIM=domain.getDim()
396               f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
397              
398               v_p=v_p*f
399               v_s=v_s*f
400    
401              
402              
403               if u0 == None:
404                  u0=Vector(0.,Solution(domain))
405               else:
406                  u0=interpolate(p0, Solution(domain ))
407                  
408               if v0 == None:
409                  v0=Vector(0.,Solution(domain))
410               else:
411                  v0=interpolate(v0, Solution(domain ))
412              
413               if dt == None:
414                      dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
415                
416               super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
417              
418               self.__wavelet=wavelet
419              
420               self.__mypde=LinearPDESystem(domain)
421               if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
422               self.__mypde.setSymmetryOn()
423               self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
424               self.__source_tag=source_tag
425              
426    
427               if DIM ==2 :
428              source_vector= [source_vector[0],source_vector[2]]
429    
430          
431           self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
432               self.__r.setTaggedValue(self.__source_tag, source_vector)
433    
434              
435               self.c33=v_p**2 * rho
436           self.c44=v_s**2 * rho
437           self.c11=(1+2*eps) * self.c33
438           self.c66=(1+2*gamma) * self.c33
439           self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
440           self.c12=self.c11-2*self.c66
441          
442            def  _getAcceleration(self, t, u):
443                 """
444                 returns the acceleraton for time t and solution u at time t
445                 """
446                 strain=symmetric(grad(u))
447                 sigma=self.__mypde.getCoefficient('X')
448                
449                 if self.__mypde.getDim() == 3:
450            e11=strain[0,0]
451            e22=strain[1,1]
452            e33=strain[2,2]
453                
454            sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
455            sigma[1,1]=self.c12*e11+self.c11*u22+self.c13*e33
456            sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
457    
458            s=self.c44*(strain[1,0]+strain[0,1])*0.5
459            sigma[0,1]=s
460            sigma[1,0]=s
461    
462            s=self.c44*(strain[2,1]+strain[1,2])*0.5
463            sigma[1,2]=s
464            sigma[2,1]=s
465                
466            s=self.c44*(strain[2,0]+strain[0,2])*0.5
467            sigma[0,2]=s
468            sigma[2,0]=s
469                 else:
470            e11=strain[0,0]
471            e22=strain[1,1]
472                
473            sigma[0,0]=self.c11*e11+self.c13*e22
474            sigma[1,1]=self.c13*e11+self.c33*e22
475                
476            s=self.c44*(strain[1,0]+strain[0,1])*0.5
477            sigma[0,1]=s
478            sigma[1,0]=s
479                
480                 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
481                 return self.__mypde.getSolution()
482    
483                
484    

Legend:
Removed from v.4561  
changed lines
  Added in v.4562

  ViewVC Help
Powered by ViewVC 1.1.26