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

Revision 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 1 month ago) by sshaw
File MIME type: text/x-python
File size: 25496 byte(s)
```python3ified things, replaced mixed whitespace and xrange calls
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2013 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland 17 http://www.uq.edu.au 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave', 'HTIWave' ] 24 25 26 from math import pi 27 import numpy as np 28 import sys 29 import time 30 from esys.escript import * 31 import esys.escript.unitsSI as U 32 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem 33 34 35 36 37 class Wavelet(object): 38 """ 39 place holder for source wavelet 40 """ 41 pass 42 43 class Ricker(Wavelet): 44 """ 45 The Ricker Wavelet w=f(t) 46 """ 47 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 50 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 55 self.__f=f_dom 56 self.__f_max=sqrt(7)*f_dom 57 self.__s=pi*self.__f 58 if t_dom == None: 59 t_dom=sqrt(drop)/self.__s 60 self.__t0=t_dom 61 62 def getCenter(self): 63 """ 64 return value of wavelet center 65 """ 66 return self.__t0 67 68 def getTimeScale(self): 69 """ 70 returns the time scale which is the inverse of the largest freqence with a significant 71 spectral component. 72 """ 73 return 1/self.__f_max 74 75 def getValue(self, t): 76 """ 77 get value of wavelet at time t 78 """ 79 e2=(self.__s*(t-self.__t0))**2 80 return (1-2*e2)*exp(-e2) 81 82 def getAcceleration(self, t): 83 """ 84 get the acceleration f''(t) at time t 85 """ 86 e2=(self.__s*(t-self.__t0))**2 87 return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2) 88 89 90 class SimpleSEGYWriter(object): 91 """ 92 as simple writer for 2D and 3D seimic lines in particular for synthetic data 93 94 Typical usage: 95 96 from esys.escript import unitsSI as U 97 sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec) 98 while n < 10: 99 sw.addRecord([i*2., i*0.67, i**2, -i*7]) 100 sw.write('example.segy') 101 102 :note: the writer uses `obsy` 103 """ 104 COORDINATE_SCALE = 1000. 105 def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"): 106 """ 107 initalize writer 108 109 :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 110 are given 111 :param source: coordinates of the source (in meters). For the 2D case a single floats is given, for the 3D case a coordinate tuples 112 :param sampling_interval: sample rate in seconds 113 :param text: a text for the header file (e.g a description) 114 """ 115 if isinstance(source, float) or isinstance(source, int) : 116 DIM=1 117 source = (source, 0.) 118 elif hasattr(source, '__len__'): 119 DIM=len(source) 120 if DIM == 1: 121 source = (source[0], 0.) 122 elif DIM==2: 123 source = (source[0], source[1] ) 124 else: 125 raise ValueError("Only 1D or 2D arrays accepted.") 126 else: 127 raise TypeError("illegal source type.") 128 if receiver_group== None: 129 if DIM == 1: 130 receiver_group=[0.] 131 else: 132 receiver_group=[(0.,0.)] 133 134 if isinstance(receiver_group, float) or isinstance(receiver_group, int) : 135 if DIM == 1: 136 rg = [ (receiver_group, 0. ) ] 137 else: 138 raise TypeError("illegal receiver_group type.") 139 elif hasattr(receiver_group, '__len__'): 140 if DIM == 1: 141 rg = [(c,0) for c in receiver_group] 142 else: 143 rg = [(c[0],c[1]) for c in receiver_group] 144 else: 145 raise TypeError("illegal receiver_group type.") 146 147 self.__source=source 148 self.__receiver_group=rg 149 self.__text=text 150 self.__trace = [ [] for i in self.__receiver_group ] 151 self.__sampling_interval = sampling_interval 152 153 def addRecord(self, record): 154 """ 155 adds a record to the traces. a time difference of sample_interval between two records is assumed. 156 The record mast be a list of as many values as given receivers or a float if a single receiver is used. 157 158 :param record: list of tracks to be added to the record. 159 """ 160 if not len(self.__receiver_group) == len(record): 161 raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record))) 162 if len(self.__receiver_group) == 1: 163 if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int): 164 self.__trace[0].append(record) 165 else: 166 self.__trace[0].append(record[0]) 167 else: 168 for i in range(len(record)): 169 self.__trace[i].append(record[i]) 170 171 def getSamplingInterval(self): 172 """ 173 returns the sampling interval in seconds. 174 """ 175 return self.__sampling_interval 176 177 def write(self, filename): 178 """ 179 writes to segy file 180 181 :param filename: file name 182 :note: the function uses the `obspy` module. 183 """ 184 from obspy import Trace, Stream, UTCDateTime 185 from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader 186 from obspy.core import AttribDict 187 188 stream=Stream() 189 190 for i in range(len(self.__receiver_group)): 191 trace = Trace(data=np.array(self.__trace[i], dtype='float32')) 192 # Attributes in trace.stats will overwrite everything in 193 # trace.stats.segy.trace_header (in Hz) 194 trace.stats.sampling_rate = 1./self.getSamplingInterval() 195 #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0) 196 if not hasattr(trace.stats, 'segy.trace_header'): 197 trace.stats.segy = {} 198 trace.stats.segy.trace_header = SEGYTraceHeader() 199 trace.stats.segy.trace_header.trace_identification_code=1 200 trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1 201 trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = int(self.COORDINATE_SCALE) 202 trace.stats.segy.trace_header.coordinate_units=1 203 trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE) 204 trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE) 205 trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE) 206 trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE) 207 208 # Add trace to stream 209 stream.append(trace) 210 # A SEGY file has file wide headers. This can be attached to the stream 211 # object. If these are not set, they will be autocreated with default 212 # values. 213 stream.stats = AttribDict() 214 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime()) 215 stream.stats.binary_file_header = SEGYBinaryFileHeader() 216 217 if getMPIRankWorld()<1: 218 stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder) 219 220 class WaveBase(object): 221 """ 222 Base for wave propagation using the Verlet scheme. 223 224 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0 225 226 with a given acceleration force as function of time. 227 228 a_n=A(t_{n-1}) 229 v_n=v_{n-1} + dt * a_n 230 u_n=u_{n-1} + dt * v_n 231 232 233 """ 234 def __init__(self, dt, u0, v0, t0=0.): 235 """ 236 set up the wave base 237 238 :param dt: time step size (need to be sufficiently small) 239 :param u0: initial value 240 :param v0: initial velocity 241 :param t0: initial time 242 """ 243 self.u=u0 244 self.v=v0 245 self.t=t0 246 self.t_last=t0 247 self.__dt=dt 248 249 def getTimeStepSize(self): 250 return self.__dt 251 252 def _getAcceleration(self, t, u): 253 """ 254 returns the acceleraton for time t and solution u at time t 255 : note: this function needs to be overwritten by a particular wave problem 256 """ 257 pass 258 259 def update(self, t): 260 """ 261 returns the solution for the next time marker t which needs to greater than the time marker from the 262 previous call. 263 """ 264 if not self.t_last <= t: 265 raise ValueError("You can not move backward in time.") 266 267 dt = self.getTimeStepSize() 268 # apply Verlet scheme until 269 while self.t < t: 270 a=self._getAcceleration(self.t, self.u) 271 self.v += dt * a 272 self.u += dt * self.v 273 self.t += dt 274 275 # now we work backwards 276 self.t_last=t 277 return t, self.u + self.v * (t-self.t) 278 279 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 282 and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone' 283 toward each boundary except the top of the domain. 284 285 :param x: location of points in the domain 286 :type x: `Data` 287 :param absorption_zone: thickness of the aborption zone 288 :param absorption_cut: value of decay function on domain boundary 289 :return: function on 'x' which is one in the iterior and decays to almost zero over a margin 290 toward the boundary. 291 """ 292 dom=x.getDomain() 293 bb=boundingBox(dom) 294 DIM=dom.getDim() 295 decay=-log(absorption_cut)/absorption_zone**2 296 f=1 297 for i in range(DIM): 298 x_i=x[i] 299 x_l=x_i-(bb[i][0]+absorption_zone) 300 m_l=whereNegative(x_l) 301 f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 ) 302 if top_absorbation or not DIM-1 == i: 303 x_r=(bb[i][1]-absorption_zone)-x_i 304 m_r=whereNegative(x_r) 305 f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 ) 306 return f 307 308 class SonicWave(WaveBase): 309 """ 310 Solving the sonic wave equation 311 312 p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p. 313 314 f(t) is wavelet acting at a point source term at positon s 315 """ 316 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): 317 """ 318 initialize the sonic wave solver 319 320 :param domain: domain of the problem 321 :type domain: `Doamin` 322 :param v_p: p-velocity field 323 :type v_p: `Scalar` 324 :param wavelet: wavelet to describe the time evolution of source term 325 :type wavelet: `Wavelet` 326 :param source_tag: tag of the source location 327 :type source_tag: 'str' or 'int' 328 :param dt: time step size. If not present a suitable time step size is calculated. 329 :param p0: initial solution. If not present zero is used. 330 :param p0_t: initial solution change rate. If not present zero is used. 331 :param absorption_zone: thickness of absorption zone 332 :param absorption_cut: boundary value of absorption decay factor 333 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion. 334 """ 335 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut) 336 v_p=v_p*f 337 338 if p0 == None: 339 p0=Scalar(0.,Solution(domain)) 340 else: 341 p0=interpolate(p0, Solution(domain )) 342 343 if p0_t == None: 344 p0_t=Scalar(0.,Solution(domain)) 345 else: 346 p0_t=interpolate(p0_t, Solution(domain )) 347 348 if dt == None: 349 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale()) 350 351 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.) 352 353 self.__wavelet=wavelet 354 self.__mypde=LinearSinglePDE(domain) 355 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING) 356 self.__mypde.setSymmetryOn() 357 self.__mypde.setValue(D=1.) 358 self.__source_tag=source_tag 359 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain())) 360 self.__vp2=v_p**2 361 362 363 def _getAcceleration(self, t, u): 364 """ 365 returns the acceleraton for time t and solution u at time t 366 """ 367 self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t)) 368 self.__mypde.setValue(X=-self.__vp2*grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r) 369 return self.__mypde.getSolution() 370 371 372 373 class VTIWave(WaveBase): 374 """ 375 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., 380 dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True): 381 """ 382 initialize the VTI wave solver 383 384 :param domain: domain of the problem 385 :type domain: `Doamin` 386 :param v_p: vertical p-velocity field 387 :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 391 :type wavelet: `Wavelet` 392 :param source_tag: tag of the source location 393 :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. 400 :param u0: initial solution. If not present zero is used. 401 :param v0: initial solution change rate. If not present zero is used. 402 :param absorption_zone: thickness of absorption zone 403 :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. 405 """ 406 DIM=domain.getDim() 407 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut) 408 409 v_p=v_p*f 410 v_s=v_s*f 411 412 413 414 if u0 == None: 415 u0=Vector(0.,Solution(domain)) 416 else: 417 u0=interpolate(p0, Solution(domain )) 418 419 if v0 == None: 420 v0=Vector(0.,Solution(domain)) 421 else: 422 v0=interpolate(v0, Solution(domain )) 423 424 if dt == None: 425 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale()) 426 427 super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.) 428 429 self.__wavelet=wavelet 430 431 self.__mypde=LinearPDESystem(domain) 432 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING) 433 self.__mypde.setSymmetryOn() 434 self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X')) 435 self.__source_tag=source_tag 436 437 438 if DIM ==2 : 439 source_vector= [source_vector[0],source_vector[2]] 440 441 442 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain())) 443 self.__r.setTaggedValue(self.__source_tag, source_vector) 444 445 446 self.c33=v_p**2 * rho 447 self.c44=v_s**2 * rho 448 self.c11=(1+2*eps) * self.c33 449 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 451 self.c12=self.c11-2*self.c66 452 453 def _getAcceleration(self, t, u): 454 """ 455 returns the acceleraton for time t and solution u at time t 456 """ 457 du = grad(u) 458 sigma=self.__mypde.getCoefficient('X') 459 460 if self.__mypde.getDim() == 3: 461 e11=du[0,0] 462 e22=du[1,1] 463 e33=du[2,2] 464 465 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 467 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33 468 469 s=self.c44*(du[2,1]+du[1,2]) 470 sigma[1,2]=s 471 sigma[2,1]=s 472 473 s=self.c44*(du[2,0]+du[0,2]) 474 sigma[0,2]=s 475 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: 483 e11=du[0,0] 484 e22=du[1,1] 485 486 sigma[0,0]=self.c11*e11+self.c13*e22 487 sigma[1,1]=self.c13*e11+self.c33*e22 488 489 s=self.c44*(du[1,0]+du[0,1]) 490 sigma[0,1]=s 491 sigma[1,0]=s 492 493 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t)) 494 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 """ 578 du = grad(u) 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

 ViewVC Help Powered by ViewVC 1.1.26