/[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 4624 - (show annotations)
Sun Jan 19 23:03:59 2014 UTC (5 years, 4 months ago) by sshaw
File MIME type: text/x-python
File size: 25132 byte(s)
removed some leftover debug prints
1
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 __all__ = ['SimpleSEGYWriter', 'Ricker', 'WaveBase', 'SonicWave', 'VTIWave',
24 'HTIWave', 'createAbsorbtionLayerFunction' ]
25
26
27 from math import pi
28 import numpy as np
29 import sys
30 import time
31 from esys.escript import *
32 import esys.escript.unitsSI as U
33 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem, VTIWavePDE
34
35 class Wavelet(object):
36 """
37 place holder for source wavelet
38 """
39 pass
40
41 class Ricker(Wavelet):
42 """
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
50 :note: maximum frequence is about 2 x the dominant frequence.
51 """
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
60 def getCenter(self):
61 """
62 return value of wavelet center
63 """
64 return self.__t0
65
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
73 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 class SimpleSEGYWriter(object):
88 """
89 as simple writer for 2D and 3D seimic lines in particular for synthetic data
90
91 Typical usage:
92
93 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
99 :note: the writer uses `obspy`
100 """
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
106 :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 if receiver_group== None:
126 if DIM == 1:
127 receiver_group=[0.]
128 else:
129 receiver_group=[(0.,0.)]
130
131 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
144 self.__source=source
145 self.__receiver_group=rg
146 self.__text=text
147 self.__trace = [ [] for i in self.__receiver_group ]
148 self.__sampling_interval = sampling_interval
149
150 def addRecord(self, record):
151 """
152 adds a record to the traces. a time difference of sample_interval between two records is assumed.
153 The record mast be a list of as many values as given receivers or a float if a single receiver is used.
154
155 :param record: list of tracks to be added to the record.
156 """
157 if not len(self.__receiver_group) == len(record):
158 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 for i in range(len(record)):
166 self.__trace[i].append(record[i])
167
168 def getSamplingInterval(self):
169 """
170 returns the sampling interval in seconds.
171 """
172 return self.__sampling_interval
173
174 def write(self, filename):
175 """
176 writes to segy file
177
178 :param filename: file name
179 :note: the function uses the `obspy` module.
180 """
181 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
190 stream=Stream()
191
192 for i in range(len(self.__receiver_group)):
193 trace = Trace(data=np.array(self.__trace[i], dtype='float32'))
194 # 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 trace.stats.segy.trace_header.scalar_to_be_applied_to_all_coordinates = -int(self.COORDINATE_SCALE)
204 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 stream.stats.textual_file_header = 'C.. '+self.__text+'\nC.. with esys.escript.downunder r%s\nC.. %s'%(getVersion(),time.asctime())
217 stream.stats.binary_file_header = SEGYBinaryFileHeader()
218
219 if getMPIRankWorld()<1:
220 stream.write(filename, format="SEGY", data_encoding=1,byteorder=sys.byteorder)
221
222 class WaveBase(object):
223 """
224 Base for wave propagation using the Verlet scheme.
225
226 u_tt = A(t,u), u(t=t0)=u0, u_t(t=t0)=v0
227
228 with a given acceleration force as function of time.
229
230 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
234
235 """
236 def __init__(self, dt, u0, v0, t0=0.):
237 """
238 set up the wave base
239
240 :param dt: time step size (need to be sufficiently small)
241 :param u0: initial value
242 :param v0: initial velocity
243 :param t0: initial time
244 """
245 self.u=u0
246 self.v=v0
247 self.t=t0
248 self.t_last=t0
249 self.__dt=dt
250
251 def getTimeStepSize(self):
252 return self.__dt
253
254 def _getAcceleration(self, t, u):
255 """
256 returns the acceleraton for time t and solution u at time t
257 : note: this function needs to be overwritten by a particular wave problem
258 """
259 pass
260
261 def update(self, t):
262 """
263 returns the solution for the next time marker t which needs to greater than the time marker from the
264 previous call.
265 """
266 if not self.t_last <= t:
267 raise ValueError("You can not move backward in time.")
268
269 dt = self.getTimeStepSize()
270 # apply Verlet scheme until
271 while self.t < t:
272 a=self._getAcceleration(self.t, self.u)
273 self.v += dt * a
274 self.u += dt * self.v
275 self.t += dt
276
277 # now we work backwards
278 self.t_last=t
279 return t, self.u + self.v * (t-self.t)
280
281 def createAbsorbtionLayerFunction(x, absorption_zone=300*U.m, absorption_cut=1.e-2, top_absorbation=False):
282 """
283 creating a distribution which is one in the interior of the domain of x
284 and is falling down to the value 'absorption_cut' over a margain of thickness 'absorption_zone'
285 toward each boundary except the top of the domain.
286
287 :param x: location of points in the domain
288 :type x: `Data`
289 :param absorption_zone: thickness of the aborption zone
290 :param absorption_cut: value of decay function on domain boundary
291 :return: function on 'x' which is one in the iterior and decays to almost zero over a margin
292 toward the boundary.
293 """
294 dom=x.getDomain()
295 bb=boundingBox(dom)
296 DIM=dom.getDim()
297 decay=-log(absorption_cut)/absorption_zone**2
298 f=1
299 for i in range(DIM):
300 x_i=x[i]
301 x_l=x_i-(bb[i][0]+absorption_zone)
302 m_l=whereNegative(x_l)
303 f=f*( (exp(-decay*(x_l*m_l)**2)-1) * m_l+1 )
304 if top_absorbation or not DIM-1 == i:
305 x_r=(bb[i][1]-absorption_zone)-x_i
306 m_r=whereNegative(x_r)
307 f=f*( (exp(-decay*(x_r*m_r)**2)-1) * m_r+1 )
308 return f
309
310 class SonicWave(WaveBase):
311 """
312 Solving the sonic wave equation
313
314 p_tt = (v_p**2 * p_i)_i + f(t) * delta_s where (p-) velocity v_p.
315
316 f(t) is wavelet acting at a point source term at positon s
317 """
318 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):
319 """
320 initialize the sonic wave solver
321
322 :param domain: domain of the problem
323 :type domain: `Domain`
324 :param v_p: p-velocity field
325 :type v_p: `Scalar`
326 :param wavelet: wavelet to describe the time evolution of source term
327 :type wavelet: `Wavelet`
328 :param source_tag: tag of the source location
329 :type source_tag: 'str' or 'int'
330 :param dt: time step size. If not present a suitable time step size is calculated.
331 :param p0: initial solution. If not present zero is used.
332 :param p0_t: initial solution change rate. If not present zero is used.
333 :param absorption_zone: thickness of absorption zone
334 :param absorption_cut: boundary value of absorption decay factor
335 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
336 """
337 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
338 v_p=v_p*f
339
340 if p0 == None:
341 p0=Scalar(0.,Solution(domain))
342 else:
343 p0=interpolate(p0, Solution(domain ))
344
345 if p0_t == None:
346 p0_t=Scalar(0.,Solution(domain))
347 else:
348 p0_t=interpolate(p0_t, Solution(domain ))
349
350 if dt == None:
351 dt=min(inf((1./5.)*domain.getSize()/v_p), wavelet.getTimeScale())
352
353 super(SonicWave, self).__init__( dt, u0=p0, v0=p0_t, t0=0.)
354
355 self.__wavelet=wavelet
356 self.__mypde=LinearSinglePDE(domain)
357 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
358 self.__mypde.setSymmetryOn()
359 self.__mypde.setValue(D=1./v_p**2)
360 self.__source_tag=source_tag
361 self.__r=Scalar(0., DiracDeltaFunctions(self.__mypde.getDomain()))
362
363
364 def _getAcceleration(self, t, u):
365 """
366 returns the acceleraton for time t and solution u at time t
367 """
368 self.__r.setTaggedValue(self.__source_tag, self.__wavelet.getAcceleration(t))
369 self.__mypde.setValue(X=-grad(u,Function(self.__mypde.getDomain())), y_dirac= self.__r)
370 return self.__mypde.getSolution()
371
372
373 class VTIWave(WaveBase):
374 """
375 Solving the VTI wave equation
376
377 :note: In case of a two dimensional domain the second spatial dimenion is depth.
378 """
379 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
380 source_vector = [0.,0.,1.], eps=0., gamma=0., delta=0., rho=1.,
381 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
382 absorption_cut=1e-2, lumping=True):
383 """
384 initialize the VTI wave solver
385
386 :param domain: domain of the problem
387 :type domain: `Domain`
388 :param v_p: vertical p-velocity field
389 :type v_p: `Scalar`
390 :param v_s: vertical s-velocity field
391 :type v_s: `Scalar`
392 :param wavelet: wavelet to describe the time evolution of source term
393 :type wavelet: `Wavelet`
394 :param source_tag: tag of the source location
395 :type source_tag: 'str' or 'int'
396 :param source_vector: source orientation vector
397 :param eps: first Thompsen parameter
398 :param delta: second Thompsen parameter
399 :param gamma: third Thompsen parameter
400 :param rho: density
401 :param dt: time step size. If not present a suitable time step size is calculated.
402 :param u0: initial solution. If not present zero is used.
403 :param v0: initial solution change rate. If not present zero is used.
404 :param absorption_zone: thickness of absorption zone
405 :param absorption_cut: boundary value of absorption decay factor
406 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
407 """
408 DIM=domain.getDim()
409 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
410
411 v_p=v_p*f
412 v_s=v_s*f
413
414
415
416 if u0 == None:
417 u0=Vector(0.,Solution(domain))
418 else:
419 u0=interpolate(p0, Solution(domain ))
420
421 if v0 == None:
422 v0=Vector(0.,Solution(domain))
423 else:
424 v0=interpolate(v0, Solution(domain ))
425
426 if dt == None:
427 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
428
429 super(VTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
430
431 self.__wavelet=wavelet
432
433 self.fastAssembler = hasattr(domain, "setAssembler")
434
435 if self.fastAssembler:
436 c33=v_p**2 * rho
437 c44=v_s**2 * rho
438 c11=(1+2*eps) * c33
439 c66=(1+2*gamma) * c44
440 c13=sqrt(2*c33*(c33-c44) * delta + (c33-c44)**2)-c44
441 c12=c11-2*c66
442 self.__mypde=VTIWavePDE(domain, [c11, c12, c13, c33, c44, c66])
443 else:
444 self.__mypde=LinearPDESystem(domain)
445 self.__mypde.setValue(X=self.__mypde.createCoefficient('X'))
446
447 if lumping:
448 self.__mypde.getSolverOptions().setSolverMethod(
449 self.__mypde.getSolverOptions().HRZ_LUMPING)
450 self.__mypde.setSymmetryOn()
451 self.__mypde.setValue(D=rho*kronecker(DIM))
452 self.__source_tag=source_tag
453
454 if DIM ==2 :
455 source_vector= [source_vector[0],source_vector[2]]
456
457 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
458 self.__r.setTaggedValue(self.__source_tag, source_vector)
459
460 if not self.fastAssembler:
461 self.c33=v_p**2 * rho
462 self.c44=v_s**2 * rho
463 self.c11=(1+2*eps) * self.c33
464 self.c66=(1+2*gamma) * self.c44
465 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
466 self.c12=self.c11-2*self.c66
467
468 def _getAcceleration(self, t, u):
469 """
470 returns the acceleraton for time t and solution u at time t
471 """
472 du = grad(u)
473 if not self.fastAssembler:
474 sigma=self.__mypde.getCoefficient('X')
475 if self.__mypde.getDim() == 3:
476 e11=du[0,0]
477 e22=du[1,1]
478 e33=du[2,2]
479 sigma[0,0]=self.c11*e11+self.c12*e22+self.c13*e33
480 sigma[1,1]=self.c12*e11+self.c11*e22+self.c13*e33
481 sigma[2,2]=self.c13*(e11+e22)+self.c33*e33
482
483 s=self.c44*(du[2,1]+du[1,2])
484 sigma[1,2]=s
485 sigma[2,1]=s
486
487 s=self.c44*(du[2,0]+du[0,2])
488 sigma[0,2]=s
489 sigma[2,0]=s
490
491 s=self.c66*(du[0,1]+du[1,0])
492 sigma[0,1]=s
493 sigma[1,0]=s
494
495
496 else:
497 e11=du[0,0]
498 e22=du[1,1]
499 sigma[0,0]=self.c11*e11+self.c13*e22
500 sigma[1,1]=self.c13*e11+self.c33*e22
501 s=self.c44*(du[1,0]+du[0,1])
502 sigma[0,1]=s
503 sigma[1,0]=s
504 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
505
506 else:
507 self.__mypde.setValue(du=du, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
508
509 return self.__mypde.getSolution()
510
511
512 class HTIWave(WaveBase):
513 """
514 Solving the HTI wave equation (along the x_0 axis)
515
516 :note: In case of a two dimensional domain a horizontal domain is considered, i.e. the depth component is dropped.
517 """
518
519 def __init__(self, domain, v_p, v_s, wavelet, source_tag,
520 source_vector = [1.,0.,0.], eps=0., gamma=0., delta=0., rho=1.,
521 dt=None, u0=None, v0=None, absorption_zone=300*U.m,
522 absorption_cut=1e-2, lumping=True):
523 """
524 initialize the VTI wave solver
525
526 :param domain: domain of the problem
527 :type domain: `Domain`
528 :param v_p: vertical p-velocity field
529 :type v_p: `Scalar`
530 :param v_s: vertical s-velocity field
531 :type v_s: `Scalar`
532 :param wavelet: wavelet to describe the time evolution of source term
533 :type wavelet: `Wavelet`
534 :param source_tag: tag of the source location
535 :type source_tag: 'str' or 'int'
536 :param source_vector: source orientation vector
537 :param eps: first Thompsen parameter
538 :param delta: second Thompsen parameter
539 :param gamma: third Thompsen parameter
540 :param rho: density
541 :param dt: time step size. If not present a suitable time step size is calculated.
542 :param u0: initial solution. If not present zero is used.
543 :param v0: initial solution change rate. If not present zero is used.
544 :param absorption_zone: thickness of absorption zone
545 :param absorption_cut: boundary value of absorption decay factor
546 :param lumping: if True mass matrix lumping is being used. This is accelerates the computing but introduces some diffusion.
547 """
548 DIM=domain.getDim()
549 f=createAbsorbtionLayerFunction(Function(domain).getX(), absorption_zone, absorption_cut)
550
551 v_p=v_p*f
552 v_s=v_s*f
553
554 if u0 == None:
555 u0=Vector(0.,Solution(domain))
556 else:
557 u0=interpolate(p0, Solution(domain ))
558
559 if v0 == None:
560 v0=Vector(0.,Solution(domain))
561 else:
562 v0=interpolate(v0, Solution(domain ))
563
564 if dt == None:
565 dt=min((1./5.)*min(inf(domain.getSize()/v_p), inf(domain.getSize()/v_s)), wavelet.getTimeScale())
566
567 super(HTIWave, self).__init__( dt, u0=u0, v0=v0, t0=0.)
568
569 self.__wavelet=wavelet
570
571 self.__mypde=LinearPDESystem(domain)
572 if lumping: self.__mypde.getSolverOptions().setSolverMethod(self.__mypde.getSolverOptions().HRZ_LUMPING)
573 self.__mypde.setSymmetryOn()
574 self.__mypde.setValue(D=rho*kronecker(DIM), X=self.__mypde.createCoefficient('X'))
575 self.__source_tag=source_tag
576
577 if DIM ==2 :
578 source_vector= [source_vector[0],source_vector[2]]
579
580 self.__r=Vector(0, DiracDeltaFunctions(self.__mypde.getDomain()))
581 self.__r.setTaggedValue(self.__source_tag, source_vector)
582
583 self.c33=v_p**2 * rho
584 self.c44=v_s**2 * rho
585 self.c11=(1+2*eps) * self.c33
586 self.c66=(1+2*gamma) * self.c44
587 self.c13=sqrt(2*self.c33*(self.c33-self.c44) * delta + (self.c33-self.c44)**2)-self.c44
588 self.c23=self.c33-2*self.c66
589
590 def _getAcceleration(self, t, u):
591 """
592 returns the acceleraton for time t and solution u at time t
593 """
594 du = grad(u)
595 sigma=self.__mypde.getCoefficient('X')
596
597 if self.__mypde.getDim() == 3:
598 e11=du[0,0]
599 e22=du[1,1]
600 e33=du[2,2]
601
602 sigma[0,0]=self.c11*e11+self.c13*(e22+e33)
603 sigma[1,1]=self.c13*e11+self.c33*e22+self.c23*e33
604 sigma[2,2]=self.c13*e11+self.c23*e22+self.c33*e33
605
606 s=self.c44*(du[2,1]+du[1,2])
607 sigma[1,2]=s
608 sigma[2,1]=s
609
610 s=self.c66*(du[2,0]+du[0,2])
611 sigma[0,2]=s
612 sigma[2,0]=s
613
614 s=self.c66*(du[0,1]+du[1,0])
615 sigma[0,1]=s
616 sigma[1,0]=s
617
618 else:
619 e11=du[0,0]
620 e33=du[1,1]
621 sigma[0,0]=self.c11*e11+self.c13*e33
622 sigma[1,1]=self.c13*e11+self.c33*e33
623
624 s=self.c66*(du[1,0]+du[0,1])
625 sigma[0,1]=s
626 sigma[1,0]=s
627
628 self.__mypde.setValue(X=-sigma, y_dirac= self.__r * self.__wavelet.getAcceleration(t))
629 return self.__mypde.getSolution()
630

  ViewVC Help
Powered by ViewVC 1.1.26