/[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 4562 - (show annotations)
Wed Dec 4 05:41:10 2013 UTC (5 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 19600 byte(s)
VTI wave added  
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' ]
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 drop=18
53 self.__f=f_dom
54 self.__f_max=sqrt(7)*f_dom
55 self.__s=pi*self.__f
56 if t_dom == None:
57 t_dom=sqrt(drop)/self.__s
58 self.__t0=t_dom
59
60 def getCenter(self):
61 """
62 return value of wavelet center
63 """
64 return self.__t0
65
66 def getTimeScale(self):
67 """
68 returns the time scale which is the inverse of the largest freqence with a significant
69 spectral component.
70 """
71 return 1/self.__f_max
72
73 def getValue(self, t):
74 """
75 get value of wavelet at time t
76 """
77 e2=(self.__s*(t-self.__t0))**2
78 return (1-2*e2)*exp(-e2)
79
80 def getAcceleration(self, t):
81 """
82 get the acceleration f''(t) at time t
83 """
84 e2=(self.__s*(t-self.__t0))**2
85 return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)
86
87
88 class SimpleSEGYWriter(object):
89 """
90 as simple writer for 2D and 3D seimic lines in particular for synthetic data
91
92 Typical usage:
93
94 from esys.escript import unitsSI as U
95 sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)
96 while n < 10:
97 sw.addRecord([i*2., i*0.67, i**2, -i*7])
98 sw.write('example.segy')
99
100 :note: the writer uses `obsy`
101 """
102 COORDINATE_SCALE = 1000.
103 def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):
104 """
105 initalize writer
106
107 :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
108 are given
109 :param source: coordinates of the source (in meters). For the 2D case a single floats is given, for the 3D case a coordinate tuples
110 :param sampling_interval: sample rate in seconds
111 :param text: a text for the header file (e.g a description)
112 """
113 if isinstance(source, float) or isinstance(source, int) :
114 DIM=1
115 source = (source, 0.)
116 elif hasattr(source, '__len__'):
117 DIM=len(source)
118 if DIM == 1:
119 source = (source[0], 0.)
120 elif DIM==2:
121 source = (source[0], source[1] )
122 else:
123 raise ValueError("Only 1D or 2D arrays accepted.")
124 else:
125 raise TypeError("illegal source type.")
126 if receiver_group== None:
127 if DIM == 1:
128 receiver_group=[0.]
129 else:
130 receiver_group=[(0.,0.)]
131
132 if isinstance(receiver_group, float) or isinstance(receiver_group, int) :
133 if DIM == 1:
134 rg = [ (receiver_group, 0. ) ]
135 else:
136 raise TypeError("illegal receiver_group type.")
137 elif hasattr(receiver_group, '__len__'):
138 if DIM == 1:
139 rg = [(c,0) for c in receiver_group]
140 else:
141 rg = [(c[0],c[1]) for c in receiver_group]
142 else:
143 raise TypeError("illegal receiver_group type.")
144
145 self.__source=source
146 self.__receiver_group=rg
147 self.__text=text
148 self.__trace = [ [] for i in self.__receiver_group ]
149 self.__sampling_interval = sampling_interval
150
151 def addRecord(self, record):
152 """
153 adds a record to the traces. a time difference of sample_interval between two records is assumed.
154 The record mast be a list of as many values as given receivers or a float if a single receiver is used.
155
156 :param record: list of tracks to be added to the record.
157 """
158 if not len(self.__receiver_group) == len(record):
159 raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))
160 if len(self.__receiver_group) == 1:
161 if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
162 self.__trace[0].append(record)
163 else:
164 self.__trace[0].append(record[0])
165 else:
166 for i in xrange(len(record)):
167 self.__trace[i].append(record[i])
168
169 def getSamplingInterval(self):
170 """
171 returns the sampling interval in seconds.
172 """
173 return self.__sampling_interval
174
175 def write(self, filename):
176 """
177 writes to segy file
178
179 :param filename: file name
180 :note: the function uses the `obspy` module.
181 """
182 from obspy import Trace, Stream, UTCDateTime
183 from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
184 from obspy.core import AttribDict
185
186 stream=Stream()
187
188 for i in xrange(len(self.__receiver_group)):
189 trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
190 # Attributes in trace.stats will overwrite everything in
191 # trace.stats.segy.trace_header (in Hz)
192 trace.stats.sampling_rate = 1./self.getSamplingInterval()
193 #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0)
194 if not hasattr(trace.stats, 'segy.trace_header'):
195 trace.stats.segy = {}
196 trace.stats.segy.trace_header = SEGYTraceHeader()
197 trace.stats.segy.trace_header.trace_identification_code=1
198 trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
199 trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = int(self.COORDINATE_SCALE)
200 trace.stats.segy.trace_header.coordinate_units=1
201 trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE)
202 trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE)
203 trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE)
204 trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE)
205
206 # Add trace to stream
207 stream.append(trace)
208 # A SEGY file has file wide headers. This can be attached to the stream
209 # object. If these are not set, they will be autocreated with default
210 # values.
211 stream.stats = AttribDict()
212 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
213 stream.stats.binary_file_header = SEGYBinaryFileHeader()
214
215 if getMPIRankWorld()<1:
216 stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
217
218 class WaveBase(object):
219 """
220 Base for wave propagation using the Verlet scheme.
221
222 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
223
224 with a given acceleration force as function of time.
225
226 a_n=A(t_{n-1})
227 v_n=v_{n-1} + dt * a_n
228 u_n=u_{n-1} + dt * v_n
229
230
231 """
232 def __init__(self, dt, u0, v0, t0=0.):
233 """
234 set up the wave base
235
236 :param dt: time step size (need to be sufficiently small)
237 :param u0: initial value
238 :param v0: initial velocity
239 :param t0: initial time
240 """
241 self.u=u0
242 self.v=v0
243 self.t=t0
244 self.t_last=t0
245 self.__dt=dt
246
247 def getTimeStepSize(self):
248 return self.__dt
249
250 def _getAcceleration(self, t, u):
251 """
252 returns the acceleraton for time t and solution u at time t
253 : note: this function needs to be overwritten by a particular wave problem
254 """
255 pass
256
257 def update(self, t):
258 """
259 returns the solution for the next time marker t which needs to greater than the time marker from the
260 previous call.
261 """
262 if not self.t_last <= t:
263 raise ValueError("You can not move backward in time.")
264
265 dt = self.getTimeStepSize()
266 # apply Verlet scheme until
267 while self.t < t:
268 a=self._getAcceleration(self.t, self.u)
269 self.v += dt * a
270 self.u += dt * self.v
271 self.t += dt
272
273 # now we work backwards
274 self.t_last=t
275 return t, self.u + self.v * (t-self.t)
276
277 def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2):
278 """
279 creating a distribution which is one in the interior of the domain of x
280 and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
281 toward each boundary except the top of the domain.
282
283 :param x: location of points in the domain
284 :type x: `Data`
285 :param absorption_zone: thickness of the aborption zone
286 :param absorption_cut: value of decay function on domain boundary
287 :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
288 toward the boundary.
289 """
290 dom=x.getDomain()
291 bb=boundingBox(dom)
292 DIM=dom.getDim()
293 decay=-log(absorption_cut)/absorption_zone**2
294 f=1
295 for i in xrange(DIM):
296 x_i=x[i]
297 x_l=x_i-(bb[i][0]+absorption_zone)
298 m_l=whereNegative(x_l)
299 f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
300 if not DIM-1 == i:
301 x_r=(bb[i][1]-absorption_zone)-x_i
302 m_r=whereNegative(x_r)
303 f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
304 return f
305
306 class SonicWave(WaveBase):
307 """
308 Solving the sonic wave equation
309
310 p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.
311
312 f(t) is wavelet acting at a point source term at positon s
313 """
314 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):
315 """
316 initialize the sonic wave solver
317
318 :param domain: domain of the problem
319 :type domain: `Doamin`
320 :param v_p: p-velocity field
321 :type v_p: `Scalar`
322 :param wavelet: wavelet to describe the time evolution of source term
323 :type wavelet: `Wavelet`
324 :param source_tag: tag of the source location
325 :type source_tag: 'str' or 'int'
326 :param dt: time step size. If not present a suitable time step size is calculated.
327 :param p0: initial solution. If not present zero is used.
328 :param p0_t: initial solution change rate. If not present zero is used.
329 :param absorption_zone: thickness of absorption zone
330 :param absorption_cut: boundary value of absorption decay factor
331 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
332 """
333 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
334 v_p=v_p*f
335
336 if p0 == None:
337 p0=Scalar(0.,Solution(domain))
338 else:
339 p0=interpolate(p0, Solution(domain ))
340
341 if p0_t == None:
342 p0_t=Scalar(0.,Solution(domain))
343 else:
344 p0_t=interpolate(p0_t, Solution(domain ))
345
346 if dt == None:
347 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
348
349 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
350
351 self.__wavelet=wavelet
352 self.__mypde=LinearSinglePDE(domain)
353 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
354 self.__mypde.setSymmetryOn()
355 self.__mypde.setValue(D=1.)
356 self.__source_tag=source_tag
357 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
358 self.__vp2=v_p**2
359
360
361 def _getAcceleration(self, t, u):
362 """
363 returns the acceleraton for time t and solution u at time t
364 """
365 self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t))
366 self.__mypde.setValue(X=-self.__vp2*grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
367 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

  ViewVC Help
Powered by ViewVC 1.1.26