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

  ViewVC Help
Powered by ViewVC 1.1.26