/[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 4691 - (show annotations)
Wed Feb 19 05:23:17 2014 UTC (5 years, 7 months ago) by sshaw
File MIME type: text/x-python
File size: 35122 byte(s)
added support for HTI waves to the fast wave assembler
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 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 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave', 'HTIWave', 'createAbsorbtionLayerFunction', 'SonicHTIWave' , "TTIWave"]
25
26
27 from math import pi
28 import numpy as np
29 import sys
30 import time
31 from esys.escript import *
32 import esys.escript.unitsSI as U
33 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem, WavePDE
34
35 class Wavelet(object):
36 """
37 place holder for source wavelet
38 """
39 pass
40
41 class Ricker(Wavelet):
42 """
43 The Ricker Wavelet w=f(t)
44 """
45 def __init__(self, f_dom=40, t_dom=None):
46 """
47 set up Ricker wavelet wih dominant frequence f_dom and center at time t_dom. If t_dom
48 is not given an estimate for suitable t_dom is calculated so f(0)~0.
49
50 :note: maximum frequence is about 2 x the dominant frequence.
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 class SimpleSEGYWriter(object):
88 """
89 as simple writer for 2D and 3D seimic lines in particular for synthetic data
90
91 Typical usage:
92
93 from esys.escript import unitsSI as U
94 sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)
95 while n < 10:
96 sw.addRecord([i*2., i*0.67, i**2, -i*7])
97 sw.write('example.segy')
98
99 :note: the writer uses `obspy`
100 """
101 COORDINATE_SCALE = 1000.
102 def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):
103 """
104 initalize writer
105
106 :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
107 are given
108 :param source: coordinates of the source (in meters). For the 2D case a single floats is given, for the 3D case a coordinate tuples
109 :param sampling_interval: sample rate in seconds
110 :param text: a text for the header file (e.g a description)
111 """
112 if isinstance(source, float) or isinstance(source, int) :
113 DIM=1
114 source = (source, 0.)
115 elif hasattr(source, '__len__'):
116 DIM=len(source)
117 if DIM == 1:
118 source = (source[0], 0.)
119 elif DIM==2:
120 source = (source[0], source[1] )
121 else:
122 raise ValueError("Only 1D or 2D arrays accepted.")
123 else:
124 raise TypeError("illegal source type.")
125 if receiver_group== None:
126 if DIM == 1:
127 receiver_group=[0.]
128 else:
129 receiver_group=[(0.,0.)]
130
131 if isinstance(receiver_group, float) or isinstance(receiver_group, int) :
132 if DIM == 1:
133 rg = [ (receiver_group, 0. ) ]
134 else:
135 raise TypeError("illegal receiver_group type.")
136 elif hasattr(receiver_group, '__len__'):
137 if DIM == 1:
138 rg = [(c,0) for c in receiver_group]
139 else:
140 rg = [(c[0],c[1]) for c in receiver_group]
141 else:
142 raise TypeError("illegal receiver_group type.")
143
144 self.__source=source
145 self.__receiver_group=rg
146 self.__text=text
147 self.__trace = [ [] for i in self.__receiver_group ]
148 self.__sampling_interval = sampling_interval
149
150 def addRecord(self, record):
151 """
152 adds a record to the traces. a time difference of sample_interval between two records is assumed.
153 The record mast be a list of as many values as given receivers or a float if a single receiver is used.
154
155 :param record: list of tracks to be added to the record.
156 """
157 if not len(self.__receiver_group) == len(record):
158 raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))
159 if len(self.__receiver_group) == 1:
160 if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
161 self.__trace[0].append(record)
162 else:
163 self.__trace[0].append(record[0])
164 else:
165 for i in range(len(record)):
166 self.__trace[i].append(record[i])
167
168 def getSamplingInterval(self):
169 """
170 returns the sampling interval in seconds.
171 """
172 return self.__sampling_interval
173
174 def write(self, filename):
175 """
176 writes to segy file
177
178 :param filename: file name
179 :note: the function uses the `obspy` module.
180 """
181 try:
182 from obspy import Trace, Stream, UTCDateTime
183 from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
184 from obspy.core import AttribDict
185 except ImportError as e:
186 raise RuntimeError("This feature (SimpleSEGYWriter.write())"+\
187 " depends on obspy, which is not installed, see "+\
188 "https://github.com/obspy/obspy for install guide")
189
190 stream=Stream()
191
192 for i in range(len(self.__receiver_group)):
193 trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
194 # Attributes in trace.stats will overwrite everything in
195 # trace.stats.segy.trace_header (in Hz)
196 trace.stats.sampling_rate = 1./self.getSamplingInterval()
197 #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0)
198 if not hasattr(trace.stats, 'segy.trace_header'):
199 trace.stats.segy = {}
200 trace.stats.segy.trace_header = SEGYTraceHeader()
201 trace.stats.segy.trace_header.trace_identification_code=1
202 trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
203 trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = -int(self.COORDINATE_SCALE)
204 trace.stats.segy.trace_header.coordinate_units=1
205 trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE)
206 trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE)
207 trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE)
208 trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE)
209
210 # Add trace to stream
211 stream.append(trace)
212 # A SEGY file has file wide headers. This can be attached to the stream
213 # object. If these are not set, they will be autocreated with default
214 # values.
215 stream.stats = AttribDict()
216 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
217 stream.stats.binary_file_header = SEGYBinaryFileHeader()
218
219 if getMPIRankWorld()<1:
220 stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
221
222 class WaveBase(object):
223 """
224 Base for wave propagation using the Verlet scheme.
225
226 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
227
228 with a given acceleration force as function of time.
229
230 a_n=A(t_{n-1})
231 v_n=v_{n-1} + dt * a_n
232 u_n=u_{n-1} + dt * v_n
233
234
235 """
236 def __init__(self, dt, u0, v0, t0=0.):
237 """
238 set up the wave base
239
240 :param dt: time step size (need to be sufficiently small)
241 :param u0: initial value
242 :param v0: initial velocity
243 :param t0: initial time
244 """
245 self.u=u0
246 self.v=v0
247 self.t=t0
248 self.n=0
249 self.t_last=t0
250 self.__dt=dt
251
252 def getTimeStepSize(self):
253 return self.__dt
254
255 def _getAcceleration(self, t, u):
256 """
257 returns the acceleraton for time t and solution u at time t
258 : note: this function needs to be overwritten by a particular wave problem
259 """
260 pass
261
262 def update(self, t):
263 """
264 returns the solution for the next time marker t which needs to greater than the time marker from the
265 previous call.
266 """
267 if not self.t_last <= t:
268 raise ValueError("You can not move backward in time.")
269
270 dt = self.getTimeStepSize()
271 # apply Verlet scheme until
272 while self.t < t:
273 a=self._getAcceleration(self.t, self.u)
274 self.v += dt * a
275 self.u += dt * self.v
276 self.t += dt
277 self.n += 1
278
279 # now we work backwards
280 self.t_last=t
281 return t, self.u + self.v * (t-self.t)
282
283 def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
284 """
285 creating a distribution which is one in the interior of the domain of x
286 and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
287 toward each boundary except the top of the domain.
288
289 :param x: location of points in the domain
290 :type x: `Data`
291 :param absorption_zone: thickness of the aborption zone
292 :param absorption_cut: value of decay function on domain boundary
293 :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
294 toward the boundary.
295 """
296 dom=x.getDomain()
297 bb=boundingBox(dom)
298 DIM=dom.getDim()
299 decay=-log(absorption_cut)/absorption_zone**2
300 f=1
301 for i in range(DIM):
302 x_i=x[i]
303 x_l=x_i-(bb[i][0]+absorption_zone)
304 m_l=whereNegative(x_l)
305 f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
306 if top_absorbation or not DIM-1 == i:
307 x_r=(bb[i][1]-absorption_zone)-x_i
308 m_r=whereNegative(x_r)
309 f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
310 return f
311
312 class SonicWave(WaveBase):
313 """
314 Solving the sonic wave equation
315
316 p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.
317
318 f(t) is wavelet acting at a point source term at positon s
319 """
320 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):
321 """
322 initialize the sonic wave solver
323
324 :param domain: domain of the problem
325 :type domain: `Domain`
326 :param v_p: p-velocity field
327 :type v_p: `Scalar`
328 :param wavelet: wavelet to describe the time evolution of source term
329 :type wavelet: `Wavelet`
330 :param source_tag: tag of the source location
331 :type source_tag: 'str' or 'int'
332 :param dt: time step size. If not present a suitable time step size is calculated.
333 :param p0: initial solution. If not present zero is used.
334 :param p0_t: initial solution change rate. If not present zero is used.
335 :param absorption_zone: thickness of absorption zone
336 :param absorption_cut: boundary value of absorption decay factor
337 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
338 """
339 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
340 v_p=v_p*f
341
342 if p0 == None:
343 p0=Scalar(0.,Solution(domain))
344 else:
345 p0=interpolate(p0, Solution(domain ))
346
347 if p0_t == None:
348 p0_t=Scalar(0.,Solution(domain))
349 else:
350 p0_t=interpolate(p0_t, Solution(domain ))
351
352 if dt == None:
353 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
354
355 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
356
357 self.__wavelet=wavelet
358 self.__mypde=LinearSinglePDE(domain)
359 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
360 self.__mypde.setSymmetryOn()
361 self.__mypde.setValue(D=1./v_p**2)
362 self.__source_tag=source_tag
363 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
364
365
366 def _getAcceleration(self, t, u):
367 """
368 returns the acceleraton for time t and solution u at time t
369 """
370 self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t))
371 self.__mypde.setValue(X=-grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
372 return self.__mypde.getSolution()
373
374
375 class VTIWave(WaveBase):
376 """
377 Solving the VTI wave equation
378
379 :note: In case of a two dimensional domain the second spatial dimenion is depth.
380 """
381 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
382 source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
383 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
384 absorption_cut=1e-2, lumping=True):
385 """
386 initialize the VTI wave solver
387
388 :param domain: domain of the problem
389 :type domain: `Domain`
390 :param v_p: vertical p-velocity field
391 :type v_p: `Scalar`
392 :param v_s: vertical s-velocity field
393 :type v_s: `Scalar`
394 :param wavelet: wavelet to describe the time evolution of source term
395 :type wavelet: `Wavelet`
396 :param source_tag: tag of the source location
397 :type source_tag: 'str' or 'int'
398 :param source_vector: source orientation vector
399 :param eps: first Thompsen parameter
400 :param delta: second Thompsen parameter
401 :param gamma: third Thompsen parameter
402 :param rho: density
403 :param dt: time step size. If not present a suitable time step size is calculated.
404 :param u0: initial solution. If not present zero is used.
405 :param v0: initial solution change rate. If not present zero is used.
406 :param absorption_zone: thickness of absorption zone
407 :param absorption_cut: boundary value of absorption decay factor
408 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
409 """
410 DIM=domain.getDim()
411 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
412
413 v_p=v_p*f
414 v_s=v_s*f
415
416
417
418 if u0 == None:
419 u0=Vector(0.,Solution(domain))
420 else:
421 u0=interpolate(p0, Solution(domain ))
422
423 if v0 == None:
424 v0=Vector(0.,Solution(domain))
425 else:
426 v0=interpolate(v0, Solution(domain ))
427
428 if dt == None:
429 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
430
431 super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
432
433 self.__wavelet=wavelet
434
435 self.fastAssembler = hasattr(domain, "setAssembler")
436 self.c33=v_p**2 * rho
437 self.c44=v_s**2 * rho
438 self.c11=(1+2*eps) * self.c33
439 self.c66=(1+2*gamma) * self.c44
440 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
441 self.c12=self.c11-2*self.c66
442
443 if self.fastAssembler:
444 self.__mypde=WavePDE(domain, [("c11", self.c11),
445 ("c12", self.c12), ("c13", self.c13), ("c33", self.c33),
446 ("c44", self.c44), ("c66", self.c66)])
447 else:
448 self.__mypde=LinearPDESystem(domain)
449 self.__mypde.setValue(X=self.__mypde.createCoefficient('X'))
450
451 if lumping:
452 self.__mypde.getSolverOptions().setSolverMethod(
453 self.__mypde.getSolverOptions().HRZ_LUMPING)
454 self.__mypde.setSymmetryOn()
455 self.__mypde.setValue(D=rho*kronecker(DIM))
456 self.__source_tag=source_tag
457
458 if DIM ==2 :
459 source_vector= [source_vector[0],source_vector[2]]
460
461 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
462 self.__r.setTaggedValue(self.__source_tag, source_vector)
463
464 def _getAcceleration(self, t, u):
465 """
466 returns the acceleraton for time t and solution u at time t
467 """
468 du = grad(u)
469 if not self.fastAssembler:
470 sigma=self.__mypde.getCoefficient('X')
471 if self.__mypde.getDim() == 3:
472 e11=du[0,0]
473 e22=du[1,1]
474 e33=du[2,2]
475 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
476 sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
477 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
478
479 s=self.c44*(du[2,1]+du[1,2])
480 sigma[1,2]=s
481 sigma[2,1]=s
482
483 s=self.c44*(du[2,0]+du[0,2])
484 sigma[0,2]=s
485 sigma[2,0]=s
486
487 s=self.c66*(du[0,1]+du[1,0])
488 sigma[0,1]=s
489 sigma[1,0]=s
490
491
492 else:
493 e11=du[0,0]
494 e22=du[1,1]
495 sigma[0,0]=self.c11*e11+self.c13*e22
496 sigma[1,1]=self.c13*e11+self.c33*e22
497 s=self.c44*(du[1,0]+du[0,1])
498 sigma[0,1]=s
499 sigma[1,0]=s
500 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
501
502 else:
503 self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
504
505 return self.__mypde.getSolution()
506
507
508 class HTIWave(WaveBase):
509 """
510 Solving the HTI wave equation (along the x_0 axis)
511
512 :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
513 """
514
515 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
516 source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
517 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
518 absorption_cut=1e-2, lumping=True):
519 """
520 initialize the VTI wave solver
521
522 :param domain: domain of the problem
523 :type domain: `Domain`
524 :param v_p: vertical p-velocity field
525 :type v_p: `Scalar`
526 :param v_s: vertical s-velocity field
527 :type v_s: `Scalar`
528 :param wavelet: wavelet to describe the time evolution of source term
529 :type wavelet: `Wavelet`
530 :param source_tag: tag of the source location
531 :type source_tag: 'str' or 'int'
532 :param source_vector: source orientation vector
533 :param eps: first Thompsen parameter
534 :param delta: second Thompsen parameter
535 :param gamma: third Thompsen parameter
536 :param rho: density
537 :param dt: time step size. If not present a suitable time step size is calculated.
538 :param u0: initial solution. If not present zero is used.
539 :param v0: initial solution change rate. If not present zero is used.
540 :param absorption_zone: thickness of absorption zone
541 :param absorption_cut: boundary value of absorption decay factor
542 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
543 """
544 DIM=domain.getDim()
545 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
546
547 v_p=v_p*f
548 v_s=v_s*f
549
550 if u0 == None:
551 u0=Vector(0.,Solution(domain))
552 else:
553 u0=interpolate(p0, Solution(domain ))
554
555 if v0 == None:
556 v0=Vector(0.,Solution(domain))
557 else:
558 v0=interpolate(v0, Solution(domain ))
559
560 if dt == None:
561 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
562
563 super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
564
565 self.__wavelet=wavelet
566
567 self.fastAssembler = hasattr(domain, "setAssembler")
568 self.c33=v_p**2 * rho
569 self.c44=v_s**2 * rho
570 self.c11=(1+2*eps) * self.c33
571 self.c66=(1+2*gamma) * self.c44
572 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
573 self.c23=self.c33-2*self.c66
574
575 if self.fastAssembler:
576 self.__mypde=WavePDE(domain, [("c11", self.c11),
577 ("c23", self.c23), ("c13", self.c13), ("c33", self.c33),
578 ("c44", self.c44), ("c66", self.c66)])
579 else:
580 self.__mypde=LinearPDESystem(domain)
581 self.__mypde.setValue(X=self.__mypde.createCoefficient('X'))
582
583 if lumping:
584 self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
585 self.__mypde.setSymmetryOn()
586 self.__mypde.setValue(D=rho*kronecker(DIM))
587 self.__source_tag=source_tag
588
589 if DIM ==2 :
590 source_vector= [source_vector[0],source_vector[2]]
591
592 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
593 self.__r.setTaggedValue(self.__source_tag, source_vector)
594
595
596
597 def _getAcceleration(self, t, u):
598 """
599 returns the acceleraton for time t and solution u at time t
600 """
601 du = grad(u)
602 if self.fastAssembler:
603 self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
604 else:
605 sigma=self.__mypde.getCoefficient('X')
606
607 if self.__mypde.getDim() == 3:
608 e11=du[0,0]
609 e22=du[1,1]
610 e33=du[2,2]
611
612 sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
613 sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
614 sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
615
616 s=self.c44*(du[2,1]+du[1,2])
617 sigma[1,2]=s
618 sigma[2,1]=s
619
620 s=self.c66*(du[2,0]+du[0,2])
621 sigma[0,2]=s
622 sigma[2,0]=s
623
624 s=self.c66*(du[0,1]+du[1,0])
625 sigma[0,1]=s
626 sigma[1,0]=s
627
628 else:
629 e11=du[0,0]
630 e22=du[1,1]
631 sigma[0,0]=self.c11*e11+self.c13*e22
632 sigma[1,1]=self.c13*e11+self.c33*e22
633
634 s=self.c66*(du[1,0]+du[0,1])
635 sigma[0,1]=s
636 sigma[1,0]=s
637 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
638
639 return self.__mypde.getSolution()
640
641 class TTIWave(WaveBase):
642 """
643 Solving the 2D TTI wave equation with
644
645 sigma_xx= c11*e_xx + c13*e_zz + c15*e_xz
646 sigma_zz= c13*e_xx + c33*e_zz + c35*e_xz
647 sigma_xz= c15*e_xx + c35*e_zz + c55*e_xz
648
649 the coefficicients c11, c13, etc are calculated from the tompsen parameters eps and delta and the tilt theta
650
651 :note: currently only the 2D case is supported.
652 """
653
654 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
655 source_vector = [0.,1.], eps=0., delta=0., theta=0., rho=1.,
656 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
657 absorption_cut=1e-2, lumping=True):
658 """
659 initialize the TTI wave solver
660
661 :param domain: domain of the problem
662 :type domain: `Domain`
663 :param v_p: vertical p-velocity field
664 :type v_p: `Scalar`
665 :param v_s: vertical s-velocity field
666 :type v_s: `Scalar`
667 :param wavelet: wavelet to describe the time evolution of source term
668 :type wavelet: `Wavelet`
669 :param source_tag: tag of the source location
670 :type source_tag: 'str' or 'int'
671 :param source_vector: source orientation vector
672 :param eps: first Thompsen parameter
673 :param delta: second Thompsen parameter
674 :param theta: tilting (in Rad)
675 :param rho: density
676 :param dt: time step size. If not present a suitable time step size is calculated.
677 :param u0: initial solution. If not present zero is used.
678 :param v0: initial solution change rate. If not present zero is used.
679 :param absorption_zone: thickness of absorption zone
680 :param absorption_cut: boundary value of absorption decay factor
681 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
682 """
683 DIM=domain.getDim()
684 if not DIM == 2:
685 raise ValueError("Only 2D is supported.")
686 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
687
688 v_p=v_p*f
689 v_s=v_s*f
690
691 if u0 == None:
692 u0=Vector(0.,Solution(domain))
693 else:
694 u0=interpolate(p0, Solution(domain ))
695
696 if v0 == None:
697 v0=Vector(0.,Solution(domain))
698 else:
699 v0=interpolate(v0, Solution(domain ))
700
701 if dt == None:
702 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
703
704 super(TTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
705
706 self.__wavelet=wavelet
707
708 self.__mypde=LinearPDESystem(domain)
709 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
710 self.__mypde.setSymmetryOn()
711 self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
712 self.__source_tag=source_tag
713
714 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
715 self.__r.setTaggedValue(self.__source_tag, source_vector)
716
717 c0_33=v_p**2 * rho
718 c0_66=v_s**2 * rho
719 c0_11=(1+2*eps) * c0_33
720 c0_13=sqrt(2*c0_33*(c0_33-c0_66) * delta + (c0_33-c0_66)**2)-c0_66
721
722 self.c11= c0_11*cos(theta)**4 - 2*c0_13*cos(theta)**4 + 2*c0_13*cos(theta)**2 + c0_33*sin(theta)**4 - 4*c0_66*cos(theta)**4 + 4*c0_66*cos(theta)**2
723 self.c13= -c0_11*cos(theta)**4 + c0_11*cos(theta)**2 + c0_13*sin(theta)**4 + c0_13*cos(theta)**4 - c0_33*cos(theta)**4 + c0_33*cos(theta)**2 + 4*c0_66*cos(theta)**4 - 4*c0_66*cos(theta)**2
724 self.c16= (-2*c0_11*cos(theta)**2 - 4*c0_13*sin(theta)**2 + 2*c0_13 + 2*c0_33*sin(theta)**2 - 8*c0_66*sin(theta)**2 + 4*c0_66)*sin(theta)*cos(theta)/2
725 self.c33= c0_11*sin(theta)**4 - 2*c0_13*cos(theta)**4 + 2*c0_13*cos(theta)**2 + c0_33*cos(theta)**4 - 4*c0_66*cos(theta)**4 + 4*c0_66*cos(theta)**2
726 self.c36= (2*c0_11*cos(theta)**2 - 2*c0_11 + 4*c0_13*sin(theta)**2 - 2*c0_13 + 2*c0_33*cos(theta)**2 + 8*c0_66*sin(theta)**2 - 4*c0_66)*sin(theta)*cos(theta)/2
727 self.c66= -c0_11*cos(theta)**4 + c0_11*cos(theta)**2 + 2*c0_13*cos(theta)**4 - 2*c0_13*cos(theta)**2 - c0_33*cos(theta)**4 + c0_33*cos(theta)**2 + c0_66*sin(theta)**4 + 3*c0_66*cos(theta)**4 - 2*c0_66*cos(theta)**2
728
729 def _getAcceleration(self, t, u):
730 """
731 returns the acceleraton for time t and solution u at time t
732 """
733 du = grad(u)
734 sigma=self.__mypde.getCoefficient('X')
735
736 e_xx=du[0,0]
737 e_zz=du[1,1]
738 e_xz=du[0,1]+du[1,0]
739
740
741 sigma[0,0]= self.c11 * e_xx + self.c13 * e_zz + self.c16 * e_xz
742 sigma[1,1]= self.c13 * e_xx + self.c33 * e_zz + self.c36 * e_xz
743 sigma_xz = self.c16 * e_xx + self.c36 * e_zz + self.c66 * e_xz
744
745 sigma[0,1]=sigma_xz
746 sigma[1,0]=sigma_xz
747
748 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
749 return self.__mypde.getSolution()
750
751 class SonicHTIWave(WaveBase):
752 """
753 Solving the HTI wave equation (along the x_0 axis) with azimuth (rotation around verticle axis)
754 under the assumption of zero shear wave velocities
755 The unknowns are the transversal (along x_0) and vertial stress (Q, P)
756
757 :note: In case of a two dimensional domain the second spatial dimenion is depth.
758 """
759 def __init__(self, domain, v_p, wavelet, source_tag, source_vector = [1.,0.], eps=0., delta=0., azimuth=0.,
760 dt=None, p0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
761 """
762 initialize the HTI wave solver
763
764 :param domain: domain of the problem
765 :type domain: `Doamin`
766 :param v_p: vertical p-velocity field
767 :type v_p: `Scalar`
768 :param v_s: vertical s-velocity field
769 :type v_s: `Scalar`
770 :param wavelet: wavelet to describe the time evolution of source term
771 :type wavelet: `Wavelet`
772 :param source_tag: tag of the source location
773 :type source_tag: 'str' or 'int'
774 :param source_vector: source orientation vector
775 :param eps: first Thompsen parameter
776 :param azimuth: azimuth (rotation around verticle axis)
777 :param gamma: third Thompsen parameter
778 :param rho: density
779 :param dt: time step size. If not present a suitable time step size is calculated.
780 :param p0: initial solution (Q(t=0), P(t=0)). If not present zero is used.
781 :param v0: initial solution change rate. If not present zero is used.
782 :param absorption_zone: thickness of absorption zone
783 :param absorption_cut: boundary value of absorption decay factor
784 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
785 """
786 DIM=domain.getDim()
787 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
788
789 self.v2_p=v_p**2
790 self.v2_t=self.v2_p*sqrt(1+2*delta)
791 self.v2_n=self.v2_p*(1+2*eps)
792
793 if p0 == None:
794 p0=Data(0.,(2,),Solution(domain))
795 else:
796 p0=interpolate(p0, Solution(domain ))
797
798 if v0 == None:
799 v0=Data(0.,(2,),Solution(domain))
800 else:
801 v0=interpolate(v0, Solution(domain ))
802
803 if dt == None:
804 dt=min(min(inf(domain.getSize()/sqrt(self.v2_p)), inf(domain.getSize()/sqrt(self.v2_t)), inf(domain.getSize()/sqrt(self.v2_n))) , wavelet.getTimeScale())*0.2
805
806 super(SonicHTIWave, self).__init__( dt, u0=p0, v0=v0, t0=0.)
807
808 self.__wavelet=wavelet
809
810 self.__mypde=LinearPDESystem(domain)
811 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
812 self.__mypde.setSymmetryOn()
813 self.__mypde.setValue(D=kronecker(2), X=self.__mypde.createCoefficient('X'))
814 self.__source_tag=source_tag
815
816
817 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
818 self.__r.setTaggedValue(self.__source_tag, source_vector)
819
820 def _getAcceleration(self, t, u):
821 """
822 returns the acceleraton for time t and solution u at time t
823 """
824 dQ = grad(u[0])[0]
825 dP = grad(u[1])[1:]
826 sigma=self.__mypde.getCoefficient('X')
827
828 sigma[0,0] = self.v2_n*dQ
829 sigma[0,1:] = self.v2_t*dP
830 sigma[1,0] = self.v2_t*dQ
831 sigma[1,1:] = self.v2_p*dP
832
833 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
834 return self.__mypde.getSolution()

  ViewVC Help
Powered by ViewVC 1.1.26