[cig-commits] r22200 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Jun 10 06:04:10 PDT 2013
Author: xie.zhinan
Date: 2013-06-10 06:04:10 -0700 (Mon, 10 Jun 2013)
New Revision: 22200
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
the first edition of adjoint simulation with PML
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -59,6 +59,29 @@
! returns elastic mass matrix
if( PML_CONDITIONS ) then
+ if(ACOUSTIC_SIMULATION)then
+ allocate(rmass_elastic_interface(nglob),stat=ier)
+ if(ier /= 0) stop 'error allocating array rmass'
+ rmass_elastic_interface(:) = 0._CUSTOM_REAL
+ do ispec=1,nspec
+ if( ispec_is_elastic(ispec) ) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_elastic_interface(iglob) = rmass_elastic_interface(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)) )
+ else
+ rmass_elastic_interface(iglob) = rmass_elastic_interface(iglob) + &
+ jacobianl * weight * rhostore(i,j,k,ispec)
+ endif
+ enddo;enddo;enddo
+ endif
+ enddo
+ endif
+
call create_mass_matrices_pml_elastic(nspec,ibool)
else
do ispec=1,nspec
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -324,6 +324,11 @@
if(PML_CONDITIONS)then
if( ELASTIC_SIMULATION ) then
deallocate(rmass)
+ if(PML_CONDITIONS)then
+ if(ACOUSTIC_SIMULATION)then
+ write(IOUT)rmass_elastic_interface
+ endif
+ endif
endif
if( ACOUSTIC_SIMULATION) then
deallocate(rmass_acoustic)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -208,8 +208,9 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass,rmass_acoustic,&
rmass_solid_poroelastic,rmass_fluid_poroelastic
-! mass matrix for interface !ZN
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface !ZN
+! mass matrix for interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
! mass matrix contributions
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -112,6 +112,11 @@
! elastic
if( ELASTIC_SIMULATION ) then
write(IOUT) rmass
+ if(PML_CONDITIONS)then
+ if(ACOUSTIC_SIMULATION)then
+ write(IOUT)rmass_elastic_interface
+ endif
+ endif
if( APPROXIMATE_OCEAN_LOAD) then
write(IOUT) rmass_ocean_load
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -34,7 +34,8 @@
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
- potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
+ potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
+ SIMULATION_TYPE,backward_simulation,accel_interface)
! returns the updated pressure array: potential_dot_dot_acoustic
@@ -47,7 +48,10 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic,&
potential_dot_dot_acoustic_interface
- real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: rmemory_coupling_ac_el_displ
+ integer :: SIMULATION_TYPE
+ logical :: backward_simulation
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel_interface
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: rmemory_coupling_ac_el_displ
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -95,15 +99,34 @@
! elastic displacement on global point
if(PML_CONDITIONS)then
- if(is_CPML(ispec))then
- ispec_CPML = spec_to_CPML(ispec)
- call pml_compute_memory_variables_acoustic_elastic(ispec_CPML,iface,iglob,i,j,k,&
- displ_x,displ_y,displ_z,displ,veloc,&
- num_coupling_ac_el_faces,rmemory_coupling_ac_el_displ)
+ if(.not. backward_simulation)then
+ if(is_CPML(ispec))then
+ if(SIMULATION_TYPE == 1)then
+ ispec_CPML = spec_to_CPML(ispec)
+ call pml_compute_memory_variables_acoustic_elastic(ispec_CPML,iface,iglob,i,j,k,&
+ displ_x,displ_y,displ_z,displ,veloc,&
+ num_coupling_ac_el_faces,rmemory_coupling_ac_el_displ)
+ endif
+
+ if(SIMULATION_TYPE == 3)then
+ displ_x = -accel_interface(1,iglob)
+ displ_y = -accel_interface(1,iglob)
+ displ_z = -accel_interface(1,iglob)
+ endif
+
+ else
+ displ_x = displ(1,iglob)
+ displ_y = displ(2,iglob)
+ displ_z = displ(3,iglob)
+ endif
else
- displ_x = displ(1,iglob)
- displ_y = displ(2,iglob)
- displ_z = displ(3,iglob)
+ if(is_CPML(ispec))then
+! left blank, since no operation needed
+ else
+ displ_x = displ(1,iglob)
+ displ_y = displ(2,iglob)
+ displ_z = displ(3,iglob)
+ endif
endif
else
displ_x = displ(1,iglob)
@@ -134,8 +157,13 @@
! it also means you have to calculate and update this here first before
! calculating the coupling on the elastic side for the acceleration...
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + jacobianw*displ_n
- potential_dot_dot_acoustic_interface(iglob) = potential_dot_dot_acoustic_interface(iglob) &
+
+ if(PML_CONDITIONS)then
+ if(is_CPML(ispec))then
+ potential_dot_dot_acoustic_interface(iglob) = potential_dot_dot_acoustic_interface(iglob) &
+ jacobianw*displ_n
+ endif
+ endif
enddo ! igll
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -33,19 +33,23 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
+ PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface,&
+ SIMULATION_TYPE,backward_simulation,accel_interface,&
+ rmemory_coupling_el_ac_potential,spec_to_CPML, &
+ potential_acoustic,potential_dot_acoustic)
! returns the updated acceleration array: accel
implicit none
include 'constants.h'
- integer :: NSPEC_AB,NGLOB_AB
+ integer :: NSPEC_AB,NGLOB_AB,SIMULATION_TYPE
+ logical :: backward_simulation
! displacement and pressure
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic,&
- potential_dot_dot_acoustic_interface
+ potential_acoustic,potential_dot_acoustic
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -69,8 +73,14 @@
integer :: i,j,k
! CPML
+ integer :: ispec_CPML
+ integer :: spec_to_CPML(NSPEC_AB)
logical :: PML_CONDITIONS
logical :: is_CPML(NSPEC_AB)
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel_interface
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: rmemory_coupling_el_ac_potential
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
+
! loops on all coupling faces
do iface = 1,num_coupling_ac_el_faces
@@ -93,12 +103,29 @@
iglob = ibool(i,j,k,ispec)
! acoustic pressure on global point
- if(PML_CONDITIONS)then
- if(is_CPML(ispec))then
- pressure = - potential_dot_dot_acoustic_interface(iglob)
- else
+ if(PML_CONDITIONS)then
+ if(.not. backward_simulation)then
+ if(is_CPML(ispec))then
+ if(SIMULATION_TYPE == 1)then
+ pressure = - potential_dot_dot_acoustic_interface(iglob)
+ endif
+
+ if(SIMULATION_TYPE == 3)then
+ ispec_CPML = spec_to_CPML(ispec)
+ call pml_compute_memory_variables_elastic_acoustic(ispec_CPML,iface,iglob,i,j,k,&
+ pressure,potential_acoustic,potential_dot_acoustic,&
+ num_coupling_ac_el_faces,rmemory_coupling_el_ac_potential)
+ endif
+ else
pressure = - potential_dot_dot_acoustic(iglob)
- endif
+ endif
+ else
+ if(is_CPML(ispec))then
+! left blank, since no operation needed
+ else
+ pressure = - potential_dot_dot_acoustic(iglob)
+ endif
+ endif
else
pressure = - potential_dot_dot_acoustic(iglob)
endif
@@ -126,6 +153,12 @@
accel(2,iglob) = accel(2,iglob) + jacobianw*ny*pressure
accel(3,iglob) = accel(3,iglob) + jacobianw*nz*pressure
+ if(SIMULATION_TYPE == 3)then
+ accel_interface(1,iglob) = accel_interface(1,iglob) + jacobianw*nx*pressure
+ accel_interface(2,iglob) = accel_interface(2,iglob) + jacobianw*ny*pressure
+ accel_interface(3,iglob) = accel_interface(3,iglob) + jacobianw*nz*pressure
+ endif
+
enddo ! igll
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -72,6 +72,12 @@
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces,ispec_is_acoustic)
+ if(PML_CONDITIONS)then
+ if(ELASTIC_SIMULATION ) then
+ potential_dot_dot_acoustic_interface = 0.0
+ endif
+ endif
+
! distinguishes two runs: for elements on MPI interfaces, and elements within the partitions
do iphase=1,2
@@ -130,7 +136,8 @@
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
- potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
+ potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
+ SIMULATION_TYPE,.false.,accel_interface)
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
@@ -143,7 +150,8 @@
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
- potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
+ potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
+ SIMULATION_TYPE,.false.,accel_interface)
endif
endif
endif
@@ -171,7 +179,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_src,dt,t0, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -186,7 +194,6 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
my_neighbours_ext_mesh, &
request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
-
else
! waits for send/receive requests to be completed and assembles values
@@ -195,7 +202,6 @@
max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
-
endif !phase_is_inner
enddo
@@ -401,7 +407,8 @@
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
- potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ)
+ potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
+ SIMULATION_TYPE,.true.,accel_interface)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -84,7 +84,7 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,.false.)
+ phase_ispec_inner_elastic,.false.,accel_interface,ACOUSTIC_SIMULATION)
endif
@@ -116,7 +116,12 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
+ PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface,&
+ SIMULATION_TYPE,.false.,accel_interface,&
+ rmemory_coupling_el_ac_potential,spec_to_CPML,&
+ potential_acoustic,potential_dot_acoustic)
+
+
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adoint definition: pressure^\dagger=potential^\dagger
@@ -127,7 +132,11 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
+ PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface,&
+ SIMULATION_TYPE,.false.,accel_interface,&
+ rmemory_coupling_el_ac_potential,spec_to_CPML,&
+ potential_acoustic,potential_dot_acoustic)
+
endif
endif ! num_coupling_ac_el_faces
@@ -200,6 +209,16 @@
accel(2,:) = accel(2,:)*rmassy(:)
accel(3,:) = accel(3,:)*rmassz(:)
+ if(SIMULATION_TYPE == 3)then
+ if(PML_CONDITIONS)then
+ if(ACOUSTIC_SIMULATION)then
+ accel_interface(1,:) = accel_interface(1,:)*rmass_elastic_interface(:)
+ accel_interface(2,:) = accel_interface(2,:)*rmass_elastic_interface(:)
+ accel_interface(3,:) = accel_interface(3,:)*rmass_elastic_interface(:)
+ endif
+ endif
+ endif
+
! updates acceleration with ocean load term
if(APPROXIMATE_OCEAN_LOAD) then
call compute_coupling_ocean(NSPEC_AB,NGLOB_AB, &
@@ -230,6 +249,12 @@
veloc(:,iglob) = 0.0
displ(:,iglob) = 0.0
+ if(SIMULATION_TYPE ==3)then
+ if(ACOUSTIC_SIMULATION)then
+ accel_interface(:,iglob) = 0.0
+ endif
+ endif
+
enddo
endif ! ispec_is_elastic
!!! endif
@@ -331,7 +356,7 @@
b_dsdx_top,b_dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,.true.)
+ phase_ispec_inner_elastic,.true.,accel_interface,ACOUSTIC_SIMULATION)
endif
@@ -360,7 +385,11 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface)
+ PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface,&
+ SIMULATION_TYPE,.true.,accel_interface,&
+ rmemory_coupling_el_ac_potential,spec_to_CPML,&
+ potential_acoustic,potential_dot_acoustic)
+
endif ! num_coupling_ac_el_faces
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -50,7 +50,7 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,backward_simulation)
+ phase_ispec_inner_elastic,backward_simulation,accel_interface,ACOUSTIC_SIMULATION)
use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS
use pml_par, only: NSPEC_CPML,is_CPML, spec_to_CPML, accel_elastic_CPML, &
@@ -73,6 +73,7 @@
implicit none
integer :: NSPEC_AB,NGLOB_AB
+ logical :: ACOUSTIC_SIMULATION
! displacement, velocity and acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
@@ -139,6 +140,7 @@
! C-PML absorbing boundary conditions
logical :: PML_CONDITIONS
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel_interface
! CPML adjoint
logical :: backward_simulation
@@ -803,6 +805,11 @@
! do not merge this second line with the first using an ".and." statement
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if(is_CPML(ispec)) then
+ if(SIMULATION_TYPE == 3)then
+ if(ACOUSTIC_SIMULATION)then
+ accel_interface(:,iglob) = accel(:,iglob)
+ endif
+ endif
accel(1,iglob) = accel(1,iglob) - accel_elastic_CPML(1,i,j,k)
accel(2,iglob) = accel(2,iglob) - accel_elastic_CPML(2,i,j,k)
accel(3,iglob) = accel(3,iglob) - accel_elastic_CPML(3,i,j,k)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -817,7 +817,7 @@
if (SIMULATION_TYPE == 3 .and. .NOT. GPU_MODE) then
! acoustic backward fields
if( ACOUSTIC_SIMULATION ) then
- if(PML_CONDITIONS)then !ZN
+ if(PML_CONDITIONS)then
if(nglob_interface_PML_acoustic > 0)then
call read_potential_on_pml_interface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic,&
nglob_interface_PML_acoustic,b_PML_potential,b_reclen_PML_potential)
@@ -833,7 +833,7 @@
! elastic backward fields
if( ELASTIC_SIMULATION ) then
- if(PML_CONDITIONS)then !ZN
+ if(PML_CONDITIONS)then
if(nglob_interface_PML_elastic > 0)then
call read_field_on_pml_interface(b_accel,b_veloc,b_displ,nglob_interface_PML_elastic,&
b_PML_field,b_reclen_PML_field)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -182,6 +182,13 @@
if(ier /= 0) stop 'error allocating rmemory_coupling_ac_el_displ array'
endif
+ if(SIMULATION_TYPE == 3)then
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ allocate(rmemory_coupling_el_ac_potential(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_coupling_el_ac_potential array'
+ endif
+ endif
+
spec_to_CPML(:) = 0
CPML_type(:) = 0
@@ -256,6 +263,12 @@
rmemory_coupling_ac_el_displ(:,:,:,:,:,:) = 0._CUSTOM_REAL
endif
+ if(SIMULATION_TYPE == 3)then
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ rmemory_coupling_el_ac_potential(:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
+ endif
+
! fields on PML interface will be reconstructed for adjoint simulations
! using snapshot files of wavefields
if( PML_CONDITIONS ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -426,8 +426,8 @@
! A3 = temp_A3 + (it+0.0)*deltat*temp_A4 + ((it+0.0)*deltat)**2*temp_A5
! A4 = -temp_A4 -2.0*(it+0.0)*deltat*temp_A5
! A5 = temp_A5
-!ZN the full experssion of A3,A4,A5 are given by above equation, here we use reduced
-!ZN exprssion of A3,A4,A5 in order to stabilized the code.
+!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
+!!! exprssion of A3,A4,A5 in order to stabilized the code.
A3 = temp_A3
A4 = 0.0
@@ -786,8 +786,8 @@
! A4 = -temp_A4 -2.0*it*deltat*temp_A5
! A5 = temp_A5
-!ZN the full experssion of A3,A4,A5 are given by above equation, here we use reduced
-!ZN exprssion of A3,A4,A5 in order to stabilized the code.
+!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
+!!! exprssion of A3,A4,A5 in order to stabilized the code.
A3 = temp_A3
A4 = 0.0
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -2586,3 +2586,255 @@
end subroutine pml_compute_memory_variables_acoustic_elastic
+!
+!=====================================================================
+!
+subroutine pml_compute_memory_variables_elastic_acoustic(ispec_CPML,iface,iglob,i,j,k,&
+ pressure,potential_acoustic,potential_dot_acoustic,&
+ num_coupling_ac_el_faces,rmemory_coupling_el_ac_potential)
+ ! calculates C-PML elastic memory variables and computes stress sigma
+
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+
+ use specfem_par, only: NGLOB_AB,it,deltat
+ use pml_par,only : CPML_regions,k_store_x,k_store_y,k_store_z,d_store_x,d_store_y,d_store_z,alpha_store
+ use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,&
+ CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
+
+ implicit none
+
+ integer, intent(in) :: ispec_CPML,iface,iglob,num_coupling_ac_el_faces
+ real(kind=CUSTOM_REAL) :: pressure
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: &
+ rmemory_coupling_el_ac_potential
+
+ ! local parameters
+ integer :: i,j,k
+ real(kind=CUSTOM_REAL) :: bb,coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2
+ real(kind=CUSTOM_REAL) :: A12,A13,A14
+
+
+ if( CPML_regions(ispec_CPML) == CPML_X_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- X-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML)
+ A14 = 0.d0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
+ !------------------------------------------------------------------------------
+ !---------------------------- Y-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_z
+ A12 = k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_y(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- Z-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_z
+ A12 = 1.d0
+ A13 = 0.d0
+ A14 = 0.d0
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = 0.d0
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XY-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+ A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = coef0_2 * rmemory_coupling_el_ac_potential(i,j,k,iface,2) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * it*deltat * coef1_2 &
+ + (potential_acoustic(iglob)) * it*deltat * coef2_2
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XZ-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- YZ-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_z
+ A12 = k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_y(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XYZ-corner C-PML --------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+ A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_el_ac_potential(i,j,k,iface,1) = coef0_1 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * coef1_1 &
+ + (potential_acoustic(iglob)) * coef2_1
+ rmemory_coupling_el_ac_potential(i,j,k,iface,2) = coef0_2 * rmemory_coupling_el_ac_potential(i,j,k,iface,2) &
+ + (potential_acoustic(iglob) + deltat * potential_dot_acoustic(iglob)) * it*deltat * coef1_2 &
+ + (potential_acoustic(iglob)) * it*deltat * coef2_2
+
+ pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ else
+ stop 'wrong PML flag in PML memory variable calculation routine'
+ endif
+
+end subroutine pml_compute_memory_variables_elastic_acoustic
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -108,6 +108,9 @@
! C-PML contribution to update displacement on elastic/acoustic interface
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: rmemory_coupling_ac_el_displ
+ ! C-PML contribution to update displacement on elastic/acoustic interface
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_coupling_el_ac_potential
+
! --------------------------------------------------------------------------------------------
! for adjoint tomography
! need the array stored the points on interface between PML and interior computational domain
@@ -115,7 +118,7 @@
integer :: nglob_interface_PML_acoustic,nglob_interface_PML_elastic
integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
- integer :: b_reclen_PML_field,b_reclen_PML_potential !ZN
+ integer :: b_reclen_PML_field,b_reclen_PML_potential
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_PML_field,b_PML_potential
end module pml_par
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -302,6 +302,17 @@
rmassy(:) = 1._CUSTOM_REAL / rmassy(:)
rmassz(:) = 1._CUSTOM_REAL / rmassz(:)
+ if(PML_CONDITIONS)then
+ if(ACOUSTIC_SIMULATION)then
+ call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_elastic_interface, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh)
+ where(rmass_elastic_interface <= 0._CUSTOM_REAL) rmass_elastic_interface = 1._CUSTOM_REAL
+ rmass_elastic_interface(:) = 1._CUSTOM_REAL / rmass_elastic_interface(:)
+ endif
+ endif
+
! ocean load
if(APPROXIMATE_OCEAN_LOAD ) then
call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_ocean_load, &
@@ -1292,14 +1303,14 @@
SAVE_FORWARD, &
COMPUTE_AND_STORE_STRAIN, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
-!ZN epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, &
+!!! epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, &
epsilondev_xz,epsilondev_yz, &
ATTENUATION, &
size(R_xx), &
R_xx,R_yy,R_xy,R_xz,R_yz, &
-!ZN R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, &
+!!! R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, &
one_minus_sum_beta,factor_common, &
-!ZN one_minus_sum_beta_kappa,factor_commonkappa, &
+!!! one_minus_sum_beta_kappa,factor_commonkappa, &
alphaval,betaval,gammaval, &
APPROXIMATE_OCEAN_LOAD,rmass_ocean_load, &
NOISE_TOMOGRAPHY, &
@@ -1319,12 +1330,12 @@
COMPUTE_AND_STORE_STRAIN, &
epsilon_trace_over_3, &
b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
-!ZN b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
+!!! b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, &
b_epsilondev_xz,b_epsilondev_yz, &
b_epsilon_trace_over_3, &
ATTENUATION,size(R_xx), &
b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
-!ZN b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
+!!! b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, &
b_alphaval,b_betaval,b_gammaval, &
APPROXIMATE_HESS_KL)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -139,6 +139,18 @@
! allocates mass matrix
allocate(rmass(NGLOB_AB),stat=ier)
+ if(PML_CONDITIONS)then !ZN
+ if(ACOUSTIC_SIMULATION)then !ZN
+ allocate(rmass_elastic_interface(NGLOB_AB),stat=ier) !ZN
+ if( ier /= 0 ) stop 'error allocating array rmass_elastic_interface' !ZN
+ rmass_elastic_interface(:) = 0._CUSTOM_REAL
+ if(SIMULATION_TYPE == 3)then
+ allocate(accel_interface(NDIM,NGLOB_AB),stat=ier) !ZN
+ if( ier /= 0 ) stop 'error allocating array accel_interface'
+ accel_interface(:,:) = 0._CUSTOM_REAL
+ endif
+ endif !ZN
+ endif !ZN
if( ier /= 0 ) stop 'error allocating array rmass'
! initializes mass matrix contributions
allocate(rmassx(NGLOB_AB), &
@@ -214,6 +226,12 @@
! reads mass matrices
read(27,iostat=ier) rmass
if( ier /= 0 ) stop 'error reading in array rmass'
+ if(PML_CONDITIONS)then !ZN
+ if(ACOUSTIC_SIMULATION)then !ZN
+ read(27,iostat=ier) rmass_elastic_interface !ZN
+ if( ier /= 0 ) stop 'error reading in array rmass_elastic_interface' !ZN
+ endif !ZN
+ endif !ZN
if( APPROXIMATE_OCEAN_LOAD ) then
! ocean mass matrix
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-06-09 21:04:06 UTC (rev 22199)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-06-10 13:04:10 UTC (rev 22200)
@@ -316,6 +316,10 @@
! mass matrix
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
+! PML on fluid-solid interface
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
+
! Stacey
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
More information about the CIG-COMMITS
mailing list