/[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 5197 - (show annotations)
Wed Oct 8 04:41:15 2014 UTC (4 years, 7 months ago) by sshaw
File MIME type: text/x-python
File size: 35326 byte(s)
allows absorption zones of zero size in seismic models
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, SolverOptions
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 Sets up a Ricker wavelet wih dominant frequence `f_dom` and
48 center at time `t_dom`. If `t_dom` is not given an estimate
49 for suitable `t_dom` is calculated so f(0)~0.
50
51 :note: maximum frequence is about 2 x the dominant frequence.
52 """
53 drop=18
54 self.__f=f_dom
55 self.__f_max=sqrt(7)*f_dom
56 self.__s=pi*self.__f
57 if t_dom == None:
58 t_dom=sqrt(drop)/self.__s
59 self.__t0=t_dom
60
61 def getCenter(self):
62 """
63 Return value of wavelet center
64 """
65 return self.__t0
66
67 def getTimeScale(self):
68 """
69 Returns the time scale which is the inverse of the largest
70 frequence with a significant spectral component.
71 """
72 return 1/self.__f_max
73
74 def getValue(self, t):
75 """
76 get value of wavelet at time `t`
77 """
78 e2=(self.__s*(t-self.__t0))**2
79 return (1-2*e2)*exp(-e2)
80
81 def getAcceleration(self, t):
82 """
83 get the acceleration f''(t) at time `t`
84 """
85 e2=(self.__s*(t-self.__t0))**2
86 return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)
87
88 class SimpleSEGYWriter(object):
89 """
90 A simple writer for 2D and 3D seismic 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
98 ` sw.addRecord([i*2., i*0.67, i**2, -i*7])`
99
100 `sw.write('example.segy')`
101
102 :note: the writer uses `obspy`
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
185 try:
186 from obspy import Trace, Stream, UTCDateTime
187 from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
188 from obspy.core import AttribDict
189 except ImportError as e:
190 raise RuntimeError("This feature (SimpleSEGYWriter.write())"+\
191 " depends on obspy, which is not installed, see "+\
192 "https://github.com/obspy/obspy for install guide")
193
194 if getMPISizeWorld() > 1:
195 raise RuntimeError("Writing segy files with multiple ranks is"+\
196 " not yet supported.")
197
198 stream=Stream()
199
200 for i in range(len(self.__receiver_group)):
201 trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
202 # Attributes in trace.stats will overwrite everything in
203 # trace.stats.segy.trace_header (in Hz)
204 trace.stats.sampling_rate = 1./self.getSamplingInterval()
205 #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0)
206 if not hasattr(trace.stats, 'segy.trace_header'):
207 trace.stats.segy = {}
208 trace.stats.segy.trace_header = SEGYTraceHeader()
209 trace.stats.segy.trace_header.trace_identification_code=1
210 trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
211 trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = -int(self.COORDINATE_SCALE)
212 trace.stats.segy.trace_header.coordinate_units=1
213 trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE)
214 trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE)
215 trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE)
216 trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE)
217
218 # Add trace to stream
219 stream.append(trace)
220 # A SEGY file has file wide headers. This can be attached to the stream
221 # object. If these are not set, they will be autocreated with default
222 # values.
223 stream.stats = AttribDict()
224 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
225 stream.stats.binary_file_header = SEGYBinaryFileHeader()
226
227 if getMPIRankWorld()<1:
228 stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
229
230 class WaveBase(object):
231 """
232 Base for wave propagation using the Verlet scheme.
233
234 ``u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0``
235
236 with a given acceleration force as function of time.
237
238 a_n=A(t_{n-1})
239 v_n=v_{n-1} + dt * a_n
240 u_n=u_{n-1} + dt * v_n
241
242
243 """
244 def __init__(self, dt, u0, v0, t0=0.):
245 """
246 set up the wave base
247
248 :param dt: time step size (need to be sufficiently small)
249 :param u0: initial value
250 :param v0: initial velocity
251 :param t0: initial time
252 """
253 self.u=u0
254 self.v=v0
255 self.t=t0
256 self.n=0
257 self.t_last=t0
258 self.__dt=dt
259
260 def getTimeStepSize(self):
261 return self.__dt
262
263 def _getAcceleration(self, t, u):
264 """
265 returns the acceleraton for time `t` and solution `u` at time `t`
266 : note: this function needs to be overwritten by a particular wave problem
267 """
268 pass
269
270 def update(self, t):
271 """
272 returns the solution for the next time marker t which needs to greater than the time marker from the
273 previous call.
274 """
275 if not self.t_last <= t:
276 raise ValueError("You can not move backward in time.")
277
278 dt = self.getTimeStepSize()
279 # apply Verlet scheme until
280 while self.t < t:
281 a=self._getAcceleration(self.t, self.u)
282 self.v += dt * a
283 self.u += dt * self.v
284 self.t += dt
285 self.n += 1
286
287 # now we work backwards
288 self.t_last=t
289 return t, self.u + self.v * (t-self.t)
290
291 def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
292 """
293 Creates a distribution which is one in the interior of the domain of `x`
294 and is falling down to the value 'absorption_cut' over a margin of thickness 'absorption_zone'
295 toward each boundary except the top of the domain.
296
297 :param x: location of points in the domain
298 :type x: `Data`
299 :param absorption_zone: thickness of the aborption zone
300 :param absorption_cut: value of decay function on domain boundary
301 :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
302 toward the boundary.
303 """
304 if absorption_zone is None or absorption_zone == 0:
305 return 1
306
307 dom=x.getDomain()
308 bb=boundingBox(dom)
309 DIM=dom.getDim()
310 decay=-log(absorption_cut)/absorption_zone**2
311 f=1
312 for i in range(DIM):
313 x_i=x[i]
314 x_l=x_i-(bb[i][0]+absorption_zone)
315 m_l=whereNegative(x_l)
316 f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
317 if top_absorbation or not DIM-1 == i:
318 x_r=(bb[i][1]-absorption_zone)-x_i
319 m_r=whereNegative(x_r)
320 f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
321 return f
322
323 class SonicWave(WaveBase):
324 """
325 Solving the sonic wave equation
326
327 `p_tt = (v_p**2 * p_i)_i + f(t) * delta_s` where (p-) velocity v_p.
328
329 f(t) is wavelet acting at a point source term at positon s
330 """
331 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):
332 """
333 initialize the sonic wave solver
334
335 :param domain: domain of the problem
336 :type domain: `Domain`
337 :param v_p: p-velocity field
338 :type v_p: `Scalar`
339 :param wavelet: wavelet to describe the time evolution of source term
340 :type wavelet: `Wavelet`
341 :param source_tag: tag of the source location
342 :type source_tag: 'str' or 'int'
343 :param dt: time step size. If not present a suitable time step size is calculated.
344 :param p0: initial solution. If not present zero is used.
345 :param p0_t: initial solution change rate. If not present zero is used.
346 :param absorption_zone: thickness of absorption zone
347 :param absorption_cut: boundary value of absorption decay factor
348 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
349 """
350 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
351 v_p=v_p*f
352
353 if p0 == None:
354 p0=Scalar(0.,Solution(domain))
355 else:
356 p0=interpolate(p0, Solution(domain ))
357
358 if p0_t == None:
359 p0_t=Scalar(0.,Solution(domain))
360 else:
361 p0_t=interpolate(p0_t, Solution(domain ))
362
363 if dt == None:
364 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
365
366 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
367
368 self.__wavelet=wavelet
369 self.__mypde=LinearSinglePDE(domain)
370 if lumping: self.__mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
371 self.__mypde.setSymmetryOn()
372 self.__mypde.setValue(D=1./v_p**2)
373 self.__source_tag=source_tag
374 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
375
376
377 def _getAcceleration(self, t, u):
378 """
379 returns the acceleraton for time t and solution u at time t
380 """
381 self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getValue(t))
382 self.__mypde.setValue(X=-grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
383 return self.__mypde.getSolution()
384
385
386 class VTIWave(WaveBase):
387 """
388 Solving the VTI wave equation
389
390 :note: In case of a two dimensional domain the second spatial dimenion is depth.
391 """
392 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
393 source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
394 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
395 absorption_cut=1e-2, lumping=True):
396 """
397 initialize the VTI wave solver
398
399 :param domain: domain of the problem
400 :type domain: `Domain`
401 :param v_p: vertical p-velocity field
402 :type v_p: `Scalar`
403 :param v_s: vertical s-velocity field
404 :type v_s: `Scalar`
405 :param wavelet: wavelet to describe the time evolution of source term
406 :type wavelet: `Wavelet`
407 :param source_tag: tag of the source location
408 :type source_tag: 'str' or 'int'
409 :param source_vector: source orientation vector
410 :param eps: first Thompsen parameter
411 :param delta: second Thompsen parameter
412 :param gamma: third Thompsen parameter
413 :param rho: density
414 :param dt: time step size. If not present a suitable time step size is calculated.
415 :param u0: initial solution. If not present zero is used.
416 :param v0: initial solution change rate. If not present zero is used.
417 :param absorption_zone: thickness of absorption zone
418 :param absorption_cut: boundary value of absorption decay factor
419 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
420 """
421 DIM=domain.getDim()
422 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
423
424 v_p=v_p*f
425 v_s=v_s*f
426
427
428
429 if u0 == None:
430 u0=Vector(0.,Solution(domain))
431 else:
432 u0=interpolate(p0, Solution(domain ))
433
434 if v0 == None:
435 v0=Vector(0.,Solution(domain))
436 else:
437 v0=interpolate(v0, Solution(domain ))
438
439 if dt == None:
440 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
441
442 super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
443
444 self.__wavelet=wavelet
445
446 self.fastAssembler = hasattr(domain, "setAssembler")
447 self.c33=v_p**2 * rho
448 self.c44=v_s**2 * rho
449 self.c11=(1+2*eps) * self.c33
450 self.c66=(1+2*gamma) * self.c44
451 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
452 self.c12=self.c11-2*self.c66
453
454 if self.fastAssembler:
455 self.__mypde=WavePDE(domain, [("c11", self.c11),
456 ("c12", self.c12), ("c13", self.c13), ("c33", self.c33),
457 ("c44", self.c44), ("c66", self.c66)])
458 else:
459 self.__mypde=LinearPDESystem(domain)
460 self.__mypde.setValue(X=self.__mypde.createCoefficient('X'))
461
462 if lumping:
463 self.__mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
464 self.__mypde.setSymmetryOn()
465 self.__mypde.setValue(D=rho*kronecker(DIM))
466 self.__source_tag=source_tag
467
468 if DIM ==2 :
469 source_vector= [source_vector[0],source_vector[2]]
470
471 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
472 self.__r.setTaggedValue(self.__source_tag, source_vector)
473
474 def _getAcceleration(self, t, u):
475 """
476 returns the acceleraton for time `t` and solution `u` at time `t`
477 """
478 du = grad(u)
479 if not self.fastAssembler:
480 sigma=self.__mypde.getCoefficient('X')
481 if self.__mypde.getDim() == 3:
482 e11=du[0,0]
483 e22=du[1,1]
484 e33=du[2,2]
485 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
486 sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
487 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
488
489 s=self.c44*(du[2,1]+du[1,2])
490 sigma[1,2]=s
491 sigma[2,1]=s
492
493 s=self.c44*(du[2,0]+du[0,2])
494 sigma[0,2]=s
495 sigma[2,0]=s
496
497 s=self.c66*(du[0,1]+du[1,0])
498 sigma[0,1]=s
499 sigma[1,0]=s
500
501
502 else:
503 e11=du[0,0]
504 e22=du[1,1]
505 sigma[0,0]=self.c11*e11+self.c13*e22
506 sigma[1,1]=self.c13*e11+self.c33*e22
507 s=self.c44*(du[1,0]+du[0,1])
508 sigma[0,1]=s
509 sigma[1,0]=s
510 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getValue(t))
511
512 else:
513 self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getValue(t))
514
515 return self.__mypde.getSolution()
516
517
518 class HTIWave(WaveBase):
519 """
520 Solving the HTI wave equation (along the x_0 axis)
521
522 :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
523 """
524
525 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
526 source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
527 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
528 absorption_cut=1e-2, lumping=True):
529 """
530 initialize the VTI wave solver
531
532 :param domain: domain of the problem
533 :type domain: `Domain`
534 :param v_p: vertical p-velocity field
535 :type v_p: `Scalar`
536 :param v_s: vertical s-velocity field
537 :type v_s: `Scalar`
538 :param wavelet: wavelet to describe the time evolution of source term
539 :type wavelet: `Wavelet`
540 :param source_tag: tag of the source location
541 :type source_tag: 'str' or 'int'
542 :param source_vector: source orientation vector
543 :param eps: first Thompsen parameter
544 :param delta: second Thompsen parameter
545 :param gamma: third Thompsen parameter
546 :param rho: density
547 :param dt: time step size. If not present a suitable time step size is calculated.
548 :param u0: initial solution. If not present zero is used.
549 :param v0: initial solution change rate. If not present zero is used.
550 :param absorption_zone: thickness of absorption zone
551 :param absorption_cut: boundary value of absorption decay factor
552 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
553 """
554 DIM=domain.getDim()
555 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
556
557 v_p=v_p*f
558 v_s=v_s*f
559
560 if u0 == None:
561 u0=Vector(0.,Solution(domain))
562 else:
563 u0=interpolate(p0, Solution(domain ))
564
565 if v0 == None:
566 v0=Vector(0.,Solution(domain))
567 else:
568 v0=interpolate(v0, Solution(domain ))
569
570 if dt == None:
571 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
572
573 super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
574
575 self.__wavelet=wavelet
576
577 self.fastAssembler = hasattr(domain, "setAssembler")
578 self.c33=v_p**2 * rho
579 self.c44=v_s**2 * rho
580 self.c11=(1+2*eps) * self.c33
581 self.c66=(1+2*gamma) * self.c44
582 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
583 self.c23=self.c33-2*self.c66
584
585 if self.fastAssembler:
586 self.__mypde=WavePDE(domain, [("c11", self.c11),
587 ("c23", self.c23), ("c13", self.c13), ("c33", self.c33),
588 ("c44", self.c44), ("c66", self.c66)])
589 else:
590 self.__mypde=LinearPDESystem(domain)
591 self.__mypde.setValue(X=self.__mypde.createCoefficient('X'))
592
593 if lumping:
594 self.__mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
595 self.__mypde.setSymmetryOn()
596 self.__mypde.setValue(D=rho*kronecker(DIM))
597 self.__source_tag=source_tag
598
599 if DIM ==2 :
600 source_vector= [source_vector[0],source_vector[2]]
601
602 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
603 self.__r.setTaggedValue(self.__source_tag, source_vector)
604
605
606
607 def _getAcceleration(self, t, u):
608 """
609 returns the acceleraton for time `t` and solution `u` at time `t`
610 """
611 du = grad(u)
612 if self.fastAssembler:
613 self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getValue(t))
614 else:
615 sigma=self.__mypde.getCoefficient('X')
616
617 if self.__mypde.getDim() == 3:
618 e11=du[0,0]
619 e22=du[1,1]
620 e33=du[2,2]
621
622 sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
623 sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
624 sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
625
626 s=self.c44*(du[2,1]+du[1,2])
627 sigma[1,2]=s
628 sigma[2,1]=s
629
630 s=self.c66*(du[2,0]+du[0,2])
631 sigma[0,2]=s
632 sigma[2,0]=s
633
634 s=self.c66*(du[0,1]+du[1,0])
635 sigma[0,1]=s
636 sigma[1,0]=s
637
638 else:
639 e11=du[0,0]
640 e22=du[1,1]
641 sigma[0,0]=self.c11*e11+self.c13*e22
642 sigma[1,1]=self.c13*e11+self.c33*e22
643
644 s=self.c66*(du[1,0]+du[0,1])
645 sigma[0,1]=s
646 sigma[1,0]=s
647 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getValue(t))
648
649 return self.__mypde.getSolution()
650
651 class TTIWave(WaveBase):
652 """
653 Solving the 2D TTI wave equation with
654
655 `sigma_xx= c11*e_xx + c13*e_zz + c15*e_xz`
656 `sigma_zz= c13*e_xx + c33*e_zz + c35*e_xz`
657 `sigma_xz= c15*e_xx + c35*e_zz + c55*e_xz`
658
659 the coefficients `c11`, `c13`, etc are calculated from the tompsen parameters `eps`, `delta` and the tilt `theta`
660
661 :note: currently only the 2D case is supported.
662 """
663
664 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
665 source_vector = [0.,1.], eps=0., delta=0., theta=0., rho=1.,
666 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
667 absorption_cut=1e-2, lumping=True):
668 """
669 initialize the TTI wave solver
670
671 :param domain: domain of the problem
672 :type domain: `Domain`
673 :param v_p: vertical p-velocity field
674 :type v_p: `Scalar`
675 :param v_s: vertical s-velocity field
676 :type v_s: `Scalar`
677 :param wavelet: wavelet to describe the time evolution of source term
678 :type wavelet: `Wavelet`
679 :param source_tag: tag of the source location
680 :type source_tag: 'str' or 'int'
681 :param source_vector: source orientation vector
682 :param eps: first Thompsen parameter
683 :param delta: second Thompsen parameter
684 :param theta: tilting (in Rad)
685 :param rho: density
686 :param dt: time step size. If not present a suitable time step size is calculated.
687 :param u0: initial solution. If not present zero is used.
688 :param v0: initial solution change rate. If not present zero is used.
689 :param absorption_zone: thickness of absorption zone
690 :param absorption_cut: boundary value of absorption decay factor
691 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
692 """
693 DIM=domain.getDim()
694 if not DIM == 2:
695 raise ValueError("Only 2D is supported.")
696 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
697
698 v_p=v_p*f
699 v_s=v_s*f
700
701 if u0 == None:
702 u0=Vector(0.,Solution(domain))
703 else:
704 u0=interpolate(p0, Solution(domain ))
705
706 if v0 == None:
707 v0=Vector(0.,Solution(domain))
708 else:
709 v0=interpolate(v0, Solution(domain ))
710
711 if dt == None:
712 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
713
714 super(TTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
715
716 self.__wavelet=wavelet
717
718 self.__mypde=LinearPDESystem(domain)
719 if lumping: self.__mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
720 self.__mypde.setSymmetryOn()
721 self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
722 self.__source_tag=source_tag
723
724 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
725 self.__r.setTaggedValue(self.__source_tag, source_vector)
726
727 c0_33=v_p**2 * rho
728 c0_66=v_s**2 * rho
729 c0_11=(1+2*eps) * c0_33
730 c0_13=sqrt(2*c0_33*(c0_33-c0_66) * delta + (c0_33-c0_66)**2)-c0_66
731
732 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
733 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
734 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
735 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
736 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
737 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
738
739 def _getAcceleration(self, t, u):
740 """
741 returns the acceleraton for time `t` and solution `u` at time `t`
742 """
743 du = grad(u)
744 sigma=self.__mypde.getCoefficient('X')
745
746 e_xx=du[0,0]
747 e_zz=du[1,1]
748 e_xz=du[0,1]+du[1,0]
749
750
751 sigma[0,0]= self.c11 * e_xx + self.c13 * e_zz + self.c16 * e_xz
752 sigma[1,1]= self.c13 * e_xx + self.c33 * e_zz + self.c36 * e_xz
753 sigma_xz = self.c16 * e_xx + self.c36 * e_zz + self.c66 * e_xz
754
755 sigma[0,1]=sigma_xz
756 sigma[1,0]=sigma_xz
757
758 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getValue(t))
759 return self.__mypde.getSolution()
760
761 class SonicHTIWave(WaveBase):
762 """
763 Solving the HTI wave equation (along the x_0 axis) with azimuth (rotation around verticle axis)
764 under the assumption of zero shear wave velocities
765 The unknowns are the transversal (along x_0) and vertial stress (Q, P)
766
767 :note: In case of a two dimensional domain the second spatial dimenion is depth.
768 """
769 def __init__(self, domain, v_p, wavelet, source_tag, source_vector = [1.,0.], eps=0., delta=0., azimuth=0.,
770 dt=None, p0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
771 """
772 initialize the HTI wave solver
773
774 :param domain: domain of the problem
775 :type domain: `Doamin`
776 :param v_p: vertical p-velocity field
777 :type v_p: `Scalar`
778 :param v_s: vertical s-velocity field
779 :type v_s: `Scalar`
780 :param wavelet: wavelet to describe the time evolution of source term
781 :type wavelet: `Wavelet`
782 :param source_tag: tag of the source location
783 :type source_tag: 'str' or 'int'
784 :param source_vector: source orientation vector
785 :param eps: first Thompsen parameter
786 :param azimuth: azimuth (rotation around verticle axis)
787 :param gamma: third Thompsen parameter
788 :param rho: density
789 :param dt: time step size. If not present a suitable time step size is calculated.
790 :param p0: initial solution (Q(t=0), P(t=0)). If not present zero is used.
791 :param v0: initial solution change rate. If not present zero is used.
792 :param absorption_zone: thickness of absorption zone
793 :param absorption_cut: boundary value of absorption decay factor
794 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
795 """
796 DIM=domain.getDim()
797 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
798
799 self.v2_p=v_p**2
800 self.v2_t=self.v2_p*sqrt(1+2*delta)
801 self.v2_n=self.v2_p*(1+2*eps)
802
803 if p0 == None:
804 p0=Data(0.,(2,),Solution(domain))
805 else:
806 p0=interpolate(p0, Solution(domain ))
807
808 if v0 == None:
809 v0=Data(0.,(2,),Solution(domain))
810 else:
811 v0=interpolate(v0, Solution(domain ))
812
813 if dt == None:
814 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
815
816 super(SonicHTIWave, self).__init__( dt, u0=p0, v0=v0, t0=0.)
817
818 self.__wavelet=wavelet
819
820 self.__mypde=LinearPDESystem(domain)
821 if lumping: self.__mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
822 self.__mypde.setSymmetryOn()
823 self.__mypde.setValue(D=kronecker(2), X=self.__mypde.createCoefficient('X'))
824 self.__source_tag=source_tag
825
826
827 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
828 self.__r.setTaggedValue(self.__source_tag, source_vector)
829
830 def _getAcceleration(self, t, u):
831 """
832 returns the acceleraton for time `t` and solution `u` at time `t`
833 """
834 dQ = grad(u[0])[0]
835 dP = grad(u[1])[1:]
836 sigma=self.__mypde.getCoefficient('X')
837
838 sigma[0,0] = self.v2_n*dQ
839 sigma[0,1:] = self.v2_t*dP
840 sigma[1,0] = self.v2_t*dQ
841 sigma[1,1:] = self.v2_p*dP
842
843 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getValue(t))
844 return self.__mypde.getSolution()

  ViewVC Help
Powered by ViewVC 1.1.26