same sources + different order = different results

[joao-bapdm]

Doing a simulation with the same sources leads to different results depending on the order sources are added.

Using the 2d example from the example/ directory:

from pywave import *
import numpy as np

# shape of the grid
shape = (512, 512)

# spacing
spacing = (15.0, 15.0)

# propagation time
time = 1500

# Velocity model
vel = np.zeros(shape, dtype=np.float32)
vel[:] = 1500.0
vel[200:] = 2000.0
velModel = Model(ndarray=vel)

# Compiler
compiler = Compiler(cc="gcc", cflags="-O3 -shared")

# domain extension (damping + spatial order halo)
extension = BoundaryProcedures(
    nbl=50,
    boundary_condition=(("NN", "NN"), ("NN", "NN")),
    damping_polynomial_degree=3,
    alpha=0.0001,
)

# Wavelet
wavelet = Wavelet(frequency=5.0)

# Source
source = Source(kws_half_width=1, wavelet=wavelet)
source.add(position=(270, 240))
source.add(position=(255.5, 255.5))

# receivers
receivers = Receiver(kws_half_width=1)
for i in range(512):
    receivers.add(position=(255.5, i))

setup = Setup(
    velocity_model=velModel,
    sources=source,
    receivers=receivers,
    boundary_config=extension,
    spacing=spacing,
    propagation_time=time,
    jumps=0,
    compiler=compiler,
    space_order=2,
)

solver = AcousticSolver(setup=setup)

wavefields, rec, exec_time = solver.forward()

print(wavefields.shape)

if len(wavefields.shape) > 2:
    count = 0
    for wavefield in wavefields:
        plot_wavefield(wavefield, file_name="arq-" + str(count))
        count += 1
else:
    plot_wavefield(wavefields)

print("Forward execution time: %f seconds" % exec_time)

plot_shotrecord(rec)

swapping lines 34 and 35 leads to different shot records:

Adding first (255.5, 255.5), then (270, 240):

Adding first (270, 240), then (255.5, 255.5):

Unfortunately the bug does not come from my latest pull request, as I generated these images in commit 0b1b157.

[krober10nd]

It appears at first glance to be that the two wavefields are being saved in the same image.

[joao-bapdm]

I've made some kind of mistake, both orders give the same result now. So I wasn't adding the second source somehow..

Could someone just confirm swapping the order doesn't change anything?

[krober10nd]

Sure, but first how were you calling add_sources that created this issue?

[krober10nd]

This is a good test @jaimesouza run several simulations and then reverse the order of the sources and compare the results.

[joao-bapdm]

# Source
source = Source(kws_half_width=1, wavelet=wavelet)
source.add(position=(255.5, 255.5))
source.add(position=(270, 240))

[krober10nd]

And what's the way that didn't produce the error?

[joao-bapdm]

And what's the way that didn't produce the error?

# Source
source = Source(kws_half_width=1, wavelet=wavelet)
source.add(position=(255.5, 255.5))
source.add(position=(270, 240))

Not sure what I did before. I would like someone to run the same code just swapping the sources orders to confirm it's working, then I'll close this issue.

[krober10nd]

Yep I got a different answer just switching the source add calls

[jaimesouza]

After a careful look, I found the problem. It is a bug in the C implementation. I am gonna fix it.

[jaimesouza]

It is fixed. @joao-bapdm and @krober10nd , can you try the latest commit, please?

[krober10nd]

The example that previously did not work now appears to be working. Thanks.

[joao-bapdm]

Was getting different results in 0b1b157, now the order doesn't matter anymore.