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

  ViewVC Help
Powered by ViewVC 1.1.26