Question regarding implementing a source

[Costaki33]

Hi Alan, I'm trying to implement a source for the elastic wave solver. However, I am not sure what the units the elastic class/function is using.

Looking at the documentation:

source_amplitudes_y – A Tensor with dimensions [shot, source, time] containing time samples of the source wavelets for sources oriented in the first spatial dimension.

source_amplitudes_x – A Tensor containing source wavelet time samples for sources oriented in the second spatial dimension.

source_amplitudes_p – A Tensor containing source wavelet time samples for pressure sources.

_y and _x seem to be in meters? I am not so sure about _p.

Could you please provide more clarification?

Best,
Constantinos

[ar4]

Hi Constantinos,

One of the improvements planned for the next release is to make the description of what exactly the sources and receiver amplitudes represent clearer. As part of this, I have realised that there is a scaling factor missing from the recently added pressure sources that I will correct in the next release.

The directional body force sources are added to the velocity fields: v_new = v_old + dt * b * f (where b is the buoyancy and f is the user-supplied source amplitude). I think the units are therefore (m/s) * (1/s) * (kg / m^3) = kg / (s^2 * m^2) = N / m^3.

The pressure sources are added to the normal stress fields: sigma_ii_new = sigma_ii_old - dt * p / 2 (in the current release). I think the units are therefore (P) * (1/s) = kg / (s^3 * m) = N / (s * m^2).

It should be possible to construct arbitrary sources using only an array of directional body force sources, but you have to be careful with the grid locations of the different components because of the staggered grid.

[Costaki33]

Hi Alan,

Thanks for the response. I think I am still a bit confused. So as of right now, is it safe to assume that the units of the pressure sources are of Pa per second? What are the units used when using the ricker function? And in the sigma_ii_new = sigma_ii_old - dt * p / 2 , is p = pressure? Is the equation divided by 2 to have the change in stress to be equally distributed across the stress field in both the x and y directions?

I am essentially trying to use DW to create a ricker source that has a amplitude of 1 Pa, with a period of 0.01 to simulate a dilated stress source.

Do you have any advice on the conversion factor or the implementation of the ricker source with these units?

[ar4]

Hi Constantinos,

I'm sorry for the confusion. No one has asked about the precise details of the sources, or how to create one with a specific amplitude, before, so it is an area that is a bit lacking in Deepwave at the moment. I am glad that you are asking now, so that this part of Deepwave becomes more rigorous.

The scaling factor that I mentioned in my previous comment relates to the fact that the pressure sources and receivers in the elastic propagator do not seem to obey source-receiver reciprocity. I haven't looked into it in depth yet (I'm still working on finishing the 1D and 3D propagators), but the current implementation seems to be the standard one, so I might not change it (although I will remove the division by two that is currently applied to the pressure source amplitudes, as I don't think that is standard).

Yes, the pressure source is a pressure rate with units of Pa / s. The source amplitudes that you specify are added directly to the normal stresses in the update equation, after multiplication by -dt (and division by two in the current implementation, but this will be removed in the next release).

Using a wavelet that is the time derivative of the source wavelet you want should therefore work. When I tried it, however, the amplitude was half of what it should be when using the 3D propagator and only one third of what it should be when using the 2D propagator. I'm not sure yet where these factors come from. When I apply a scaling factor to compensate for this, the result looks quite good (the recorded pressure is the receiver recording for a receiver at the same location as the source):

In the current (2D-only) release, there is the additional factor of a half applied to the source amplitudes, so you would need to multiply the source amplitudes by six.

[ar4]

I realised that the recorded pressure is probably smaller than the integral over the injected pressure rate because the pressure is propagating away from the source, so I think it is behaving as expected.

Is that the kind of source that you wanted? Pressure sources (and other types of sources) can also be implemented using an array of directional body forces. As I mentioned before, that requires some care with the locations because of the staggered grid, but I would be happy to help you implement that if you wish.

[Costaki33]

Okay this is really helpful, thanks Alan. And thanks for being on top of this and updating the documentation so we can have more clarification of what is going on mathematically.

And yes, I believe that's the source I am looking for. So if I wanted to have a source with an amplitude of 1 Pa with a period of 0.01 at a dt = 0.0025s, would the working implementation be:
`
dt = 2.5e-4 # 0.00025 s
nt = 2040 # same as original cy 2000
ts1, ts2 = 0.0, 2.5e-3
f0 = 1.0 / 0.01 # 100 Hz
peak_time = 0.5 * (ts1 + ts2) # 0.00125 s (Ricker center)
t = torch.arange(nt, dtype=torch.float32, device=device) * dt # Time axis

amplitude_pa = 1.0
f0 = 1.0 / 0.01 # 1/T = f
peak_time = 0.5 * (ts1+ts2) # 1.25e-3 s # Single Ricker: amplitude=1 Pa, period=0.01 -> f0=100 Hz, peak centered in [ts1,ts2]

s_ricker = deepwave.wavelets.ricker(
freq=f0,
length=nt,
dt=dt,
peak_time=peak_time + 0.005, # Shift by 5ms to account for mtime offset
dtype=torch.float32,
).to(device)

s_derivative = torch.zeros_like(s_ricker)
s_derivative[1:-1] = (s_ricker[2:] - s_ricker[:-2]) / (2*dt)
s_derivative[0] = (s_ricker[1] - s_ricker[0]) / dt
s_derivative[-1] = (s_ricker[-1] - s_ricker[-2]) / dt

p = 6.0 * s_derivative # [nt], units: Pa/s
print(f"Final source peak (Pa/s): {p.abs().max().item():.1f}")
`

[Costaki33]

Another potential workaround that I found (might be incorrect - could use your guidance to make sure it's correct) is:

`dt = 2.5e-4 # 0.00025 s
nt = 2040 # same as original cy 2000
ts1, ts2 = 0.0, 2.5e-3
f0 = 1.0 / 0.01 # 100 Hz
peak_time = 0.5 * (ts1 + ts2) # 0.00125 s (Ricker center)
t = torch.arange(nt, dtype=torch.float32, device=device) * dt

amplitude_pa = 1.0
f0 = 1.0 / 0.01 # 1/T = f
peak_time = 0.5 * (ts1+ts2) # 1.25e-3 s

s_ricker = deepwave.wavelets.ricker(
freq=f0,
length=nt,
dt=dt,
peak_time=peak_time + 0.005, # Shift by 5ms to account for mtime offset
dtype=torch.float32,
).to(device)

s_derivative = torch.zeros_like(s_ricker)
s_derivative[1:-1] = (s_ricker[2:] - s_ricker[:-2]) / (2*dt)
s_derivative[0] = (s_ricker[1] - s_ricker[0]) / dt
s_derivative[-1] = (s_ricker[-1] - s_ricker[-2]) / dt

p0 = s_derivative # Pa/s up to an unknown scale
p0_cpu = p0.detach().cpu().numpy()
dt_float = float(dt)
sigma_test = -0.5 * dt_float * np.cumsum(p0_cpu)
sigma_peak = np.abs(sigma_test).max()
print("Unscaled stress peak (Pa):", sigma_peak)

amplitude_pa = 1.0
k = amplitude_pa / sigma_peak
print("Scale factor k:", k)

p = -k * s_derivative # [nt], units Pa/s # negative because we want a dilational stress source
source_amplitudes_p = p.view(1, 1, -1)

p_cpu = p.detach().cpu().numpy()
sigma_dw = -0.5 * dt_float * np.cumsum(p_cpu)
print("New stress peak (Pa):", np.abs(sigma_dw).max())

Unscaled stress peak (Pa): 0.5582149
Scale factor k: 1.7914247628370796
New stress peak (Pa): 0.9999999
`

[ar4]

I think there are two ways of correctly implementing the type of source that you want. The first is to use the time derivative of the pressure source you want as the amplitudes for the pressure rate ("p") type source. The second is to use dipole body force sources in each spatial dimension. I demonstrate both below, and confirm that they match the expected far-field pressure. I use the (not yet released) 3D propagator as the analytical solution is easier to calculate in 3D. In v0.0.23 you would need to multiply the source amplitudes by 2 to compensate for the division that has since been removed.

Also, I was overly optimistic in predicting that the next version would be released in the middle of the month. It is working, but still needs more polishing. I'm sorry for the delay.

import torch
import matplotlib.pyplot as plt
import deepwave

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
nz = 11
ny = 11
nx = 31
dx = 1.0
nt = 100
dt = 0.001
freq = 100.0
t_0 = 1.2 / freq

vp = 2000.0  # P-wave velocity (m/s)
vs = 1200.0  # S-wave velocity (m/s)
rho = 2200.0  # Density (kg/m^3

vp_tensor = torch.full((nz, ny, nx), vp, device=device)
vs_tensor = torch.full((nz, ny, nx), vs, device=device)
rho_tensor = torch.full((nz, ny, nx), rho, device=device)

lamb, mu, buoyancy = deepwave.common.vpvsrho_to_lambmubuoyancy(vp_tensor, vs_tensor, rho_tensor)

dsr = nx - 2  # Source-receiver distance (in cells)
receiver_locations_p = torch.tensor([[[nz // 2, ny // 2, 1 + dsr]]], dtype=torch.long, device=device)
ricker = deepwave.wavelets.ricker(freq, nt, dt, t_0).reshape(1, 1, -1).to(device)

def _shift(f: torch.Tensor, t_shift: float) -> torch.Tensor:
    """Shift by a fractional amount in time."""
    nt = f.shape[-1]
    w = torch.fft.rfftfreq(nt, 1.0)
    fw = torch.fft.rfft(f, dim=-1)
    fw_shifted = fw * torch.exp(-1j * 2 * torch.pi * w * t_shift)
    return torch.fft.irfft(fw_shifted, nt, dim=-1)

r = dsr * dx  # Source-receiver distance in metres
ricker_tt = torch.nn.functional.pad((ricker[..., 2:] - 2*ricker[..., 1:-1] + ricker[..., :-2]) / dt**2, (1, 1))
expected = _shift(ricker_tt * (vp**2-4/3*vs**2)/vp**4 / (4 * torch.pi * r) * dx**3, r / vp / dt)

# 1: Pressure rate source
# The rate of pressure injection is specified in units of Pa / s.
# We use the time derivative of the desired pressure source as the rate.

source_locations_p = torch.tensor([[[nz // 2, ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_p = torch.gradient(ricker, spacing=dt, dim=-1)[0]

out_1 = deepwave.elastic(
    lamb, mu, buoyancy,
    grid_spacing=dx,
    dt=dt,
    source_amplitudes_p=source_amplitudes_p,
    source_locations_p=source_locations_p,
    receiver_locations_p=receiver_locations_p,
    pml_freq=freq,
)[-4].cpu()

# 2: Dipole body force sources
# Dipole body forces in each dimension are centred on the desired pressure source location.
# Body forces are injected on the corresponding velocity component.
# Because of the staggered grid (see https://ausargeo.com/deepwave/elastic), the positive
# dipole component is in the same grid cell as the pressure source location, while the
# negative dipole components are in the grid cell before it. This will create a dipole
# centred on the normal stress component location.
# Directional body forces are specified in force/volume.

source_locations_z = torch.tensor([[[nz // 2 - 1, ny // 2, 1],
                                    [nz // 2, ny // 2, 1]]], dtype=torch.long, device=device)
source_locations_y = torch.tensor([[[nz // 2, ny // 2 - 1, 1],
                                    [nz // 2, ny // 2, 1]]], dtype=torch.long, device=device)
source_locations_x = torch.tensor([[[nz // 2, ny // 2, 1 - 1],
                                    [nz // 2, ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_dipole = torch.cat([-ricker / dx, ricker / dx], dim=-2)

out_2 = deepwave.elastic(
    lamb, mu, buoyancy,
    grid_spacing=dx,
    dt=dt,
    source_amplitudes_z=source_amplitudes_dipole,
    source_amplitudes_y=source_amplitudes_dipole,
    source_amplitudes_x=source_amplitudes_dipole,
    source_locations_z=source_locations_z,
    source_locations_y=source_locations_y,
    source_locations_x=source_locations_x,
    receiver_locations_p=receiver_locations_p,
    pml_freq=freq,
)[-4].cpu()

plt.plot(expected.flatten(), linewidth=8)
plt.plot(out_1.flatten(), linewidth=4)
plt.plot(out_2.flatten())

[Costaki33]

Hi Alan, thanks for the code, this is really helpful. And no worries about the 3D implementation; the community and I are patient with the eventual release and very much looking forward to it when you're ready! I implemented this code in 2D and I got some slightly different results. I had to make some slight modifications to make the code work in the current code version:

import torch
import deepwave
import matplotlib.pyplot as plt
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

ny = 11
nx = 31
dx = 1.0
nt = 100
dt = 0.001
freq = 100.0 # Hz
t_0 = 1.2 / freq # peak time -> want peak to be at 1.2s, so divide by frequency to get time for peak_time
t = torch.arange(nt, dtype=torch.long, device=device) * dt

vp = 2000.0 # P-wave velocity (m/s)
vs = 1200.0 # S-wave velocity (m/s)
rho = 2200.0 # Density (kg/m^3

vp_tensor = torch.full((ny, nx), vp, device=device)
vs_tensor = torch.full((ny, nx), vs, device=device)
rho_tensor = torch.full((ny, nx), rho, device=device)

lamb, mu, buoyancy = deepwave.common.vpvsrho_to_lambmubuoyancy(vp_tensor, vs_tensor, rho_tensor)

dsr = nx - 2 # Source-receiver distance (in cells)
receiver_locations_p = torch.tensor([[[ny // 2, 1 + dsr]]], dtype=torch.long, device=device)
ricker = deepwave.wavelets.ricker(freq, nt, dt, t_0).reshape(1, 1, -1).to(device)

def _shift(f: torch.Tensor, t_shift: float) -> torch.Tensor:
"""
Shift by a fractional amount in time.
Transforms the input time-series wavelet into a frequency domain,
rotates each frequency component by a phase proportional to the desired time shift, and
transforms it back into the time domain
"""
nt = f.shape[-1]
w = torch.fft.rfftfreq(nt, 1.0).to(device) # gives the frequency grid for the real-values FFT, w = cycles per sample
fw = torch.fft.rfft(f, dim=-1) # compute real-valued FFT along time axis, fw = complex spectrum of f
fw_shifted = fw * torch.exp(-1j * 2 * torch.pi * w * t_shift) #
return torch.fft.irfft(fw_shifted, nt, dim=-1)

r = dsr * dx # Source-receiver distance in metres
ricker_tt = torch.nn.functional.pad((ricker[..., 2:] - 2*ricker[..., 1:-1] + ricker[..., :-2]) / dt**2, (1, 1)) # second time derivative of the ricker pulse

expected = _shift(
ricker_tt
* (vp2 - 4.0/3.0 * vs2) / vp4
/ (4.0 * torch.pi * r) # 3D-like decay; treat as a scale factor in 2D
* dx2, # area element instead of volume
r / vp / dt
)

source_locations_p = torch.tensor([[[ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_p = torch.gradient(ricker, spacing=dt, dim=-1)[0]

out_1 = deepwave.elastic(
lamb, mu, buoyancy,
grid_spacing=dx,
dt=dt,
source_amplitudes_p=source_amplitudes_p,
source_locations_p=source_locations_p,
receiver_locations_p=receiver_locations_p,
pml_freq=freq,
)

vy1, vx1, sigmayy1, sigmaxy1, sigmaxx1, m_vyy1, m_vyx1, m_vxy1, m_vxx1, m_sigmayyy1, m_sigmaxyy1, m_sigmaxyx1, m_sigmaxxx1, rec_p1, rec_y1, rec_x1 = out_1

source_locations_y = torch.tensor([[[ny // 2 - 1, 1],
[ny // 2, 1]]], dtype=torch.long, device=device)
source_locations_x = torch.tensor([[[ny // 2, 1 - 1],
[ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_dipole = torch.cat([-ricker / dx, ricker / dx], dim=-2)

out_2 = deepwave.elastic(
lamb, mu, buoyancy,
grid_spacing=dx,
dt=dt,
source_amplitudes_y=source_amplitudes_dipole,
source_amplitudes_x=source_amplitudes_dipole,
source_locations_y=source_locations_y,
source_locations_x=source_locations_x,
receiver_locations_p=receiver_locations_p,
pml_freq=freq,
)

vy2, vx2, sigmayy2, sigmaxy2, sigmaxx2, m_vyy2, m_vyx2, m_vxy2, m_vxx2, m_sigmayyy2, m_sigmaxyy2, m_sigmaxyx2, m_sigmaxxx2, rec_p2, rec_y2, rec_x2, = out_2

e = expected[0, 0].cpu().numpy() # analytic
u1 = rec_p1[0, 0].cpu().numpy() # pressure-rate
u2 = rec_p2[0, 0].cpu().numpy() # dipole

plt.figure(figsize=(6,4))
plt.plot(e, label="expected", linewidth=8)
plt.plot(u1, label="pressure-rate", linewidth=4)
plt.plot(u2, label="dipole")
plt.legend()
plt.grid(True)
plt.show()

Do you have any advice on the _shift() with its implementation in 2D? I feel that is why my results are different from yours.

[ar4]

Hi Constantinos,

2D wave propagation changes the waveform. Calculating the effect of this on the analytic solution is complicated, but here is a modification of the example to use an approximation of the expected solution.

import torch
import matplotlib.pyplot as plt
import scipy
import math
import deepwave

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
ny = 11
nx = 31
dx = 0.5
nt = 100
dt = 0.001
freq = 100.0
t_0 = 1.2 / freq

vp = 2000.0  # P-wave velocity (m/s)
vs = 1200.0  # S-wave velocity (m/s)
rho = 2200.0  # Density (kg/m^3

vp_tensor = torch.full((ny, nx), vp, device=device)
vs_tensor = torch.full((ny, nx), vs, device=device)
rho_tensor = torch.full((ny, nx), rho, device=device)

lamb, mu, buoyancy = deepwave.common.vpvsrho_to_lambmubuoyancy(vp_tensor, vs_tensor, rho_tensor)

dsr = nx - 2  # Source-receiver distance (in cells)
receiver_locations_p = torch.tensor([[[ny // 2, 1 + dsr]]], dtype=torch.long, device=device)
ricker = deepwave.wavelets.ricker(freq, nt, dt, t_0).reshape(1, 1, -1).to(device)

r = dsr * dx  # Source-receiver distance in metres

# Estimate expected response in 2D medium
w = torch.fft.rfftfreq(nt, dt)
fw = torch.fft.rfft(ricker, dim=-1)
omega = 2 * math.pi * w
h = scipy.special.hankel1(0, -omega.numpy() * r / vp)
g = 1j / 4 * torch.tensor(h)
g[0] = 0
S = -g * fw * dx**2 * omega**2
s = torch.fft.irfft(S, nt, dim=-1)
expected = s * (vp**2-vs**2)/vp**4

# 1: Pressure rate source
# The rate of pressure injection is specified in units of Pa / s.
# We use the time derivative of the desired pressure source as the rate.

source_locations_p = torch.tensor([[[ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_p = torch.gradient(ricker, spacing=dt, dim=-1)[0]

out_1 = deepwave.elastic(
    lamb, mu, buoyancy,
    grid_spacing=dx,
    dt=dt,
    source_amplitudes_p=source_amplitudes_p,
    source_locations_p=source_locations_p,
    receiver_locations_p=receiver_locations_p,
    pml_freq=freq,
)[-3].cpu()

# 2: Dipole body force sources
# Dipole body forces in each dimension are centred on the desired pressure source location.
# Body forces are injected on the corresponding velocity component.
# Because of the staggered grid (see https://ausargeo.com/deepwave/elastic), the positive
# dipole component is in the same grid cell as the pressure source location, while the
# negative dipole components are in the grid cell before it. This will create a dipole
# centred on the normal stress component location.
# Directional body forces are specified in force/volume.

source_locations_y = torch.tensor([[[ny // 2 - 1, 1],
                                    [ny // 2, 1]]], dtype=torch.long, device=device)
source_locations_x = torch.tensor([[[ny // 2, 1 - 1],
                                    [ny // 2, 1]]], dtype=torch.long, device=device)
source_amplitudes_dipole = torch.cat([-ricker / dx, ricker / dx], dim=-2)

out_2 = deepwave.elastic(
    lamb, mu, buoyancy,
    grid_spacing=dx,
    dt=dt,
    source_amplitudes_y=source_amplitudes_dipole,
    source_amplitudes_x=source_amplitudes_dipole,
    source_locations_y=source_locations_y,
    source_locations_x=source_locations_x,
    receiver_locations_p=receiver_locations_p,
    pml_freq=freq,
)[-3].cpu()

plt.plot(expected.flatten(), linewidth=8)
plt.plot(out_1.flatten(), linewidth=4)
plt.plot(out_2.flatten())

[Costaki33]

Hi Alan,

I see, thanks so much for this. I am happy that my pressure-rate and dipole implementations I had were replicated by you. You're right that the implementation in 2D is quite complicated with using Green's function, but I'm glad that with your help I am on the right track.

I'm sorry to ask again; I took your code had it run on my machine and it produced slightly different results:

I'm not sure if I am doing anything wrong, I didn't change any parameter values or code. What versions of CUDA, NumPy, SciPy, and DeepWave are you using? It may be the solvers in the backend that are causing these discrepancies. I am using:

PyTorch version: 2.7.1+cu126
Scipy version: 1.15.3
Device: cuda
Max out_1: 0.0005880516255274415, Max out_2: 0.0013677392853423953
Max expected: 0.0013705615419894457

[ar4]

Hi Constantinos,

In v0.0.23 I unfortunately had an extra factor of 0.5 multiplied by the pressure sources (corrected in the upcoming v0.0.24), so for the pressure source method above you will need to multiply the source amplitudes by 2.

[Costaki33]

Ah! Got it! Thanks Alan!

Thanks for taking the time for answering my questions, I really appreciate it. My team at the Texas Seismological Network are trying to use DW as our solver for our FWI solution so this has been really helpful for comparing our existing solver's results with those from DW - which are really good.

[Costaki33]

Hey Alan,

Congrats on the pre-release! My team and I are looking forward to testing it against our solvers to make DW better.

I have a question again:

I am trying to ensure that the receiver waveform outputs from 2D

DW matches the existing solution we have using another solver, WAVE3D (is 2D capable). We have a reference wavelet from a sample problem using WAVE3D, which we are plugging into DW using the same geometry as our reference model. In theory, the receivers in the model geometry should receiver from both DW and WAVE roughly the exact waveforms.

Our input Ricker wavelet looks correct:

However, we are seeing discrepancies at the receiver locations. For example, our source is located at (200, 190) in our 400x400 grid. When we plot receiver data at, for instance (350, 365), Our arrivals are matching in time but we are seeing a drift of the energy of the wavelet in time. I attached some sample images:

We also plotted the Fast Fourier Transform of both the reference and outputted waveforms to compare the frequency contents:

We are really not sure why we are seeing a difference between the two results. We suspect that there is a source space shift issue and a frequency shift issue but not sure how to debug this.

Similarly, we are noticing a scaling factor issue between the two results for them to be comparable (we have to utilize a scale object in our plotting). Scaling factor values can be has high as x50e6 when very close to the source (210, 190), and decreasing "linearly" as we move away from the source (x21e6 at receiver location (250, 190), x15e6 at receiver location (290, 190), x12e6 at receiver location (350, 190), and x11e6 at (110, 365)). We believe this may be because of a time dilation issue but we're not sure at the moment where in the source code that could be coming from.

Could you please help us decipher why there is a discrepancy between our two results even when using the same source as input into the elastic solver? I've attached the source CSV file and supporting code to use the below code that I have built to run DW on our reference data.

We think this is really important for validating DW's results and would really appreciate your feedback, as this is beyond our expertise.

If I can be any clearer on anything, please let me know. I am using v0.0.23 of DW.

Thanks,
Constantinos

import os
import math
import time
import torch
import scipy
import deepwave
import matplotlib
matplotlib.use("Agg")   # headless-safe backend
import numpy as np
import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
from deepwave.wavelets import ricker
import matplotlib.animation as animation
from deepwave.common import vpvsrho_to_lambmubuoyancy
from wave3d_parser import parse_histories, parse_config

param_file =  # './seis_ai/seis_ai_m.txt'
seis_csv   =  # './seis_ai/seis_ai.csv'

%matplotlib widget
dtype  = torch.float32
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
torch.set_printoptions(threshold=10000)
np.set_printoptions(threshold=10000)


# ---------------------------------------------------------
# Load Wave3D data and config (needed before we define source_1d)
# ---------------------------------------------------------
config     = parse_config(filepath=param_file)
v_waveform = parse_histories(csv_file=seis_csv)   # shape (nt, n_cols)

time = v_waveform[:, 0]
n_cols = v_waveform.shape[1]
n_receivers = (n_cols - 2) // 2  # 2 columns per receiver (y,x)

print(f"Wave3D: nt={v_waveform.shape[0]}, n_receivers={n_receivers}")

# ---------------------------------------------------------
# Model Parameters (400 x 400, homogeneous)
# ---------------------------------------------------------
ny, nx = config.grid.ny, config.grid.nx
sy, sx = config.source.j, config.source.i  # shot / logical center
dx = config.grid.dx          # m
nt = config.cycles
dt = config.timestep       # s
pml_freq = 1.0 / config.source.period       # Hz
ts1, ts2 = config.source.ts1, config.source.ts2
peak_time = (ts1 + ts2)           # 0.00125 s
t_0 = peak_time
pml_width = config.boundary.pml_thickness

# Velocity Model
vp  = config.material.cp   # P-wave velocity (m/s)
vs  = config.material.cs   # S-wave velocity (m/s)
rho = config.material.density   # Density (kg/m^3)

vp_tensor  = torch.full((ny, nx), vp,  device=device)
vs_tensor  = torch.full((ny, nx), vs,  device=device)
rho_tensor = torch.full((ny, nx), rho, device=device)

lamb, mu, buoyancy = deepwave.common.vpvsrho_to_lambmubuoyancy(
    vp_tensor, vs_tensor, rho_tensor
)

# ---------------------------------------------------------
# Receiver locations
# ---------------------------------------------------------
OFFSET = 0   # change to -1 if your (i,j) start at 1

def to_locations(coords):
    arr = torch.tensor(coords, dtype=torch.long, device=device)
    if OFFSET != 0:
        arr = arr + OFFSET  # apply -1 if converting 1-based -> 0-based
    # add shot dimension
    return arr.unsqueeze(0)   # [1, N, 2]

# y-velocity receiver coordinates (note: flip order to y, x from x, y)
y_coords = [hist.yvel[1::-1] for hist in config.history_points]

# x-velocity receiver coordinates (note: flip order to y, x from x, y)
x_coords = [hist.xvel[1::-1] for hist in config.history_points]

receiver_locations_y = to_locations(y_coords)
receiver_locations_x = to_locations(x_coords)
receiver_locations_p = receiver_locations_y

if n_receivers != len(y_coords):
    raise ValueError(
        f"CSV has {n_receivers} receivers but y_coords has {len(y_coords)}. "
        "They must match in order and count."
    )

# ---------------------------------------------------------
# Spatial regularization kernel
# ---------------------------------------------------------
def regularize_point_source(
    source_1d: torch.Tensor,
    center: tuple[int, int],
    ny: int,
    nx: int,
    *,
    sigma_grid: float = 1.5,
    method: str = "gaussian",
    radius_scale: float = 1.5,
    preserve_integral: bool = True,
    device=None,
) -> torch.Tensor:
    """
    Take a 1D time series 'source_1d' and spread it in space around `center`
    on a (ny, nx) grid using a 2D kernel (Gaussian or B-spline).

    Returns
    -------
    sourcing_field : (T, ny, nx) tensor on `device`
        Time-varying source field suitable for converting into Deepwave
        source_amplitudes_p / source_locations_p.
    """
    if device is None:
        device = source_1d.device

    T = int(source_1d.shape[0])
    cy, cx = center

    # (T,H,W) tensor with a pure point source at (cy,cx)
    sourcing_field = torch.zeros(T, ny, nx, device=device, dtype=torch.float32)
    sourcing_field[:, cy, cx] = source_1d.to(device).to(torch.float32)

    # --- spatial kernel only ---
    method = method.lower()
    if method == "gaussian":
        kernel_radius = int(math.ceil(4 * sigma_grid))
        kernel_size   = 2 * kernel_radius + 1

        y, x = torch.meshgrid(
            torch.arange(kernel_size, device=device, dtype=torch.float32) - kernel_radius,
            torch.arange(kernel_size, device=device, dtype=torch.float32) - kernel_radius,
            indexing="ij",
        )
        kernel = torch.exp(-(x**2 + y**2) / (2.0 * sigma_grid**2))

    elif method == "bspline":
        kernel_radius = int(math.ceil(2 * radius_scale))
        kernel_size   = 2 * kernel_radius + 1

        y, x = torch.meshgrid(
            torch.arange(kernel_size, device=device, dtype=torch.float32) - kernel_radius,
            torch.arange(kernel_size, device=device, dtype=torch.float32) - kernel_radius,
            indexing="ij",
        )
        r = torch.sqrt(x**2 + y**2)
        r_norm = r / radius_scale

        kernel = torch.zeros_like(r)
        region1 = r_norm < 1.0
        region2 = (r_norm >= 1.0) & (r_norm < 2.0)

        q1 = r_norm[region1]
        kernel[region1] = (4.0 - 6.0 * q1**2 + 3.0 * q1**3)

        q2 = r_norm[region2]
        kernel[region2] = (2.0 - q2) ** 3
    else:
        raise ValueError(f"Unknown kernel method {method!r}")

    if preserve_integral:
        kernel = kernel / kernel.sum()

    kernel = kernel.view(1, 1, kernel.shape[0], kernel.shape[1])

    field_channels = sourcing_field.unsqueeze(1)  # (T,1,H,W)
    kernel_h, kernel_w = kernel.shape[-2], kernel.shape[-1]
    padding = (kernel_h // 2, kernel_w // 2)

    smoothed = torch.nn.functional.conv2d(
        field_channels,
        kernel,
        padding=padding,
        stride=1,
    )  # (T,1,H,W)

    smoothed = smoothed.squeeze(1)  # (T,H,W)
    return smoothed

# ---------------------------------------------------------
# 1D temporal source: USE REFERENCE WAVELET FROM Wave3D
# ---------------------------------------------------------
# Choose which receiver/component to treat as "reference wavelet"
src_irec       = 1      # receiver index in the CSV to use as wavelet
src_component  = "x"    # "x" or "y"

if src_irec < 0 or src_irec >= n_receivers:
    raise IndexError(f"src_irec {src_irec} out of range [0, {n_receivers-1}]")

col_base_src = 2 + 2 * src_irec
if src_component == "y":
    src_col = col_base_src
elif src_component == "x":
    src_col = col_base_src + 1
else:
    raise ValueError("src_component must be 'x' or 'y'")

# Reference source wavelet from Wave3D
source_1d_np = v_waveform[:nt, src_col]
source_1d    = torch.from_numpy(source_1d_np).to(device).to(torch.float32)

# Optional: remove any DC offset (this is a shift, not a scaling)
source_1d = source_1d - source_1d.mean()
source_1d_centered_np = source_1d.detach().cpu().numpy()

print("Using reference wavelet from CSV:")
print(f"  src_irec     = {src_irec}")
print(f"  src_component= {src_component}")
print(f"  source_1d peak abs (centered) = {float(source_1d.abs().max()):.6e}")

# Plot the reference wavelet AND the actual Deepwave source wavelet
plt.figure(figsize=(8,4))
plt.plot(time[:nt], source_1d_np, linewidth=2, label="Reference wavelet (CSV)")
plt.plot(time[:nt], source_1d_centered_np, linewidth=2, linestyle="--",
         label="Wavelet used as Deepwave source (mean-removed)")
plt.xlabel("Time (s)")
plt.ylabel("Velocity")
plt.title(f"Reference wavelet vs Deepwave source (rec {src_irec}, {src_component})")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

# ---------------------------------------------------------
# 3D regularized source field (T, ny, nx)
# ---------------------------------------------------------
kernel_type = 'gaussian'
sigma_grid  = 2.0
radius_scale = 0.1   # very tight; you can adjust

sourcing_field = regularize_point_source(
    source_1d=source_1d,
    center=(sy, sx),         # (cy, cx) in index space
    ny=ny,
    nx=nx,
    sigma_grid=sigma_grid,
    method=kernel_type,
    radius_scale=radius_scale,
    preserve_integral=True,
    device=device,
)

# ---------------------------------------------------------
# Energy match: unsmoothed vs smoothed (L2 in space–time)
# ---------------------------------------------------------
dx_phys = dx
dt_phys = dt

sourcing_field_ref = torch.zeros_like(sourcing_field)
sourcing_field_ref[:, sy, sx] = source_1d.to(sourcing_field_ref.device)

E_ref = (sourcing_field_ref ** 2).sum() * dx_phys * dx_phys * dt_phys
E_dw  = (sourcing_field      ** 2).sum() * dx_phys * dx_phys * dt_phys

if E_dw <= 0:
    raise RuntimeError("DW source field has zero or negative energy; check construction.")

scale_energy = torch.sqrt(E_ref / E_dw)
sourcing_field = sourcing_field * scale_energy

print(f"E_ref (unsmoothed): {E_ref.item():.6e}")
print(f"E_dw  (smoothed)  : {E_dw.item():.6e}")
print(f"Energy scale applied to smoothed field: {scale_energy.item():.6e}")
print(f"Smoothed field peak AFTER scaling: {float(sourcing_field.abs().max()):.3e}")

# ---------------------------------------------------------
# Convert regularized field to Deepwave sparse source
# ---------------------------------------------------------
nonzero_mask   = (sourcing_field.abs().sum(dim=0) > 1e-20)   # (ny,nx)
source_indices = torch.nonzero(nonzero_mask, as_tuple=False) # (n_src,2)[y,x]
n_src          = source_indices.shape[0]
print(f"Number of active source points: {n_src}")

source_amplitudes_list = []
for i in range(n_src):
    y_idx, x_idx = source_indices[i]
    trace_i = sourcing_field[:, y_idx, x_idx]  # (nt,)
    source_amplitudes_list.append(trace_i)

source_amplitudes_p = torch.stack(source_amplitudes_list, dim=0)  # (n_src, nt)
source_amplitudes_p = source_amplitudes_p.unsqueeze(0)            # (1, n_src, nt)

source_locations_p = source_indices.unsqueeze(0)                  # (1, n_src, 2)

print("source_amplitudes_p shape:", source_amplitudes_p.shape)
print("source_locations_p shape :", source_locations_p.shape)
print("Max amplitude per source :", source_amplitudes_p[0].abs().max(dim=1)[0][:5])

# ---------------------------------------------------------
# Collect snapshots during propagation
# ---------------------------------------------------------
snapshots = []                   # will store CPU tensors [ny, nx]
steps     = []                   # deepwave internal step index (integer)
snapshot_stride    = 1
callback_frequency = 1

def forward_cb(state):
    vy_t = state.get_wavefield("vy_0")[0]   # [ny, nx] on device
    vx_t = state.get_wavefield("vx_0")[0]
    if state.step % snapshot_stride == 0:
        vabs = torch.sqrt(vy_t*vy_t + vx_t*vx_t).detach().to("cpu").clone()
        snapshots.append(vabs)
        steps.append(int(state.step))

(
    vy_f1, vx_f1,
    sigmayy_f1, sigmaxy_f1, sigmaxx_f1,
    m_vyy_f1, m_vyx_f1, m_vxy_f1, m_vxx_f1,
    m_sigmayyy_f1, m_sigmaxyy_f1, m_sigmaxyx_f1, m_sigmaxxx_f1,
    rec_p1, rec_vy1, rec_vx1,
) = deepwave.elastic(
    lamb, mu, buoyancy,
    grid_spacing=dx,
    dt=dt,
    accuracy=4,
    source_amplitudes_p=source_amplitudes_p,
    source_locations_p=source_locations_p,
    receiver_locations_p=receiver_locations_p,
    receiver_locations_x=receiver_locations_x,
    receiver_locations_y=receiver_locations_y,
    pml_width=pml_width,
    pml_freq=pml_freq,
    forward_callback=forward_cb,
    callback_frequency=callback_frequency,
)



def get_trace(irec:int, component:str='y', N:int=2000):
     # Wave3D reference: columns [time, maxvel, y0, x0, y1, x1, ...]
    col_base = 2 + (2 * irec)
    if component == "y":
        col_expected = col_base
    else:
        col_expected = col_base + 1
    return v_waveform[:N, col_expected]

# ---------------------------------------------------------
# Plot receiver waveforms: Deepwave vs Wave3D
# ---------------------------------------------------------
def plot_receiver_waveform_velocity(
    irec: int,
    component: str,
    N: int,
    time: np.ndarray,
    rec_vx1: torch.Tensor,
    rec_vy1: torch.Tensor,
    x_coords: list[tuple[int,int]],
    y_coords: list[tuple[int,int]],
    kernel_type: str,
    sigma_grid: float,
    receiver_locations_p: torch.Tensor,
):
    """
    Compare Wave3D vs Deepwave at receiver index `irec`.

    component: 'x' or 'y'
    N: number of time samples to plot
    """

    if irec < 0 or irec >= receiver_locations_p.shape[1]:
        raise IndexError(f"irec {irec} out of range [0, {receiver_locations_p.shape[1]-1}]")

    if component not in ("x", "y"):
        raise ValueError("component must be 'x' or 'y'")

    # --- choose Deepwave array / coord list ---
    if component == "x":
        rec_dw = rec_vx1  # [1, nrec_x, nt]
        coords_list = x_coords
        vel_label = "vx"
    else:
        rec_dw = rec_vy1
        coords_list = y_coords
        vel_label = "vy"

    # Deepwave trace
    rec_dw_np = rec_dw[0, irec, :N].detach().cpu().numpy()

    expected_np = get_trace(irec, component=component, N=N)
    t_axis = time[:N]

    # For plotting only: rescale DW to match peak of reference
    rec_dw_tensor = rec_dw[0, irec, :N].detach()
    expected_tensor = torch.from_numpy(expected_np).to(rec_dw_tensor.device)
    if rec_dw_tensor.abs().max() > 0:
        scale = expected_tensor.abs().max() / rec_dw_tensor.abs().max()
    else:
        scale = 1.0

    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(1,1,1)
    ax.plot(t_axis, expected_np, linewidth=2, label="Wave3D ref")
    ax.plot(t_axis, rec_dw_np * float(scale), linewidth=2, label="Deepwave (scaled)")
    ax.set_ylabel(f"Velocity ({vel_label}): scale={scale:.1e}")
    ax.set_xlabel("Time (s)")
    ax.set_title(
        f"Trace at grid location {coords_list[irec]} "
        f"({vel_label}), kernel={kernel_type}, σ={sigma_grid}"
    )
    ax.grid(True, alpha=0.3)
    ax.legend()
    fig.tight_layout()
    plt.show()

dw_waveform = rec_vx1[0, irec]
gt_waveform = torch.from_numpy(get_trace(irec, component='y', N=dw_waveform.shape[0])).to(device)
def get_fft(wf):
    fft = torch.fft.rfft(wf)
    spectrum = torch.abs(fft)
    frequencies = torch.fft.rfftfreq(wf.shape[0], d=dt)
    return frequencies, spectrum

gt_freq, gt_sepc = get_fft(gt_waveform)
dw_freq, dw_spec = get_fft(dw_waveform)

scale = gt_sepc.max() / dw_spec.max()

# Plotting the magnitude spectrum
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1,1,1)
ax.plot(gt_freq.cpu().numpy(), gt_sepc.cpu().numpy(), linewidth=2, label="Wave3D ref")
ax.plot(dw_freq.cpu().numpy(), (dw_spec*scale).cpu().numpy(), linewidth=2, label="Deepwave")
ax.set_ylabel(f"Magnitude (scale={scale:.1e})")
ax.set_xlabel("Frequency (Hz)")
ax.set_title(f"Trace at grid location yvel {y_coords[irec]}")
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()

seis_ai_m.txt
seis_ai.csv
wave3d_parser.py

[ar4]

Hi Constantinos,

This is great. It's a big help for you to do these comparisons and validations. Thank you for doing it, and thank you also for the clear code that I was able to run to reproduce the results.

I don't have any immediate explanation for the difference, other than suggesting testing it with a smaller grid spacing to see whether numerical dispersion might be the cause. There seems to be quite a few people eagerly waiting for the release of the 3D propagator, so I propose that I work on finishing that and then I will dedicate some time to investigating this in more depth.

Thank you again.

[Costaki33]

Hi Alan,

Great, that sounds like a plan!