[cig-commits] r22666 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Jul 24 11:09:53 PDT 2013
Author: xie.zhinan
Date: 2013-07-24 11:09:53 -0700 (Wed, 24 Jul 2013)
New Revision: 22666
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/specfem2D.F90
Log:
add the mpi support for using PML in adjoint tomography
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-07-24 17:11:03 UTC (rev 22665)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-07-24 18:09:53 UTC (rev 22666)
@@ -61,7 +61,7 @@
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
rmemory_potential_acoust_LDDRK,alpha_LDDRK,beta_LDDRK, &
rmemory_acoustic_dux_dx_LDDRK,rmemory_acoustic_dux_dz_LDDRK,&
- deltat,PML_BOUNDARY_CONDITIONS)
+ deltat,PML_BOUNDARY_CONDITIONS,STACEY_BOUNDARY_CONDITIONS)
! compute forces for the acoustic elements
@@ -147,7 +147,7 @@
double precision :: deltat
real(kind=CUSTOM_REAL) :: A0, A1, A2, A3, A4, A5, A6, A7, A8
- logical :: PML_BOUNDARY_CONDITIONS
+ logical :: PML_BOUNDARY_CONDITIONS,STACEY_BOUNDARY_CONDITIONS
!coefficients and memory variables when using CPML with LDDRK
integer :: stage_time_scheme,i_stage
@@ -750,7 +750,7 @@
! for Stacey paraxial absorbing conditions (more precisely: Sommerfeld in the case of a fluid) we implement them here
- if(.not. PML_BOUNDARY_CONDITIONS .and. anyabs) then
+ if(STACEY_BOUNDARY_CONDITIONS) then
do ispecabs=1,nelemabs
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-07-24 17:11:03 UTC (rev 22665)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-07-24 18:09:53 UTC (rev 22666)
@@ -69,7 +69,7 @@
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation)
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation,STACEY_BOUNDARY_CONDITIONS)
! compute forces for the elastic elements
@@ -92,7 +92,8 @@
integer, dimension(nelemabs) :: ib_top
integer :: stage_time_scheme,i_stage,nadj_rec_local
- logical :: anyabs,assign_external_model,initialfield,ATTENUATION_VISCOELASTIC_SOLID,add_Bielak_conditions
+ logical :: anyabs,assign_external_model,initialfield,ATTENUATION_VISCOELASTIC_SOLID,add_Bielak_conditions,&
+ STACEY_BOUNDARY_CONDITIONS
logical :: ADD_SPRING_TO_STACEY
real(kind=CUSTOM_REAL) :: x_center_spring,z_center_spring
@@ -1390,8 +1391,7 @@
!
!--- absorbing boundaries
!
-! if(anyabs .and. .not. PML_BOUNDARY_CONDITIONS .and. backward_simulation) then
- if(anyabs .and. .not. PML_BOUNDARY_CONDITIONS) then
+ if(STACEY_BOUNDARY_CONDITIONS) then
count_left=1
count_right=1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-07-24 17:11:03 UTC (rev 22665)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-07-24 18:09:53 UTC (rev 22666)
@@ -571,7 +571,7 @@
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
output_grid_ASCII,output_grid_Gnuplot,ATTENUATION_VISCOELASTIC_SOLID,output_postscript_snapshot,output_color_image, &
plot_lowerleft_corner_only,add_Bielak_conditions,output_energy,READ_EXTERNAL_SEP_FILE, &
- output_wavefield_dumps,use_binary_for_wavefield_dumps,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE
+ output_wavefield_dumps,use_binary_for_wavefield_dumps,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,STACEY_BOUNDARY_CONDITIONS
double precision :: ROTATE_PML_ANGLE
@@ -1098,7 +1098,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'
-
!! DK DK Dec 2011: add a small crack (discontinuity) in the medium manually
npgeo_ori = npgeo
if(ADD_A_SMALL_CRACK_IN_THE_MEDIUM) npgeo = npgeo + NB_POINTS_TO_ADD_TO_NPGEO
@@ -1595,7 +1594,13 @@
nspec_left,nspec_right,nspec_bottom,nspec_top, &
ib_right,ib_left,ib_bottom,ib_top)
+ if(anyabs .and. (.not. PML_BOUNDARY_CONDITIONS))then
+ STACEY_BOUNDARY_CONDITIONS = .true.
+ else
+ STACEY_BOUNDARY_CONDITIONS = .false.
+ endif
+
if( anyabs ) then
! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
if(ipass == 1) then
@@ -5197,7 +5202,7 @@
rmemory_acoustic_dux_dx,rmemory_acoustic_dux_dz,&
rmemory_potential_acoust_LDDRK,alpha_LDDRK,beta_LDDRK, &
rmemory_acoustic_dux_dx_LDDRK,rmemory_acoustic_dux_dz_LDDRK,&
- deltat,PML_BOUNDARY_CONDITIONS)
+ deltat,PML_BOUNDARY_CONDITIONS,STACEY_BOUNDARY_CONDITIONS)
if( SIMULATION_TYPE == 3 ) then
@@ -5217,7 +5222,6 @@
if(PML_BOUNDARY_CONDITIONS)then
do i = 1, nglob_interface
- b_potential_dot_dot_acoustic(point_interface(i)) = pml_interface_history_potential_dot_dot(i,NSTEP-it+1)
b_potential_dot_acoustic(point_interface(i)) = pml_interface_history_potential_dot(i,NSTEP-it+1)
b_potential_acoustic(point_interface(i)) = pml_interface_history_potential(i,NSTEP-it+1)
enddo
@@ -5243,14 +5247,13 @@
rmemory_potential_acoust_LDDRK,alpha_LDDRK,beta_LDDRK, &
rmemory_acoustic_dux_dx_LDDRK,rmemory_acoustic_dux_dz_LDDRK,&
! deltat,PML_BOUNDARY_CONDITIONS)
- deltat,.false.)
+ deltat,.false.,STACEY_BOUNDARY_CONDITIONS)
if(PML_BOUNDARY_CONDITIONS)then
do ispec = 1,nspec
do i = 1, NGLLX
do j = 1, NGLLZ
if(.not. elastic(ispec) .and. .not. poroelastic(ispec) .and. is_pml(ispec))then
- b_potential_dot_dot_acoustic(ibool(i,j,ispec)) = 0.
b_potential_dot_acoustic(ibool(i,j,ispec)) = 0.
b_potential_acoustic(ibool(i,j,ispec)) = 0.
endif
@@ -5261,7 +5264,6 @@
if(PML_BOUNDARY_CONDITIONS)then
do i = 1, nglob_interface
- b_potential_dot_dot_acoustic(point_interface(i)) = pml_interface_history_potential_dot_dot(i,NSTEP-it+1)
b_potential_dot_acoustic(point_interface(i)) = pml_interface_history_potential_dot(i,NSTEP-it+1)
b_potential_acoustic(point_interface(i)) = pml_interface_history_potential(i,NSTEP-it+1)
enddo
@@ -5306,16 +5308,6 @@
endif
endif ! if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
- if(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)then
- if(any_acoustic .and. nglob_interface > 0)then
- do i = 1, nglob_interface
- write(72)potential_dot_dot_acoustic(point_interface(i)),&
- potential_dot_acoustic(point_interface(i)),&
- potential_acoustic(point_interface(i))
- enddo
- endif
- endif
-
endif ! end of test if any acoustic element
! *********************************************************
@@ -5631,6 +5623,26 @@
#endif
+ if(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)then
+ if(any_acoustic .and. nglob_interface > 0)then
+ do i = 1, nglob_interface
+ write(72)potential_dot_dot_acoustic(point_interface(i)),&
+ potential_dot_acoustic(point_interface(i)),&
+ potential_acoustic(point_interface(i))
+ enddo
+ endif
+ endif
+
+ if(SIMULATION_TYPE == 3)then
+ if(any_acoustic .and. nglob_interface > 0)then
+ if(PML_BOUNDARY_CONDITIONS)then
+ do i = 1, nglob_interface
+ b_potential_dot_dot_acoustic(point_interface(i)) = pml_interface_history_potential_dot_dot(i,NSTEP-it+1)
+ enddo
+ endif
+ endif
+ endif
+
! ************************************************************************************
! ************* multiply by the inverse of the mass matrix and update velocity
! ************************************************************************************
@@ -5773,7 +5785,7 @@
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.)
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.,STACEY_BOUNDARY_CONDITIONS)
if(SIMULATION_TYPE == 3)then
if(PML_BOUNDARY_CONDITIONS)then
@@ -5781,7 +5793,6 @@
do i = 1, NGLLX
do j = 1, NGLLZ
if(elastic(ispec) .and. is_pml(ispec))then
- b_accel_elastic(:,ibool(i,j,ispec)) = 0.
b_veloc_elastic(:,ibool(i,j,ispec)) = 0.
b_displ_elastic(:,ibool(i,j,ispec)) = 0.
endif
@@ -5792,9 +5803,6 @@
if(PML_BOUNDARY_CONDITIONS)then
do i = 1, nglob_interface
- b_accel_elastic(1,point_interface(i)) = pml_interface_history_accel(1,i,NSTEP-it+1)
- b_accel_elastic(2,point_interface(i)) = pml_interface_history_accel(2,i,NSTEP-it+1)
- b_accel_elastic(3,point_interface(i)) = pml_interface_history_accel(3,i,NSTEP-it+1)
b_veloc_elastic(1,point_interface(i)) = pml_interface_history_veloc(1,i,NSTEP-it+1)
b_veloc_elastic(2,point_interface(i)) = pml_interface_history_veloc(2,i,NSTEP-it+1)
b_veloc_elastic(3,point_interface(i)) = pml_interface_history_veloc(3,i,NSTEP-it+1)
@@ -5835,13 +5843,12 @@
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
!ZN PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
- .false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
+ .false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.,STACEY_BOUNDARY_CONDITIONS)
if(PML_BOUNDARY_CONDITIONS)then
do ispec = 1,nspec
do i = 1, NGLLX
do j = 1, NGLLZ
if(elastic(ispec) .and. is_pml(ispec))then
- b_accel_elastic(:,ibool(i,j,ispec)) = 0.
b_veloc_elastic(:,ibool(i,j,ispec)) = 0.
b_displ_elastic(:,ibool(i,j,ispec)) = 0.
endif
@@ -5852,9 +5859,6 @@
if(PML_BOUNDARY_CONDITIONS)then
do i = 1, nglob_interface
- b_accel_elastic(1,point_interface(i)) = pml_interface_history_accel(1,i,NSTEP-it+1)
- b_accel_elastic(2,point_interface(i)) = pml_interface_history_accel(2,i,NSTEP-it+1)
- b_accel_elastic(3,point_interface(i)) = pml_interface_history_accel(3,i,NSTEP-it+1)
b_veloc_elastic(1,point_interface(i)) = pml_interface_history_veloc(1,i,NSTEP-it+1)
b_veloc_elastic(2,point_interface(i)) = pml_interface_history_veloc(2,i,NSTEP-it+1)
b_veloc_elastic(3,point_interface(i)) = pml_interface_history_veloc(3,i,NSTEP-it+1)
@@ -5945,19 +5949,6 @@
endif ! if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)
- if(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)then
- if(any_elastic .and. nglob_interface > 0)then
- do i = 1, nglob_interface
- write(71)accel_elastic(1,point_interface(i)),accel_elastic(2,point_interface(i)),&
- accel_elastic(3,point_interface(i)),&
- veloc_elastic(1,point_interface(i)),veloc_elastic(2,point_interface(i)),&
- veloc_elastic(3,point_interface(i)),&
- displ_elastic(1,point_interface(i)),displ_elastic(2,point_interface(i)),&
- displ_elastic(3,point_interface(i))
- enddo
- endif
- endif
-
endif !if(any_elastic)
! *********************************************************
@@ -6477,7 +6468,30 @@
endif
#endif
+ if(PML_BOUNDARY_CONDITIONS .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1)then
+ if(any_elastic .and. nglob_interface > 0)then
+ do i = 1, nglob_interface
+ write(71)accel_elastic(1,point_interface(i)),accel_elastic(2,point_interface(i)),&
+ accel_elastic(3,point_interface(i)),&
+ veloc_elastic(1,point_interface(i)),veloc_elastic(2,point_interface(i)),&
+ veloc_elastic(3,point_interface(i)),&
+ displ_elastic(1,point_interface(i)),displ_elastic(2,point_interface(i)),&
+ displ_elastic(3,point_interface(i))
+ enddo
+ endif
+ endif
+ if(SIMULATION_TYPE == 3)then
+ if(PML_BOUNDARY_CONDITIONS)then
+ do i = 1, nglob_interface
+ b_accel_elastic(1,point_interface(i)) = pml_interface_history_accel(1,i,NSTEP-it+1)
+ b_accel_elastic(2,point_interface(i)) = pml_interface_history_accel(2,i,NSTEP-it+1)
+ b_accel_elastic(3,point_interface(i)) = pml_interface_history_accel(3,i,NSTEP-it+1)
+ enddo
+ endif
+ endif
+
+
! ************************************************************************************
! ************* multiply by the inverse of the mass matrix and update velocity
! ************************************************************************************
More information about the CIG-COMMITS
mailing list