Mimatch between incident and total field leading to scattering artifacts around sources and within PML

[whajjali]

Describe the bug
When trying to obtain the scattered pressure field by subtracting the total and incident fields computed using the simulate_wave_propagation function in time_varying.py using FourierSeries sound speed field for the background, there's a mismatch between both fields around the source location and within the PML which leads to scattering artifacts.

To Reproduce
To illustrate this behavior, I modified the single_source_simulation function in the FWI.ipynb notebook provided in the documentation to return the full pressure field in the domain. The incident pressure field p_inc was computed using a homogeneous sound speed field of 1480 m/s, and the total field p_tot was computed using the provided sound speed field with the 3 circle reflectors. The scattered pressure field is then computed as p_scat = p_tot - p_inc. As seen in the first animation (p_scattered.mp4) attached, there's a mismatch around the source and within the PML leading to the scattering artifacts away from the circle reflectors.

Expected behavior
The incident and total fields should match exactly until the propagating total field hits the reflectors.

Desktop (please complete the following information):

• OS: Ubuntu
• Version 22.04

Potential Fix
I think that this issue is caused by having different reference sound speed values used in the PSTD solver for the incident and total fields. Specifically, the function currently computes the reference sound speed used to compute the reference wavenumber as c_ref = functional(medium.sound_speed)(jnp.amax) in both the simulate_wave_propagation and ongrid_wave_prop_params functions. As a result, the c_ref values are different for the case of a homogeneous background and a medium with the 3 circle reflectors leading to the mismatch when going from the spectral to the spatial domain using the fft and ifft calls. A potential fix would be to use the same reference sound speed value for both computations (potential provided as input to the simulate_wave_propagation function). I tried that and was able to get rid of the artifacts as seen in the second animation attached (p_scattered_fixed_cref.mp4)

p_scattered.mp4

p_scattered_fixed_cref.mp4

[astanziola]

Hi! Really sorry for getting back to you very late on those issues and thanks for such a nice issue report!

Yes, you are right, the user should have a way to control the reference sound speed in case that the k-space operators and PML are not explicitly carried over from one simulation to the other one.
I totally agree with your solution too, do you want to draft a PR on this if you already have this solved? Otherwise I can work on in in the next days, it shouldn't take too long

[astanziola]

This has been fixed in the #221 branch.

It turns out that there were a lot of loose ends that I had to fix, to make sure tha this was not just another dirty patch.

The time simulation code now has a TimeWavePropagationSettings class that can be used to set the reference sound speed, among other things. For the latter, it allows the user to generally define a function that takes the Medium object and returns the sound speed reference value.

For your use case, one can do:

c_ref_value = 1540.0
c_ref_function = lambda m: c_ref_value 
simulation_settings = TimeWavePropagationSettings(c_ref=c_ref_function)

I know it is a bit of a convoluted way of defining a constant reference value, but in this way it can support different, dynamical ways of defining this value. For example, the default method uses

simulation_settings = TimeWavePropagationSettings(c_ref = lambda m: m.max_sound_speed)

If you could give me the code to run your example, I can transform it into a test to make sure that your use case is covered.

[whajjali]

Hi! Thanks again for the fixes above and the thorough explanation. I managed to modify my case using the example you provided above and it now works well.

[astanziola]

This should be fixed in the main :)