/[escript]/trunk/downunder/py_src/seismic.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/seismic.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4599 - (hide annotations)
Sun Dec 15 23:51:13 2013 UTC (5 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 25794 byte(s)
merged 3.4.1 release changes into trunk
1 gross 4538
2     ##############################################################################
3     #
4     # Copyright (c) 2003-2013 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 since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15    
16     __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17     http://www.uq.edu.au
18     Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23 gross 4573 __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave', 'HTIWave' ]
24 gross 4538
25    
26 gross 4541 from math import pi
27 gross 4538 import numpy as np
28     import sys
29     import time
30     from esys.escript import *
31     import esys.escript.unitsSI as U
32 gross 4562 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
33 gross 4538
34 gross 4544
35    
36    
37 gross 4540 class Wavelet(object):
38 sshaw 4556 """
39     place holder for source wavelet
40     """
41     pass
42    
43 gross 4540 class Ricker(Wavelet):
44 sshaw 4556 """
45     The Ricker Wavelet w=f(t)
46     """
47     def __init__(self, f_dom=40, t_dom=None):
48     """
49     set up Ricker wavelet wih dominant frequence f_dom and center at time t_dom. If t_dom
50     is not given an estimate for suitable t_dom is calculated so f(0)~0.
51 gross 4573
52     :note: maximum frequence is about 2 x the dominant frequence.
53 sshaw 4556 """
54     drop=18
55     self.__f=f_dom
56     self.__f_max=sqrt(7)*f_dom
57     self.__s=pi*self.__f
58     if t_dom == None:
59     t_dom=sqrt(drop)/self.__s
60     self.__t0=t_dom
61    
62     def getCenter(self):
63     """
64     return value of wavelet center
65     """
66     return self.__t0
67 gross 4540
68     def getTimeScale(self):
69     """
70     returns the time scale which is the inverse of the largest freqence with a significant
71     spectral component.
72     """
73 sshaw 4556 return 1/self.__f_max
74 gross 4538
75 sshaw 4556 def getValue(self, t):
76     """
77     get value of wavelet at time t
78     """
79     e2=(self.__s*(t-self.__t0))**2
80     return (1-2*e2)*exp(-e2)
81 gross 4540
82 sshaw 4556 def getAcceleration(self, t):
83     """
84     get the acceleration f''(t) at time t
85     """
86     e2=(self.__s*(t-self.__t0))**2
87 gross 4540 return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)
88    
89    
90 gross 4538 class SimpleSEGYWriter(object):
91 sshaw 4556 """
92     as simple writer for 2D and 3D seimic lines in particular for synthetic data
93 gross 4538
94 sshaw 4556 Typical usage:
95 gross 4538
96 sshaw 4556 from esys.escript import unitsSI as U
97     sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)
98     while n < 10:
99     sw.addRecord([i*2., i*0.67, i**2, -i*7])
100     sw.write('example.segy')
101 gross 4538
102 sshaw 4599 :note: the writer uses `obspy`
103 sshaw 4556 """
104     COORDINATE_SCALE = 1000.
105     def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):
106     """
107     initalize writer
108 gross 4538
109 sshaw 4556 :param receiver_group: list of receiver coordinates (in meters). For the 2D case a list of floats is given, for the 3D case a list of coordinate tuples
110     are given
111     :param source: coordinates of the source (in meters). For the 2D case a single floats is given, for the 3D case a coordinate tuples
112     :param sampling_interval: sample rate in seconds
113     :param text: a text for the header file (e.g a description)
114     """
115     if isinstance(source, float) or isinstance(source, int) :
116     DIM=1
117     source = (source, 0.)
118     elif hasattr(source, '__len__'):
119     DIM=len(source)
120     if DIM == 1:
121     source = (source[0], 0.)
122     elif DIM==2:
123     source = (source[0], source[1] )
124     else:
125     raise ValueError("Only 1D or 2D arrays accepted.")
126     else:
127     raise TypeError("illegal source type.")
128     if receiver_group== None:
129     if DIM == 1:
130     receiver_group=[0.]
131     else:
132     receiver_group=[(0.,0.)]
133 gross 4538
134 sshaw 4556 if isinstance(receiver_group, float) or isinstance(receiver_group, int) :
135     if DIM == 1:
136     rg = [ (receiver_group, 0. ) ]
137     else:
138     raise TypeError("illegal receiver_group type.")
139     elif hasattr(receiver_group, '__len__'):
140     if DIM == 1:
141     rg = [(c,0) for c in receiver_group]
142     else:
143     rg = [(c[0],c[1]) for c in receiver_group]
144     else:
145     raise TypeError("illegal receiver_group type.")
146 gross 4538
147 sshaw 4556 self.__source=source
148     self.__receiver_group=rg
149     self.__text=text
150 gross 4538 self.__trace = [ [] for i in self.__receiver_group ]
151     self.__sampling_interval = sampling_interval
152    
153 sshaw 4556 def addRecord(self, record):
154     """
155     adds a record to the traces. a time difference of sample_interval between two records is assumed.
156     The record mast be a list of as many values as given receivers or a float if a single receiver is used.
157 gross 4538
158     :param record: list of tracks to be added to the record.
159 sshaw 4556 """
160     if not len(self.__receiver_group) == len(record):
161     raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))
162     if len(self.__receiver_group) == 1:
163     if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
164     self.__trace[0].append(record)
165     else:
166     self.__trace[0].append(record[0])
167     else:
168 sshaw 4576 for i in range(len(record)):
169 sshaw 4556 self.__trace[i].append(record[i])
170 gross 4538
171 sshaw 4556 def getSamplingInterval(self):
172     """
173     returns the sampling interval in seconds.
174     """
175 gross 4538 return self.__sampling_interval
176    
177 sshaw 4556 def write(self, filename):
178 gross 4538 """
179     writes to segy file
180    
181     :param filename: file name
182     :note: the function uses the `obspy` module.
183     """
184 sshaw 4599 try:
185     from obspy import Trace, Stream, UTCDateTime
186     from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
187     from obspy.core import AttribDict
188     except ImportError as e:
189     raise RuntimeError("This feature (SimpleSEGYWriter.write())"+\
190     " depends on obspy, which is not installed, see "+\
191     "https://github.com/obspy/obspy for install guide")
192 gross 4538
193 sshaw 4556 stream=Stream()
194    
195 sshaw 4576 for i in range(len(self.__receiver_group)):
196 sshaw 4556 trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
197 gross 4538 # Attributes in trace.stats will overwrite everything in
198     # trace.stats.segy.trace_header (in Hz)
199     trace.stats.sampling_rate = 1./self.getSamplingInterval()
200     #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0)
201     if not hasattr(trace.stats, 'segy.trace_header'):
202     trace.stats.segy = {}
203     trace.stats.segy.trace_header = SEGYTraceHeader()
204     trace.stats.segy.trace_header.trace_identification_code=1
205     trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
206     trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = int(self.COORDINATE_SCALE)
207     trace.stats.segy.trace_header.coordinate_units=1
208     trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE)
209     trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE)
210     trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE)
211     trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE)
212    
213     # Add trace to stream
214     stream.append(trace)
215     # A SEGY file has file wide headers. This can be attached to the stream
216     # object. If these are not set, they will be autocreated with default
217     # values.
218     stream.stats = AttribDict()
219 gross 4549 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
220 gross 4538 stream.stats.binary_file_header = SEGYBinaryFileHeader()
221    
222     if getMPIRankWorld()<1:
223     stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
224 gross 4545
225     class WaveBase(object):
226     """
227     Base for wave propagation using the Verlet scheme.
228 gross 4538
229 gross 4545 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
230 gross 4538
231 gross 4545 with a given acceleration force as function of time.
232    
233     a_n=A(t_{n-1})
234     v_n=v_{n-1} + dt * a_n
235     u_n=u_{n-1} + dt * v_n
236    
237    
238     """
239     def __init__(self, dt, u0, v0, t0=0.):
240     """
241     set up the wave base
242    
243     :param dt: time step size (need to be sufficiently small)
244     :param u0: initial value
245     :param v0: initial velocity
246     :param t0: initial time
247     """
248     self.u=u0
249     self.v=v0
250     self.t=t0
251     self.t_last=t0
252     self.__dt=dt
253    
254     def getTimeStepSize(self):
255 sshaw 4556 return self.__dt
256    
257 gross 4545 def _getAcceleration(self, t, u):
258     """
259     returns the acceleraton for time t and solution u at time t
260     : note: this function needs to be overwritten by a particular wave problem
261     """
262     pass
263    
264     def update(self, t):
265     """
266     returns the solution for the next time marker t which needs to greater than the time marker from the
267     previous call.
268     """
269     if not self.t_last <= t:
270     raise ValueError("You can not move backward in time.")
271    
272     dt = self.getTimeStepSize()
273     # apply Verlet scheme until
274     while self.t < t:
275     a=self._getAcceleration(self.t, self.u)
276 sshaw 4556 self.v += dt * a
277     self.u += dt * self.v
278 gross 4545 self.t += dt
279    
280     # now we work backwards
281     self.t_last=t
282     return t, self.u + self.v * (t-self.t)
283    
284 gross 4573 def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
285 gross 4545 """
286     creating a distribution which is one in the interior of the domain of x
287     and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
288     toward each boundary except the top of the domain.
289    
290     :param x: location of points in the domain
291     :type x: `Data`
292     :param absorption_zone: thickness of the aborption zone
293     :param absorption_cut: value of decay function on domain boundary
294     :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
295     toward the boundary.
296     """
297     dom=x.getDomain()
298     bb=boundingBox(dom)
299 gross 4546 DIM=dom.getDim()
300 gross 4545 decay=-log(absorption_cut)/absorption_zone**2
301     f=1
302 sshaw 4576 for i in range(DIM):
303 gross 4545 x_i=x[i]
304     x_l=x_i-(bb[i][0]+absorption_zone)
305 sshaw 4556 m_l=whereNegative(x_l)
306     f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
307 gross 4573 if top_absorbation or not DIM-1 == i:
308 sshaw 4556 x_r=(bb[i][1]-absorption_zone)-x_i
309     m_r=whereNegative(x_r)
310     f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
311 gross 4546 return f
312 gross 4545
313     class SonicWave(WaveBase):
314 sshaw 4556 """
315     Solving the sonic wave equation
316    
317     p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.
318    
319     f(t) is wavelet acting at a point source term at positon s
320     """
321     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):
322 gross 4545 """
323     initialize the sonic wave solver
324    
325     :param domain: domain of the problem
326     :type domain: `Doamin`
327     :param v_p: p-velocity field
328     :type v_p: `Scalar`
329     :param wavelet: wavelet to describe the time evolution of source term
330     :type wavelet: `Wavelet`
331     :param source_tag: tag of the source location
332     :type source_tag: 'str' or 'int'
333     :param dt: time step size. If not present a suitable time step size is calculated.
334     :param p0: initial solution. If not present zero is used.
335     :param p0_t: initial solution change rate. If not present zero is used.
336     :param absorption_zone: thickness of absorption zone
337     :param absorption_cut: boundary value of absorption decay factor
338 gross 4549 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
339 gross 4545 """
340     f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
341     v_p=v_p*f
342    
343 sshaw 4556 if p0 == None:
344     p0=Scalar(0.,Solution(domain))
345     else:
346     p0=interpolate(p0, Solution(domain ))
347    
348     if p0_t == None:
349     p0_t=Scalar(0.,Solution(domain))
350     else:
351     p0_t=interpolate(p0_t, Solution(domain ))
352    
353     if dt == None:
354 gross 4545 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
355    
356 gross 4546 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
357 gross 4545
358     self.__wavelet=wavelet
359 gross 4562 self.__mypde=LinearSinglePDE(domain)
360 gross 4549 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
361     self.__mypde.setSymmetryOn()
362 gross 4545 self.__mypde.setValue(D=1.)
363 sshaw 4556 self.__source_tag=source_tag
364 gross 4545 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
365 sshaw 4556 self.__vp2=v_p**2
366 gross 4545
367 gross 4546
368     def _getAcceleration(self, t, u):
369 gross 4545 """
370     returns the acceleraton for time t and solution u at time t
371     """
372     self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t))
373 gross 4549 self.__mypde.setValue(X=-self.__vp2*grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
374 gross 4545 return self.__mypde.getSolution()
375    
376    
377    
378 gross 4562 class VTIWave(WaveBase):
379     """
380     Solving the VTI wave equation
381 gross 4573
382     :note: In case of a two dimensional domain the second spatial dimenion is depth.
383 gross 4562 """
384     def __init__(self, domain, v_p, v_s, wavelet, source_tag, source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
385     dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
386     """
387     initialize the VTI wave solver
388    
389     :param domain: domain of the problem
390     :type domain: `Doamin`
391 gross 4573 :param v_p: vertical p-velocity field
392     :type v_p: `Scalar`
393     :param v_s: vertical s-velocity field
394     :type v_s: `Scalar`
395 gross 4562 :param wavelet: wavelet to describe the time evolution of source term
396     :type wavelet: `Wavelet`
397     :param source_tag: tag of the source location
398 gross 4573 :type source_tag: 'str' or 'int'
399     :param source_vector: source orientation vector
400     :param eps: first Thompsen parameter
401     :param delta: second Thompsen parameter
402     :param gamma: third Thompsen parameter
403     :param rho: density
404 gross 4562 :param dt: time step size. If not present a suitable time step size is calculated.
405 gross 4573 :param u0: initial solution. If not present zero is used.
406     :param v0: initial solution change rate. If not present zero is used.
407 gross 4562 :param absorption_zone: thickness of absorption zone
408     :param absorption_cut: boundary value of absorption decay factor
409     :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
410     """
411     DIM=domain.getDim()
412     f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
413    
414     v_p=v_p*f
415     v_s=v_s*f
416    
417    
418    
419     if u0 == None:
420     u0=Vector(0.,Solution(domain))
421     else:
422     u0=interpolate(p0, Solution(domain ))
423    
424     if v0 == None:
425     v0=Vector(0.,Solution(domain))
426     else:
427     v0=interpolate(v0, Solution(domain ))
428    
429     if dt == None:
430     dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
431    
432     super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
433    
434     self.__wavelet=wavelet
435    
436     self.__mypde=LinearPDESystem(domain)
437     if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
438     self.__mypde.setSymmetryOn()
439     self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
440     self.__source_tag=source_tag
441    
442    
443     if DIM ==2 :
444 sshaw 4576 source_vector= [source_vector[0],source_vector[2]]
445 gross 4562
446 sshaw 4576
447     self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
448 gross 4562 self.__r.setTaggedValue(self.__source_tag, source_vector)
449    
450    
451     self.c33=v_p**2 * rho
452 sshaw 4576 self.c44=v_s**2 * rho
453     self.c11=(1+2*eps) * self.c33
454     self.c66=(1+2*gamma) * self.c44
455     self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
456     self.c12=self.c11-2*self.c66
457    
458 gross 4562 def _getAcceleration(self, t, u):
459     """
460     returns the acceleraton for time t and solution u at time t
461     """
462 gross 4563 du = grad(u)
463 gross 4562 sigma=self.__mypde.getCoefficient('X')
464    
465     if self.__mypde.getDim() == 3:
466 sshaw 4576 e11=du[0,0]
467     e22=du[1,1]
468     e33=du[2,2]
469 gross 4562
470 sshaw 4576 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
471     sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
472     sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
473 gross 4562
474 sshaw 4576 s=self.c44*(du[2,1]+du[1,2])
475     sigma[1,2]=s
476     sigma[2,1]=s
477 gross 4573
478 sshaw 4576 s=self.c44*(du[2,0]+du[0,2])
479     sigma[0,2]=s
480     sigma[2,0]=s
481    
482     s=self.c66*(du[0,1]+du[1,0])
483     sigma[0,1]=s
484     sigma[1,0]=s
485    
486 gross 4573
487 gross 4562 else:
488 sshaw 4576 e11=du[0,0]
489     e22=du[1,1]
490 gross 4562
491 sshaw 4576 sigma[0,0]=self.c11*e11+self.c13*e22
492     sigma[1,1]=self.c13*e11+self.c33*e22
493 gross 4562
494 sshaw 4576 s=self.c44*(du[1,0]+du[0,1])
495     sigma[0,1]=s
496     sigma[1,0]=s
497 gross 4562
498     self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
499     return self.__mypde.getSolution()
500    
501 gross 4573 class HTIWave(WaveBase):
502     """
503     Solving the HTI wave equation (along the x_0 axis)
504    
505     :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
506     """
507     def __init__(self, domain, v_p, v_s, wavelet, source_tag, source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
508     dt=None, u0=None, v0=None, absorption_zone=300*U.m, absorption_cut=1e-2, lumping=True):
509     """
510     initialize the VTI wave solver
511    
512     :param domain: domain of the problem
513     :type domain: `Doamin`
514     :param v_p: vertical p-velocity field
515     :type v_p: `Scalar`
516     :param v_s: vertical s-velocity field
517     :type v_s: `Scalar`
518     :param wavelet: wavelet to describe the time evolution of source term
519     :type wavelet: `Wavelet`
520     :param source_tag: tag of the source location
521     :type source_tag: 'str' or 'int'
522     :param source_vector: source orientation vector
523     :param eps: first Thompsen parameter
524     :param delta: second Thompsen parameter
525     :param gamma: third Thompsen parameter
526     :param rho: density
527     :param dt: time step size. If not present a suitable time step size is calculated.
528     :param u0: initial solution. If not present zero is used.
529     :param v0: initial solution change rate. If not present zero is used.
530     :param absorption_zone: thickness of absorption zone
531     :param absorption_cut: boundary value of absorption decay factor
532     :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
533     """
534     DIM=domain.getDim()
535     f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut, top_absorbation=True)
536    
537     v_p=v_p*f
538     v_s=v_s*f
539    
540     if u0 == None:
541     u0=Vector(0.,Solution(domain))
542     else:
543     u0=interpolate(p0, Solution(domain ))
544    
545     if v0 == None:
546     v0=Vector(0.,Solution(domain))
547     else:
548     v0=interpolate(v0, Solution(domain ))
549    
550     if dt == None:
551     dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
552 gross 4562
553 gross 4573 super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
554    
555     self.__wavelet=wavelet
556    
557     self.__mypde=LinearPDESystem(domain)
558     if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
559     self.__mypde.setSymmetryOn()
560     self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
561     self.__source_tag=source_tag
562    
563 gross 4562
564 gross 4573 if DIM ==2 :
565 sshaw 4576 source_vector= [source_vector[0],source_vector[1]]
566 gross 4573
567 sshaw 4576
568     self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
569 gross 4573 self.__r.setTaggedValue(self.__source_tag, source_vector)
570    
571    
572     self.c33=v_p**2 * rho
573 sshaw 4576 self.c44=v_s**2 * rho
574     self.c11=(1+2*eps) * self.c33
575     self.c66=(1+2*gamma) * self.c44
576     self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
577     self.c23=self.c33-2*self.c66
578    
579 gross 4573 def _getAcceleration(self, t, u):
580     """
581     returns the acceleraton for time t and solution u at time t
582     """
583     du = grad(u)
584     sigma=self.__mypde.getCoefficient('X')
585    
586     if self.__mypde.getDim() == 3:
587 sshaw 4576 e11=du[0,0]
588     e22=du[1,1]
589     e33=du[2,2]
590 gross 4573
591 sshaw 4576 sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
592     sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
593     sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
594 gross 4573
595 sshaw 4576 s=self.c44*(du[2,1]+du[1,2])
596     sigma[1,2]=s
597     sigma[2,1]=s
598 gross 4573
599 sshaw 4576 s=self.c66*(du[2,0]+du[0,2])
600     sigma[0,2]=s
601     sigma[2,0]=s
602    
603     s=self.c66*(du[0,1]+du[1,0])
604     sigma[0,1]=s
605     sigma[1,0]=s
606 gross 4573
607     else:
608 sshaw 4576 e11=du[0,0]
609     e22=du[1,1]
610 gross 4573
611 sshaw 4576 sigma[0,0]=self.c11*e11+self.c13*e22
612     sigma[1,1]=self.c13*e11+self.c33*e22
613    
614     s=self.c66*(du[0,1]+du[1,0])
615     sigma[0,1]=s
616     sigma[1,0]=s
617 gross 4573
618     self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
619     return self.__mypde.getSolution()
620    

  ViewVC Help
Powered by ViewVC 1.1.26