restructuring/simplifying absorbing boundary conditions for adjoint simulations

[EtienneBachmann]

Hi all,

With the arrival of the option "UNDO_ATTENUATION_AND_OR_PML", the previous mechanism implemented for the management of the PML in reconstructing the forward wavefield became obsolete. This was based on saving at each time step the attenuation value on each GLL point of the layer, which is probably more demanding in terms of memory-I/O than the new feature. Thus, to make the code easier we should suppress this option, and set a warning saying to use UNDO_ATTENUATION_AND_OR_PML=.true. when SAVE_FORWARD and PML are true.

I would like also to mimic the 3D code in the way the compute_forces_visco* is written, ie removing the compute_forces_viscoelastic_backward that is currently used. Aside code duplication, this routine is also wrong, it lets imagine that a reverted attenuation is possible using e1, e11, and e13. (even b_e1 etc would be wrong). Since the memory variable update is based on the current and previous wavefield value, on a reverse time scheme it does not make sense to expect reverting the equations, which would require the value of the wavefield in the future to update correctly the memory variables. Thus, reverting attenuation directly is not only numerically unstable, it is also impossible because of the scheme we use.

Best regards, Etienne

[komatits]

Hi Etienne,

Very good idea. Option UNDO_ATTENUATION_AND_OR_PML is very convenient to time-revert PML, no need to keep older ways. (that is why I renamed the flag from UNDO_ATTENUATION to UNDO_ATTENUATION_AND_OR_PML a few months ago).

Thanks, Best, Dimitri.

On 05/19/2018 06:45 PM, EtienneBachmann wrote:

Hi all,

With the arrival of the option "UNDO_ATTENUATION_AND_OR_PML", the previous mechanism implemented for the management of the PML in reconstructing the forward wavefield became obsolete. This was based on saving at each time step the attenuation value on each GLL point of the layer, which is probably more demanding in terms of memory-I/O than the new feature. Thus, to make the code easier we should suppress this option, and set a warning saying to use UNDO_ATTENUATION_AND_OR_PML=.true. when SAVE_FORWARD and PML are true.

I would like also to mimic the 3D code in the way the compute_forces_visco* is written, ie removing the compute_forces_viscoelastic_backward that is currently used. Aside code duplication, this routine is also wrong, it lets imagine that a reverted attenuation is possible using e1, e11, and e13. (even b_e1 etc would be wrong). Since the memory variable update is based on the current and previous wavefield value, on a reverse time scheme it does not make sense to expect reverting the equations, which would require the value of the wavefield in the future to update correctly the memory variables. Thus, reverting attenuation directly is not only numerically unstable, it is also impossible because of the scheme we use.

Best regards, Etienne

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/geodynamics/specfem2d/issues/955, or mute the thread https://github.com/notifications/unsubscribe-auth/AFjDKduGSRVsqlOscQX_Ur2Xcd5B1uxaks5t0Ew7gaJpZM4UFxfe.

-- Dimitri Komatitsch, CNRS Research Director (DR CNRS) Laboratory of Mechanics and Acoustics, Marseille, France http://komatitsch.free.fr

[EtienneBachmann]

Hi Dimitri,

Perfect, I will work on it now. I was checking that actually the current UNDO_ATTENUATION is still using this old mechanism, PML are read and not reconstructed. As I imagined, PML files are around 100 times bigger than the ones of UNDO_ATTENUATION, and still I work on a small test case. I will implement PML the correct way, with reconstruction and not reading.

Best, Etienne

[EtienneBachmann]

Hi Dimitri,

Regarding the UNDO_ATTENUATION (for attenuation, not pml), it seems there is an error : the formerly called "displ_elastic_old" (now dux_dxl_old, duz_dzl_old & dux_dzl_plus_duz_dxl_old) is not saved in the undo_att database. It means that when you change of subset, at the first iteration of the next subset it is the unupdated displ_elastic_old that is used to update e1,e11,e13, which is incorrect. I checked, the error is also present in specfem3D globe. Since the error arrives only during one timestep, this has hopefully few impact. (and memory variables are reset at each subset so there is no cumulative effect). Dimitri, in the UNDO_ATTENUATION test you did, did you notice any small discrepancy between the forward and the reconstructed forward wavefield?

Thanks, Best,

Etienne

[EtienneBachmann]

I made some experiments, and indeed there is a (small) impact on the computed kernels because the reconstructed forward wavefield is incorrect. To get the reference kernel with attenuation, I chose NT_DUMP_ATTENUATION = NSTEP, which leads to a perfect reconstruction, since the beginning reference is 0 everywhere. Now, I reran the simulation with NT_DUMP_ATTENUATION = NSTEP/10. The image kernel_diff_current_way shows the difference between the ref kernel and the other kernel, which is non zero. This difference does not seem to be significant (around 1e-5 of the magnitude of the true kernel) but can depend on the Q chosen. Solving this bug is easy, I just saved the dux_dxl_old, duz_dzl_old & dux_dzl_plus_duz_dxl_old variables in the undo database. I reproduced the experiment, the difference is shown in "kernel_diff_new_way", and the difference is now truly 0.

I will add this to specfem2D, but this should be done in specfem3d globe too.

Best, Etienne

[komatits]

Hi Etienne,

Thanks! Very useful, because indeed after our 2016 paper on undoing attenuation was published ( http://komatitsch.free.fr/preprints/GJI_undo_attenuation_2016.pdf ) we noticed that Figure 3(e) should be exactly zero, but it is not. We kept wondering for a while what the small bug was, well I guess you have just found it :-)

Thus I confirm that it should be present in 3D_GLOBE as well, since that paper was based on it.

Thanks again, Best wishes, Dimitri.

On 05/21/2018 06:30 PM, EtienneBachmann wrote:

I made some experiments, and indeed there is a (small) impact on the computed kernels because the reconstructed forward wavefield is incorrect. To get the reference kernel with attenuation, I chose NT_DUMP_ATTENUATION = NSTEP, which leads to a perfect reconstruction, since the beginning reference is 0 everywhere. Now, I reran the simulation with NT_DUMP_ATTENUATION = NSTEP/10. The image kernel_diff_current_way shows the difference between the ref kernel and the other kernel, which is non zero. This difference does not seem to be significant (around 1e-5 of the magnitude of the true kernel) but can depend on the Q chosen. Solving this bug is easy, I just saved the dux_dxl_old, duz_dzl_old & dux_dzl_plus_duz_dxl_old variables in the undo database. I reproduced the experiment, the difference is shown in "kernel_diff_new_way", and the difference is now truly 0.

I will add this to specfem2D, but this should be done in specfem3d globe too.

Best, Etienne

kernel_rho_true https://user-images.githubusercontent.com/8615806/40318120-9b893d6a-5cf1-11e8-96f7-5d233c7c48f5.png kernel_diff_current_way https://user-images.githubusercontent.com/8615806/40317954-0ca42ee8-5cf1-11e8-92da-2c0ea0fe2b1a.png kernel_diff_new_way https://user-images.githubusercontent.com/8615806/40317955-0cb998f0-5cf1-11e8-9f78-afaddef141a5.png

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/geodynamics/specfem2d/issues/955#issuecomment-390707777, or mute the thread https://github.com/notifications/unsubscribe-auth/AFjDKSKiUwPGJJkvDF37xKDF4NyABgzEks5t0uuxgaJpZM4UFxfe.

-- Dimitri Komatitsch, CNRS Research Director (DR CNRS) Laboratory of Mechanics and Acoustics, Marseille, France http://komatitsch.free.fr

[EtienneBachmann]

Hi Dimitri,

Perfect! I am working on a commit to fix it all.

Best Etienne

[komatits]

Hi Etienne,

Perfect, thank you so much!

Best, Dimitri.

On 05/22/2018 06:27 PM, EtienneBachmann wrote:

Hi Dimitri,

Perfect! I am working on a commit to fix it all.

Best Etienne

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/geodynamics/specfem2d/issues/955#issuecomment-391054746, or mute the thread https://github.com/notifications/unsubscribe-auth/AFjDKa-rxk2gMNI2uuNQl4sl7MK8R6Evks5t1DxegaJpZM4UFxfe.

-- Dimitri Komatitsch, CNRS Research Director (DR CNRS) Laboratory of Mechanics and Acoustics, Marseille, France http://komatitsch.free.fr