[cig-commits] r22715 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Aug 19 00:43:55 PDT 2013
Author: xie.zhinan
Date: 2013-08-19 00:43:54 -0700 (Mon, 19 Aug 2013)
New Revision: 22715
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
fix one bug in LDDRK implementation for CFS-UPML
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-08-18 18:50:20 UTC (rev 22714)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-08-19 07:43:54 UTC (rev 22715)
@@ -424,7 +424,7 @@
PML_duz_dxl_new = 0._CUSTOM_REAL
PML_duz_dzl_new = 0._CUSTOM_REAL
if(stage_time_scheme == 6) then
- displ_elastic_new = displ_elastic + c_LDDRK(i_stage) * deltat * veloc_elastic
+ displ_elastic_new = displ_elastic !+ c_LDDRK(i_stage) * deltat * veloc_elastic
else
displ_elastic_new = displ_elastic + deltat * veloc_elastic
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-08-18 18:50:20 UTC (rev 22714)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-08-19 07:43:54 UTC (rev 22715)
@@ -57,14 +57,14 @@
ibegin_edge1,iend_edge1,ibegin_edge3,iend_edge3, &
ibegin_edge4,iend_edge4,ibegin_edge2,iend_edge2, &
rmass_inverse_elastic_three,&
- nelemabs,vsext,xix,xiz,gammaz,gammax &
+ nelemabs,vsext,xix,xiz,gammaz,gammax, &
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
- ,coord &
+ coord, &
#endif
- ,K_x_store,K_z_store,is_PML,&
- d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,region_CPML, & !which_PML_elem, &
- nspec_PML,spec_to_PML)
+ K_x_store,K_z_store,is_PML,&
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,region_CPML, &
+ nspec_PML,spec_to_PML,time_stepping_scheme)
! builds the global mass matrix
@@ -93,9 +93,8 @@
! inverse mass matrices
integer :: nglob_elastic
-
real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,&
- rmass_inverse_elastic_three
+ rmass_inverse_elastic_three
integer :: nglob_acoustic
real(kind=CUSTOM_REAL), dimension(nglob_acoustic) :: rmass_inverse_acoustic
@@ -134,8 +133,8 @@
K_x_store,K_z_store,d_x_store,d_z_store
logical, dimension(nspec) :: is_PML
logical :: PML_BOUNDARY_CONDITIONS,this_element_has_PML
-! logical, dimension(4,nspec) :: which_PML_elem
integer, dimension(nspec) :: region_CPML
+
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
double precision, dimension(NDIM,nglob_elastic), intent(in) :: coord
@@ -143,6 +142,9 @@
#endif
!! DK DK added this for Guenneau, March 2012
+!! time scheme
+ integer :: time_stepping_scheme
+
! initialize mass matrix
if(any_elastic) rmass_inverse_elastic_one(:) = 0._CUSTOM_REAL
if(any_elastic) rmass_inverse_elastic_three(:) = 0._CUSTOM_REAL
@@ -167,7 +169,6 @@
endif
if( poroelastic(ispec) ) then
-
! material is poroelastic
rhol_s = density(1,kmato(ispec))
@@ -187,85 +188,107 @@
! for elastic medium
else if( elastic(ispec) ) then
- this_element_has_PML = .false.
- if(PML_BOUNDARY_CONDITIONS .and. size(is_PML) > 1) then
+ this_element_has_PML = .false.
+ if(PML_BOUNDARY_CONDITIONS .and. size(is_PML) > 1) then
! do not merge this condition with the above line because array is_PML() sometimes has a dummy size of 1
- if (is_PML(ispec)) this_element_has_PML = .true.
- endif
+ if (is_PML(ispec)) this_element_has_PML = .true.
+ endif
- if (this_element_has_PML) then
+ if(this_element_has_PML) then
+ ispec_PML=spec_to_PML(ispec)
+ if(time_stepping_scheme == 1)then
+ if(region_CPML(ispec) == CPML_X_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ else if (region_CPML(ispec) == CPML_XY_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML)&
+ + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)+&
+ d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ else if(region_CPML(ispec) == CPML_Y_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
+ + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ endif
+ else
+ if(region_CPML(ispec) == CPML_X_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML))
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ else if (region_CPML(ispec) == CPML_XY_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML))
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ else if(region_CPML(ispec) == CPML_Y_ONLY) then
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML))
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ endif
+ endif
+ else
- ispec_PML=spec_to_PML(ispec)
- if (region_CPML(ispec) == CPML_X_ONLY) then
- rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
- + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
- rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
- else if (region_CPML(ispec) == CPML_XY_ONLY) then
- rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML)&
- + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)+&
- d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
- rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
- else if(region_CPML(ispec) == CPML_Y_ONLY) then
- rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
- + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
- rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
- endif
-
- else
-
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
include "include_mass_Guenneau.f90"
#endif
!! DK DK added this for Guenneau, March 2012
- rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
-
+ rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
#ifdef USE_GUENNEAU
endif
#endif
- rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+ endif
- endif
-
! for acoustic medium
else
- this_element_has_PML = .false.
- if(PML_BOUNDARY_CONDITIONS .and. size(is_PML) > 1) then
+ this_element_has_PML = .false.
+ if(PML_BOUNDARY_CONDITIONS .and. size(is_PML) > 1) then
! do not merge this condition with the above line because array is_PML() sometimes has a dummy size of 1
- if (is_PML(ispec)) this_element_has_PML = .true.
- endif
+ if (is_PML(ispec)) this_element_has_PML = .true.
+ endif
- if (this_element_has_PML) then
+ if(this_element_has_PML) then
- ispec_PML=spec_to_PML(ispec)
- if (region_CPML(ispec) == CPML_X_ONLY) then
- rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
- + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
- else if (region_CPML(ispec) == CPML_XY_ONLY) then
- rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML)&
- + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)&
- +d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
- else if(region_CPML(ispec) == CPML_Y_ONLY) then
- rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
- + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
- endif
-
- else
- rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
- + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
- endif
-
+ ispec_PML=spec_to_PML(ispec)
+ if(time_stepping_scheme == 1)then
+ if(region_CPML(ispec) == CPML_X_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ + d_x_store(i,j,ispec_PML) * deltat / 2.d0)
+ else if (region_CPML(ispec) == CPML_XY_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML) &
+ + (d_x_store(i,j,ispec_PML)*k_z_store(i,j,ispec_PML)&
+ + d_z_store(i,j,ispec_PML)*k_x_store(i,j,ispec_PML)) * deltat / 2.d0)
+ else if(region_CPML(ispec) == CPML_Y_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML)&
+ + d_z_store(i,j,ispec_PML)* deltat / 2.d0)
+ endif
+ else
+ if(region_CPML(ispec) == CPML_X_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML))
+ else if (region_CPML(ispec) == CPML_XY_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML) * K_z_store(i,j,ispec_PML))
+ else if(region_CPML(ispec) == CPML_Y_ONLY) then
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(i,j,ispec_PML))
+ endif
+ endif
+ else
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
+ endif
endif
-
enddo
enddo
enddo ! of do ispec = 1,nspec
@@ -278,7 +301,7 @@
!--- DK and Zhinan Xie: one per component of the wave field i.e. one per spatial dimension.
!--- DK and Zhinan Xie: This was also suggested by Jean-Paul Ampuero in 2003.
!
- if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs) then
+ if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs .and. time_stepping_scheme /= 1) then
count_left=1
count_right=1
count_bottom=1
@@ -296,11 +319,8 @@
if(codeabs(IEDGE4,ispecabs)) then
i = 1
-
do j = 1,NGLLZ
-
iglob = ibool(i,j,ispec)
-
! external velocity model
if(assign_external_model) then
cpl = vpext(i,j,ispec)
@@ -514,7 +534,7 @@
endif ! end of absorbing boundaries
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs) then
+ if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs .and. time_stepping_scheme /= 1) then
do ispecabs=1,nelemabs
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-08-18 18:50:20 UTC (rev 22714)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-08-19 07:43:54 UTC (rev 22715)
@@ -584,7 +584,7 @@
! reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
Rcoef = 0.001d0
- ALPHA_MAX_PML = PI*f0_temp ! from Festa and Vilotte
+ ALPHA_MAX_PML = 0.25d0*PI*f0_temp ! from Festa and Vilotte
! check that NPOWER is okay
if(NPOWER < 1) stop 'NPOWER must be greater than 1'
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-08-18 18:50:20 UTC (rev 22714)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-08-19 07:43:54 UTC (rev 22715)
@@ -3251,14 +3251,14 @@
ibegin_edge1,iend_edge1,ibegin_edge3,iend_edge3, &
ibegin_edge4,iend_edge4,ibegin_edge2,iend_edge2, &
rmass_inverse_elastic_three,&
- nelemabs,vsext,xix,xiz,gammaz,gammax &
+ nelemabs,vsext,xix,xiz,gammaz,gammax, &
!! DK DK added this for Guenneau, March 2012
#ifdef USE_GUENNEAU
- ,coord &
+ coord, &
#endif
- ,K_x_store,K_z_store,is_PML,&
+ K_x_store,K_z_store,is_PML,&
d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,region_CPML, &
- nspec_PML,spec_to_PML)
+ nspec_PML,spec_to_PML,time_stepping_scheme)
#ifdef USE_MPI
if ( nproc > 1 ) then
More information about the CIG-COMMITS
mailing list