search

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
398 stars 225 forks source link

Runing both PML and attenuation causes error results #356

Closed specfem3d-zhang-ksu closed 9 years ago

specfem3d-zhang-ksu commented 9 years ago

I am running the example of homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides by turning ATTENUATION false or true in the Par_file and change Qmu in the nummaterial_velocity_file.

The result with ATTENUATION=.true. is totally wrong. Is it a problem that the example is only setup for CPML or there is a bug in the code?

By examining the source code, I think that there is a bug is more probable. For example, in the subroutine compute_forces_viscoelastic_Dev_5p, there is code segment ! use first order Taylor expansion of displacement for local storage of stresses ! at this current time step, to fix attenuation in a consistent way if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = deltat_veloc(1,iglob) dummyy_loc_att(INDEX_IJK) = deltat_veloc(2,iglob) dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob) ENDDO_LOOP_IJK else if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then ! do not merge this second line with the first using an ".and." statement ! because array is_CPML() is unallocated when PML_CONDITIONS is false if (is_CPML(ispec)) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = displ_old(1,iglob) dummyy_loc_att(INDEX_IJK) = displ_old(2,iglob) dummyz_loc_att(INDEX_IJK) = displ_old(3,iglob) ENDDO_LOOP_IJK endif endif

which seams to indicate we can not have both CPML and attenuation.

Has anyone experienced this issue?

komatits commented 9 years ago

Hi,

Yes, PML in the 3D code is currently for acoustic/elastic media only (i.e. viscoelasticity cannot be turned on).

Let me cc my colleague Zhinan Xie because he is currently working on the code and thus he could probably remove that constraint (note that if so PML will not be perfectly matched strictly speaking because it is written without attenuation, thus the viscoelastic medium will be in contact with an elastic PML; but that works fine in practice if Q factors are not too small (typically not below 50 or so).

PS for Zhinan: Zhang Ksu is right, there is no reason to stop the code in that case, we should just run with an elastic PML in contact with the viscoelastic medium, as we do in the 2D code. Thus the code in his email below should be slightly modified (and the PML arrays should be allocated even if attenuation is on I guess).

Thanks, Dimitri.

On 01/20/2015 12:00 AM, specfem3d-zhang-ksu wrote:

I am running the example of homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides by turning ATTENUATION false or true in the Par_file and change Qmu in the nummaterial_velocity_file.

The result with ATTENUATION=.true. is totally wrong. Is it a problem that the example is only setup for CPML or there is a bug in the code?

By examining the source code, I think that there is a bug is more probable. For example, in the subroutine compute_forces_viscoelastic_Dev_5p, there is code segment ! use first order Taylor expansion of displacement for local storage of stresses ! at this current time step, to fix attenuation in a consistent way if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = deltat/veloc(1,iglob) dummyy_loc_att(INDEX_IJK) = deltat/veloc(2,iglob) dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob) ENDDO_LOOP_IJK else if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then ! do not merge this second line with the first using an ".and." statement ! because array is_CPML() is unallocated when PML_CONDITIONS is false if (is_CPML(ispec)) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = displ_old(1,iglob) dummyy_loc_att(INDEX_IJK) = displ_old(2,iglob) dummyz_loc_att(INDEX_IJK) = displ_old(3,iglob) ENDDO_LOOP_IJK endif endif

which seams to indicate we can not have both CPML and attenuation.

Has anyone experienced this issue?

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/356.

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

specfem3d-zhang-ksu commented 9 years ago

Thanks, Dimitri. I have replaced the code segment by the following one, 

 ! if both PML and attenuation are used, the PML region must be elastic
 if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then
    ! do not merge this second line with the first using an ".and." statement
    ! because array is_CPML() is unallocated when PML_CONDITIONS is false
    if (is_CPML(ispec)) then
       ! PML region
       ispec_CPML = spec_to_CPML(ispec)
       DO_LOOP_IJK

         dummyx_loc_att(INDEX_IJK) = PML_displ_old(1,INDEX_IJK,ispec_CPML)
         dummyy_loc_att(INDEX_IJK) = PML_displ_old(2,INDEX_IJK,ispec_CPML)
         dummyz_loc_att(INDEX_IJK) = PML_displ_old(3,INDEX_IJK,ispec_CPML)

       ENDDO_LOOP_IJK
    else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
       ! attenuation region
       DO_LOOP_IJK
          iglob = ibool(INDEX_IJK,ispec)
          dummyx_loc_att(INDEX_IJK) = deltat*veloc(1,iglob)
          dummyy_loc_att(INDEX_IJK) = deltat*veloc(2,iglob)
          dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)

       ENDDO_LOOP_IJK
    endif
 else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
    ! only attenuation
    ! use first order Taylor expansion of displacement for local storage of stresses
    ! at this current time step, to fix attenuation in a consistent way
    DO_LOOP_IJK
       iglob = ibool(INDEX_IJK,ispec)
       dummyx_loc_att(INDEX_IJK) = deltat*veloc(1,iglob)
       dummyy_loc_att(INDEX_IJK) = deltat*veloc(2,iglob)
       dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)
    ENDDO_LOOP_IJK

 endif

And similarly the other places, the results are much better, but I have a strong spike at the source location all the time, no idea what is going on.

komatits commented 9 years ago

OK, let me cc Zhinan, he is currently cleaning the 3D PML code and thus he can probably fix that as well.

Thanks for reporting the bug.

Best wishes, Dimitri.

On 01/21/2015 12:11 AM, specfem3d-zhang-ksu wrote:

Thanks, Dimitri. I have replaced the code segment by the following one,

| ! if both PML and attenuation are used, the PML region must be elastic if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then ! do not merge this second line with the first using an ".and." statement ! because array is_CPML() is unallocated when PML_CONDITIONS is false if (is_CPML(ispec)) then ! PML region ispec_CPML = spec_to_CPML(ispec) DO_LOOP_IJK

      dummyx_loc_att(INDEX_IJK) = PML_displ_old(1,INDEX_IJK,ispec_CPML)
      dummyy_loc_att(INDEX_IJK) = PML_displ_old(2,INDEX_IJK,ispec_CPML)
      dummyz_loc_att(INDEX_IJK) = PML_displ_old(3,INDEX_IJK,ispec_CPML)

    ENDDO_LOOP_IJK
 else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
    ! attenuation region
    DO_LOOP_IJK
       iglob = ibool(INDEX_IJK,ispec)
       dummyx_loc_att(INDEX_IJK) = deltat*veloc(1,iglob)
       dummyy_loc_att(INDEX_IJK) = deltat*veloc(2,iglob)
       dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob)

    ENDDO_LOOP_IJK
 endif

else if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then ! only attenuation ! use first order Taylor expansion of displacement for local storage of stresses ! at this current time step, to fix attenuation in a consistent way DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = deltat_veloc(1,iglob) dummyy_loc_att(INDEX_IJK) = deltat_veloc(2,iglob) dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob) ENDDO_LOOP_IJK

endif |

And similarly the other places, the results are much better, but I have a strong spike at the source location all the time, no idea what is going on.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/356#issuecomment-70754088.

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

komatits commented 9 years ago

Hi Dimitri, hi all,

Then I will add that into the list in specfem3d for me. Yes, you are right, that we should make that work.

Thank you so much. Best regards Zhinan

At 2015-01-21 06:58:21, "Dimitri Komatitsch" komatitsch@lma.cnrs-mrs.fr wrote:

Hi,

Yes, PML in the 3D code is currently for acoustic/elastic media only (i.e. viscoelasticity cannot be turned on).

Let me cc my colleague Zhinan Xie because he is currently working on the code and thus he could probably remove that constraint (note that if so PML will not be perfectly matched strictly speaking because it is written without attenuation, thus the viscoelastic medium will be in contact with an elastic PML; but that works fine in practice if Q factors are not too small (typically not below 50 or so).

PS for Zhinan: Zhang Ksu is right, there is no reason to stop the code in that case, we should just run with an elastic PML in contact with the viscoelastic medium, as we do in the 2D code. Thus the code in his email below should be slightly modified (and the PML arrays should be allocated even if attenuation is on I guess).

Thanks, Dimitri.

On 01/20/2015 12:00 AM, specfem3d-zhang-ksu wrote:

I am running the example of homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides by turning ATTENUATION false or true in the Par_file and change Qmu in the nummaterial_velocity_file.

The result with ATTENUATION=.true. is totally wrong. Is it a problem that the example is only setup for CPML or there is a bug in the code?

By examining the source code, I think that there is a bug is more probable. For example, in the subroutine compute_forces_viscoelastic_Dev_5p, there is code segment ! use first order Taylor expansion of displacement for local storage of stresses ! at this current time step, to fix attenuation in a consistent way if (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = deltat/veloc(1,iglob) dummyy_loc_att(INDEX_IJK) = deltat/veloc(2,iglob) dummyz_loc_att(INDEX_IJK) = deltat*veloc(3,iglob) ENDDO_LOOP_IJK else if (PML_CONDITIONS .and. (.not. backward_simulation) .and. NSPEC_CPML > 0) then ! do not merge this second line with the first using an ".and." statement ! because array is_CPML() is unallocated when PML_CONDITIONS is false if (is_CPML(ispec)) then DO_LOOP_IJK iglob = ibool(INDEX_IJK,ispec) dummyx_loc_att(INDEX_IJK) = displ_old(1,iglob) dummyy_loc_att(INDEX_IJK) = displ_old(2,iglob) dummyz_loc_att(INDEX_IJK) = displ_old(3,iglob) ENDDO_LOOP_IJK endif endif

which seams to indicate we can not have both CPML and attenuation.

Has anyone experienced this issue?

¡ª Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/356.

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

komatits commented 9 years ago

Fixed by Zhinan.