/[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 4549 - (hide annotations)
Fri Nov 29 05:24:40 2013 UTC (5 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 13765 byte(s)
argument for lumping added.
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 4546 __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave' ]
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 4546 from esys.escript.linearPDEs import LinearPDE
33 gross 4538
34 gross 4544
35    
36    
37 gross 4540 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     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     """
52 gross 4541 drop=18
53 gross 4540 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 gross 4541
60 gross 4540 def getCenter(self):
61     """
62     return value of wavelet center
63     """
64 gross 4541 return self.__t0
65 gross 4540
66     def getTimeScale(self):
67     """
68     returns the time scale which is the inverse of the largest freqence with a significant
69     spectral component.
70     """
71     return 1/self.__f_max
72 gross 4538
73 gross 4540 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    
80     def getAcceleration(self, t):
81     """
82     get the acceleration f''(t) at time t
83     """
84     e2=(self.__s*(t-self.__t0))**2
85     return 2*self.__s**2*(-4*e2**2 + 12*e2 - 3)*exp(-e2)
86    
87    
88 gross 4538 class SimpleSEGYWriter(object):
89     """
90     as simple writer for 2D and 3D seimic lines in particular for synthetic data
91    
92     Typical usage:
93    
94     from esys.escript import unitsSI as U
95     sw=SimpleSEGYWriter([0.,100*U.m,200*U,m,300.], source=200*U.m, sampling_interval=4*U.msec)
96     while n < 10:
97     sw.addRecord([i*2., i*0.67, i**2, -i*7])
98     sw.write('example.segy')
99    
100     :note: the writer uses `obsy`
101     """
102     COORDINATE_SCALE = 1000.
103     def __init__(self, receiver_group=None, source=0., sampling_interval=4*U.msec, text="some seimic data"):
104     """
105     initalize writer
106    
107     :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
108     are given
109     :param source: coordinates of the source (in meters). For the 2D case a single floats is given, for the 3D case a coordinate tuples
110     :param sampling_interval: sample rate in seconds
111     :param text: a text for the header file (e.g a description)
112     """
113     if isinstance(source, float) or isinstance(source, int) :
114     DIM=1
115     source = (source, 0.)
116     elif hasattr(source, '__len__'):
117     DIM=len(source)
118     if DIM == 1:
119     source = (source[0], 0.)
120     elif DIM==2:
121     source = (source[0], source[1] )
122     else:
123     raise ValueError("Only 1D or 2D arrays accepted.")
124     else:
125     raise TypeError("illegal source type.")
126     if receiver_group== None:
127     if DIM == 1:
128     receiver_group=[0.]
129     else:
130     receiver_group=[(0.,0.)]
131    
132     if isinstance(receiver_group, float) or isinstance(receiver_group, int) :
133     if DIM == 1:
134     rg = [ (receiver_group, 0. ) ]
135     else:
136     raise TypeError("illegal receiver_group type.")
137     elif hasattr(receiver_group, '__len__'):
138     if DIM == 1:
139     rg = [(c,0) for c in receiver_group]
140     else:
141     rg = [(c[0],c[1]) for c in receiver_group]
142     else:
143     raise TypeError("illegal receiver_group type.")
144    
145     self.__source=source
146     self.__receiver_group=rg
147     self.__text=text
148     self.__trace = [ [] for i in self.__receiver_group ]
149     self.__sampling_interval = sampling_interval
150    
151     def addRecord(self, record):
152     """
153     adds a record to the traces. a time difference of sample_interval between two records is assumed.
154     The record mast be a list of as many values as given receivers or a float if a single receiver is used.
155    
156     :param record: list of tracks to be added to the record.
157     """
158     if not len(self.__receiver_group) == len(record):
159     raise ValueError("expected number of records is %s but %s given."%(len(self.__receiver_group), len(record)))
160     if len(self.__receiver_group) == 1:
161     if isinstance(self.__receiver_group, float) or isinstance(self.__receiver_group, int):
162     self.__trace[0].append(record)
163     else:
164     self.__trace[0].append(record[0])
165     else:
166     for i in xrange(len(record)):
167     self.__trace[i].append(record[i])
168    
169     def getSamplingInterval(self):
170     """
171     returns the sampling interval in seconds.
172     """
173     return self.__sampling_interval
174    
175     def write(self, filename):
176     """
177     writes to segy file
178    
179     :param filename: file name
180     :note: the function uses the `obspy` module.
181     """
182     from obspy import Trace, Stream, UTCDateTime
183     from obspy.segy.segy import SEGYTraceHeader, SEGYBinaryFileHeader
184     from obspy.core import AttribDict
185    
186     stream=Stream()
187    
188     for i in xrange(len(self.__receiver_group)):
189     trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
190     # Attributes in trace.stats will overwrite everything in
191     # trace.stats.segy.trace_header (in Hz)
192     trace.stats.sampling_rate = 1./self.getSamplingInterval()
193     #trace.stats.starttime = UTCDateTime(2011,11,11,0,0,0)
194     if not hasattr(trace.stats, 'segy.trace_header'):
195     trace.stats.segy = {}
196     trace.stats.segy.trace_header = SEGYTraceHeader()
197     trace.stats.segy.trace_header.trace_identification_code=1
198     trace.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
199     trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = int(self.COORDINATE_SCALE)
200     trace.stats.segy.trace_header.coordinate_units=1
201     trace.stats.segy.trace_header.source_coordinate_x=int(self.__source[0] * self.COORDINATE_SCALE)
202     trace.stats.segy.trace_header.source_coordinate_y=int(self.__source[1] * self.COORDINATE_SCALE)
203     trace.stats.segy.trace_header.group_coordinate_x=int(self.__receiver_group[i][0] * self.COORDINATE_SCALE)
204     trace.stats.segy.trace_header.group_coordinate_y=int(self.__receiver_group[i][1] * self.COORDINATE_SCALE)
205    
206     # Add trace to stream
207     stream.append(trace)
208     # A SEGY file has file wide headers. This can be attached to the stream
209     # object. If these are not set, they will be autocreated with default
210     # values.
211     stream.stats = AttribDict()
212 gross 4549 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
213 gross 4538 stream.stats.binary_file_header = SEGYBinaryFileHeader()
214    
215     if getMPIRankWorld()<1:
216     stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
217 gross 4545
218     class WaveBase(object):
219     """
220     Base for wave propagation using the Verlet scheme.
221 gross 4538
222 gross 4545 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
223 gross 4538
224 gross 4545 with a given acceleration force as function of time.
225    
226     a_n=A(t_{n-1})
227     v_n=v_{n-1} + dt * a_n
228     u_n=u_{n-1} + dt * v_n
229    
230    
231     """
232     def __init__(self, dt, u0, v0, t0=0.):
233     """
234     set up the wave base
235    
236     :param dt: time step size (need to be sufficiently small)
237     :param u0: initial value
238     :param v0: initial velocity
239     :param t0: initial time
240     """
241     self.u=u0
242     self.v=v0
243     self.t=t0
244     self.t_last=t0
245     self.__dt=dt
246    
247     def getTimeStepSize(self):
248     return self.__dt
249    
250     def _getAcceleration(self, t, u):
251     """
252     returns the acceleraton for time t and solution u at time t
253     : note: this function needs to be overwritten by a particular wave problem
254     """
255     pass
256    
257     def update(self, t):
258     """
259     returns the solution for the next time marker t which needs to greater than the time marker from the
260     previous call.
261     """
262     if not self.t_last <= t:
263     raise ValueError("You can not move backward in time.")
264    
265     dt = self.getTimeStepSize()
266     # apply Verlet scheme until
267     while self.t < t:
268     a=self._getAcceleration(self.t, self.u)
269     self.v += dt * a
270     self.u += dt * self.v
271     self.t += dt
272    
273     # now we work backwards
274     self.t_last=t
275     return t, self.u + self.v * (t-self.t)
276    
277     def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2):
278     """
279     creating a distribution which is one in the interior of the domain of x
280     and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
281     toward each boundary except the top of the domain.
282    
283     :param x: location of points in the domain
284     :type x: `Data`
285     :param absorption_zone: thickness of the aborption zone
286     :param absorption_cut: value of decay function on domain boundary
287     :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
288     toward the boundary.
289     """
290     dom=x.getDomain()
291     bb=boundingBox(dom)
292 gross 4546 DIM=dom.getDim()
293 gross 4545 decay=-log(absorption_cut)/absorption_zone**2
294     f=1
295     for i in xrange(DIM):
296     x_i=x[i]
297     x_l=x_i-(bb[i][0]+absorption_zone)
298     m_l=whereNegative(x_l)
299     f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
300     if not DIM-1 == i:
301     x_r=(bb[i][1]-absorption_zone)-x_i
302     m_r=whereNegative(x_r)
303     f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
304 gross 4546 return f
305 gross 4545
306     class SonicWave(WaveBase):
307     """
308     Solving the sonic wave equation
309    
310     p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.
311    
312     f(t) is wavelet acting at a point source term at positon s
313     """
314 gross 4549 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):
315 gross 4545 """
316     initialize the sonic wave solver
317    
318     :param domain: domain of the problem
319     :type domain: `Doamin`
320     :param v_p: p-velocity field
321     :type v_p: `Scalar`
322     :param wavelet: wavelet to describe the time evolution of source term
323     :type wavelet: `Wavelet`
324     :param source_tag: tag of the source location
325     :type source_tag: 'str' or 'int'
326     :param dt: time step size. If not present a suitable time step size is calculated.
327     :param p0: initial solution. If not present zero is used.
328     :param p0_t: initial solution change rate. If not present zero is used.
329     :param absorption_zone: thickness of absorption zone
330     :param absorption_cut: boundary value of absorption decay factor
331 gross 4549 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
332 gross 4545 """
333     f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
334     v_p=v_p*f
335    
336     if p0 == None:
337     p0=Scalar(0.,Solution(domain))
338     else:
339 gross 4546 p0=interpolate(p0, Solution(domain ))
340 gross 4545
341     if p0_t == None:
342     p0_t=Scalar(0.,Solution(domain))
343     else:
344 gross 4546 p0_t=interpolate(p0_t, Solution(domain ))
345 gross 4545
346     if dt == None:
347     dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
348    
349 gross 4546 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
350 gross 4545
351     self.__wavelet=wavelet
352     self.__mypde=LinearPDE(domain)
353 gross 4549 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
354     self.__mypde.setSymmetryOn()
355 gross 4545 self.__mypde.setValue(D=1.)
356     self.__source_tag=source_tag
357     self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
358     self.__vp2=v_p**2
359    
360 gross 4546
361     def _getAcceleration(self, t, u):
362 gross 4545 """
363     returns the acceleraton for time t and solution u at time t
364     """
365     self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t))
366 gross 4549 self.__mypde.setValue(X=-self.__vp2*grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
367 gross 4545 return self.__mypde.getSolution()
368    
369    
370    

  ViewVC Help
Powered by ViewVC 1.1.26