/[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 5665 - (show annotations)
Fri Jun 19 02:08:01 2015 UTC (4 years ago) by sshaw
File MIME type: text/x-python
File size: 36399 byte(s)
cleaning up seismic.py whitespace and fixing spelling mistakes in function names and params, the existing ones have been left with warnings about them being deprecated (createAbsorptionLayerFunction)

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

  ViewVC Help
Powered by ViewVC 1.1.26