- Email: [email protected]
Extending the stability limit of explicit scheme with spatial filtering for solving wave equations
Journal Pre-proof Extending the stability limit of explicit scheme with spatial filtering for solving wave equations
Yingjie Gao, Jinhai Zhang, Zhenxing Yao
PII:
S0021-9991(19)30537-6
DOI:
https://doi.org/10.1016/j.jcp.2019.07.051
Reference:
YJCPH 8853
To appear in:
Journal of Computational Physics
Received date:
24 September 2018
Revised date:
27 May 2019
Accepted date:
22 July 2019
Please cite this article as: Y. Gao et al., Extending the stability limit of explicit scheme with spatial filtering for solving wave equations, J. Comput. Phys. (2019), doi: https://doi.org/10.1016/j.jcp.2019.07.051.
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.
Highlights • • • • •
We extend the stability limit of the explicit scheme using spatial filtering. Filtering out the high-wavenumber noise before it boosts can guarantee the stability. We further reduce the time-dispersion error by inverse time dispersion transform. The waveform is highly accurate compared with the reference waveforms. The time step is much larger than the maximum time step allowed by the CFL stability limit.
1
Extending the stability limit of explicit scheme with spatial
2
filtering for solving wave equations
3 4
Yingjie Gaoa, b , Jinhai Zhanga, b, *, and Zhenxing Yaoa, b
5
a Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics,
6
Chinese Academy of Sciences, Beijing, China
7
b Institutions of Earth Science, Chinese Academy of Sciences, Beijing, China
8
Running Head: Extending stability limit by spatial filtering
9 10 11
Submitted to Journal of Computational Physics on Sep 24, 2018
12
Revised manuscript submitted on Jan 24, 2019
13
Revised manuscript submitted on May 27, 2019
14
*Corresponding author: Dr. Jinhai Zhang
15
Chinese Academy of Sciences
16
Institute of Geology and Geophysics
17
Beijing 100029
18
China.
19
E-mail: [email protected]
20
Phone number: +86-10-82998013 1
21 22
Extending the stability limit of explicit scheme with spatial
23
filtering for solving wave equations
24
ABSTRACT
25
The explicit finite-difference scheme is widely used for solving the wave equation;
26
however, its time step is strictly constrained by the Courant-Friedrichs-Lewy (CFL)
27
stability limit. We apply spatial filtering to reduce the high-wavenumber noise caused by
28
using a larger time step beyond the CFL stability limit. We further incorporate the
29
forward and inverse time dispersion transforms to reduce the time-dispersion error. This
30
method allows us to use a time step larger than the maximum time step allowed by the
31
CFL stability limit.
32 33
.
34 35 36 37 38 39 40
Key words: explicit scheme; spatial filtering; CFL stability limit; time-dispersion error; large time step.
41 2
42 43
1. Introduction
44
The explicit scheme, i.e., iteratively updating the wavefields along the temporal
45
direction in the time domain, is the most popular methods for simulating seismic
46
wavefields because of its simplicity in numerical implementation [1]. However, its
47
maximum time step is strictly restricted by the Courant-Friedrichs-Lewy (CFL) stability
48
limit [2], which leads to a heavy computational burden in the presence of small-scale
49
structures and high-velocity targets. We have to adopt a small enough time step to
50
guarantee that the updated wavefields after each time step do not exceed the spatial range
51
of one cell; otherwise, we could encounter stability problems. The CFL stability limit
52
leads to a great computational burden, especially for imaging fine-structure targets (e.g.,
53
channels and fractures) using a small spatial interval. For example, when the spatial
54
interval is reduced by half, the computational cost for a 3D model would become 16
55
times that for the original spatial interval. We should develop new methods to overcome
56
the stability limit with an acceptable increase of computational cost when imaging
57
fine-structure targets.
58
In the electromagnetic wave simulation field, Li et al. [3] extended the CFL stability
59
limit by perturbing the unstable eigenvalues to be stable, which is called eigenvalue
60
perturbation [4]. Gao et al. [5] cooperated the forward and inverse time dispersion
61
transform methods [6-8] into the eigenvalue perturbation to solve the wave equation in
3
62
the seismic wave simulation field. As a preprocessing algorithm, the eigenvalue
63
perturbation does not increase the amount of computation during the iteration process.
64
In fluid dynamics, Ecer et al. [9] found that a larger time step beyond the CFL
65
stability limit would introduce unstable high-wavenumber wavefields; thus, one can filter
66
out the unstable high-wavenumber wavefields to enlarge the time step of an explicit
67
scheme. Sarris [10] proposed a spatial filtering algorithm to deal with the unstable
68
problems caused by a large time step when solving Maxwell’s equations. Using a
69
low-pass filter to the unstable components distributed in high-wavenumbers, the time step
70
of an explicit scheme can be successfully extended beyond the CFL stability limit
71
[11-13].
72
In this paper, we extend the CFL stability limit of the explicit scheme to solve the
73
wave equation. Proper low-pass spatial filtering [10] is used to reduce the slightly
74
boosted amplitude in each time step. The presented scheme allows us to apply a much
75
larger time step beyond the CFL stability limit. The spatial filtering can be perfectly
76
incorporated into the pseudospectral method and its computational cost can be neglected
77
compared with that of the Fourier transforms. To reduce the time-dispersion error caused
78
by using a large time step beyond the CFL stability limit, we incorporate the forward time
79
dispersion-transform (FTDT) method and the inverse time dispersion-transform (ITDT)
80
method [6-8]. The maximum time step in our test is around 2-3 times larger than the CFL
81
stability limit allowed. The memory demand is almost the same compared with that of the
4
82
classic pseudospectral method and the generated wavefield is shown to be highly
83
accurate.
84 85 86
2. Methodology We start with the 2D scalar wave equation
1 ∂ 2u ∂ 2u ∂ 2u = + , c 2 ∂t 2 ∂x 2 ∂z 2
87
(1)
88
where u ( x, z, t ) represents the particle displacement at position
( x, z )
89
c ( x, z ) is the velocity. We use the pseudospectral (PS) method for the spatial
90
discretization [14] since its error mainly arises from the temporal discretization, which is
91
helpful for identifying the errors caused by using larger time steps without considering
92
the space-dispersion error. The PS method for equation (1) reads
and time t, and
93
∂ 2u = c 2 F − ª¬ − k 2 F + ( u ) º¼ , ∂t 2
94
where k 2 = k x2 + k z2 , with kx and kz being the wavenumbers along the x- and z-
95
directions, respectively; F+ and F− are 2D forward and inverse Fourier transforms,
96
respectively.
97 98
{
}
(2)
Based on the Taylor expansion, a general time discretization scheme for the temporal derivative in equation (1) can be obtained by the following formula [15, 16] j
99
J Δt 2 j − 2 § ∂ 2 u n · u n +1 − 2u n + u n −1 = 2 ¦ ¨ 2 ¸ , Δt 2 j =1 ( 2 j ) ! © ∂t ¹
5
(3)
100
where 2J is the order of the temporal finite-difference (FD) scheme, and u n represents
101
the wavefield at time nΔt , and Δt is the time step. In practice, we only use a limited
102
number of terms on the right-hand side of equation (3) by the Lax-Wendroff methods,
103
which use spatial derivatives to replace high-order temporal derivatives [15, 17] 2 n u n +1 − 2u n + u n −1 ∂ 2u n 2§∂ u c = + ¨ 2 Δt 2 ∂z 2 © ∂x
104
j
J · Δt 2 j − 2 § ∂ 2 u n · 2 + ¸ ¦ ¨ 2 ¸ , j = 2 ( 2 j ) ! © ∂t ¹ ¹
(4)
105
and the higher-order temporal derivatives are obtained by the following recursive formula
106
[17] j
2 n § ∂ 2u n · ∂ 2u n · § ∂ 2u n · 2§∂ u c = + ¨ 2 ¸ ¨ 2 ¸¨ ¸ ∂ z 2 ¹ © ∂t 2 ¹ © ∂t ¹ © ∂x
107
j −1
.
(5)
108
The second-order (2nd-order) temporal FD scheme (i.e., J=1) and the fourth-order
109
(4th-order) temporal FD scheme (i.e., J=2) are commonly used in numerical simulations
110
[15-17]. The CFL stability limit for the 2D PS method with the 2nd-order temporal FD scheme
111 112
(i.e., PS (2nd-order)) and the 4th-order temporal FD scheme(i.e., PS (4th-order)) [14, 17]
( cmax k Δt )
113 114
115
2
≤1
4
(6)
and
( c k Δt ) 0 ≤ max 4
2
( c k Δt ) − max 48
4
≤ 1ˈ
(7)
116
respectively; and cmax is the maximum velocity. According to the CFL stability limit, the
117
time step Δ t must be limited by Δt ≤ 2 ( cmax kmax ) and Δt ≤ 2 3 ( cmax kmax ) for PS
118
(2nd-order) and PS (4th-order) schemes, respectively. Using a velocity c=4000 m/s and a 6
119
spatial grid interval ǻx=ǻz=h=10 m, we can obtain the maximum time step Δtmax ,
120
which corresponds to the maximum wavenumber kmax = 2π h (with k xmax =k zmax = π h )
121
for the 2D case, as shown in Table 1.
122
Fig. 1 shows the critical curves of the CFL stability limit for equation (6) and
123
equation (7). Obviously, the CFL stability limit is to make the time step not exceed the
124
maximum time step Δtmax , which corresponds to the k max , as shown by the shaded area
125
(or cyan area) in Fig. 1. For a given time step Δt that is larger than Δtmax , the
126
corresponding wavenumber kmax can be obtained from the critical curve. For the PS
127
(2nd-order)
128
kmax = 2 ( cmax Δ t ) and kmax = 2 3 ( cmax Δt ) , respectively.
and
the
PS
(4th-order)
schemes,
kmax
can
be
expressed
as
129
We can see that the wavefield components with a wavenumber smaller than kmax
130
are still sufficient for the CFL stability limit, as shown in Table 1. The kmax gradually
131
decreases with an increasing time step, which means that there is a possibility that we
132
could use a larger time step by limiting the range of wavenumbers using a proper
133
low-pass filter defined as [10, 13]
134
°1, F (k ) = ® °¯0,
for k x2 + k z2 ≤ kmax , otherwise,
(8)
135
where kmax is the corner frequency of the low-pass filter. The spatial filtering algorithm
136
is implemented by the following stages per time step [10]:
7
137
n+1 (1) After wavefield u has been updated at time
( n + 1) Δt
using equation (4), we
138
transform the wavefield to the wavenumber domain using a 2D Fourier transform
139
U n+1 =F + ( u n+1 ) .
140 141 142 143
n+1 =F ( k ) ⋅U n+1 in the wavenumber (2) Compute the spatially filtered wavefield U domain, where ⋅ is a Hadamard product. (3) Apply inverse 2D Fourier transform to obtain the filtered wavefield in the space
(
)
domain un+1 =F − U n+1 , and un +1 is taken as the input wavefield for the next iteration.
144
The spatial filtering algorithm can be incorporated perfectly into the implementation
145
of the PS method, since the PS method is implemented by shuttling the wavefield
146
frequently between the space and wavenumber domains. We can directly apply the spatial
147
filtering in the wavenumber domain together with the operation of wavenumber −k 2 .
148
Therefore, the computational cost of the PS method after applying the spatial filtering for
149
each time step could be almost perfectly retained, since the time consumed by the
150
Hadamard product associated with the spatial filtering can be neglected compared with
151
that of the Fourier transforms.
152 153
3. Numerical experiments
154
We perform numerical experiments on a homogeneous model. Two schemes are
155
tested: the PS (2nd-order) and PS (4th-order) schemes. The wave velocity is c=4000 m/s.
156
The spatial grid interval is ǻx=ǻz=10 m, and the grid number is 401×401. The source is a
157
Ricker wavelet with a dominant frequency of 15 Hz (shown in Appendix A) [18], which 8
158
is located at x=0 and z=0 m. The maximum time steps for the PS (2nd-order) and PS
159
(4th-order) schemes are 1.125 and 1.949 ms, respectively. We tested several large time
160
steps (ǻt=3, 5, 7, 9, and 11 ms) for each scheme. Note that all these time steps are beyond
161
the CFL stability limit. As a preparation, we filtered the source time functions at the
162
corresponding time steps by the FTDT method before the simulation [10]. The result
163
obtained by the PS (4th-order) scheme with a time step ǻt =1 ms was taken as a
164
theoretical reference to examine the accuracy of the waveforms obtained by different
165
time steps.
166
Fig. 2a1 and 2a3 show the snapshots at 333 ms using ǻt =3 ms for the PS (2nd-order)
167
scheme and the PS (4th-order) scheme, respectively. To illustrate the unstable phenomena
168
in the snapshot, the spatial filtering algorithm is stopped after 300 ms. Obviously, many
169
unstable components appear in the snapshot at 333 ms. Fig. 2b1 and 2b3 show the
170
corresponding wavenumber snapshots of Fig. 2a1 and 2a3, respectively. The unstable
171
components are mainly distributed in the high-wavenumber regions. As shown in Fig.
172
2a2, 2b2, 2a4, and 2b4, the spatial filtering can remove the unstable components at high
173
wavenumbers while retaining the effective stable wavefield at low-wavenumbers. Fig. 2c
174
and 2d show the profiles along the red dashed lines in Fig. 2a1 and 2b1, respectively. Fig.
175
2d clearly shows that the unstable components are distributed in a high-wavenumber
176
region, which can be clearly distinguished from the effective low-wavenumber wavefield
177
components when the time step is not much larger than the upper-limit time step Δtmax .
9
178
Fig. 3a shows the waveforms obtained by different schemes using different time
179
steps, recorded at x=í800 and z=í800 m. Fig. 3b shows the waveforms after applying the
180
ITDT method. Both the PS (2nd-order) scheme and the PS (4th-order) scheme after
181
applying the ITDT method are fairly accurate up to time step ǻt = 7 ms (even ǻt = 9 ms
182
for the PS (4th-order) scheme). This indicates that the combination of the spatial filtering
183
and the ITDT method could provide highly accurate simulation results even though a
184
much larger time step beyond the CFL stability limit is used. The PS (2nd-order) scheme
185
using ǻt 9 ms and the PS (4th-order) scheme using ǻt = 11 ms still have some residual
186
errors even after applying the ITDT method. At these cases, the corner frequency of the
187
low-pass filter kmax is smaller than the maximum wavenumber of the effective stable
188
wavefield. For example, kmax is 0.078729 for the PS (4th-order) scheme using ǻt = 11
189
ms (as shown in Table 1), while the maximum wavenumber of the effective stable
190
wavefield is k § 0.08 (as shown in Fig. 2d). This means that we could not separate the
191
effective wavefield and the unstable components using spatial filtering since they are
192
completely mixed into each other.
193
Based on the above constant velocity model, we perform two additional sets of
194
numerical experiments using two other wavelets (Klauder wavelet and Ormsby wavelet,
195
as shown in Appendix A) [18], respectively. The mid-frequency of both wavelets is 15 Hz.
196
We tested several large time steps (ǻt=3, 5, 7, 9, 11, and 13 ms). Fig. 4a and Fig. 4b show
197
the waveforms after applying the ITDT method to the wavefields obtained using the
198
Klauder wavelet and the Ormsby wavelet, respectively. The wavefoms are both stable and 10
199
accurate when using time steps beyond the CFL stability limit up to ǻt=11 ms, and the
200
result obtained using the Ormsby wavelet still has high accuracy even using time step
201
ǻt=13ms. This is due to the effective frequency band of the Ormsby wavelet is the
202
narrowest (the maximum value of the effective frequency is the smallest, as shown in Fig.
203
A1b) among these three wavelets, which means a smaller corner frequency kmax (i.e. a
204
larger time step ǻt) can be used without worrying about the effective wave field would be
205
filtered out.We verify the feasibility of the spatial filtering algorithm with the ITDT
206
methods by a modified Marmousi model, as shown in Fig. 5a. The spatial grid interval is
207
ǻx=ǻz=10 m, and the grid number is 737×751. The source time function is a Ricker
208
wavelet with a dominant frequency of 15 Hz, and the source is located at x=3200 and
209
z=3000 m. We record a trace at x=5600 and z=5000 m for different time steps, as shown
210
in Fig. 5b. Fig. 5c shows the local zoom in of the waveforms in the box of Fig. 5b.
211
Obviously, both the PS (2nd-order) scheme and the PS (4th-order) scheme, using the
212
spatial filtering with the ITDT methods, can obtain acceptable results using a much larger
213
time step than the maximum time steps allowed by the CFL stability limit. However, the
214
residual errors become evident at much larger time steps at ǻt = 7 ms.
215 216
4. Discussion
217
As shown in Fig. 2b1 and 2b3, we could easily separate unstable wavefield
218
components and the effective wavefield by setting a proper filter when the actual time
219
step is not much larger than the maximum time step allowed by the CFL stability limit. 11
220
For a model that has a strong velocity contrast between the maximum and minimum
221
velocities (i.e., heterogeneity is arbitrary), the wavefield corresponding to the low
222
velocity region would be distributed in higher wavenumbers, even larger than the corner
223
frequency. In other words, the unstable wavefield components are mixed with the
224
effective wavefield components when the velocity contrast is too strong; in such cases,
225
any attempt to filter the unstable components out would ruin the effective wavefields.
226
Therefore, we map the velocity of the original Marmousi model onto 2500 to 4000 m/s to
227
illustrate the feasibility of the proposed method on moderate variation of the velocity.
228
Future works are still needed to work on strong variation of the velocity.
229
With the spatial filtering algorithm, we can now use a much larger time step beyond
230
the CFL stability limit; however, the time step cannot be infinitely large. A larger time
231
step corresponds to a smaller corner frequency of the low-pass filter, which may filter out
232
the effective wavefield components. For example, with c=4000 m/s and h=10 m, the
233
corresponding corner frequency of ǻt = 9 ms for the PS (2nd-order) scheme is 1/18;
234
unfortunately, the low-pass filter at 1/16 would ruin the effective wavefield components,
235
as shown in Fig. 2d.
236
Another issue needs to emphasize is the Nyquist limit in temporal discretization. For
237
example, with c=4000 m/s and h=100 m, the upper-limit time step Δtmax for the PS
238
(4nd-order) scheme is 19.49 ms. The maximum frequency f max of the Ricker wavelet
239
with a dominant frequency of 15 Hz is about 45 Hz (as shown in Fig. A1b), and the time
12
240
step cannot be larger than 11.11 ms according to the Nyquist limit Δt < 1 ( 2 f max ) . In
241
other words, the time step reaches the Nyquist sampling limit first, while it has not
242
reached the upper limit of the stability condition yet. Fortunately, Gao et al. [5]
243
demonstrated that the time dispersion transform methods can guarantee the simulation
244
accuracy using a time step even towards the Nyquist limit. Therefore, in this case, one
245
only needs to apply the time dispersion transforms to ensure the simulation accuracy of
246
the simulation using a large time step simulation, and it is not necessary to consider using
247
spatial filtering method to extend the CFL stability limit.
248
To illustrate the unstable wavefield components, we do not apply the spatial filtering
249
between 300 and 333 ms, as shown in Figs. 2a1 and 2a3. In other words, the spatial
250
filtering is not used for continuous ten time steps, and the instability of the high
251
wavenumber can still be filtered out, as shown in Figs. 2a2 and 2a4. This indicates that it
252
is not necessary to apply the spatial filtering for each time step. We can save some
253
computational cost at several other time steps when the spatial discretization is purely in
254
the spatial domain.
255 256
5. Conclusion
257
We introduce the spatial filtering algorithm to guarantee the stability of the explicit
258
scheme for solving the wave equation with a much larger time step beyond the
259
Courant-Friedrichs-Lewy (CFL) stability limit. The computational cost is kept low since
260
the spatial filtering can be perfectly incorporated into the pseudospectral method. 13
261
Furthermore, we apply the forward time dispersion transform method and the inverse
262
time dispersion transform method to reduce the time-dispersion error caused by the large
263
time step. Both theoretical analyses and numerical experiments show that the proposed
264
method can use a larger time step than that allowed by the CFL stability limit, which is
265
helpful for saving a large number of iteration times without suffering from the time
266
dispersion and the stability problems. The proposed method retains the superiority of the
267
conventional explicit scheme by avoiding solving matrices but breaks through its major
268
limitation on the time step.
269 270
Acknowledgements
271
We are grateful to Erik Koene for helpful discussion, derivation of source loading for
272
the PS (4th-order) scheme, and the corresponding MATLAB code. This research was
273
supported by the National Natural Science Foundation of China (41704063) and by the
274
National Major Project of China (2017ZX05008-007). YG was also supported by the
275
General Financial Grant from the China Postdoctoral science foundation (2017M610980).
276
JZ thanks the supports from the Strategic Pioneer Program on Space Science, Chinese
277
Academy of Sciences (XDA15011700) and the Foundation for Excellent Member of the
278
Youth Innovation Promotion Association, Chinese Academy of Sciences (2016).
279 280
References 14
281 282 283 284 285 286 287 288 289 290 291 292
[1] J.T. Etgen, M.J. O'Brien, Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial, Geophysics, 72 (2007) SM223–SM230. [2] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen differenzengleichungen der mathematischen physik, Mathematische Annalen (In German), 100 (1928) 32–74. [3] X. Li, C.D. Sarris, P. Triverio, Overcoming the FDTD stability limit via model order reduction and eigenvalue perturbation, in:
Microwave Symposium, 2014, pp. 1–3.
[4] X. Li, Model order reduction and stability enforcement of finite-difference time-domain equations beyond the CFL limit, in, University of Toronto (Canada), 2014. [5] Y. Gao, J. Zhang, Z. Yao, Removing the stability limit of the explicit finite-difference scheme with eigenvalue perturbation, Geophysics, 83 (2018) A93–A98. [6] M. Wang, S. Xu, Finite-difference time dispersion transforms for wave propagation, Geophysics, 80 (2015) WD19–WD25.
293
[7] Y. Gao, J. Zhang, Z. Yao, Third-order symplectic integration method with inverse time
294
dispersion transform for long-term simulation, Journal of Computational Physics, 314 (2016)
295
436–449.
296 297
[8] E.F.M. Koene, J.O.A. Robertsson, F. Broggini, F. Andersson, liminating time dispersion from seismic wave modelling, Geophysical Journal International, 213 (2018) 169–180.
298
[9] A. Ecer, N. Gopalaswamy, H. Akay, Y. Chien, Digital filtering techniques for parallel
299
computation of explicit schemes, International Journal of Computational Fluid Dynamics, 13
300
(2000) 211–222.
15
301 302
[10] C.D. Sarris, Extending the stability limit of the FDTD method with spatial filtering, IEEE Microwave and Wireless Components Letters, 21 (2011) 176–178.
303
[11] C. Chang, C.D. Sarris, A spatial filter-enabled high-resolution subgridding scheme for stable
304
FDTD modeling of multiscale geometries, in: international microwave symposium, 2011,
305
pp. 1–4.
306 307
[12] C. Chang, C.D. Sarris, A three-dimensional spatially filtered FDTD with controllable stability beyond the courant limit, in:
Microwave Symposium Digest, 2012, pp. 1–3.
308
[13] C. Chang, C.D. Sarris, A spatially filtered finite-difference time-domain scheme with
309
controllable stability beyond the CFL limit: theory and applications, IEEE Transactions on
310
Microwave Theory & Techniques, 61 (2013) 351–359.
311 312 313 314 315 316 317 318 319 320
[14] D.D. Kosloff, E. Baysal, Forward modeling by a Fourier method, Geophysics, 47 (1982) 1402–1412. [15] M.A. Dablain, The application of high-order differencing to the scalar wave equation, Geophysics, 51 (1986) 54–66. [16] J.T. Etgen, Accurate wave equation modeling, Stanford Exploration Project, 60 (1988) 131–147. [17] J.B. Chen, High-order time discretizations in seismic modeling, Geophysics, 72 (2007) SM115–SM122. [18] H.Ryan, Ricker, Ormsby, Klauder, Butterworth—a choice of wavelets. CSEG Recorder, 19 (1994) 8-9.
321 16
322
APPENDIX A: THREE TYPES OF WAVELETS
323 324
Ricker wavelet can be expressed as [18] 2 −π Ricker ( t ) = ª1 − 2π 2 f 2 ( t − t0 ) º e ¬ ¼
325
2
f 2 ( t − t0 )
2
,
(A-1)
326
where f =15 Hz is the dominant frequency, and t0 = 0.15 s is the time delay of the
327
wavelet peak.
328 329
Klauder wavelet can be expressed as [18]
° sin ªπ kt (T − t ) º¼ −2iπ f0 (t −t0 ) °½ Klauder ( t ) = real ® ¬ e ¾ π kt ¯° ¿° ,
(A-2)
330
where T =7 is the duration of the input signal. k = ( f 2 − f 2 ) T is the rate of change
331
with time, with f1 =5 Hz being the terminal low frequency, and f2 =25 Hz being the
332
terminal high frequency. f 0 = ( f 2 +f 2 ) T is the mid-frequency of bandwidth.
333
334
Ormsby wavelet can be expressed as [18] ° (π f 4 )2 kt sinc 2 ª f 4 (T − t ) º (π f3 )2 kt sinc 2 ª f3 (T − t ) º ½° ¬ ¼− ¬ ¼ − Ormsby ( t ) = ® ¾ (π f 4 − π f 3 ) (π f 4 − π f3 ) ¯° ¿° ° (π f 2 )2 kt sinc 2 ª f 2 (T − t ) º (π f1 )2 kt sinc 2 ª f1 (T − t ) º ½° ¬ ¼− ¬ ¼ , ® ¾ π f π f π f π f − − ( 2 ( 2 1) 1) °¯ °¿
(A-3)
335
where f1 =5 Hz is the low-cut frequency, f2 =10 Hz is the low-pass frequency, f 3 =20
336
Hz is the high-pass frequency, and f4 =25 Hz is the high-cut frequency.
337 338
The waveforms of the above three wavelets and spectrograms are shown in Fig. A1a and Fig. A1b, respectively.
339 17
340 341
FIGURE CAPTIONS
342
Fig. 1. The CFL stability limit of the pseudospectral method with the 2nd-order temporal
343
FD scheme (PS (2nd-order)) and the 4th-order temporal FD scheme (PS (4th-order)). The
344
red line is the CFL stability limit for the PS (2nd-order) scheme, and the blue line is that
345
for the PS (4th-order) scheme. The corner frequency of the low-pass filter is indicated as
346
kmax .
347 348
Fig. 2. Snapshots at 333 ms obtained using a time step ǻt = 3 ms and the corresponding
349
spectral amplitudes in the wavenumber domain. (a1) and (a3) are the snapshots for the PS
350
(2nd-order) scheme and the PS (4th-order) scheme, respectively. The spatial filtering
351
algorithm stops between 300 to 333 ms. (a2) and (a4) are the spatially filtered snapshots
352
of (a1) and (a3), respectively; (b1)~(b4) are the corresponding wavenumber snapshots of
353
(a1)~(a4), where the blue arches in (b1) and (b3) are the corresponding corner frequency
354
of the low-pass filter for the PS (2nd-order) scheme and the PS (4th-order) scheme,
355
respectively, with a time step ǻt =3 ms and h = 10 m. (c) and (d) are the profiles of (a1)
356
and (b1) along the red dashed line, respectively.
357 358
Fig. 3. Waveforms recorded at x = í800 and z = í800 m using different schemes and
359
different time steps. The waveform obtained by the PS (4th-order) scheme using ǻt = 1 ms
360
is taken as the theoretical reference. The waveforms in (a) are obtained using spatial
361
filtering only, while the waveforms in (b) are obtained using the spatial filtering with the
362
ITDT method.
363 18
364
Fig. 4. Waveforms obtained using the spatial filtering with the ITDT method. The
365
waveforms in (a) are obtained using the Klauder wavelet, while the waveforms in (b) are
366
are obtained using the Ormsby wavelet.
367 368 369
Fig. 5. Comparison of waveforms using different schemes and time steps. (a) Modified
370
Marmousi model, where c(x,z) ranges from 2500 to 4000 m/s. (b) Waveforms recorded at
371
x = 5600 and z =5000 m. (c) Local zoom in of the rectangle in (b).
372 373
Fig. A1. Schematic diagram of Ricker wavelet, Klauder wavelet, and Ormsby wavelet. (a)
374
Waveforms in time domain. (b) Spectrum amplitude in frequency domain.
375
19
376 377
Fig. 1. The CFL stability limit of the pseudospectral method with the 2nd-order temporal
378
FD scheme (PS (2nd-order)) and the 4th-order temporal FD scheme (PS (4th-order)). The
379
red line is the CFL stability limit for the PS (2nd-order) scheme, and the blue line is that
380
for the PS (4th-order) scheme. The corner frequency of the low-pass filter is indicated as
381
kmax .
382 383
20
384
385 386
Fig. 2. Snapshots at 333 ms obtained using a time step ǻt = 3 ms and the corresponding
387
spectral amplitudes in the wavenumber domain. (a1) and (a3) are the snapshots for the PS
388
(2nd-order) scheme and the PS (4th-order) scheme, respectively. The spatial filtering
389
algorithm stops between 300 to 333 ms. (a2) and (a4) are the spatially filtered snapshots
390
of (a1) and (a3), respectively; (b1)~(b4) are the corresponding wavenumber snapshots of
391
(a1)~(a4), where the blue arches in (b1) and (b3) are the corresponding corner frequency 21
392
of the low-pass filter for the PS (2nd-order) scheme and the PS (4th-order) scheme,
393
respectively, with a time step ǻt =3 ms and h = 10 m. (c) and (d) are the profiles of (a1)
394
and (b1) along the red dashed line, respectively.
395 396
22
397
398 399
Fig. 3. Waveforms recorded at x = í800 and z = í800 m using different schemes and
400
different time steps. The waveform obtained by the PS (4th-order) scheme using ǻt = 1 ms
401
is taken as the theoretical reference. The waveforms in (a) are obtained using spatial
402
filtering only, while the waveforms in (b) are obtained using the spatial filtering with the
403
ITDT method.
404 405
23
406 407
Fig. 4. Waveforms obtained using the spatial filtering with the ITDT method. The
408
waveforms in (a) are obtained using the Klauder wavelet, while the waveforms in (b) are
409
are obtained using the Ormsby wavelet.
410 411
24
412
413 414
Fig. 5. Comparison of waveforms using different schemes and time steps. (a) Modified
415
Marmousi model, where c(x,z) ranges from 2500 to 4000 m/s. (b) Waveforms recorded at
416
x = 5600 and z =5000 m. (c) Local zoom in of the rectangle in (b).
417 418
25
419 420
Fig. A1. Schematic diagram of Ricker wavelet, Klauder wavelet, and Ormsby wavelet. (a)
421
Waveforms in time domain. (b) Spectrum amplitude in frequency domain.
422
26
423 424
TABLE CAPTIONS
425
Table 1. The CFL stability limit of the pseudospectral method (for the 2nd-order temporal
426
FD scheme and the 4th-order temporal FD scheme) and the corresponding corner
427
frequency of the low-pass filter using different time steps*.
428 429
*This table is generated using the velocity c=4000 m/s and the spatial grid interval h=10
430
m.
27
We declare no conflict interests.
Yingjie Gao, Jinhai Zhang, Zhenxing Yao
PII:
S0021-9991(19)30537-6
DOI:
https://doi.org/10.1016/j.jcp.2019.07.051
Reference:
YJCPH 8853
To appear in:
Journal of Computational Physics
Received date:
24 September 2018
Revised date:
27 May 2019
Accepted date:
22 July 2019
Please cite this article as: Y. Gao et al., Extending the stability limit of explicit scheme with spatial filtering for solving wave equations, J. Comput. Phys. (2019), doi: https://doi.org/10.1016/j.jcp.2019.07.051.
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.
Highlights • • • • •
We extend the stability limit of the explicit scheme using spatial filtering. Filtering out the high-wavenumber noise before it boosts can guarantee the stability. We further reduce the time-dispersion error by inverse time dispersion transform. The waveform is highly accurate compared with the reference waveforms. The time step is much larger than the maximum time step allowed by the CFL stability limit.
1
Extending the stability limit of explicit scheme with spatial
2
filtering for solving wave equations
3 4
Yingjie Gaoa, b , Jinhai Zhanga, b, *, and Zhenxing Yaoa, b
5
a Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics,
6
Chinese Academy of Sciences, Beijing, China
7
b Institutions of Earth Science, Chinese Academy of Sciences, Beijing, China
8
Running Head: Extending stability limit by spatial filtering
9 10 11
Submitted to Journal of Computational Physics on Sep 24, 2018
12
Revised manuscript submitted on Jan 24, 2019
13
Revised manuscript submitted on May 27, 2019
14
*Corresponding author: Dr. Jinhai Zhang
15
Chinese Academy of Sciences
16
Institute of Geology and Geophysics
17
Beijing 100029
18
China.
19
E-mail: [email protected]
20
Phone number: +86-10-82998013 1
21 22
Extending the stability limit of explicit scheme with spatial
23
filtering for solving wave equations
24
ABSTRACT
25
The explicit finite-difference scheme is widely used for solving the wave equation;
26
however, its time step is strictly constrained by the Courant-Friedrichs-Lewy (CFL)
27
stability limit. We apply spatial filtering to reduce the high-wavenumber noise caused by
28
using a larger time step beyond the CFL stability limit. We further incorporate the
29
forward and inverse time dispersion transforms to reduce the time-dispersion error. This
30
method allows us to use a time step larger than the maximum time step allowed by the
31
CFL stability limit.
32 33
.
34 35 36 37 38 39 40
Key words: explicit scheme; spatial filtering; CFL stability limit; time-dispersion error; large time step.
41 2
42 43
1. Introduction
44
The explicit scheme, i.e., iteratively updating the wavefields along the temporal
45
direction in the time domain, is the most popular methods for simulating seismic
46
wavefields because of its simplicity in numerical implementation [1]. However, its
47
maximum time step is strictly restricted by the Courant-Friedrichs-Lewy (CFL) stability
48
limit [2], which leads to a heavy computational burden in the presence of small-scale
49
structures and high-velocity targets. We have to adopt a small enough time step to
50
guarantee that the updated wavefields after each time step do not exceed the spatial range
51
of one cell; otherwise, we could encounter stability problems. The CFL stability limit
52
leads to a great computational burden, especially for imaging fine-structure targets (e.g.,
53
channels and fractures) using a small spatial interval. For example, when the spatial
54
interval is reduced by half, the computational cost for a 3D model would become 16
55
times that for the original spatial interval. We should develop new methods to overcome
56
the stability limit with an acceptable increase of computational cost when imaging
57
fine-structure targets.
58
In the electromagnetic wave simulation field, Li et al. [3] extended the CFL stability
59
limit by perturbing the unstable eigenvalues to be stable, which is called eigenvalue
60
perturbation [4]. Gao et al. [5] cooperated the forward and inverse time dispersion
61
transform methods [6-8] into the eigenvalue perturbation to solve the wave equation in
3
62
the seismic wave simulation field. As a preprocessing algorithm, the eigenvalue
63
perturbation does not increase the amount of computation during the iteration process.
64
In fluid dynamics, Ecer et al. [9] found that a larger time step beyond the CFL
65
stability limit would introduce unstable high-wavenumber wavefields; thus, one can filter
66
out the unstable high-wavenumber wavefields to enlarge the time step of an explicit
67
scheme. Sarris [10] proposed a spatial filtering algorithm to deal with the unstable
68
problems caused by a large time step when solving Maxwell’s equations. Using a
69
low-pass filter to the unstable components distributed in high-wavenumbers, the time step
70
of an explicit scheme can be successfully extended beyond the CFL stability limit
71
[11-13].
72
In this paper, we extend the CFL stability limit of the explicit scheme to solve the
73
wave equation. Proper low-pass spatial filtering [10] is used to reduce the slightly
74
boosted amplitude in each time step. The presented scheme allows us to apply a much
75
larger time step beyond the CFL stability limit. The spatial filtering can be perfectly
76
incorporated into the pseudospectral method and its computational cost can be neglected
77
compared with that of the Fourier transforms. To reduce the time-dispersion error caused
78
by using a large time step beyond the CFL stability limit, we incorporate the forward time
79
dispersion-transform (FTDT) method and the inverse time dispersion-transform (ITDT)
80
method [6-8]. The maximum time step in our test is around 2-3 times larger than the CFL
81
stability limit allowed. The memory demand is almost the same compared with that of the
4
82
classic pseudospectral method and the generated wavefield is shown to be highly
83
accurate.
84 85 86
2. Methodology We start with the 2D scalar wave equation
1 ∂ 2u ∂ 2u ∂ 2u = + , c 2 ∂t 2 ∂x 2 ∂z 2
87
(1)
88
where u ( x, z, t ) represents the particle displacement at position
( x, z )
89
c ( x, z ) is the velocity. We use the pseudospectral (PS) method for the spatial
90
discretization [14] since its error mainly arises from the temporal discretization, which is
91
helpful for identifying the errors caused by using larger time steps without considering
92
the space-dispersion error. The PS method for equation (1) reads
and time t, and
93
∂ 2u = c 2 F − ª¬ − k 2 F + ( u ) º¼ , ∂t 2
94
where k 2 = k x2 + k z2 , with kx and kz being the wavenumbers along the x- and z-
95
directions, respectively; F+ and F− are 2D forward and inverse Fourier transforms,
96
respectively.
97 98
{
}
(2)
Based on the Taylor expansion, a general time discretization scheme for the temporal derivative in equation (1) can be obtained by the following formula [15, 16] j
99
J Δt 2 j − 2 § ∂ 2 u n · u n +1 − 2u n + u n −1 = 2 ¦ ¨ 2 ¸ , Δt 2 j =1 ( 2 j ) ! © ∂t ¹
5
(3)
100
where 2J is the order of the temporal finite-difference (FD) scheme, and u n represents
101
the wavefield at time nΔt , and Δt is the time step. In practice, we only use a limited
102
number of terms on the right-hand side of equation (3) by the Lax-Wendroff methods,
103
which use spatial derivatives to replace high-order temporal derivatives [15, 17] 2 n u n +1 − 2u n + u n −1 ∂ 2u n 2§∂ u c = + ¨ 2 Δt 2 ∂z 2 © ∂x
104
j
J · Δt 2 j − 2 § ∂ 2 u n · 2 + ¸ ¦ ¨ 2 ¸ , j = 2 ( 2 j ) ! © ∂t ¹ ¹
(4)
105
and the higher-order temporal derivatives are obtained by the following recursive formula
106
[17] j
2 n § ∂ 2u n · ∂ 2u n · § ∂ 2u n · 2§∂ u c = + ¨ 2 ¸ ¨ 2 ¸¨ ¸ ∂ z 2 ¹ © ∂t 2 ¹ © ∂t ¹ © ∂x
107
j −1
.
(5)
108
The second-order (2nd-order) temporal FD scheme (i.e., J=1) and the fourth-order
109
(4th-order) temporal FD scheme (i.e., J=2) are commonly used in numerical simulations
110
[15-17]. The CFL stability limit for the 2D PS method with the 2nd-order temporal FD scheme
111 112
(i.e., PS (2nd-order)) and the 4th-order temporal FD scheme(i.e., PS (4th-order)) [14, 17]
( cmax k Δt )
113 114
115
2
≤1
4
(6)
and
( c k Δt ) 0 ≤ max 4
2
( c k Δt ) − max 48
4
≤ 1ˈ
(7)
116
respectively; and cmax is the maximum velocity. According to the CFL stability limit, the
117
time step Δ t must be limited by Δt ≤ 2 ( cmax kmax ) and Δt ≤ 2 3 ( cmax kmax ) for PS
118
(2nd-order) and PS (4th-order) schemes, respectively. Using a velocity c=4000 m/s and a 6
119
spatial grid interval ǻx=ǻz=h=10 m, we can obtain the maximum time step Δtmax ,
120
which corresponds to the maximum wavenumber kmax = 2π h (with k xmax =k zmax = π h )
121
for the 2D case, as shown in Table 1.
122
Fig. 1 shows the critical curves of the CFL stability limit for equation (6) and
123
equation (7). Obviously, the CFL stability limit is to make the time step not exceed the
124
maximum time step Δtmax , which corresponds to the k max , as shown by the shaded area
125
(or cyan area) in Fig. 1. For a given time step Δt that is larger than Δtmax , the
126
corresponding wavenumber kmax can be obtained from the critical curve. For the PS
127
(2nd-order)
128
kmax = 2 ( cmax Δ t ) and kmax = 2 3 ( cmax Δt ) , respectively.
and
the
PS
(4th-order)
schemes,
kmax
can
be
expressed
as
129
We can see that the wavefield components with a wavenumber smaller than kmax
130
are still sufficient for the CFL stability limit, as shown in Table 1. The kmax gradually
131
decreases with an increasing time step, which means that there is a possibility that we
132
could use a larger time step by limiting the range of wavenumbers using a proper
133
low-pass filter defined as [10, 13]
134
°1, F (k ) = ® °¯0,
for k x2 + k z2 ≤ kmax , otherwise,
(8)
135
where kmax is the corner frequency of the low-pass filter. The spatial filtering algorithm
136
is implemented by the following stages per time step [10]:
7
137
n+1 (1) After wavefield u has been updated at time
( n + 1) Δt
using equation (4), we
138
transform the wavefield to the wavenumber domain using a 2D Fourier transform
139
U n+1 =F + ( u n+1 ) .
140 141 142 143
n+1 =F ( k ) ⋅U n+1 in the wavenumber (2) Compute the spatially filtered wavefield U domain, where ⋅ is a Hadamard product. (3) Apply inverse 2D Fourier transform to obtain the filtered wavefield in the space
(
)
domain un+1 =F − U n+1 , and un +1 is taken as the input wavefield for the next iteration.
144
The spatial filtering algorithm can be incorporated perfectly into the implementation
145
of the PS method, since the PS method is implemented by shuttling the wavefield
146
frequently between the space and wavenumber domains. We can directly apply the spatial
147
filtering in the wavenumber domain together with the operation of wavenumber −k 2 .
148
Therefore, the computational cost of the PS method after applying the spatial filtering for
149
each time step could be almost perfectly retained, since the time consumed by the
150
Hadamard product associated with the spatial filtering can be neglected compared with
151
that of the Fourier transforms.
152 153
3. Numerical experiments
154
We perform numerical experiments on a homogeneous model. Two schemes are
155
tested: the PS (2nd-order) and PS (4th-order) schemes. The wave velocity is c=4000 m/s.
156
The spatial grid interval is ǻx=ǻz=10 m, and the grid number is 401×401. The source is a
157
Ricker wavelet with a dominant frequency of 15 Hz (shown in Appendix A) [18], which 8
158
is located at x=0 and z=0 m. The maximum time steps for the PS (2nd-order) and PS
159
(4th-order) schemes are 1.125 and 1.949 ms, respectively. We tested several large time
160
steps (ǻt=3, 5, 7, 9, and 11 ms) for each scheme. Note that all these time steps are beyond
161
the CFL stability limit. As a preparation, we filtered the source time functions at the
162
corresponding time steps by the FTDT method before the simulation [10]. The result
163
obtained by the PS (4th-order) scheme with a time step ǻt =1 ms was taken as a
164
theoretical reference to examine the accuracy of the waveforms obtained by different
165
time steps.
166
Fig. 2a1 and 2a3 show the snapshots at 333 ms using ǻt =3 ms for the PS (2nd-order)
167
scheme and the PS (4th-order) scheme, respectively. To illustrate the unstable phenomena
168
in the snapshot, the spatial filtering algorithm is stopped after 300 ms. Obviously, many
169
unstable components appear in the snapshot at 333 ms. Fig. 2b1 and 2b3 show the
170
corresponding wavenumber snapshots of Fig. 2a1 and 2a3, respectively. The unstable
171
components are mainly distributed in the high-wavenumber regions. As shown in Fig.
172
2a2, 2b2, 2a4, and 2b4, the spatial filtering can remove the unstable components at high
173
wavenumbers while retaining the effective stable wavefield at low-wavenumbers. Fig. 2c
174
and 2d show the profiles along the red dashed lines in Fig. 2a1 and 2b1, respectively. Fig.
175
2d clearly shows that the unstable components are distributed in a high-wavenumber
176
region, which can be clearly distinguished from the effective low-wavenumber wavefield
177
components when the time step is not much larger than the upper-limit time step Δtmax .
9
178
Fig. 3a shows the waveforms obtained by different schemes using different time
179
steps, recorded at x=í800 and z=í800 m. Fig. 3b shows the waveforms after applying the
180
ITDT method. Both the PS (2nd-order) scheme and the PS (4th-order) scheme after
181
applying the ITDT method are fairly accurate up to time step ǻt = 7 ms (even ǻt = 9 ms
182
for the PS (4th-order) scheme). This indicates that the combination of the spatial filtering
183
and the ITDT method could provide highly accurate simulation results even though a
184
much larger time step beyond the CFL stability limit is used. The PS (2nd-order) scheme
185
using ǻt 9 ms and the PS (4th-order) scheme using ǻt = 11 ms still have some residual
186
errors even after applying the ITDT method. At these cases, the corner frequency of the
187
low-pass filter kmax is smaller than the maximum wavenumber of the effective stable
188
wavefield. For example, kmax is 0.078729 for the PS (4th-order) scheme using ǻt = 11
189
ms (as shown in Table 1), while the maximum wavenumber of the effective stable
190
wavefield is k § 0.08 (as shown in Fig. 2d). This means that we could not separate the
191
effective wavefield and the unstable components using spatial filtering since they are
192
completely mixed into each other.
193
Based on the above constant velocity model, we perform two additional sets of
194
numerical experiments using two other wavelets (Klauder wavelet and Ormsby wavelet,
195
as shown in Appendix A) [18], respectively. The mid-frequency of both wavelets is 15 Hz.
196
We tested several large time steps (ǻt=3, 5, 7, 9, 11, and 13 ms). Fig. 4a and Fig. 4b show
197
the waveforms after applying the ITDT method to the wavefields obtained using the
198
Klauder wavelet and the Ormsby wavelet, respectively. The wavefoms are both stable and 10
199
accurate when using time steps beyond the CFL stability limit up to ǻt=11 ms, and the
200
result obtained using the Ormsby wavelet still has high accuracy even using time step
201
ǻt=13ms. This is due to the effective frequency band of the Ormsby wavelet is the
202
narrowest (the maximum value of the effective frequency is the smallest, as shown in Fig.
203
A1b) among these three wavelets, which means a smaller corner frequency kmax (i.e. a
204
larger time step ǻt) can be used without worrying about the effective wave field would be
205
filtered out.We verify the feasibility of the spatial filtering algorithm with the ITDT
206
methods by a modified Marmousi model, as shown in Fig. 5a. The spatial grid interval is
207
ǻx=ǻz=10 m, and the grid number is 737×751. The source time function is a Ricker
208
wavelet with a dominant frequency of 15 Hz, and the source is located at x=3200 and
209
z=3000 m. We record a trace at x=5600 and z=5000 m for different time steps, as shown
210
in Fig. 5b. Fig. 5c shows the local zoom in of the waveforms in the box of Fig. 5b.
211
Obviously, both the PS (2nd-order) scheme and the PS (4th-order) scheme, using the
212
spatial filtering with the ITDT methods, can obtain acceptable results using a much larger
213
time step than the maximum time steps allowed by the CFL stability limit. However, the
214
residual errors become evident at much larger time steps at ǻt = 7 ms.
215 216
4. Discussion
217
As shown in Fig. 2b1 and 2b3, we could easily separate unstable wavefield
218
components and the effective wavefield by setting a proper filter when the actual time
219
step is not much larger than the maximum time step allowed by the CFL stability limit. 11
220
For a model that has a strong velocity contrast between the maximum and minimum
221
velocities (i.e., heterogeneity is arbitrary), the wavefield corresponding to the low
222
velocity region would be distributed in higher wavenumbers, even larger than the corner
223
frequency. In other words, the unstable wavefield components are mixed with the
224
effective wavefield components when the velocity contrast is too strong; in such cases,
225
any attempt to filter the unstable components out would ruin the effective wavefields.
226
Therefore, we map the velocity of the original Marmousi model onto 2500 to 4000 m/s to
227
illustrate the feasibility of the proposed method on moderate variation of the velocity.
228
Future works are still needed to work on strong variation of the velocity.
229
With the spatial filtering algorithm, we can now use a much larger time step beyond
230
the CFL stability limit; however, the time step cannot be infinitely large. A larger time
231
step corresponds to a smaller corner frequency of the low-pass filter, which may filter out
232
the effective wavefield components. For example, with c=4000 m/s and h=10 m, the
233
corresponding corner frequency of ǻt = 9 ms for the PS (2nd-order) scheme is 1/18;
234
unfortunately, the low-pass filter at 1/16 would ruin the effective wavefield components,
235
as shown in Fig. 2d.
236
Another issue needs to emphasize is the Nyquist limit in temporal discretization. For
237
example, with c=4000 m/s and h=100 m, the upper-limit time step Δtmax for the PS
238
(4nd-order) scheme is 19.49 ms. The maximum frequency f max of the Ricker wavelet
239
with a dominant frequency of 15 Hz is about 45 Hz (as shown in Fig. A1b), and the time
12
240
step cannot be larger than 11.11 ms according to the Nyquist limit Δt < 1 ( 2 f max ) . In
241
other words, the time step reaches the Nyquist sampling limit first, while it has not
242
reached the upper limit of the stability condition yet. Fortunately, Gao et al. [5]
243
demonstrated that the time dispersion transform methods can guarantee the simulation
244
accuracy using a time step even towards the Nyquist limit. Therefore, in this case, one
245
only needs to apply the time dispersion transforms to ensure the simulation accuracy of
246
the simulation using a large time step simulation, and it is not necessary to consider using
247
spatial filtering method to extend the CFL stability limit.
248
To illustrate the unstable wavefield components, we do not apply the spatial filtering
249
between 300 and 333 ms, as shown in Figs. 2a1 and 2a3. In other words, the spatial
250
filtering is not used for continuous ten time steps, and the instability of the high
251
wavenumber can still be filtered out, as shown in Figs. 2a2 and 2a4. This indicates that it
252
is not necessary to apply the spatial filtering for each time step. We can save some
253
computational cost at several other time steps when the spatial discretization is purely in
254
the spatial domain.
255 256
5. Conclusion
257
We introduce the spatial filtering algorithm to guarantee the stability of the explicit
258
scheme for solving the wave equation with a much larger time step beyond the
259
Courant-Friedrichs-Lewy (CFL) stability limit. The computational cost is kept low since
260
the spatial filtering can be perfectly incorporated into the pseudospectral method. 13
261
Furthermore, we apply the forward time dispersion transform method and the inverse
262
time dispersion transform method to reduce the time-dispersion error caused by the large
263
time step. Both theoretical analyses and numerical experiments show that the proposed
264
method can use a larger time step than that allowed by the CFL stability limit, which is
265
helpful for saving a large number of iteration times without suffering from the time
266
dispersion and the stability problems. The proposed method retains the superiority of the
267
conventional explicit scheme by avoiding solving matrices but breaks through its major
268
limitation on the time step.
269 270
Acknowledgements
271
We are grateful to Erik Koene for helpful discussion, derivation of source loading for
272
the PS (4th-order) scheme, and the corresponding MATLAB code. This research was
273
supported by the National Natural Science Foundation of China (41704063) and by the
274
National Major Project of China (2017ZX05008-007). YG was also supported by the
275
General Financial Grant from the China Postdoctoral science foundation (2017M610980).
276
JZ thanks the supports from the Strategic Pioneer Program on Space Science, Chinese
277
Academy of Sciences (XDA15011700) and the Foundation for Excellent Member of the
278
Youth Innovation Promotion Association, Chinese Academy of Sciences (2016).
279 280
References 14
281 282 283 284 285 286 287 288 289 290 291 292
[1] J.T. Etgen, M.J. O'Brien, Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial, Geophysics, 72 (2007) SM223–SM230. [2] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen differenzengleichungen der mathematischen physik, Mathematische Annalen (In German), 100 (1928) 32–74. [3] X. Li, C.D. Sarris, P. Triverio, Overcoming the FDTD stability limit via model order reduction and eigenvalue perturbation, in:
Microwave Symposium, 2014, pp. 1–3.
[4] X. Li, Model order reduction and stability enforcement of finite-difference time-domain equations beyond the CFL limit, in, University of Toronto (Canada), 2014. [5] Y. Gao, J. Zhang, Z. Yao, Removing the stability limit of the explicit finite-difference scheme with eigenvalue perturbation, Geophysics, 83 (2018) A93–A98. [6] M. Wang, S. Xu, Finite-difference time dispersion transforms for wave propagation, Geophysics, 80 (2015) WD19–WD25.
293
[7] Y. Gao, J. Zhang, Z. Yao, Third-order symplectic integration method with inverse time
294
dispersion transform for long-term simulation, Journal of Computational Physics, 314 (2016)
295
436–449.
296 297
[8] E.F.M. Koene, J.O.A. Robertsson, F. Broggini, F. Andersson, liminating time dispersion from seismic wave modelling, Geophysical Journal International, 213 (2018) 169–180.
298
[9] A. Ecer, N. Gopalaswamy, H. Akay, Y. Chien, Digital filtering techniques for parallel
299
computation of explicit schemes, International Journal of Computational Fluid Dynamics, 13
300
(2000) 211–222.
15
301 302
[10] C.D. Sarris, Extending the stability limit of the FDTD method with spatial filtering, IEEE Microwave and Wireless Components Letters, 21 (2011) 176–178.
303
[11] C. Chang, C.D. Sarris, A spatial filter-enabled high-resolution subgridding scheme for stable
304
FDTD modeling of multiscale geometries, in: international microwave symposium, 2011,
305
pp. 1–4.
306 307
[12] C. Chang, C.D. Sarris, A three-dimensional spatially filtered FDTD with controllable stability beyond the courant limit, in:
Microwave Symposium Digest, 2012, pp. 1–3.
308
[13] C. Chang, C.D. Sarris, A spatially filtered finite-difference time-domain scheme with
309
controllable stability beyond the CFL limit: theory and applications, IEEE Transactions on
310
Microwave Theory & Techniques, 61 (2013) 351–359.
311 312 313 314 315 316 317 318 319 320
[14] D.D. Kosloff, E. Baysal, Forward modeling by a Fourier method, Geophysics, 47 (1982) 1402–1412. [15] M.A. Dablain, The application of high-order differencing to the scalar wave equation, Geophysics, 51 (1986) 54–66. [16] J.T. Etgen, Accurate wave equation modeling, Stanford Exploration Project, 60 (1988) 131–147. [17] J.B. Chen, High-order time discretizations in seismic modeling, Geophysics, 72 (2007) SM115–SM122. [18] H.Ryan, Ricker, Ormsby, Klauder, Butterworth—a choice of wavelets. CSEG Recorder, 19 (1994) 8-9.
321 16
322
APPENDIX A: THREE TYPES OF WAVELETS
323 324
Ricker wavelet can be expressed as [18] 2 −π Ricker ( t ) = ª1 − 2π 2 f 2 ( t − t0 ) º e ¬ ¼
325
2
f 2 ( t − t0 )
2
,
(A-1)
326
where f =15 Hz is the dominant frequency, and t0 = 0.15 s is the time delay of the
327
wavelet peak.
328 329
Klauder wavelet can be expressed as [18]
° sin ªπ kt (T − t ) º¼ −2iπ f0 (t −t0 ) °½ Klauder ( t ) = real ® ¬ e ¾ π kt ¯° ¿° ,
(A-2)
330
where T =7 is the duration of the input signal. k = ( f 2 − f 2 ) T is the rate of change
331
with time, with f1 =5 Hz being the terminal low frequency, and f2 =25 Hz being the
332
terminal high frequency. f 0 = ( f 2 +f 2 ) T is the mid-frequency of bandwidth.
333
334
Ormsby wavelet can be expressed as [18] ° (π f 4 )2 kt sinc 2 ª f 4 (T − t ) º (π f3 )2 kt sinc 2 ª f3 (T − t ) º ½° ¬ ¼− ¬ ¼ − Ormsby ( t ) = ® ¾ (π f 4 − π f 3 ) (π f 4 − π f3 ) ¯° ¿° ° (π f 2 )2 kt sinc 2 ª f 2 (T − t ) º (π f1 )2 kt sinc 2 ª f1 (T − t ) º ½° ¬ ¼− ¬ ¼ , ® ¾ π f π f π f π f − − ( 2 ( 2 1) 1) °¯ °¿
(A-3)
335
where f1 =5 Hz is the low-cut frequency, f2 =10 Hz is the low-pass frequency, f 3 =20
336
Hz is the high-pass frequency, and f4 =25 Hz is the high-cut frequency.
337 338
The waveforms of the above three wavelets and spectrograms are shown in Fig. A1a and Fig. A1b, respectively.
339 17
340 341
FIGURE CAPTIONS
342
Fig. 1. The CFL stability limit of the pseudospectral method with the 2nd-order temporal
343
FD scheme (PS (2nd-order)) and the 4th-order temporal FD scheme (PS (4th-order)). The
344
red line is the CFL stability limit for the PS (2nd-order) scheme, and the blue line is that
345
for the PS (4th-order) scheme. The corner frequency of the low-pass filter is indicated as
346
kmax .
347 348
Fig. 2. Snapshots at 333 ms obtained using a time step ǻt = 3 ms and the corresponding
349
spectral amplitudes in the wavenumber domain. (a1) and (a3) are the snapshots for the PS
350
(2nd-order) scheme and the PS (4th-order) scheme, respectively. The spatial filtering
351
algorithm stops between 300 to 333 ms. (a2) and (a4) are the spatially filtered snapshots
352
of (a1) and (a3), respectively; (b1)~(b4) are the corresponding wavenumber snapshots of
353
(a1)~(a4), where the blue arches in (b1) and (b3) are the corresponding corner frequency
354
of the low-pass filter for the PS (2nd-order) scheme and the PS (4th-order) scheme,
355
respectively, with a time step ǻt =3 ms and h = 10 m. (c) and (d) are the profiles of (a1)
356
and (b1) along the red dashed line, respectively.
357 358
Fig. 3. Waveforms recorded at x = í800 and z = í800 m using different schemes and
359
different time steps. The waveform obtained by the PS (4th-order) scheme using ǻt = 1 ms
360
is taken as the theoretical reference. The waveforms in (a) are obtained using spatial
361
filtering only, while the waveforms in (b) are obtained using the spatial filtering with the
362
ITDT method.
363 18
364
Fig. 4. Waveforms obtained using the spatial filtering with the ITDT method. The
365
waveforms in (a) are obtained using the Klauder wavelet, while the waveforms in (b) are
366
are obtained using the Ormsby wavelet.
367 368 369
Fig. 5. Comparison of waveforms using different schemes and time steps. (a) Modified
370
Marmousi model, where c(x,z) ranges from 2500 to 4000 m/s. (b) Waveforms recorded at
371
x = 5600 and z =5000 m. (c) Local zoom in of the rectangle in (b).
372 373
Fig. A1. Schematic diagram of Ricker wavelet, Klauder wavelet, and Ormsby wavelet. (a)
374
Waveforms in time domain. (b) Spectrum amplitude in frequency domain.
375
19
376 377
Fig. 1. The CFL stability limit of the pseudospectral method with the 2nd-order temporal
378
FD scheme (PS (2nd-order)) and the 4th-order temporal FD scheme (PS (4th-order)). The
379
red line is the CFL stability limit for the PS (2nd-order) scheme, and the blue line is that
380
for the PS (4th-order) scheme. The corner frequency of the low-pass filter is indicated as
381
kmax .
382 383
20
384
385 386
Fig. 2. Snapshots at 333 ms obtained using a time step ǻt = 3 ms and the corresponding
387
spectral amplitudes in the wavenumber domain. (a1) and (a3) are the snapshots for the PS
388
(2nd-order) scheme and the PS (4th-order) scheme, respectively. The spatial filtering
389
algorithm stops between 300 to 333 ms. (a2) and (a4) are the spatially filtered snapshots
390
of (a1) and (a3), respectively; (b1)~(b4) are the corresponding wavenumber snapshots of
391
(a1)~(a4), where the blue arches in (b1) and (b3) are the corresponding corner frequency 21
392
of the low-pass filter for the PS (2nd-order) scheme and the PS (4th-order) scheme,
393
respectively, with a time step ǻt =3 ms and h = 10 m. (c) and (d) are the profiles of (a1)
394
and (b1) along the red dashed line, respectively.
395 396
22
397
398 399
Fig. 3. Waveforms recorded at x = í800 and z = í800 m using different schemes and
400
different time steps. The waveform obtained by the PS (4th-order) scheme using ǻt = 1 ms
401
is taken as the theoretical reference. The waveforms in (a) are obtained using spatial
402
filtering only, while the waveforms in (b) are obtained using the spatial filtering with the
403
ITDT method.
404 405
23
406 407
Fig. 4. Waveforms obtained using the spatial filtering with the ITDT method. The
408
waveforms in (a) are obtained using the Klauder wavelet, while the waveforms in (b) are
409
are obtained using the Ormsby wavelet.
410 411
24
412
413 414
Fig. 5. Comparison of waveforms using different schemes and time steps. (a) Modified
415
Marmousi model, where c(x,z) ranges from 2500 to 4000 m/s. (b) Waveforms recorded at
416
x = 5600 and z =5000 m. (c) Local zoom in of the rectangle in (b).
417 418
25
419 420
Fig. A1. Schematic diagram of Ricker wavelet, Klauder wavelet, and Ormsby wavelet. (a)
421
Waveforms in time domain. (b) Spectrum amplitude in frequency domain.
422
26
423 424
TABLE CAPTIONS
425
Table 1. The CFL stability limit of the pseudospectral method (for the 2nd-order temporal
426
FD scheme and the 4th-order temporal FD scheme) and the corresponding corner
427
frequency of the low-pass filter using different time steps*.
428 429
*This table is generated using the velocity c=4000 m/s and the spatial grid interval h=10
430
m.
27
We declare no conflict interests.