# Std import block
import time
import numpy as np
import matplotlib.pyplot as plt
from pysit import *
from pysit.gallery import horizontal_reflector
if __name__ == '__main__':
# Setup
# Define Domain
pmlx = PML(0.1, 100)
pmlz = PML(0.1, 100)
x_config = (0.1, 1.0, pmlx, pmlx)
z_config = (0.1, 0.8, pmlz, pmlz)
d = RectangularDomain(x_config, z_config)
m = CartesianMesh(d, 91, 71)
# Generate true wave speed
C, C0, m, d = horizontal_reflector(m)
# Set up shots
zmin = d.z.lbound
zmax = d.z.rbound
zpos = zmin + (1./9.)*zmax
shots = equispaced_acquisition(m,
RickerWavelet(10.0),
sources=1,
source_depth=zpos,
source_kwargs={},
receivers='max',
receiver_depth=zpos,
receiver_kwargs={},
)
# Define and configure the wave solver
trange = (0.0,3.0)
solver = ConstantDensityAcousticWave(m,
spatial_accuracy_order=2,
trange=trange,
kernel_implementation='cpp')
# Generate synthetic Seismic data
tt = time.time()
wavefields = []
base_model = solver.ModelParameters(m,{'C': C})
generate_seismic_data(shots, solver, base_model, wavefields=wavefields)
print('Data generation: {0}s'.format(time.time()-tt))
objective = TemporalLeastSquares(solver)
# Define the inversion algorithm
invalg = LBFGS(objective)
initial_value = solver.ModelParameters(m,{'C': C0})
# Execute inversion algorithm
print('Running LBFGS...')
tt = time.time()
nsteps = 5
status_configuration = {'value_frequency' : 1,
'residual_frequency' : 1,
'residual_length_frequency' : 1,
'objective_frequency' : 1,
'step_frequency' : 1,
'step_length_frequency' : 1,
'gradient_frequency' : 1,
'gradient_length_frequency' : 1,
'run_time_frequency' : 1,
'alpha_frequency' : 1,
}
# line_search = ('constant', 1e-16)
line_search = 'backtrack'
result = invalg(shots, initial_value, nsteps,
line_search=line_search,
status_configuration=status_configuration, verbose=True)
print('...run time: {0}s'.format(time.time()-tt))
obj_vals = np.array([v for k,v in list(invalg.objective_history.items())])
plt.figure()
plt.semilogy(obj_vals)
clim = C.min(),C.max()
# Do something to visualize the results
plt.figure()
plt.subplot(3,1,1)
vis.plot(C0, m, clim=clim)
plt.title('Initial Model')
plt.subplot(3,1,2)
vis.plot(C, m, clim=clim)
plt.title('True Model')
plt.subplot(3,1,3)
vis.plot(result.C, m, clim=clim)
plt.title('Reconstruction')
plt.show()