Question about the source

[yulongovo]

I would like to ask: when performing elastic wave forward modeling, we usually set one source for each shot. Does Deepwave support simulating data for a single shot with multiple sources activated simultaneously?

[ar4]

Dear @yulongovo

It is nice to hear from you again. I once again apologise for the delay in releasing the 3D propagator. I have been doing experiments with different implementations to try to improve the runtime performance, but I have now decided to focus on releasing it soon and will leave further performance optimisation until the subsequent release.

Regarding your question, you can indeed have multiple sources active simultaneously. The second dimension of the source amplitudes and locations Tensors is the "number of sources per shot" dimension, which you can use to specify multiple sources per shot. In the elastic propagator, you can also have multiple different source types (directional body forces in different dimensions, and pressure sources) active simultaneously.

If you wish to have different numbers of active sources per shot, you can set the locations for inactive sources to deepwave.IGNORE_LOCATION to deactivate them.

Here is an example that shows propagating two shots simultaneously. The first shot has three directional body forces in the y dimension and one in the x dimension, while the second shot has two in the y dimension and one in the x dimension.

import torch
import deepwave
import matplotlib.pyplot as plt

ny = 100
nx = 100
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
vp = torch.ones(ny, nx, device=device) * 1500
vs = torch.ones(ny, nx, device=device) * 1000
rho = torch.ones(ny, nx, device=device) * 2200
lamb, mu, buoyancy = deepwave.common.vpvsrho_to_lambmubuoyancy(vp, vs, rho)
freq = 25
nt = 25
dt = 0.004
peak_time = 1.5 / freq

n_shots = 2
n_sources_y_per_shot = 3

source_amplitudes_y = (
    deepwave.wavelets.ricker(freq, nt, dt, peak_time)
    .reshape(1, 1, -1)
    .repeat(n_shots, n_sources_y_per_shot, 1)
    .to(device)
)
source_locations_y = (
    torch.tensor([[[0, 0], [50, 50], [98, 98]],
                  [[deepwave.IGNORE_LOCATION, deepwave.IGNORE_LOCATION], [10, 10], [20, 20]]],
                 dtype=torch.long).to(device)
)
source_amplitudes_x = (
    deepwave.wavelets.ricker(freq, nt, dt, peak_time)
    .reshape(1, 1, -1)
    .repeat(n_shots, 1, 1)
    .to(device)
)
source_locations_x = torch.tensor([[[0, 98]], [[98, 0]]], dtype=torch.long).to(device)

out = deepwave.elastic(
    lamb, mu, buoyancy, grid_spacing=4, dt=0.004,
    source_amplitudes_y=source_amplitudes_y,
    source_amplitudes_x=source_amplitudes_x,
    source_locations_y=source_locations_y,
    source_locations_x=source_locations_x,
)

_, ax = plt.subplots(1, 2, figsize=(7, 3))
ax[0].imshow(out[0][0])
ax[1].imshow(out[0][1])
plt.show()

I hope that this answers your question.

[yulongovo]

Thank for your reply. I will learn this code.I encountered another issue during forward modeling. I used a homogeneous half-space model with S-wave velocity 300 m/s, P-wave velocity 600 m/s, and density 2000 kg/m³. The model size is 98 × 31 grid points. Receivers were placed at grid points 1–97, and I set one source at the trace 2. The simulated waveform is shown in the figure. You can see that the waveform changes on traces 4 and 5, which seems abnormal. I’m not sure whether there is a problem with the forward-modeling routine.

[ar4]

Would it be possible for you to post your code?

[yulongovo]

Please wait me for a few minutes.

[yulongovo]

#!/usr/bin/env python3

-- coding: utf-8 --

"""
Created on Thu Nov 13 22:20:49 2025

@author: myl
"""

import torch
import torch.nn as nn
import numpy as np
import deepwave
import scipy.io as spio

cuda_available = torch.cuda.is_available()
device = torch.device("cuda:0" if cuda_available else "cpu")

dx = 0.5 # step interval along x/z direction
dt = 0.00025 # time interval (e.g., 3ms)
num_dims = 2 # dimension of velocity model
nz = 31 # model shape of z dimension (depth)
ny = 98 # model shape of y dimension (when num_dims>=2)
total_t = 0.4 # totle sampling time (unit: s)
nt = int(total_t / dt)+1#1001# # number of time sampling points
peak_freq=20 # central frequency
peak_source_time =total_t#1.2/peak_freq #0.07 # the time (in secs)

true_vp=torch.ones(nz,ny,device=device) *600
true_vs=torch.ones(nz,ny,device=device) *300
true_dens=torch.ones(nz,ny,device=device) *2000

sx = torch.zeros(48, 1, 2,dtype=torch.int, device=device)
sx[..., 0]=0

for iis in range(48):
sx[iis,:,1]=iis*2+1

print(sx)

rx = torch.zeros(48, 97, 2,dtype=torch.int, device=device)
rx[..., 0]=0

for ii in range(48):
for jj in range(97):
rx[ii,jj,1]=jj
print(rx)

source_amplitudes_true=(deepwave.wavelets.ricker(peak_freq, 2*nt-1, dt, peak_source_time))
source_amplitudes_ref=source_amplitudes_true.repeat(48,1,1).to(device)

ref_amplitudes= deepwave.elastic(
*deepwave.common.vpvsrho_to_lambmubuoyancy(true_vp, true_vs,true_dens),
dx, dt,
source_amplitudes_y=source_amplitudes_ref,
source_locations_y=sx,
receiver_locations_y=rx,
pml_width=[0, 20, 30, 30],
)[-2]

[ar4]

Thank you for your code, which is very clean and easy to understand.

When I plot it, the change in the waveform seems like it might be less dramatic than in your plot:

Do you still think it might indicate a problem? When the top surface is absorbing, rather than a free surface, there is no change in the waveform:

[yulongovo]

I'm very curious why there's a sudden change at this position — the waveform at the fifth trace shows a clear drop. According to my understanding, this is a homogeneous half-space model, so such a noticeable change should not appear.

[ar4]

I was mistaken in my previous comment when I said "free surface". I should have said "rigid surface", since the width of the PML at the top surface is zero but there is no layer of zero properties at the surface to turn it into a free surface. If I create a layer of zero properties at the surface:

true_vp[0].fill_(0)
true_vs[0].fill_(0)
true_dens[0].fill_(0)

then the shift you noticed seems to go away and the result looks similar to that obtained when using Deepwave's previous elastic free surface implementation (the W-AFDA method, used up to v0.0.22 before being replaced in v0.0.23 with Zeng et al.'s Improved Vacuum Method):

[yulongovo]

I’ve tried it out. After adding those lines of code, the results indeed look much more natural and reasonable. By the way, is it still possible to download the previous version of Deepwave?

[ar4]

I am glad.

Yes, you can download previous versions, such as v0.0.22, using pip install deepwave==v0.0.22.

[ar4]

Note that, if you wish to compare with the W-AFDA method, there are some things to be aware of. One is that the W-AFDA method does not require setting properties above the free surface to zero to create a free surface. Simply setting the PML width to zero is enough. Another is that the location of the free surface differs between the two methods (W-AFDA does not have the "fictitious layer" that the Improved Vacuum Method does). The final one is that I also changed the staggered grid definition slightly between the two versions, so the different components are at slightly different locations.

[yulongovo]

Thank you for reminding me!

[yulongovo]

Hello! I'm glad to see the release of the 3D simulation. However, I still have some questions regarding the simulation. As you mentioned earlier, I need to set a layer of zeros at the free surface boundary for simulation. If I have an homogeneous half-space model (vs = 300, vp = 600, den = 2000) but make the first layer all zeros, my original model will change to a model where the velocity at the first layer grid points is zero, and everywhere else it's 300, 600, and 2000. If I perform waveform inversion, will the inversion results is a model where the velocity at the first layer grid points is zero, and everywhere else is 300, 600, and 2000? Or will I be able to obtain the original true model (vs = 300, vp = 600, den = 2000)?

[ar4]

Hello and thank you for your question.

I am unfortunately not completely certain what you are asking. Would it be possible for you to explain it again?

[yulongovo]

I'm sorry for the delay in replying to your message and sorry for not making myself clear.If we want to invert a homogeneous half-space model (vs = 300, vp = 600, den = 2000,we refer to this model as Model 1),we need to simulate the waveform data corresponding to this model as our observed data. As you mentioned earlier, I need to set a layer of zeros at the free surface boundary for simulation and everywhere else it's 300, 600, and 2000(we refer to this model as Model 2). So the observed data we obtain is based on Model 2.Then, wouldn't the final model obtained from inverting the observed data be Model 2, rather than Model 1?"