Seismic PSV with discontinuity: questions and bug report

[Pbottelin]

Hi,
and thanks for the amazing work done with Fatiendo and its open distribution to the scientific community!

I am quite new is this package, as I deal mostly with seismic refraction case studies using Pygimli.
I decided to give a try with Fatiendo because I need to simulate a complet waveform, including reflected waves (not handled by pygimli if I am not mistaken).

I downloaded and tried the 2D P-SV with moho example from the cookbook (http://www.fatiando.org/dev/cookbook/seismic_wavefd_rayleigh_wave.html). It runs nicely without error.

I then tried to adapt it to my needs, i.e. active seismic surveying:

• shorter area dimensions (e.g. 500m long x 100m in depth)
• many sensors (seismic spread, e.g. 48 sensors with 2m spacing)
• an horizontal discontinuity between layer 1 (VP=300m/s, VS=200m/s, Rho=2.0g/cm3) and layer 2 (VP=3000m/s, VS=2000m/s, Rho=2.7g/cm3)
• central frequency of blast source around 10 Hz. Source located close to the surface
• simulation duration of a few seconds (e.g. 4s)

As extra steps, I store the vertical signals as obspy stream to plot the seismic section (see figure that worked hereafter)/

Unfortunately, I cannot succeed in running those simulations.
The code is shown down there.

I get several errors, mostly:
"Traceback (most recent call last):
File "C:\Users\p.bottelin\Anaconda3\envs\fatiando_prod\lib\site-packages\matplotlib\cbook_init_.py", line 387, in process
proxy(*args, **kwargs)
File "C:\Users\p.bottelin\Anaconda3\envs\fatiando_prod\lib\site-packages\matplotlib\cbook_init_.py", line 227, in call
return mtd(*args, **kwargs)
File "C:\Users\p.bottelin\Anaconda3\envs\fatiando_prod\lib\site-packages\matplotlib\animation.py", line 1026, in _start
self._init_draw()
File "C:\Users\p.bottelin\Anaconda3\envs\fatiando_prod\lib\site-packages\matplotlib\animation.py", line 1750, in _init_draw
self._draw_frame(next(self.new_frame_seq()))
File "C:\Users\p.bottelin\Anaconda3\envs\fatiando_prod\lib\site-packages\matplotlib\animation.py", line 1772, in _draw_frame
self._drawn_artists = self._func(framedata, *self._args)
File "C:/Users/p.bottelin/Desktop/fatiando/test_2layers.py", line 25, in animate
xseismogram.set_data(times[:t + 1], xcomp[0][:t + 1])
NameError: global name 'xseismogram' is not defined"

This occurs mainly when I play with the "shape" tuple and/or the layer velocity, dt, ...
Could you give me some hints/guidelines to resolve these issues?

Thanks very much,
best regards!

Pierre

********************************** CODE ************************************************
"""
Seismic: 2D finite difference simulation of elastic P and SV wave propagation
FATIANDO package
try 2020/07/22
"""

""" ******** IMPORTS ********* """
import numpy as np
from matplotlib import animation
from fatiando import gridder
from fatiando.seismic import wavefd
from fatiando.vis import mpl

import matplotlib.pyplot as plt
import obspy

""" ******** FUNCTION DEFINITION ********* """
def animate(i):
t, p, s, xcomp, zcomp = simulation.next()
mpl.title('time: %0.1f s' % (times[t]))
wavefield.set_array((background + p + s)[::-1])
xseismogram.set_data(times[:t + 1], xcomp[0][:t + 1])
zseismogram.set_data(times[:t + 1], zcomp[0][:t + 1])
return wavefield, xseismogram, zseismogram

""" ******** MAIN ********* """

""" WORLD """

Set the parameters of the finite difference grid

area = [0, 500, 0, 300] #world [xmin, xmax, zmin, zmax] in m
shape = (100, 100) #number of finite difference nodes nx, nz
#note: shape (i.e. number of nodes) should be controled by
#Lambda minimum, which is related to Vsmini and frequency max.
#its is currently fixed by the user.

""" LAYERS """

LAYER 01

density = 2000 * np.ones(shape) #kg.m-3
pvel = 300 * np.ones(shape) #m/s
svel = pvel/2 #m/s

LAYER 02

moho = shape[0]/2 #discontinuity position
density[moho:] = 2700 #kg.m-3
pvel[moho:] = 3000 #m/s
svel[moho:] = pvel[moho:]/2 #m/s

compute Lamé parameters for layers

mu = wavefd.lame_mu(svel, density)
lamb = wavefd.lame_lamb(pvel, svel, density)

""" SOURCE """

Make a wave source from a mexican hat wavelet that vibrates in the x and z

directions equaly

central_freq = 10 #Hz
ampl=1

wavefd.Source(sx, sz, area bounding the FD, shape of nodes,

#source amplitude, source central frequency, delay before source starts)
sources = wavefd.blast_source(20, 0, area, shape, ampl,
central_freq , delay=1.0/central_freq)

""" SENSORS """
nb_of_sensors = 10
sx_start = 25 #0m
dx = 10 #m
sensor_list = [[x, 0] for x in np.arange(sx_start, sx_start+nb_of_sensors*dx, dx)] # x, z coordinate of the seismometer in m

""" SIMULATION """

Get the iterator for the simulation

dt = wavefd.maxdt(area, shape, pvel.max())
duration = 4.0 #in s
maxit = int(duration / dt)
times = np.linspace(0, maxit * dt, maxit)
snapshots = int(0.5 / dt) # Plot a snapshot of the simulation every 0.5 seconds

runs simulation

simulation = wavefd.elastic_psv(lamb,
mu,
density,
area,
dt,
maxit,
sources,
sensor_list,
snapshots,
padding=50,
taper=0.01,
xz2ps=True)

""" FIGURE ANIMATION """
fig = mpl.figure(figsize=(12, 5))
background = 10 ** -5 * ((density - density.min()) / density.max())

Start with everything zero and grab the plot so that it can be updated later

wavefield = mpl.imshow(np.zeros(shape), extent=area, vmin=-10**-6,
vmax=10**-6, cmap=mpl.cm.gray_r)
mpl.points(sensor_list, '^k')
mpl.ylim(area[2:][::-1])
mpl.xlabel('x (m)')
mpl.ylabel('z (m)')
times = np.linspace(0, maxit * dt, maxit)

This function updates the plot every few timesteps

anim = animation.FuncAnimation(
fig, animate, frames = maxit/snapshots, interval=1)
mpl.show()

""" BUILD OBSPY SIGNALS """
simulation = wavefd.elastic_psv(lamb,
mu,
density,
area,
dt,
maxit,
sources,
sensor_list,
snapshot=None,
padding=50,
taper=0.01,
xz2ps=True)
t, p, s, xcomp, zcomp = simulation.next() #go to last step of iteration
st = obspy.Stream()
for idx_sensor, sensor in enumerate(sensor_list):
tr = obspy.Trace()
tr.stats.station = 'GPZ' + str(idx_sensor).zfill(02)
tr.stats.delta = np.unique(np.round(np.diff(times),4))
tr.stats.distance = sensor[0]
tr.data = zcomp[idx_sensor]
st.append(tr)
print(st) #show seismic stream

""" SHOW FIGURE SEISMIC SECTION """
display_length = 2.0 #s
st.plot(type='section',
time_down=True,
scale=10.0,
show=True, block=True)

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible.

You might also want to take a look at our Contributing Guide and Code of Conduct.