SPECFEM / specfem3d

SPECFEM3D_Cartesian simulates acoustic (fluid), elastic (solid), coupled acoustic/elastic, poroelastic or seismic wave propagation in any type of conforming mesh of hexahedra (structured or not).
GNU General Public License v3.0
389 stars 223 forks source link

Issue with Adjusting Time Step (dt) in ./EXAMPLES/noise_tomography #1666

Open liaohaiyang1534 opened 6 months ago

liaohaiyang1534 commented 6 months ago

Hello everyone,

I've been working with the ./EXAMPLES/noise_tomography example in SPECFEM3D and have encountered an issue that I'm hoping to get some insights on.

The example runs successfully in its original configuration. I started by making a modification to the S_squared file using the following command:

python NOISE_TOMOGRAPHY.py 999 0.12 1.0 15.0 NLNM

This change (where dt=0.12) still allows the example to run successfully. However, the issue arises when I adjust the dt value in the Par_file and S_squared. It works fine for dt=0.11, but errors occur when I set dt=0.10. The error specifically happens during step2.

image

In the output_solver, the recommended maximum time step is:

image

I'm feeling very perplexed by this issue, as it has been a stumbling block for me for several months now. Could anyone provide guidance on why this error occurs when dt is set to 0.10, and how I might resolve it? Is there a specific constraint or consideration I should be aware of when adjusting the dt value?

Thank you in advance for your help!

liaohaiyang1534 commented 6 months ago

can anyone give some advice? cause this error is just so weird to me.

danielpeter commented 6 months ago

the trace files are named according to their sampling rate, see: https://github.com/SPECFEM/specfem3d/wiki/B_channel_codes

in your case, when changing DT to 0.1, the sampling rate changes and thus the trace extensions become *.BXX.*, *.BXY.* and *.BXZ.* instead of the previous *.MXX.*, *.MXY.* and *.MXZ.*.

the run script run_this_example.sh only works with the MX* extensions and aborts if those trace files are not produced. in your case, you will have to change that script and modify all occurrences of the MX* extension to use BX* instead, as your output traces will have *.BX*.semd endings.

Jan-8-2024: editing the remark below:

(sorry, i was confusing above and below... having a smaller time step should always lead to a stable simulation. so, this is just a remark if exploring larger time steps for DT would have been explored.)

note that the maximum suggested time step is an estimate based on the mesh elements sizes. if you set DT above this suggested time step size, it might be that your simulation becomes unstable and blows up after a few time steps. thus, having a DT > 0.11 in this case might produce an unstable simulation if it's run for longer times or produce strong numerical artifacts - double-check the trace outputs if they still look ok.

liaohaiyang1534 commented 6 months ago

dear daniel

it works fine! thank you very much!

liaohaiyang1534 commented 6 months ago

Dear Daniel,

I've been repeatedly tackling an issue over the past two weeks since our last correspondence: the cross-correlation functions computed seem unusual, and I suspect it's due to numerical instability (numerical dispersion), as illustrated in the attached image.

image image

As time progresses, they appear increasingly unstable. I'm puzzled about why this is happening. I've attempted various troubleshooting steps, including:

  1. Altering the DT (time step).
  2. Changing the NSTEP.
  3. Modifying the noise cycle (for generating S_squared).
  4. Adjusting my model, such as increasing or decreasing the number and size of grids, etc.

I'm still going through the manual PDF but haven't achieved satisfactory results. Could you or someone else perhaps provide more specific guidance or direction on this issue? Any advice would be greatly appreciated.

danielpeter commented 5 months ago

maybe somebody else following this issue might have some ideas...?

from my perspective it's difficult to see what your setup is and what you expect as an outcome. for example, do you have more infos on how your S_squared source time function looks like? then, how would the traces look like computed for the different steps 1 & 2, and step 3 in case you are computing kernels. for the last step 3, how are the adjoint sources calculated and how do they look like? maybe compare these outputs with the example outputs to get an idea where things might break.

in general, changing DT, NSTEP and mesh grid sizes can be double-checked by just looking at a single forward simulation. any errors in this setup should become visible in the trace outputs.

for the cross-correlation functions, the S_squared and duration of the simulation plays an important role. double check that the signal from the main receiver has reached the other receiver positions.

then for the noise kernels, double check the adjoint sources. in the noise_tomography/ example, those are created by the fortran file adj_traveltime_filter.f90. in that file, a bunch of parameters have to be set, such as the adjoint station, the measurement window and filter bandwidth. again, it's just an example to show how one could approach this and you would have to figure this out for your specific case.

liaohaiyang1534 commented 5 months ago

maybe somebody else following this issue might have some ideas...?

from my perspective it's difficult to see what your setup is and what you expect as an outcome. for example, do you have more infos on how your S_squared source time function looks like? then, how would the traces look like computed for the different steps 1 & 2, and step 3 in case you are computing kernels. for the last step 3, how are the adjoint sources calculated and how do they look like? maybe compare these outputs with the example outputs to get an idea where things might break.

in general, changing DT, NSTEP and mesh grid sizes can be double-checked by just looking at a single forward simulation. any errors in this setup should become visible in the trace outputs.

for the cross-correlation functions, the S_squared and duration of the simulation plays an important role. double check that the signal from the main receiver has reached the other receiver positions.

then for the noise kernels, double check the adjoint sources. in the noise_tomography/ example, those are created by the fortran file adj_traveltime_filter.f90. in that file, a bunch of parameters have to be set, such as the adjoint station, the measurement window and filter bandwidth. again, it's just an example to show how one could approach this and you would have to figure this out for your specific case.

thanks for your suggestions. I'll need some time to work through the debugging process, so it might take a few days. Really appreciate your assistance!