[cig-commits] r19169 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Nov 9 13:57:45 PST 2011
Author: dkomati1
Date: 2011-11-09 13:57:45 -0800 (Wed, 09 Nov 2011)
New Revision: 19169
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
removed an obsolete (unused) version of compute_forces_acoustic()
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2011-11-09 21:39:45 UTC (rev 19168)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2011-11-09 21:57:45 UTC (rev 19169)
@@ -45,396 +45,6 @@
subroutine compute_forces_acoustic(nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,b_potential_dot_dot_acoustic,b_potential_acoustic, &
- density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,rhoext,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll, &
- ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
- jbegin_left,jend_left,jbegin_right,jend_right,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_acoustic_left,&
- b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
- b_absorb_acoustic_top,nspec_left,nspec_right,&
- nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top)
-
-! compute forces for the acoustic elements
-
- implicit none
-
- include "constants.h"
-
- integer :: nglob,nspec,nelemabs,numat,it,NSTEP,SIMULATION_TYPE
-
- integer :: nspec_left,nspec_right,nspec_bottom,nspec_top
- integer, dimension(nelemabs) :: ib_left
- integer, dimension(nelemabs) :: ib_right
- integer, dimension(nelemabs) :: ib_bottom
- integer, dimension(nelemabs) :: ib_top
-
- logical :: anyabs,assign_external_model
- logical :: SAVE_FORWARD
-
- integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
- integer, dimension(nspec) :: kmato
- integer, dimension(nelemabs) :: numabs,ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
- jbegin_left,jend_left,jbegin_right,jend_right
-
- logical, dimension(nspec) :: elastic,poroelastic
- logical, dimension(4,nelemabs) :: codeabs
-
- real(kind=CUSTOM_REAL), dimension(nglob) :: &
- potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
- real(kind=CUSTOM_REAL), dimension(nglob) :: &
- b_potential_dot_dot_acoustic,b_potential_acoustic
- double precision, dimension(2,numat) :: density
- double precision, dimension(4,3,numat) :: poroelastcoef
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
-
- real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_left,NSTEP) :: b_absorb_acoustic_left
- real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_right,NSTEP) :: b_absorb_acoustic_right
- real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_top,NSTEP) :: b_absorb_acoustic_top
- real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_bottom,NSTEP) :: b_absorb_acoustic_bottom
-
-! derivatives of Lagrange polynomials
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
-
-! Gauss-Lobatto-Legendre weights
- real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
- real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
-
-!---
-!--- local variables
-!---
-
- integer :: ispec,i,j,k,iglob,ispecabs,ibegin,iend,jbegin,jend
-
-! spatial derivatives
- real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,dux_dxl,dux_dzl
- real(kind=CUSTOM_REAL) :: b_dux_dxi,b_dux_dgamma,b_dux_dxl,b_dux_dzl
- real(kind=CUSTOM_REAL) :: weight,xxi,zxi,xgamma,zgamma,jacobian1D
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: b_tempx1,b_tempx2
-
-! Jacobian matrix and determinant
- real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
-
-! material properties of the elastic medium
- real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,kappal,cpl,rhol
-
- integer :: ifirstelem,ilastelem
-
- ifirstelem = 1
- ilastelem = nspec
-
-! loop over spectral elements
- do ispec = ifirstelem,ilastelem
-
-!---
-!--- acoustic spectral element
-!---
- if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
-
- rhol = density(1,kmato(ispec))
-
-! first double loop over GLL points to compute and store gradients
- do j = 1,NGLLZ
- do i = 1,NGLLX
-
-! derivative along x and along z
- dux_dxi = ZERO
- dux_dgamma = ZERO
-
- if(SIMULATION_TYPE == 2) then
- b_dux_dxi = ZERO
- b_dux_dgamma = ZERO
- endif
-
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
- do k = 1,NGLLX
- dux_dxi = dux_dxi + potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
- dux_dgamma = dux_dgamma + potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
-
- if(SIMULATION_TYPE == 2) then
- b_dux_dxi = b_dux_dxi + b_potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
- b_dux_dgamma = b_dux_dgamma + b_potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
- endif
- enddo
-
- xixl = xix(i,j,ispec)
- xizl = xiz(i,j,ispec)
- gammaxl = gammax(i,j,ispec)
- gammazl = gammaz(i,j,ispec)
-
-! derivatives of potential
- dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
- dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
-
- if(SIMULATION_TYPE == 2) then
- b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
- b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
- endif
-
- jacobianl = jacobian(i,j,ispec)
-
-! 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
- tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*dux_dxl + gammazl*dux_dzl) / rhol
-
- if(SIMULATION_TYPE == 2) then
- b_tempx1(i,j) = wzgll(j)*jacobianl*(xixl*b_dux_dxl + xizl*b_dux_dzl) /rhol
- b_tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*b_dux_dxl + gammazl*b_dux_dzl) /rhol
- endif
-
- enddo
- enddo
-
-!
-! second double-loop over GLL to compute all the terms
-!
- do j = 1,NGLLZ
- do i = 1,NGLLX
-
- iglob = ibool(i,j,ispec)
-
-! along x direction and z direction
-! and assemble the contributions
- do k = 1,NGLLX
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
- (tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
-
- if(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) - &
- (b_tempx1(k,j)*hprimewgll_xx(k,i) + b_tempx2(i,k)*hprimewgll_zz(k,j))
- endif
- enddo
-
- enddo ! second loop over the GLL points
- enddo
-
- endif ! end of test if acoustic element
-
- enddo ! end of loop over all spectral elements
-
-!
-!--- absorbing boundaries
-!
- if(anyabs) then
-
- do ispecabs=1,nelemabs
-
- ispec = numabs(ispecabs)
-
-! get elastic parameters of current spectral element
- lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
- mul_relaxed = poroelastcoef(2,1,kmato(ispec))
- kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
- rhol = density(1,kmato(ispec))
-
- cpl = sqrt(kappal/rhol)
-
-!--- left absorbing boundary
- if(codeabs(ILEFT,ispecabs)) then
-
- i = 1
-
- jbegin = jbegin_left(ispecabs)
- jend = jend_left(ispecabs)
-
- do j = jbegin,jend
-
- iglob = ibool(i,j,ispec)
-
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
-
- xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
- zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
-
- weight = jacobian1D * wzgll(j)
-
-! Sommerfeld condition if acoustic
- if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - potential_dot_acoustic(iglob)*weight/cpl/rhol
-
- if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_left(j,ib_left(ispecabs),it) = &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
- elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = &
- b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
- endif
- endif
-
- enddo
-
- endif ! end of left absorbing boundary
-
-!--- right absorbing boundary
- if(codeabs(IRIGHT,ispecabs)) then
-
- i = NGLLX
-
- jbegin = jbegin_right(ispecabs)
- jend = jend_right(ispecabs)
-
- do j = jbegin,jend
-
- iglob = ibool(i,j,ispec)
-
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
-
- xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
- zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
-
- weight = jacobian1D * wzgll(j)
-
-! Sommerfeld condition if acoustic
- if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
-
-
- if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_right(j,ib_right(ispecabs),it) = &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
- elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = &
- b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
- endif
- endif
-
- enddo
-
- endif ! end of right absorbing boundary
-
-!--- bottom absorbing boundary
- if(codeabs(IBOTTOM,ispecabs)) then
-
- j = 1
-
- ibegin = ibegin_bottom(ispecabs)
- iend = iend_bottom(ispecabs)
-
-! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(ILEFT,ispecabs)) ibegin = 2
- if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
-
- do i = ibegin,iend
-
- iglob = ibool(i,j,ispec)
-
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
-
- xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
- zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xxi**2 + zxi**2)
-
- weight = jacobian1D * wxgll(i)
-
-! Sommerfeld condition if acoustic
- if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
-
- if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),it) = &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
- elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = &
- b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
- endif
- endif
-
- enddo
-
- endif ! end of bottom absorbing boundary
-
-!--- top absorbing boundary
- if(codeabs(ITOP,ispecabs)) then
-
- j = NGLLZ
-
- ibegin = ibegin_top(ispecabs)
- iend = iend_top(ispecabs)
-
-! exclude corners to make sure there is no contradiction on the normal
- if(codeabs(ILEFT,ispecabs)) ibegin = 2
- if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
-
- do i = ibegin,iend
-
- iglob = ibool(i,j,ispec)
-
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
-
- xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
- zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xxi**2 + zxi**2)
-
- weight = jacobian1D * wxgll(i)
-
-! Sommerfeld condition if acoustic
- if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
-
- if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
- b_absorb_acoustic_top(i,ib_top(ispecabs),it) = &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
- elseif(SIMULATION_TYPE == 2) then
- b_potential_dot_dot_acoustic(iglob) = &
- b_potential_dot_dot_acoustic(iglob) - &
- b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
- endif
- endif
-
- enddo
-
- endif ! end of top absorbing boundary
-
- enddo
-
- endif ! end of absorbing boundaries
-
- end subroutine compute_forces_acoustic
-
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
- subroutine compute_forces_acoustic_2(nglob,nspec,nelemabs,numat,it,NSTEP, &
- anyabs,assign_external_model,ibool,kmato,numabs, &
- elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic, &
density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
vpext,rhoext,hprime_xx,hprimewgll_xx, &
@@ -789,5 +399,5 @@
enddo
endif ! end of absorbing boundaries
- end subroutine compute_forces_acoustic_2
+ end subroutine compute_forces_acoustic
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-11-09 21:39:45 UTC (rev 19168)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-11-09 21:57:45 UTC (rev 19169)
@@ -4082,22 +4082,7 @@
! ************* compute forces for the acoustic elements
! *********************************************************
-! call compute_forces_acoustic(nglob,nspec,nelemabs,numat,it,NSTEP, &
-! anyabs,assign_external_model,ibool,kmato,numabs, &
-! elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
-! potential_acoustic,b_potential_dot_dot_acoustic,b_potential_acoustic, &
-! density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
-! vpext,rhoext,hprime_xx,hprimewgll_xx, &
-! hprime_zz,hprimewgll_zz,wxgll,wzgll, &
-! ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
-! jbegin_left,jend_left,jbegin_right,jend_right, &
-! SIMULATION_TYPE,SAVE_FORWARD,b_absorb_acoustic_left,&
-! b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
-! b_absorb_acoustic_top,nspec_left,nspec_right,&
-! nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top)
-
-
- call compute_forces_acoustic_2(nglob,nspec,nelemabs,numat,it,NSTEP, &
+ call compute_forces_acoustic(nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic, &
@@ -4111,7 +4096,7 @@
b_absorb_acoustic_left,b_absorb_acoustic_right, &
b_absorb_acoustic_bottom,b_absorb_acoustic_top,.false.)
if( SIMULATION_TYPE == 2 ) then
- call compute_forces_acoustic_2(nglob,nspec,nelemabs,numat,it,NSTEP, &
+ call compute_forces_acoustic(nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,poroelastic,codeabs,b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
b_potential_acoustic, &
More information about the CIG-COMMITS
mailing list