[cig-commits] r20508 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Jul 9 03:14:48 PDT 2012
Author: xie.zhinan
Date: 2012-07-09 03:14:47 -0700 (Mon, 09 Jul 2012)
New Revision: 20508
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:
add mpi support for elastic simulation with CPML
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-09 07:41:56 UTC (rev 20507)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2012-07-09 10:14:47 UTC (rev 20508)
@@ -1067,9 +1067,9 @@
iPML=ibool_PML(i,j,ispec)
iglob=ibool(i,j,ispec)
- A0 = - alpha_x_store( iPML ) * d_x_store(iPML) * k_z_store(iPML)
- A1 = d_x_store(iPML) * k_z_store(iPML)
- A2 = k_x_store(iPML) * k_z_store(iPML)
+ A0 = - alpha_x_store( iPML ) * d_x_store(iPML)
+ A1 = d_x_store(iPML)
+ A2 = k_x_store(iPML)
accel_elastic_PML(1,i,j,ispec_PML) = wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( A1 * veloc_elastic(1,iglob)+ &
@@ -1110,9 +1110,9 @@
iPML=ibool_PML(i,j,ispec)
iglob=ibool(i,j,ispec)
- A0 = - alpha_z_store( iPML ) * d_z_store(iPML) * k_x_store(iPML)
- A1 = d_z_store(iPML) * k_x_store(iPML)
- A2 = k_x_store(iPML) * k_z_store(iPML)
+ A0 = - alpha_z_store( iPML ) * d_z_store(iPML)
+ A1 = d_z_store(iPML)
+ A2 = k_z_store(iPML)
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( A1 * veloc_elastic(1,iglob)+ &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-07-09 07:41:56 UTC (rev 20507)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2012-07-09 10:14:47 UTC (rev 20508)
@@ -180,12 +180,16 @@
! for elastic medium
- if(PML_BOUNDARY_CONDITIONS)then
- if (is_PML(ispec)) then
+ if (is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS) then
iPML=ibool_PML(i,j,ispec)
+! rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+! + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
+! + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
+
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
- + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
+ + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (1.0d0 &
+ (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
+
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
else
@@ -204,25 +208,9 @@
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)
-#ifdef USE_GUENNEAU
- endif
-#endif
- rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
- endif
-
-
-
else
! for acoustic medium
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-09 07:41:56 UTC (rev 20507)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2012-07-09 10:14:47 UTC (rev 20508)
@@ -43,8 +43,10 @@
subroutine pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML,ibool_PML, &
- npoin_PML,icorner_iglob,NELEM_PML_THICKNESS)
+ npoin_PML,icorner_iglob,NELEM_PML_THICKNESS,&
+ coord,myrank)
+
implicit none
include 'constants.h'
@@ -68,6 +70,9 @@
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
integer, dimension(:), allocatable :: iPML_to_iglob
logical, dimension(nspec) :: is_PML
+ integer :: myrank,ier,ispecpml
+ character(len=256) :: prname
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: coord
!!!detection of PML elements
@@ -202,6 +207,26 @@
write(IOUT,*) "number of PML spectral elements :", nspec_PML
write(IOUT,*) "number of PML spectral points :", npoin_PML
+ write(prname,230) myrank
+ open(unit=1234,file=prname,status='unknown')
+ 230 format('./OUTPUT_FILES/is_pml',i5.5)
+
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+ ispecpml=0
+ do ispec=1,nspec
+ if(is_pml(ispec))then
+ write(1234,*)myrank,'myrank'
+ write(1234,*)is_pml(ispec),coord(1,ibool(3,3,ispec)),coord(2,ibool(3,3,ispec))
+ write(1234,*)which_PML_elem(1,ispec),which_PML_elem(2,ispec),&
+ which_PML_elem(3,ispec),which_PML_elem(4,ispec)
+ ispecpml=ispecpml+1
+ endif
+ enddo
+ write(1234,*)ispecpml,'spec number of pml in myrank_',myrank
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
+
end subroutine pml_init
!
@@ -258,9 +283,6 @@
thickness_PML_x_min_right_glob,thickness_PML_x_max_right_glob,&
thickness_PML_z_min_top_glob,thickness_PML_z_max_top_glob,&
thickness_PML_x_min_left_glob,thickness_PML_x_max_left_glob
-! thickness_PML_z_bottom_glob,thickness_PML_x_right_glob,&
-! thickness_PML_z_top_glob,thickness_PML_x_left_glob
-
double precision :: xmin_glob, xmax_glob, zmin_glob, zmax_glob, vpmax_glob
#endif
@@ -413,9 +435,9 @@
d0_z_bottom = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_bottom)
d0_z_top = - (NPOWER + 1) * vpmax * log(Rcoef) / (2.d0 * thickness_PML_z_top)
- if (myrank == 0) then
+! if (myrank == 0) then
write(IOUT,*)
- write(IOUT,*) 'PML properties -------'
+ write(IOUT,*) 'PML properties -------',myrank,'myrank'
write(IOUT,*) ' Vpmax=', vpmax
write(IOUT,*) ' log(Rcoef)=',log(Rcoef)
write(IOUT,*) ' thickness_PML_z_bottom =',thickness_PML_z_bottom
@@ -426,7 +448,7 @@
write(IOUT,*) ' d0_right =', d0_x_right
write(IOUT,*) ' d0_top =', d0_z_top
write(IOUT,*) ' d0_left =', d0_x_left
- endif
+! endif
d_x = ZERO
d_z = ZERO
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-09 07:41:56 UTC (rev 20507)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2012-07-09 10:14:47 UTC (rev 20508)
@@ -1050,11 +1050,6 @@
stop 'RK and LDDRK time scheme not supported for adjoint inversion'
if(nproc /= nproc_read_from_database) stop 'must always have nproc == nproc_read_from_database'
-#ifdef USE_MPI
- if(PML_BOUNDARY_CONDITIONS)then
- stop 'PML_BOUNDARY_CONDITIONS do not ready for mpi version of the code'
- endif
-#endif
!! DK DK Dec 2011: add a small crack (discontinuity) in the medium manually
npgeo_ori = npgeo
@@ -1320,6 +1315,12 @@
stop 'PML boundary conditions do not ready for poroelastic simulation'
endif
+#ifdef USE_MPI
+ if(PML_BOUNDARY_CONDITIONS .and. (any_acoustic .or. any_poroelastic )then
+ stop 'PML_BOUNDARY_CONDITIONS do not ready for mpi version of the code'
+ endif
+#endif
+
! allocate memory variables for attenuation
if(ipass == 1) then
allocate(e1(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -2833,7 +2834,8 @@
call pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,which_PML_poin,spec_to_PML,ibool_PML, &
- npoin_PML,icorner_iglob,NELEM_PML_THICKNESS)
+ npoin_PML,icorner_iglob,NELEM_PML_THICKNESS,&
+ coord,myrank)
deallocate(icorner_iglob)
deallocate(which_PML_poin)
More information about the CIG-COMMITS
mailing list