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

  ViewVC Help
Powered by ViewVC 1.1.26