In [1]:
import numpy as np
from scipy.special import hankel2
from examples.seismic.acoustic import AcousticWaveSolver
from examples.seismic import Model, RickerSource, Receiver
from devito import clear_cache, set_log_level
from opescibench import LinePlotter
import matplotlib.pyplot as plt
from matplotlib import cm
In [2]:
# Model with fixed time step value
class ModelBench(Model):
"""
Physical model used for accuracy benchmarking.
The critical dt is made small enough to ignore
time discretization errors
"""
@property
def critical_dt(self):
"""Critical computational time step value."""
return .1
Verification¶
We compute the error between the numerical and reference solutions for varying spatial discretization order and grid spacing. We also compare the time to solution to the error for these parameters.
In [3]:
# Discretization order
orders = (2, 4, 6, 8, 10)
norder = len(orders)
# Domain sizes and gird spacing
shapes = ((201, 2.0), (161, 2.5), (101, 4.0))
dx = [2.0, 2.5, 4.0]
nshapes = len(shapes)
# Number of time steps
nt = 1501
# Time axis
time = np.linspace(0., 150., nt)
# Source peak frequency
f0 = .07
In [4]:
# Fine grid model
c0 = 1.5
model = ModelBench(vp=c0, origin=(0., 0.), spacing=(.5, .5), shape=(801, 801), nbpml=40, dtype=np.float64)
In [5]:
# Source and receiver
src = RickerSource(name='src', grid=model.grid, f0=f0, time=time)
src.coordinates.data[0, :] = np.array(model.domain_size) * .5
src.data[:] = 1e2*src.data[:] / (.5*c0)**2
# Define receiver geometry (spread across x, just below surface)
rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=1)
rec.coordinates.data[:, 0] = 260.
rec.coordinates.data[:, 1] = 260.
Reference solution for numerical convergence¶
In [6]:
solver = AcousticWaveSolver(model, source=src, receiver=rec, time_order=2, space_order=20)
ref_rec, ref_u, _ = solver.forward()
GNUCompiler: compiled /var/folders/mx/qs0dn9rx7zn6dz2zvwv7tkk00000gn/T/devito-cjqed58a/7261032ccc57ff87dc245bd04b59bcf3be82b92b.c [1.09 s] ========================================================================================= Section section_1<1499,1> with OI=0.62 computed in 0.001 s [0.06 GFlops/s] Section section_2<1499,1> with OI=0.77 computed in 0.082 s [0.00 GFlops/s] Section main<1499,861,861> with OI=2.86 computed in 6.601 s [12.12 GFlops/s, 0.17 GPts/s] =========================================================================================
Analytical solution for comparison with the reference numerical solution¶
The analytical solution of the 2D acoustic wave-equation with a source pulse is defined as:
$$ u_s(r, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ -i \pi H_0^{(2)}\left(k r \right) q(\omega) e^{i\omega t} d\omega\} $$$$ r = \sqrt{(x - x_{src})^2+(y - y_{src})^2} $$where $H_0^{(2)}$ is the Hankel function of the second kind, $F(\omega)$ is the Fourier spectrum of the source time function at angular frequencies $\omega$ and $k = \frac{\omega}{v}$ is the wavenumber.
We look at the analytical and numerical solution at a single grid point. We ensure that this grid point is on-the-grid for all discretizations analyised in the further verification.
In [7]:
# Source and receiver coordinates
sx, sz = src.coordinates.data[0, :]
rx, rz = rec.coordinates.data[0, :]
# Define an Ricker wavelet shifted to zero lag for the Fourier transform
def ricker(f, T, dt, t0):
t = np.linspace(-t0, T-t0, int(T/dt))
tt = (np.pi**2) * (f**2) * (t**2)
y = 1e2*(1.0 - 2.0 * tt) * np.exp(- tt)
return y
def analytical(nt, model, time, **kwargs):
dt = kwargs.get('dt', model.critical_dt)
# Fourier constants
nf = int(nt/2 + 1)
fnyq = 1. / (2 * dt)
df = 1.0 / time[-1]
faxis = df * np.arange(nf)
rick= ricker(f0, time[-1], dt, 1.0/f0)
# Take the FOurier transform ofthe source time-function
R = np.fft.fft(rick/(c0**2))
R = R[0:nf]
nf = len(R)
# Compute the Hankel function and multiply by the source spectrum
U_a = np.zeros((nf), dtype=complex)
for a in range(1, nf-1):
k = 2 * np.pi * faxis[a] / c0
tmp = k * np.sqrt(((rx - sx))**2 + ((rz - sz))**2)
U_a[a] = -1j * np.pi * hankel2(0.0, tmp) * R[a]
# Do inverse fft on 0:dt:T and you have analytical solution
U_t = 1.0/(2.0 * np.pi) * np.real(np.fft.ifft(U_a[:], nt))
return np.real(U_t)
In [8]:
time1 = np.linspace(0.0, 3000., 30001)
U_t = analytical(30001, model, time1, dt=time1[1] - time1[0])
U_t = U_t[0:1501]
In [9]:
# Plot wavefield and source/rec position
plt.figure()
plt.imshow(ref_u.data[1,:,:], vmin=-1., vmax=1., cmap="seismic")
plt.plot(2*sx+40, 2*sz+40, 'r*', markersize=11, label='source') # plot position of the source in model, add nbpml for correct position
plt.plot(2*rx+40, 2*rz+40, 'k^', markersize=8, label='receiver') # plot position of the receiver in model, add nbpml for correct position
plt.legend()
plt.xlabel('x position (m)')
plt.ylabel('z position (m)')
plt.savefig('wavefieldperf.pdf')
# Plot trace
plt.figure()
plt.plot(time, ref_rec.data[:, 0], '-b', label='numerical')
plt.plot(time, U_t[:], '--r', label='analytical')
plt.xlabel('time (ms)')
plt.ylabel('amplitude')
plt.legend()
plt.savefig('ref.pdf')
plt.show()
In [10]:
error_time = np.zeros(5)
error_time[0] = np.linalg.norm(U_t[:-1] - ref_rec.data[:-1, 0], 2)/np.sqrt(nt)
print(error_time[0])
0.00128222243058
Convergence in time¶
We first show the convergence of the time discretization for a fix high-order spatial discretization (20th order).
After we show that the time discretization converges in $O(dt^2)$ and therefore only contains the error in time, we will take the numerical solution for dt=.1ms as a reference for the spatial discretization analysis.
In [11]:
dt = [.1, 0.08, .075, 0.0625, .05]
nnt = [1501, 1876, 2001, 2401, 3001]
for i in range(1, 5):
# Time axis
time = np.linspace(0., 150., nnt[i])
# Source and receiver
src = RickerSource(name='src', grid=model.grid, f0=f0, time=time)
src.coordinates.data[0, :] = np.array(model.domain_size) * .5
src.data[:] = 100.0 * src.data[:] / (.5*c0)**2
# Define receiver geometry (spread across x, just below surface)
rec = Receiver(name='rec', grid=model.grid, ntime=nnt[i], npoint=1)
rec.coordinates.data[:, 0] = 260.
rec.coordinates.data[:, 1] = 260.
solver = AcousticWaveSolver(model, source=src, receiver=rec, time_order=2, space_order=20)
ref_rec1, ref_u1, _ = solver.forward(dt=dt[i])
time1 = np.linspace(0.0, 3000., 20*(nnt[i]-1) + 1)
U_t1 = analytical(20*(nnt[i]-1) + 1, model, time1, dt=time1[1] - time1[0])
U_t1 = U_t1[0:nnt[i]]
error_time[i] = np.linalg.norm(U_t1[:-1] - ref_rec1.data[:-1, 0], 2)/np.sqrt(nnt[i])
print("error for dt= %s is %s \n" % (dt[i], error_time[i]))
GNUCompiler: compiled /var/folders/mx/qs0dn9rx7zn6dz2zvwv7tkk00000gn/T/devito-cjqed58a/7261032ccc57ff87dc245bd04b59bcf3be82b92b.c [0.33 s] ========================================================================================= Section section_1<1874,1> with OI=0.62 computed in 0.002 s [0.06 GFlops/s] Section section_2<1874,1> with OI=0.77 computed in 0.084 s [0.00 GFlops/s] Section main<1874,861,861> with OI=2.86 computed in 8.354 s [11.97 GFlops/s, 0.17 GPts/s] =========================================================================================
error for dt= 0.08 is 0.000850230956588
GNUCompiler: compiled /var/folders/mx/qs0dn9rx7zn6dz2zvwv7tkk00000gn/T/devito-cjqed58a/7261032ccc57ff87dc245bd04b59bcf3be82b92b.c [0.72 s] ========================================================================================= Section section_1<1999,1> with OI=0.62 computed in 0.002 s [0.05 GFlops/s] Section section_2<1999,1> with OI=0.77 computed in 0.090 s [0.00 GFlops/s] Section main<1999,861,861> with OI=2.86 computed in 9.517 s [11.21 GFlops/s, 0.16 GPts/s] =========================================================================================
error for dt= 0.075 is 0.000755921837441
GNUCompiler: compiled /var/folders/mx/qs0dn9rx7zn6dz2zvwv7tkk00000gn/T/devito-cjqed58a/7261032ccc57ff87dc245bd04b59bcf3be82b92b.c [0.37 s] ========================================================================================= Section section_1<2399,1> with OI=0.62 computed in 0.003 s [0.05 GFlops/s] Section section_2<2399,1> with OI=0.77 computed in 0.103 s [0.00 GFlops/s] Section main<2399,861,861> with OI=2.86 computed in 13.471 s [9.51 GFlops/s, 0.13 GPts/s] =========================================================================================
error for dt= 0.0625 is 0.000542844945677
GNUCompiler: compiled /var/folders/mx/qs0dn9rx7zn6dz2zvwv7tkk00000gn/T/devito-cjqed58a/7261032ccc57ff87dc245bd04b59bcf3be82b92b.c [0.71 s] ========================================================================================= Section section_1<2999,1> with OI=0.62 computed in 0.003 s [0.05 GFlops/s] Section section_2<2999,1> with OI=0.77 computed in 0.140 s [0.00 GFlops/s] Section main<2999,861,861> with OI=2.86 computed in 13.685 s [11.70 GFlops/s, 0.16 GPts/s] =========================================================================================
error for dt= 0.05 is 0.000363658585916
In [12]:
dt = [.1, 0.08, .075, 0.0625, .05]
plt.figure(figsize=(20, 10))
theory = [t**2 for t in dt]
theory = [error_time[0]*th/theory[0] for th in theory]
plt.loglog([t for t in dt], error_time, '-ob', label=('Numerical'), linewidth=4, markersize=10)
plt.loglog([t for t in dt], theory, '-^r', label=('Theory (2nd order)'), linewidth=4, markersize=10)
for x, y, a in zip([t for t in dt], theory, [('dt = %s ms' % (t)) for t in dt]):
plt.annotate(a, xy=(x, y), xytext=(4, 2),
textcoords='offset points', size=20,
horizontalalignment='left', verticalalignment='top')
plt.xlabel("Time-step $dt$ (ms)", fontsize=20)
plt.ylabel("$|| u_{num} - u_{ana}||_2$", fontsize=20)
plt.tick_params(axis='both', which='both', labelsize=20)
plt.tight_layout()
plt.xlim((0.05, 0.1))
plt.legend(fontsize=20, ncol=4, fancybox=True, loc='best')
plt.savefig("TimeConvergence.pdf", format='pdf', facecolor='white',
orientation='landscape', bbox_inches='tight')
plt.show()
In [13]:
np.polyfit(np.log([t for t in dt]), np.log(error_time), deg=1)
Out[13]:
array([ 1.81745159, -2.47763232])
Convergence in space¶
We have a correct refernce solution we can use for space discretization analysis
In [17]:
errorl2 = np.zeros((norder, nshapes))
timing = np.zeros((norder, nshapes))
set_log_level("ERROR")
ind_o = -1
for spc in orders:
ind_o +=1
ind_spc = -1
for nn, h in shapes:
ind_spc += 1
clear_cache()
time = np.linspace(0., 150., nt)
model = ModelBench(vp=c0, origin=(0., 0.), spacing=(h, h), shape=(nn, nn),
nbpml=40, dtype=np.float64)
src = RickerSource(name='src', grid=model.grid, f0=f0, time=time)
src.coordinates.data[0, :] = 200.
src.data[:] = 100. * src.data[:] / (c0*h)**2
# Define receiver geometry
rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=1)
rec.coordinates.data[:, 0] = 260.
rec.coordinates.data[:, 1] = 260.
solver = AcousticWaveSolver(model, source=src, receiver=rec, time_order=2, space_order=spc)
loc_rec, loc_u, summary = solver.forward()
# Compare to reference solution
errorl2[ind_o, ind_spc] = np.linalg.norm(loc_rec.data[:] - ref_rec.data[:], np.inf)
timing[ind_o, ind_spc] = summary.timings['main']
print("starting space order %s with (%s, %s) grid points the error is %s for %s seconds runtime" %
(spc, nn, nn, errorl2[ind_o, ind_spc], timing[ind_o, ind_spc]))
starting space order 2 with (201, 201) grid points the error is 0.598397364492 for 0.549626 seconds runtime starting space order 2 with (161, 161) grid points the error is 0.928768514866 for 0.432005 seconds runtime starting space order 2 with (101, 101) grid points the error is 1.6639252837 for 0.325349 seconds runtime starting space order 4 with (201, 201) grid points the error is 0.035965613536 for 0.649475 seconds runtime starting space order 4 with (161, 161) grid points the error is 0.0906846693164 for 0.404651 seconds runtime starting space order 4 with (101, 101) grid points the error is 0.533946654328 for 0.234384 seconds runtime starting space order 6 with (201, 201) grid points the error is 0.00363877161591 for 0.661981 seconds runtime starting space order 6 with (161, 161) grid points the error is 0.0136802978641 for 0.4465 seconds runtime starting space order 6 with (101, 101) grid points the error is 0.194683041346 for 0.226194 seconds runtime starting space order 8 with (201, 201) grid points the error is 0.000704455837682 for 0.680788 seconds runtime starting space order 8 with (161, 161) grid points the error is 0.00303256891849 for 0.445353 seconds runtime starting space order 8 with (101, 101) grid points the error is 0.0929970041058 for 0.251205 seconds runtime starting space order 10 with (201, 201) grid points the error is 0.000334591667136 for 0.487278 seconds runtime starting space order 10 with (161, 161) grid points the error is 0.000999726295846 for 0.426763 seconds runtime starting space order 10 with (101, 101) grid points the error is 0.0516523892572 for 0.23192 seconds runtime
In [20]:
stylel = ('-^k', '-^b', '-^r', '-^g', '-^c')
plt.figure(figsize=(20, 10))
for i in range(0, 5):
plt.loglog(errorl2[i, :], timing[i, :], stylel[i], label=('order %s' % orders[i]), linewidth=4, markersize=10)
for x, y, a in zip(errorl2[i, :], timing[i, :], [('dx = %s m' % (sc)) for sc in dx]):
plt.annotate(a, xy=(x, y), xytext=(4, 2),
textcoords='offset points', size=20)
plt.xlabel("$|| u_{num} - u_{ref}||_{inf}$", fontsize=20)
plt.ylabel("Runtime (sec)", fontsize=20)
plt.tick_params(axis='both', which='both', labelsize=20)
plt.tight_layout()
plt.legend(fontsize=20, ncol=3, fancybox=True, loc='lower left')
plt.savefig("TimeAccuracy.pdf", format='pdf', facecolor='white',
orientation='landscape', bbox_inches='tight')
plt.show()
In [16]:
stylel = ('-^k', '-^b', '-^r', '-^g', '-^c')
style2 = ('--k', '--b', '--r', '--g', '--c')
plt.figure(figsize=(20, 10))
for i in range(0, 5):
theory = [k**(orders[i]) for k in dx]
theory = [errorl2[i, 2]*th/theory[2] for th in theory]
plt.loglog([sc for sc in dx], errorl2[i, :], stylel[i], label=('Numerical order %s' % orders[i]),
linewidth=4, markersize=10)
plt.loglog([sc for sc in dx], theory, style2[i], label=('Theory order %s' % orders[i]),
linewidth=4, markersize=10)
# for x, y, a in zip([sc for sc in dx], errorl2[i, :], [('dx = %s m' % (sc)) for sc in dx]):
# plt.annotate(a, xy=(x, y), xytext=(4, 2),
# textcoords='offset points', size=20,
# horizontalalignment='left', verticalalignment='top')
plt.xlabel("Grid spacing $dx$ (m)", fontsize=20)
plt.ylabel("$||u_{num} - u_{ref}||_{inf}$", fontsize=20)
plt.tick_params(axis='both', which='both', labelsize=20)
plt.tight_layout()
plt.legend(fontsize=20, ncol=4, fancybox=True, loc='lower right')
plt.xlim((2.0, 4.0))
plt.savefig("Convergence.pdf", format='pdf', facecolor='white',
orientation='landscape', bbox_inches='tight')
plt.show()
