/[escript]/trunk/downunder/py_src/seismic.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (4 years, 8 months 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 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
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 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
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