[cig-commits] r21054 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Nov 19 09:12:38 PST 2012
Author: xie.zhinan
Date: 2012-11-19 09:12:38 -0800 (Mon, 19 Nov 2012)
New Revision: 21054
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
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:
use the new flag of pml
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-11-19 17:03:52 UTC (rev 21053)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2012-11-19 17:12:38 UTC (rev 21054)
@@ -55,7 +55,7 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
b_absorb_acoustic_bottom,b_absorb_acoustic_top,IS_BACKWARD_FIELD,&
- is_PML,nspec_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
rmemory_potential_acoustic,&
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -130,7 +130,7 @@
!CPML coefficients and memory variables
integer :: nspec_PML,ispec_PML
- logical, dimension(4,nspec) :: which_PML_elem
+ integer, dimension(nspec) :: region_CPML
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: spec_to_PML
@@ -236,8 +236,7 @@
if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec))then
ispec_PML=spec_to_PML(ispec)
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
@@ -298,8 +297,8 @@
dux_dzl = PML_dux_dzl(i,j,ispec_PML) + A5 * rmemory_acoustic_dux_dz(i,j,ispec_PML)
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
!------------------------------------------------------------------------------
!---------------------------- CORNER ------------------------------------------
!------------------------------------------------------------------------------
@@ -368,7 +367,7 @@
dux_dzl = PML_dux_dzl(i,j,ispec_PML) + A5 * rmemory_acoustic_dux_dz(i,j,ispec_PML)
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
!------------------------------------------------------------------------------
!---------------------------- TOP & BOTTOM ------------------------------------
@@ -437,7 +436,6 @@
! if external density model
if(assign_external_model) rhol = rhoext(i,j,ispec)
-
! for acoustic medium
! also add GLL integration weights
tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl) / rhol
@@ -463,8 +461,7 @@
ispec_PML=spec_to_PML(ispec)
iglob=ibool(i,j,ispec)
if(stage_time_scheme == 1) then
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
@@ -487,8 +484,8 @@
rmemory_potential_acoustic(2,i,j,ispec_PML)=0.0
- elseif ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
!------------------------------------------------------------------------------
!-------------------------------- CORNER --------------------------------------
@@ -528,7 +525,7 @@
+ (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob))*(it+0.5)*deltat * coef1 &
+ potential_acoustic(iglob) *(it-0.5)*deltat * coef2
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
!------------------------------------------------------------------------------
!-------------------------------- TOP & BOTTOM --------------------------------
@@ -556,9 +553,7 @@
endif
-
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
@@ -586,8 +581,8 @@
A3 * rmemory_potential_acoustic(1,i,j,ispec_PML) + &
A4 * rmemory_potential_acoustic(2,i,j,ispec_PML))
- elseif ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
!------------------------------------------------------------------------------
!-------------------------------- CORNER --------------------------------------
!------------------------------------------------------------------------------
@@ -645,7 +640,7 @@
A3 * rmemory_potential_acoustic(1,i,j,ispec_PML) + &
A4 * rmemory_potential_acoustic(2,i,j,ispec_PML))
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
!------------------------------------------------------------------------------
!-------------------------------- TOP & BOTTOM --------------------------------
!------------------------------------------------------------------------------
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-11-19 17:03:52 UTC (rev 21053)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-11-19 17:12:38 UTC (rev 21054)
@@ -63,7 +63,7 @@
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,nadj_rec_local, &
- is_PML,nspec_PML,spec_to_PML,which_PML_elem, &
+ is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
@@ -217,7 +217,7 @@
!CPML coefficients and memory variables
integer :: nspec_PML,ispec_PML
- logical, dimension(4,nspec) :: which_PML_elem
+ integer, dimension(nspec) :: region_CPML
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: spec_to_PML
logical :: PML_BOUNDARY_CONDITIONS
@@ -416,8 +416,7 @@
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
!---------------------- A8--------------------------
A8 = - d_x_store(i,j,ispec_PML) / (k_x_store(i,j,ispec_PML) ** 2)
@@ -509,8 +508,8 @@
!------------------------------------------------------------------------------
!---------------------------- CORNER ------------------------------------------
!------------------------------------------------------------------------------
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
!---------------------- A8--------------------------
A8 = (k_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
@@ -603,7 +602,7 @@
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
!------------------------------------------------------------------------------
!---------------------------- TOP & BOTTOM ------------------------------------
!------------------------------------------------------------------------------
@@ -889,8 +888,7 @@
!------------------------------------------------------------------------------
!---------------------------- LEFT & RIGHT ------------------------------------
!------------------------------------------------------------------------------
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
bb = alpha_x_store(i,j,ispec_PML)
if(stage_time_scheme == 1) then
@@ -916,8 +914,8 @@
!------------------------------------------------------------------------------
!-------------------------------- CORNER --------------------------------------
!------------------------------------------------------------------------------
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
!------------------------------------------------------------
bb = alpha_x_store(i,j,ispec_PML)
@@ -963,7 +961,7 @@
+ displ_elastic(3,iglob)*(it-0.5)*deltat * coef2
end if
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
!------------------------------------------------------------------------------
!-------------------------------- TOP & BOTTOM --------------------------------
!------------------------------------------------------------------------------
@@ -994,8 +992,7 @@
endif
- if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
A0 = - alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML)
A1 = d_x_store(i,j,ispec_PML)
@@ -1040,8 +1037,8 @@
)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
A0 = d_x_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML) &
- alpha_x_store(i,j,ispec_PML) * d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML) &
@@ -1112,7 +1109,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!corner
- else
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) then
A0 = - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)
A1 = d_z_store(i,j,ispec_PML)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-11-19 17:03:52 UTC (rev 21053)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-11-19 17:12:38 UTC (rev 21054)
@@ -60,7 +60,7 @@
,coord &
#endif
,K_x_store,K_z_store,is_PML,&
- d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem,&
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,region_CPML, & !which_PML_elem, &
nspec_PML,spec_to_PML)
! builds the global mass matrix
@@ -130,7 +130,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
+! 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
@@ -191,20 +192,23 @@
if (this_element_has_PML) then
ispec_PML=spec_to_PML(ispec)
- if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+! if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+! .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) 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)
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+! elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+! (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) 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
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) 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)
@@ -241,18 +245,21 @@
if (this_element_has_PML) then
ispec_PML=spec_to_PML(ispec)
- if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+! if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+! .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) 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)
- elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
- (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+! elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+! (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+ elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) 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
+ elseif(region_CPML(ispec) == CPML_TOP .or. region_CPML(ispec) == CPML_BOTTOM) 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)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-11-19 17:03:52 UTC (rev 21053)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-11-19 17:12:38 UTC (rev 21054)
@@ -114,7 +114,7 @@
do i_coef=2,NELEM_PML_THICKNESS
do ispec=1,nspec
- if (.not.which_PML_elem(ibound,ispec)) then
+ if (.not. which_PML_elem(ibound,ispec)) then
do j=1,NGLLZ,NGLLZ-1
do i=1,NGLLX,NGLLX-1
iglob=ibool(i,j,ispec)
@@ -156,6 +156,75 @@
enddo ! end loop on the 4 boundaries
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do ispec=1,nspec
+ if(is_PML(ispec)) then
+
+! element is in the left cpml layer
+ if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .true.) .and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .false.) ) then
+ region_CPML(ispec) = CPML_LEFT
+! element is in the right cpml layer
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .true.) .and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .false.) ) then
+ region_CPML(ispec) = CPML_RIGHT
+! element is in the top cpml layer
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .true. ) .and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .false.) ) then
+ region_CPML(ispec) = CPML_TOP
+! element is in the bottom cpml layer
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .false.) .and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .true. ) ) then
+ region_CPML(ispec) = CPML_BOTTOM
+! element is in the left-top cpml corner
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .false. )) then
+ region_CPML(ispec) = CPML_TOP_LEFT
+! element is in the right-top cpml corner
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .false. )) then
+ region_CPML(ispec) = CPML_TOP_RIGHT
+! element is in the left-bottom cpml corner
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .true. )) then
+ region_CPML(ispec) = CPML_BOTTOM_LEFT
+! element is in the right-bottom cpml corner
+ else if( &
+ (which_PML_elem(ILEFT,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(IRIGHT,ispec) .eqv. .true. ).and. &
+ (which_PML_elem(ITOP,ispec) .eqv. .false. ).and. &
+ (which_PML_elem(IBOTTOM,ispec) .eqv. .true. )) then
+ region_CPML(ispec) = CPML_BOTTOM_RIGHT
+ else
+
+ region_CPML(ispec) = 0
+
+ endif
+ endif
+ enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
!construction of table to use less memory for absorbing coefficients
spec_to_PML=0
nspec_PML=0
@@ -183,75 +252,8 @@
endif
enddo
- do ispec=1,nspec
- if(is_PML(ispec)) then
+ write(IOUT,*) "number of PML spectral elements :", nspec_PML
-! element is in the left cpml layer
- if(region_CPML(ispec) == CPML_LEFT) then
- which_PML_elem(ILEFT,ispec) = .true.
- which_PML_elem(IRIGHT,ispec) = .false.
- which_PML_elem(ITOP,ispec) = .false.
- which_PML_elem(IBOTTOM,ispec) = .false.
-
-! element is in the right cpml layer
- else if(region_CPML(ispec) == CPML_RIGHT) then
- which_PML_elem(ILEFT,ispec) = .false.
- which_PML_elem(IRIGHT,ispec) = .true.
- which_PML_elem(ITOP,ispec) = .false.
- which_PML_elem(IBOTTOM,ispec) = .false.
-
-! element is in the top cpml layer
- else if(region_CPML(ispec) == CPML_TOP) then
- which_PML_elem(ILEFT,ispec) = .false.
- which_PML_elem(IRIGHT,ispec) = .false.
- which_PML_elem(ITOP,ispec) = .true.
- which_PML_elem(IBOTTOM,ispec) = .false.
-
-! element is in the bottom cpml layer
- else if(region_CPML(ispec) == CPML_BOTTOM) then
- which_PML_elem(ILEFT,ispec) = .false.
- which_PML_elem(IRIGHT,ispec) = .false.
- which_PML_elem(ITOP,ispec) = .false.
- which_PML_elem(IBOTTOM,ispec) = .true.
-
-! element is in the left-top cpml corner
- else if(region_CPML(ispec) == CPML_TOP_LEFT) then
- which_PML_elem(ILEFT,ispec) = .true.
- which_PML_elem(IRIGHT,ispec) = .false.
- which_PML_elem(ITOP,ispec) = .true.
- which_PML_elem(IBOTTOM,ispec) = .false.
-
-! element is in the right-top cpml corner
- else if(region_CPML(ispec) == CPML_TOP_RIGHT) then
- which_PML_elem(ILEFT,ispec) = .false.
- which_PML_elem(IRIGHT,ispec) = .true.
- which_PML_elem(ITOP,ispec) = .true.
- which_PML_elem(IBOTTOM,ispec) = .false.
-
-! element is in the left-bottom cpml corner
- else if(region_CPML(ispec) == CPML_BOTTOM_LEFT) then
- which_PML_elem(ILEFT,ispec) = .true.
- which_PML_elem(IRIGHT,ispec) = .false.
- which_PML_elem(ITOP,ispec) = .false.
- which_PML_elem(IBOTTOM,ispec) = .true.
-
-! element is in the right-bottom cpml corner
- else if(region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
- which_PML_elem(ILEFT,ispec) = .false.
- which_PML_elem(IRIGHT,ispec) = .true.
- which_PML_elem(ITOP,ispec) = .false.
- which_PML_elem(IBOTTOM,ispec) = .true.
-
- else
-
- stop 'incorrect CPML flag for a PML element'
-
- endif
- endif
- enddo
-
- write(IOUT,*) "number of PML spectral elements :", nspec_PML
-
endif
end subroutine pml_init
@@ -261,7 +263,7 @@
!
subroutine define_PML_coefficients(npoin,nspec,is_PML,ibool,coord,&
- which_PML_elem,kmato,density,poroelastcoef,numat,f0_temp,&
+ region_CPML,kmato,density,poroelastcoef,numat,f0_temp,&
myrank,&
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
nspec_PML,spec_to_PML)
@@ -278,7 +280,7 @@
double precision :: f0_temp
logical, dimension(nspec) :: is_PML
- logical, dimension(4,nspec) :: which_PML_elem
+ integer, dimension(nspec) :: region_CPML
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
@@ -378,25 +380,29 @@
do i=1,NGLLX
!!!bottom_case
- if (which_PML_elem(IBOTTOM,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_RIGHT .or. region_CPML(ispec) == CPML_BOTTOM_LEFT .or. &
+ region_CPML(ispec) == CPML_BOTTOM) then
thickness_PML_z_max_bottom=max(coord(2,ibool(i,j,ispec)),thickness_PML_z_max_bottom)
thickness_PML_z_min_bottom=min(coord(2,ibool(i,j,ispec)),thickness_PML_z_min_bottom)
endif
!!!right case
- if (which_PML_elem(IRIGHT,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_RIGHT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_RIGHT) then
thickness_PML_x_max_right=max(coord(1,ibool(i,j,ispec)),thickness_PML_x_max_right)
thickness_PML_x_min_right=min(coord(1,ibool(i,j,ispec)),thickness_PML_x_min_right)
endif
!!!top case
- if (which_PML_elem(ITOP,ispec)) then
+ if (region_CPML(ispec) == CPML_TOP_RIGHT .or. region_CPML(ispec) == CPML_TOP_LEFT .or. &
+ region_CPML(ispec) == CPML_TOP) then
thickness_PML_z_max_top=max(coord(2,ibool(i,j,ispec)),thickness_PML_z_max_top)
thickness_PML_z_min_top=min(coord(2,ibool(i,j,ispec)),thickness_PML_z_min_top)
endif
!!!left case
- if (which_PML_elem(ILEFT,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_TOP_LEFT .or. &
+ region_CPML(ispec) == CPML_LEFT) then
thickness_PML_x_max_left=max(coord(1,ibool(i,j,ispec)),thickness_PML_x_max_left)
thickness_PML_x_min_left=min(coord(1,ibool(i,j,ispec)),thickness_PML_x_min_left)
endif
@@ -513,7 +519,8 @@
zval = coord(2,iglob)
!!!! ---------- bottom edge
- if (which_PML_elem(IBOTTOM,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_RIGHT .or. region_CPML(ispec) == CPML_BOTTOM_LEFT .or. &
+ region_CPML(ispec) == CPML_BOTTOM) then
! define damping profile at the grid points
abscissa_in_PML = zoriginbottom - zval
if(abscissa_in_PML >= 0.d0) then
@@ -541,7 +548,8 @@
endif
!!!! ---------- top edge
- if (which_PML_elem(ITOP,ispec)) then
+ if (region_CPML(ispec) == CPML_TOP_RIGHT .or. region_CPML(ispec) == CPML_TOP_LEFT .or. &
+ region_CPML(ispec) == CPML_TOP) then
! define damping profile at the grid points
abscissa_in_PML = zval - zorigintop
if(abscissa_in_PML >= 0.d0) then
@@ -563,7 +571,8 @@
endif
!!!! ---------- right edge
- if (which_PML_elem(IRIGHT,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_RIGHT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
+ region_CPML(ispec) == CPML_RIGHT) then
! define damping profile at the grid points
abscissa_in_PML = xval - xoriginright
if(abscissa_in_PML >= 0.d0) then
@@ -585,7 +594,8 @@
endif
!!!! ---------- left edge
- if (which_PML_elem(ILEFT,ispec)) then
+ if (region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_TOP_LEFT .or. &
+ region_CPML(ispec) == CPML_LEFT) then
! define damping profile at the grid points
abscissa_in_PML = xoriginleft - xval
if(abscissa_in_PML >= 0.d0) then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-11-19 17:03:52 UTC (rev 21053)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-11-19 17:12:38 UTC (rev 21054)
@@ -2888,7 +2888,7 @@
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
read_external_mesh,region_CPML)
- deallocate(region_CPML)
+ deallocate(which_PML_elem)
deallocate(icorner_iglob)
if (nspec_PML==0) nspec_PML=1
@@ -2912,7 +2912,7 @@
alpha_z_store(:,:,:) = 0
call define_PML_coefficients(nglob,nspec,is_PML,ibool,coord,&
- which_PML_elem,kmato,density,poroelastcoef,numat,f0(1),&
+ region_CPML,kmato,density,poroelastcoef,numat,f0(1),&
myrank,&
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
nspec_PML,spec_to_PML)
@@ -3075,7 +3075,7 @@
,coord &
#endif
,K_x_store,K_z_store,is_PML,&
- d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem,&
+ d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,region_CPML, &
nspec_PML,spec_to_PML)
#ifdef USE_MPI
@@ -4994,7 +4994,7 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
b_absorb_acoustic_bottom,b_absorb_acoustic_top,.false.,&
- is_PML,nspec_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
rmemory_potential_acoustic,&
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -5015,7 +5015,7 @@
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
b_absorb_acoustic_bottom,b_absorb_acoustic_top,.true.,&
- is_PML,nspec_PML,spec_to_PML,which_PML_elem,&
+ is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store,&
rmemory_potential_acoustic,&
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
@@ -5512,7 +5512,7 @@
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
- is_PML,nspec_PML,spec_to_PML,which_PML_elem, &
+ is_PML,nspec_PML,spec_to_PML, region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
More information about the CIG-COMMITS
mailing list