[cig-commits] r22123 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Tue May 21 12:43:56 PDT 2013
Author: xie.zhinan
Date: 2013-05-21 12:43:56 -0700 (Tue, 21 May 2013)
New Revision: 22123
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.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/finalize_simulation.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_par.f90
Log:
commit the first edition of adjoint simulation with PML in case of elastic simulation or acoustic simulation; optimized the memory usage when using PML by removing unneeded arraies
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-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -58,8 +58,8 @@
use specfem_par_elastic
use specfem_par_poroelastic
use pml_par,only: spec_to_CPML,is_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl,&
- rmemory_potential_acoustic,rmemory_coupling_ac_el_displ
-
+ rmemory_potential_acoustic,rmemory_coupling_ac_el_displ,nglob_interface_PML_acoustic,& !ZN
+ b_PML_potential,b_reclen_PML_potential !ZN
implicit none
! local parameters
@@ -118,8 +118,7 @@
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl,&
- rmemory_potential_acoustic,potential_dot_dot_acoustic_interface)
+ .false.,potential_dot_dot_acoustic_interface)
endif
! adjoint simulations
@@ -144,8 +143,7 @@
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl,&
- rmemory_potential_acoustic,potential_dot_dot_acoustic_interface)
+ .true.,potential_dot_dot_acoustic_interface)
endif
endif
@@ -385,7 +383,8 @@
if(PML_CONDITIONS)then
do iface=1,num_abs_boundary_faces
ispec = abs_boundary_ispec(iface)
- if(ispec_is_inner(ispec) .eqv. phase_is_inner) then
+!ZN It is better to move this into do iphase=1,2 loop
+!ZN if(ispec_is_inner(ispec) .eqv. phase_is_inner) then
if(ispec_is_acoustic(ispec) .and. is_CPML(ispec) ) then
! reference gll points on boundary face
do igll = 1,NGLLSQUARE
@@ -404,7 +403,7 @@
endif
enddo
endif ! ispec_is_acoustic
- endif
+!ZN endif
enddo
endif
@@ -460,6 +459,15 @@
call acoustic_enforce_free_surf_cuda(Mesh_pointer,STACEY_INSTEAD_OF_FREE_SURFACE)
endif
+ if(PML_CONDITIONS)then !ZN
+ if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD)then
+ if(nglob_interface_PML_acoustic > 0)then
+ call save_potential_on_pml_interface(potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic,&
+ nglob_interface_PML_acoustic,b_PML_potential,b_reclen_PML_potential)
+ endif
+ endif
+ endif
+
end subroutine compute_forces_acoustic
!
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -35,8 +35,7 @@
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- rmemory_dpotential_dxl,rmemory_dpotential_dyl,rmemory_dpotential_dzl,&
- rmemory_potential_acoustic,potential_dot_dot_acoustic_interface)
+ backward_simulation,potential_dot_dot_acoustic_interface)
! computes forces for acoustic elements
!
@@ -49,7 +48,8 @@
k_store_x,k_store_y,k_store_z,d_store_x,d_store_y,d_store_z,alpha_store,&
PML_dpotential_dxl,PML_dpotential_dyl,PML_dpotential_dzl,&
PML_dpotential_dxl_new,PML_dpotential_dyl_new,PML_dpotential_dzl_new,&
- potential_dot_dot_acoustic_CPML
+ potential_dot_dot_acoustic_CPML,rmemory_dpotential_dxl,rmemory_dpotential_dyl,&
+ rmemory_dpotential_dzl,rmemory_potential_acoustic
implicit none
@@ -82,15 +82,13 @@
integer :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
integer, dimension(num_phase_ispec_acoustic,2) :: phase_ispec_inner_acoustic
-! CPML
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2) :: &
- rmemory_dpotential_dxl, rmemory_dpotential_dyl, rmemory_dpotential_dzl
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3) :: rmemory_potential_acoustic
-
! CPML fluid-solid interface
logical :: ELASTIC_SIMULATION
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
+! CPML adjoint
+ logical :: backward_simulation
+
! local variables
real(kind=CUSTOM_REAL) :: hp1,hp2,hp3
@@ -146,7 +144,7 @@
temp3l = temp3l + chi_elem(i,j,l)*hprime_zz(k,l)
enddo
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -192,7 +190,7 @@
dpotentialdzl = xizl*temp1l + etazl*temp2l + gammazl*temp3l
! stores derivatives of ux, uy and uz with respect to x, y and z
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -227,7 +225,7 @@
enddo
enddo
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -265,7 +263,7 @@
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - ( temp1l + temp2l + temp3l )
! updates potential_dot_dot_acoustic with contribution from each C-PML element
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
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-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -90,14 +90,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, &
- rmemory_dux_dxl_x, rmemory_duy_dyl_x, rmemory_duz_dzl_x, &
- rmemory_dux_dyl_x, rmemory_dux_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
- rmemory_dux_dxl_y, rmemory_duz_dzl_y, rmemory_duy_dyl_y, &
- rmemory_duy_dxl_y, rmemory_duy_dzl_y, rmemory_duz_dyl_y, rmemory_dux_dyl_y, &
- rmemory_dux_dxl_z, rmemory_duy_dyl_z, rmemory_duz_dzl_z, &
- rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duy_dzl_z, rmemory_dux_dzl_z, &
- rmemory_displ_elastic)
+ phase_ispec_inner_elastic,.false.)
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) &
@@ -134,7 +127,7 @@
rmemory_duy_dxl_y, rmemory_duy_dzl_y, rmemory_duz_dyl_y, rmemory_dux_dyl_y, &
rmemory_dux_dxl_z, rmemory_duy_dyl_z, rmemory_duz_dzl_z, &
rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duy_dzl_z, rmemory_dux_dzl_z, &
- rmemory_displ_elastic)
+ rmemory_displ_elastic,.true.)
endif
@@ -388,7 +381,8 @@
if(PML_CONDITIONS)then
do iface=1,num_abs_boundary_faces
ispec = abs_boundary_ispec(iface)
- if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+!ZN It is better to move this into do iphase=1,2 loop
+!ZN if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
! reference gll points on boundary face
do igll = 1,NGLLSQUARE
@@ -406,7 +400,7 @@
enddo
endif ! ispec_is_elastic
- endif
+!ZN endif
enddo
endif
@@ -435,6 +429,14 @@
if( APPROXIMATE_OCEAN_LOAD ) call kernel_3_b_cuda(Mesh_pointer, NGLOB_AB, deltatover2,b_deltatover2)
endif
+ if(PML_CONDITIONS)then !ZN
+ if(SIMULATION_TYPE == 1 .and. SAVE_FORWARD)then
+ if(nglob_interface_PML_elastic > 0)then
+ call save_field_on_pml_interface(displ,veloc,accel,nglob_interface_PML_elastic,&
+ b_PML_field,b_reclen_PML_field)
+ endif
+ endif
+ endif
end subroutine compute_forces_viscoelastic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -50,14 +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, &
- rmemory_dux_dxl_x, rmemory_duy_dyl_x, rmemory_duz_dzl_x, &
- rmemory_dux_dyl_x, rmemory_dux_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
- rmemory_dux_dxl_y, rmemory_duz_dzl_y, rmemory_duy_dyl_y, &
- rmemory_duy_dxl_y, rmemory_duy_dzl_y, rmemory_duz_dyl_y, rmemory_dux_dyl_y, &
- rmemory_dux_dxl_z, rmemory_duy_dyl_z, rmemory_duz_dzl_z, &
- rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duy_dzl_z, rmemory_dux_dzl_z, &
- rmemory_displ_elastic)
+ phase_ispec_inner_elastic,backward_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, &
@@ -66,7 +59,14 @@
PML_duz_dxl, PML_duz_dyl, PML_duz_dzl, &
PML_dux_dxl_new, PML_dux_dyl_new, PML_dux_dzl_new, &
PML_duy_dxl_new, PML_duy_dyl_new, PML_duy_dzl_new, &
- PML_duz_dxl_new, PML_duz_dyl_new, PML_duz_dzl_new
+ PML_duz_dxl_new, PML_duz_dyl_new, PML_duz_dzl_new, &
+ rmemory_dux_dxl_x, rmemory_duy_dyl_x, rmemory_duz_dzl_x, &
+ rmemory_dux_dyl_x, rmemory_dux_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
+ rmemory_dux_dxl_y, rmemory_duz_dzl_y, rmemory_duy_dyl_y, &
+ rmemory_duy_dxl_y, rmemory_duy_dzl_y, rmemory_duz_dyl_y, rmemory_dux_dyl_y, &
+ rmemory_dux_dxl_z, rmemory_duy_dyl_z, rmemory_duz_dzl_z, &
+ rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duy_dzl_z, rmemory_dux_dzl_z, &
+ rmemory_displ_elastic
use fault_solver_dynamic, only : Kelvin_Voigt_eta
use specfem_par, only : FULL_ATTENUATION_SOLID
@@ -139,15 +139,10 @@
! C-PML absorbing boundary conditions
logical :: PML_CONDITIONS
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2) :: &
- rmemory_dux_dxl_x, rmemory_duy_dyl_x, rmemory_duz_dzl_x, &
- rmemory_dux_dyl_x, rmemory_dux_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
- rmemory_dux_dxl_y, rmemory_duz_dzl_y, rmemory_duy_dyl_y, &
- rmemory_duy_dxl_y, rmemory_duy_dzl_y, rmemory_duz_dyl_y, rmemory_dux_dyl_y, &
- rmemory_dux_dxl_z, rmemory_duy_dyl_z, rmemory_duz_dzl_z, &
- rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duy_dzl_z, rmemory_dux_dzl_z
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3) :: rmemory_displ_elastic
+! CPML adjoint
+ logical :: backward_simulation
+
! local parameters
integer :: i_SLS,imodulo_N_SLS
integer :: ispec,iglob,ispec_p,num_elements
@@ -264,7 +259,7 @@
enddo
enddo
enddo
- else if(PML_CONDITIONS) then
+ else if(PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -350,7 +345,7 @@
tempz3_att(i,j,k) = tempz3_att(i,j,k) + dummyz_loc_att(i,j,l)*hp3
enddo
- else if(PML_CONDITIONS) then
+ else if(PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -414,7 +409,7 @@
duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
! stores derivatives of ux, uy and uz with respect to x, y and z
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -496,7 +491,7 @@
epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
- else if(PML_CONDITIONS) then
+ else if(PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -690,7 +685,7 @@
!! DK DK that when PML_CONDITIONS is on then you do not compute the tempx, tempy, tempz arrays
!! DK DK (even in non-PML elements!!), even though such arrays are needed below;
!! DK DK shouldn't there be at least a "if (is_CPML(ispec))" test as well here, or something like that?
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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(.not.is_CPML(ispec)) then
@@ -735,7 +730,7 @@
enddo
enddo
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
@@ -804,7 +799,7 @@
fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
! updates acceleration with contribution from each C-PML element
- if (PML_CONDITIONS) then
+ if (PML_CONDITIONS .and. (.not. backward_simulation)) then
! 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
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -183,58 +183,65 @@
deallocate(alpha_store)
deallocate(spec_to_CPML)
deallocate(CPML_type)
- deallocate(PML_dux_dxl)
- deallocate(PML_dux_dyl)
- deallocate(PML_dux_dzl)
- deallocate(PML_duy_dxl)
- deallocate(PML_duy_dyl)
- deallocate(PML_duy_dzl)
- deallocate(PML_duz_dxl)
- deallocate(PML_duz_dyl)
- deallocate(PML_duz_dzl)
- deallocate(PML_dux_dxl_new)
- deallocate(PML_dux_dyl_new)
- deallocate(PML_dux_dzl_new)
- deallocate(PML_duy_dxl_new)
- deallocate(PML_duy_dyl_new)
- deallocate(PML_duy_dzl_new)
- deallocate(PML_duz_dxl_new)
- deallocate(PML_duz_dyl_new)
- deallocate(PML_duz_dzl_new)
- deallocate(PML_dpotential_dxl)
- deallocate(PML_dpotential_dyl)
- deallocate(PML_dpotential_dzl)
- deallocate(PML_dpotential_dxl_new)
- deallocate(PML_dpotential_dyl_new)
- deallocate(PML_dpotential_dzl_new)
- deallocate(rmemory_dux_dxl_x)
- deallocate(rmemory_dux_dyl_x)
- deallocate(rmemory_dux_dzl_x)
- deallocate(rmemory_duy_dxl_x)
- deallocate(rmemory_duy_dyl_x)
- deallocate(rmemory_duz_dxl_x)
- deallocate(rmemory_duz_dzl_x)
- deallocate(rmemory_dux_dxl_y)
- deallocate(rmemory_dux_dyl_y)
- deallocate(rmemory_duy_dxl_y)
- deallocate(rmemory_duy_dyl_y)
- deallocate(rmemory_duy_dzl_y)
- deallocate(rmemory_duz_dyl_y)
- deallocate(rmemory_duz_dzl_y)
- deallocate(rmemory_dux_dxl_z)
- deallocate(rmemory_dux_dzl_z)
- deallocate(rmemory_duy_dyl_z)
- deallocate(rmemory_duy_dzl_z)
- deallocate(rmemory_duz_dxl_z)
- deallocate(rmemory_duz_dyl_z)
- deallocate(rmemory_duz_dzl_z)
- deallocate(rmemory_dpotential_dxl)
- deallocate(rmemory_dpotential_dyl)
- deallocate(rmemory_dpotential_dzl)
- deallocate(rmemory_displ_elastic)
- deallocate(rmemory_potential_acoustic)
- deallocate(accel_elastic_CPML)
- deallocate(potential_dot_dot_acoustic_CPML)
+
+ if( ELASTIC_SIMULATION ) then
+ deallocate(PML_dux_dxl)
+ deallocate(PML_dux_dyl)
+ deallocate(PML_dux_dzl)
+ deallocate(PML_duy_dxl)
+ deallocate(PML_duy_dyl)
+ deallocate(PML_duy_dzl)
+ deallocate(PML_duz_dxl)
+ deallocate(PML_duz_dyl)
+ deallocate(PML_duz_dzl)
+ deallocate(PML_dux_dxl_new)
+ deallocate(PML_dux_dyl_new)
+ deallocate(PML_dux_dzl_new)
+ deallocate(PML_duy_dxl_new)
+ deallocate(PML_duy_dyl_new)
+ deallocate(PML_duy_dzl_new)
+ deallocate(PML_duz_dxl_new)
+ deallocate(PML_duz_dyl_new)
+ deallocate(PML_duz_dzl_new)
+ deallocate(rmemory_dux_dxl_x)
+ deallocate(rmemory_dux_dyl_x)
+ deallocate(rmemory_dux_dzl_x)
+ deallocate(rmemory_duy_dxl_x)
+ deallocate(rmemory_duy_dyl_x)
+ deallocate(rmemory_duz_dxl_x)
+ deallocate(rmemory_duz_dzl_x)
+ deallocate(rmemory_dux_dxl_y)
+ deallocate(rmemory_dux_dyl_y)
+ deallocate(rmemory_duy_dxl_y)
+ deallocate(rmemory_duy_dyl_y)
+ deallocate(rmemory_duy_dzl_y)
+ deallocate(rmemory_duz_dyl_y)
+ deallocate(rmemory_duz_dzl_y)
+ deallocate(rmemory_dux_dxl_z)
+ deallocate(rmemory_dux_dzl_z)
+ deallocate(rmemory_duy_dyl_z)
+ deallocate(rmemory_duy_dzl_z)
+ deallocate(rmemory_duz_dxl_z)
+ deallocate(rmemory_duz_dyl_z)
+ deallocate(rmemory_duz_dzl_z)
+ deallocate(rmemory_displ_elastic)
+ deallocate(accel_elastic_CPML)
+ endif
+
+ if( ACOUSTIC_SIMULATION ) then
+ deallocate(PML_dpotential_dxl)
+ deallocate(PML_dpotential_dyl)
+ deallocate(PML_dpotential_dzl)
+ deallocate(PML_dpotential_dxl_new)
+ deallocate(PML_dpotential_dyl_new)
+ deallocate(PML_dpotential_dzl_new)
+ deallocate(rmemory_dpotential_dxl)
+ deallocate(rmemory_dpotential_dyl)
+ deallocate(rmemory_dpotential_dzl)
+ deallocate(rmemory_potential_acoustic)
+ deallocate(potential_dot_dot_acoustic_CPML)
+ endif
+
endif
deallocate(ibelm_xmin)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -729,6 +729,7 @@
use specfem_par_acoustic
use specfem_par_elastic
use specfem_par_poroelastic
+ use pml_par
implicit none
@@ -788,6 +789,12 @@
if (SIMULATION_TYPE == 3 .and. .NOT. GPU_MODE) then
! acoustic backward fields
if( ACOUSTIC_SIMULATION ) then
+ if(PML_CONDITIONS)then !ZN
+ 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)
+ endif
+ endif
b_potential_acoustic(:) = b_potential_acoustic(:) &
+ b_deltat * b_potential_dot_acoustic(:) &
+ b_deltatsqover2 * b_potential_dot_dot_acoustic(:)
@@ -798,6 +805,12 @@
! elastic backward fields
if( ELASTIC_SIMULATION ) then
+ if(PML_CONDITIONS)then !ZN
+ 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)
+ endif
+ endif
b_displ(:,:) = b_displ(:,:) + b_deltat*b_veloc(:,:) + b_deltatsqover2*b_accel(:,:)
b_veloc(:,:) = b_veloc(:,:) + b_deltatover2*b_accel(:,:)
b_accel(:,:) = 0._CUSTOM_REAL
@@ -876,13 +889,13 @@
! memory variables if attenuation
if( ATTENUATION ) then
- if(FULL_ATTENUATION_SOLID) read(27) b_R_trace !ZN
+ if(FULL_ATTENUATION_SOLID) read(27) b_R_trace
read(27) b_R_xx
read(27) b_R_yy
read(27) b_R_xy
read(27) b_R_xz
read(27) b_R_yz
- if(FULL_ATTENUATION_SOLID) read(27) b_epsilondev_trace !ZN
+ if(FULL_ATTENUATION_SOLID) read(27) b_epsilondev_trace
read(27) b_epsilondev_xx
read(27) b_epsilondev_yy
read(27) b_epsilondev_xy
@@ -952,13 +965,13 @@
! wavefields
call transfer_b_fields_to_device(NDIM*NGLOB_AB,b_displ,b_veloc,b_accel, Mesh_pointer)
endif
- if(FULL_ATTENUATION_SOLID) read(27) b_R_trace !ZN
+ if(FULL_ATTENUATION_SOLID) read(27) b_R_trace
read(27) b_R_xx
read(27) b_R_yy
read(27) b_R_xy
read(27) b_R_xz
read(27) b_R_yz
- if(FULL_ATTENUATION_SOLID) read(27) b_epsilondev_trace !ZN
+ if(FULL_ATTENUATION_SOLID) read(27) b_epsilondev_trace
read(27) b_epsilondev_xx
read(27) b_epsilondev_yy
read(27) b_epsilondev_xy
@@ -1016,13 +1029,13 @@
size(epsilondev_xx))
endif
- if(FULL_ATTENUATION_SOLID) write(27) R_trace !ZN
+ if(FULL_ATTENUATION_SOLID) write(27) R_trace
write(27) R_xx
write(27) R_yy
write(27) R_xy
write(27) R_xz
write(27) R_yz
- if(FULL_ATTENUATION_SOLID) write(27) epsilondev_trace !ZN
+ if(FULL_ATTENUATION_SOLID) write(27) epsilondev_trace
write(27) epsilondev_xx
write(27) epsilondev_yy
write(27) epsilondev_xy
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -29,7 +29,7 @@
subroutine pml_allocate_arrays()
use pml_par
- use specfem_par, only: NSPEC_AB
+ use specfem_par, only: NSPEC_AB,PML_CONDITIONS,SIMULATION_TYPE,SAVE_FORWARD,NSTEP,myrank,prname !
use constants, only: NDIM,NGLLX,NGLLY,NGLLZ
use specfem_par_acoustic, only: ACOUSTIC_SIMULATION,num_coupling_ac_el_faces
use specfem_par_elastic, only: ELASTIC_SIMULATION
@@ -37,7 +37,8 @@
implicit none
! local parameters
- integer :: ier
+ integer :: ier,b_nglob_interface_PML_elastic,b_nglob_interface_PML_acoustic
+ integer(kind=8) :: filesize
! C-PML spectral elements local indexing
allocate(spec_to_CPML(NSPEC_AB),stat=ier)
@@ -47,128 +48,133 @@
allocate(CPML_type(NSPEC_CPML),stat=ier)
if(ier /= 0) stop 'error allocating array CPML_type'
- ! stores derivatives of ux, uy and uz with respect to x, y and z
- allocate(PML_dux_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dxl array'
- allocate(PML_dux_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dyl array'
- allocate(PML_dux_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dzl array'
- allocate(PML_duy_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dxl array'
- allocate(PML_duy_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dyl array'
- allocate(PML_duy_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dzl array'
- allocate(PML_duz_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dxl array'
- allocate(PML_duz_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dyl array'
- allocate(PML_duz_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dzl array'
+ if( ELASTIC_SIMULATION) then
+ ! stores derivatives of ux, uy and uz with respect to x, y and z
+ allocate(PML_dux_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dxl array'
+ allocate(PML_dux_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dyl array'
+ allocate(PML_dux_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dzl array'
+ allocate(PML_duy_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dxl array'
+ allocate(PML_duy_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dyl array'
+ allocate(PML_duy_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dzl array'
+ allocate(PML_duz_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dxl array'
+ allocate(PML_duz_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dyl array'
+ allocate(PML_duz_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dzl array'
- allocate(PML_dux_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dxl_new array'
- allocate(PML_dux_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dyl_new array'
- allocate(PML_dux_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dux_dzl_new array'
- allocate(PML_duy_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dxl_new array'
- allocate(PML_duy_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dyl_new array'
- allocate(PML_duy_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duy_dzl_new array'
- allocate(PML_duz_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dxl_new array'
- allocate(PML_duz_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dyl_new array'
- allocate(PML_duz_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_duz_dzl_new array'
+ allocate(PML_dux_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dxl_new array'
+ allocate(PML_dux_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dyl_new array'
+ allocate(PML_dux_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dux_dzl_new array'
+ allocate(PML_duy_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dxl_new array'
+ allocate(PML_duy_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dyl_new array'
+ allocate(PML_duy_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duy_dzl_new array'
+ allocate(PML_duz_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dxl_new array'
+ allocate(PML_duz_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dyl_new array'
+ allocate(PML_duz_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_duz_dzl_new array'
- ! stores derivatives of potential with respect to x, y and z
- allocate(PML_dpotential_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
- allocate(PML_dpotential_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
- allocate(PML_dpotential_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
+ ! stores C-PML memory variables
+ allocate(rmemory_dux_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dxl_x array'
+ allocate(rmemory_dux_dyl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dyl_x array'
+ allocate(rmemory_dux_dzl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dzl_x array'
+ allocate(rmemory_duy_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dxl_x array'
+ allocate(rmemory_duy_dyl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dyl_x array'
+ allocate(rmemory_duz_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dxl_x array'
+ allocate(rmemory_duz_dzl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dzl_x array'
- allocate(PML_dpotential_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
- allocate(PML_dpotential_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
- allocate(PML_dpotential_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
+ allocate(rmemory_dux_dxl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dxl_y array'
+ allocate(rmemory_dux_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dyl_y array'
+ allocate(rmemory_duy_dxl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dxl_y array'
+ allocate(rmemory_duy_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dyl_y array'
+ allocate(rmemory_duy_dzl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dzl_y array'
+ allocate(rmemory_duz_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dyl_y array'
+ allocate(rmemory_duz_dzl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dzl_y array'
- ! stores C-PML memory variables
- allocate(rmemory_dux_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dxl_x array'
- allocate(rmemory_dux_dyl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dyl_x array'
- allocate(rmemory_dux_dzl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dzl_x array'
- allocate(rmemory_duy_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dxl_x array'
- allocate(rmemory_duy_dyl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dyl_x array'
- allocate(rmemory_duz_dxl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dxl_x array'
- allocate(rmemory_duz_dzl_x(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dzl_x array'
+ allocate(rmemory_dux_dxl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dxl_z array'
+ allocate(rmemory_dux_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dux_dzl_z array'
+ allocate(rmemory_duy_dyl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dyl_z array'
+ allocate(rmemory_duy_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duy_dzl_z array'
+ allocate(rmemory_duz_dxl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dxl_z array'
+ allocate(rmemory_duz_dyl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dyl_z array'
+ allocate(rmemory_duz_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_duz_dzl_z array'
- allocate(rmemory_dux_dxl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dxl_y array'
- allocate(rmemory_dux_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dyl_y array'
- allocate(rmemory_duy_dxl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dxl_y array'
- allocate(rmemory_duy_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dyl_y array'
- allocate(rmemory_duy_dzl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dzl_y array'
- allocate(rmemory_duz_dyl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dyl_y array'
- allocate(rmemory_duz_dzl_y(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dzl_y array'
+ ! stores C-PML memory variables needed for displacement
+ allocate(rmemory_displ_elastic(NDIM,NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_displ_elastic array'
- allocate(rmemory_dux_dxl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dxl_z array'
- allocate(rmemory_dux_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dux_dzl_z array'
- allocate(rmemory_duy_dyl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dyl_z array'
- allocate(rmemory_duy_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duy_dzl_z array'
- allocate(rmemory_duz_dxl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dxl_z array'
- allocate(rmemory_duz_dyl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dyl_z array'
- allocate(rmemory_duz_dzl_z(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_duz_dzl_z array'
+ ! stores C-PML contribution to update acceleration to the global mesh
+ allocate(accel_elastic_CPML(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating accel_elastic_CPML array'
+ endif
- allocate(rmemory_dpotential_dxl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dpotential_dxl array'
- allocate(rmemory_dpotential_dyl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dpotential_dyl array'
- allocate(rmemory_dpotential_dzl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_dpotential_dzl array'
+ if( ACOUSTIC_SIMULATION) then
+ ! stores derivatives of potential with respect to x, y and z
+ allocate(PML_dpotential_dxl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
+ allocate(PML_dpotential_dyl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
+ allocate(PML_dpotential_dzl(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl array'
- ! stores C-PML memory variables needed for displacement
- allocate(rmemory_displ_elastic(NDIM,NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_displ_elastic array'
+ allocate(PML_dpotential_dxl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
+ allocate(PML_dpotential_dyl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
+ allocate(PML_dpotential_dzl_new(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating PML_dpotential_dxl_new array'
- ! stores C-PML memory variables needed for potential
- allocate(rmemory_potential_acoustic(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_potential_acoustic array'
+ ! stores C-PML memory variables
+ allocate(rmemory_dpotential_dxl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dpotential_dxl array'
+ allocate(rmemory_dpotential_dyl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dpotential_dyl array'
+ allocate(rmemory_dpotential_dzl(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_dpotential_dzl array'
- ! stores C-PML contribution to update acceleration to the global mesh
- allocate(accel_elastic_CPML(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating accel_elastic_CPML array'
+ ! stores C-PML memory variables needed for potential
+ allocate(rmemory_potential_acoustic(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_potential_acoustic array'
- ! stores C-PML contribution to update the second derivative of the potential to the global mesh
- allocate(potential_dot_dot_acoustic_CPML(NGLLX,NGLLY,NGLLZ),stat=ier)
- if(ier /= 0) stop 'error allocating potential_dot_dot_acoustic_CPML array'
+ ! stores C-PML contribution to update the second derivative of the potential to the global mesh
+ allocate(potential_dot_dot_acoustic_CPML(NGLLX,NGLLY,NGLLZ),stat=ier)
+ if(ier /= 0) stop 'error allocating potential_dot_dot_acoustic_CPML array'
+ endif
! stores C-PML contribution on elastic/acoustic interface
if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
@@ -180,72 +186,188 @@
CPML_type(:) = 0
- PML_dux_dxl(:,:,:) = 0._CUSTOM_REAL
- PML_dux_dyl(:,:,:) = 0._CUSTOM_REAL
- PML_dux_dzl(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dxl(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dyl(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dzl(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dxl(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dyl(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dzl(:,:,:) = 0._CUSTOM_REAL
+ if( ELASTIC_SIMULATION) then
+ PML_dux_dxl(:,:,:) = 0._CUSTOM_REAL
+ PML_dux_dyl(:,:,:) = 0._CUSTOM_REAL
+ PML_dux_dzl(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dxl(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dyl(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dzl(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dxl(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dyl(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dzl(:,:,:) = 0._CUSTOM_REAL
- PML_dux_dxl_new(:,:,:) = 0._CUSTOM_REAL
- PML_dux_dyl_new(:,:,:) = 0._CUSTOM_REAL
- PML_dux_dzl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dxl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dyl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duy_dzl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dxl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dyl_new(:,:,:) = 0._CUSTOM_REAL
- PML_duz_dzl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_dux_dxl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_dux_dyl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_dux_dzl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dxl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dyl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duy_dzl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dxl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dyl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_duz_dzl_new(:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dxl(:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dyl(:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dzl(:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dyl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dzl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dyl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dzl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dxl_new(:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dyl_new(:,:,:) = 0._CUSTOM_REAL
- PML_dpotential_dzl_new(:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dxl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dxl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dzl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dzl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dyl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dzl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dyl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dxl_x(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dzl_x(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dxl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dux_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dyl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duy_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dxl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dyl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_duz_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dxl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dxl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dzl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dyl_y(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dzl_y(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_displ_elastic(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ accel_elastic_CPML(:,:,:,:) = 0._CUSTOM_REAL
+ endif
- rmemory_dux_dxl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dux_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dyl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duy_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dxl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dyl_z(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_duz_dzl_z(:,:,:,:,:) = 0._CUSTOM_REAL
+ if( ACOUSTIC_SIMULATION) then
+ PML_dpotential_dxl(:,:,:) = 0._CUSTOM_REAL
+ PML_dpotential_dyl(:,:,:) = 0._CUSTOM_REAL
+ PML_dpotential_dzl(:,:,:) = 0._CUSTOM_REAL
- rmemory_dpotential_dxl(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dpotential_dyl(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_dpotential_dzl(:,:,:,:,:) = 0._CUSTOM_REAL
+ PML_dpotential_dxl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_dpotential_dyl_new(:,:,:) = 0._CUSTOM_REAL
+ PML_dpotential_dzl_new(:,:,:) = 0._CUSTOM_REAL
- rmemory_displ_elastic(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dpotential_dxl(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dpotential_dyl(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_dpotential_dzl(:,:,:,:,:) = 0._CUSTOM_REAL
- rmemory_potential_acoustic(:,:,:,:,:) = 0._CUSTOM_REAL
+ rmemory_potential_acoustic(:,:,:,:,:) = 0._CUSTOM_REAL
+ potential_dot_dot_acoustic_CPML(:,:,:) = 0._CUSTOM_REAL
+ endif
- accel_elastic_CPML(:,:,:,:) = 0._CUSTOM_REAL
-
- potential_dot_dot_acoustic_CPML(:,:,:) = 0._CUSTOM_REAL
-
if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
rmemory_coupling_ac_el_displ(:,:,:,:,:,:) = 0._CUSTOM_REAL
endif
+! fields on PML interface will be reconstructed for adjoint simulations
+! using snapshot files of wavefields
+ if( PML_CONDITIONS ) then
+ ! opens absorbing wavefield saved/to-be-saved by forward simulations
+ if( nglob_interface_PML_elastic > 0 .and. (SIMULATION_TYPE == 3 .or. &
+ (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) ) then
+
+ b_nglob_interface_PML_elastic = nglob_interface_PML_elastic
+
+ ! elastic domains
+ if( ELASTIC_SIMULATION) then
+ ! allocates wavefield
+ allocate(b_PML_field(9,b_nglob_interface_PML_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_field'
+
+ ! size of single record
+ b_reclen_PML_field = CUSTOM_REAL * 9 * nglob_interface_PML_elastic
+
+ ! check integer size limit: size of b_reclen_PML_field must fit onto an 4-byte integer
+ if( nglob_interface_PML_elastic > 2147483646 / (CUSTOM_REAL * 9) ) then
+ print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_PML_field
+ print *,' ',CUSTOM_REAL, NDIM, 9, nglob_interface_PML_elastic
+ print*,'bit size fortran: ',bit_size(b_reclen_PML_field)
+ call exit_MPI(myrank,"error b_reclen_PML_field integer limit")
+ endif
+
+ ! total file size
+ filesize = b_reclen_PML_field
+ filesize = filesize*NSTEP
+
+ if (SIMULATION_TYPE == 3) then
+ call open_file_abs_r(0,trim(prname)//'absorb_PML_field.bin', &
+ len_trim(trim(prname)//'absorb_PML_field.bin'), &
+ filesize)
+
+ else
+ call open_file_abs_w(0,trim(prname)//'absorb_PML_field.bin', &
+ len_trim(trim(prname)//'absorb_PML_field.bin'), &
+ filesize)
+ endif
+ endif
+ else
+ ! needs dummy array
+ b_nglob_interface_PML_elastic = 1
+ if( ELASTIC_SIMULATION ) then
+ allocate(b_PML_field(9,b_nglob_interface_PML_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_field'
+ endif
+ endif
+ else ! PML_CONDITIONS
+ b_nglob_interface_PML_elastic = 1
+ if( ELASTIC_SIMULATION ) then
+ allocate(b_PML_field(9,b_nglob_interface_PML_elastic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_field'
+ endif
+ endif
+
+! acoustic domain
+
+ if( PML_CONDITIONS ) then
+ ! opens absorbing wavefield saved/to-be-saved by forward simulations
+ if( nglob_interface_PML_acoustic > 0 .and. (SIMULATION_TYPE == 3 .or. &
+ (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) ) then
+
+ b_nglob_interface_PML_acoustic = nglob_interface_PML_acoustic
+
+ ! elastic domains
+ if( ACOUSTIC_SIMULATION) then
+ ! allocates wavefield
+ allocate(b_PML_potential(3,b_nglob_interface_PML_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_potential'
+
+ ! size of single record
+ b_reclen_PML_potential = CUSTOM_REAL * nglob_interface_PML_acoustic
+
+ ! check integer size limit: size of b_reclen_PML_field must fit onto an 4-byte integer
+ if( nglob_interface_PML_acoustic > 2147483646 / (CUSTOM_REAL) ) then
+ print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_PML_potential
+ print *,' ',CUSTOM_REAL, nglob_interface_PML_acoustic
+ print*,'bit size fortran: ',bit_size(b_reclen_PML_potential)
+ call exit_MPI(myrank,"error b_reclen_PML_potential integer limit")
+ endif
+
+ ! total file size
+ filesize = b_reclen_PML_potential
+ filesize = filesize*NSTEP
+
+ if (SIMULATION_TYPE == 3) then
+ call open_file_abs_r(1,trim(prname)//'absorb_PML_potential.bin', &
+ len_trim(trim(prname)//'absorb_PML_potential.bin'), &
+ filesize)
+
+ else
+ call open_file_abs_w(1,trim(prname)//'absorb_PML_potential.bin', &
+ len_trim(trim(prname)//'absorb_PML_potential.bin'), &
+ filesize)
+ endif
+ endif
+ else
+ ! needs dummy array
+ b_nglob_interface_PML_acoustic = 1
+ if( ACOUSTIC_SIMULATION ) then
+ allocate(b_PML_potential(3,b_nglob_interface_PML_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_potential'
+ endif
+ endif
+ else ! PML_CONDITIONS
+ b_nglob_interface_PML_acoustic = 1
+ if( ACOUSTIC_SIMULATION ) then
+ allocate(b_PML_potential(3,b_nglob_interface_PML_acoustic),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_PML_potential'
+ endif
+ endif
+
end subroutine pml_allocate_arrays
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -806,3 +806,129 @@
enddo
end subroutine pml_compute_accel_contribution_acoustic
+!
+!=====================================================================
+!
+subroutine save_field_on_pml_interface(displ,veloc,accel,nglob_interface_PML_elastic,&
+ b_PML_field,b_reclen_PML_field)
+ use specfem_par, only: NGLOB_AB,it
+ use constants, only: CUSTOM_REAL,NDIM
+ integer, intent(in) :: nglob_interface_PML_elastic,b_reclen_PML_field
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB), intent(in) :: displ,veloc,accel
+ real(kind=CUSTOM_REAL), dimension(9,nglob_interface_PML_elastic) :: b_PML_field
+
+ integer :: iglob
+
+ do iglob = 1, nglob_interface_PML_elastic
+ b_PML_field(1,iglob) = displ(1,iglob)
+ b_PML_field(2,iglob) = displ(2,iglob)
+ b_PML_field(3,iglob) = displ(3,iglob)
+
+ b_PML_field(4,iglob) = veloc(1,iglob)
+ b_PML_field(5,iglob) = veloc(2,iglob)
+ b_PML_field(6,iglob) = veloc(3,iglob)
+
+ b_PML_field(7,iglob) = accel(1,iglob)
+ b_PML_field(8,iglob) = accel(2,iglob)
+ b_PML_field(9,iglob) = accel(3,iglob)
+ enddo
+
+ call write_abs(0,b_PML_field,b_reclen_PML_field,it)
+
+end subroutine save_field_on_pml_interface
+!
+!=====================================================================
+!
+subroutine read_field_on_pml_interface(b_accel,b_veloc,b_displ,nglob_interface_PML_elastic,&
+ b_PML_field,b_reclen_PML_field)
+ use specfem_par, only: NGLOB_AB,NSPEC_AB,ibool,NSTEP,it
+ use pml_par, only: NSPEC_CPML,CPML_to_spec
+ use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ
+ integer, intent(in) :: nglob_interface_PML_elastic,b_reclen_PML_field
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: b_displ,b_veloc,b_accel
+ real(kind=CUSTOM_REAL), dimension(9,nglob_interface_PML_elastic) :: b_PML_field
+
+ integer :: iglob,ispec,ispec_pml,i,j,k
+
+ do ispec_pml = 1, NSPEC_CPML
+ ispec = CPML_to_spec(ispec_pml)
+ do i = 1, NGLLX; do j = 1, NGLLY; do k = 1, NGLLZ
+ iglob = ibool(i,j,k,ispec)
+ b_displ(:,iglob) = 0._CUSTOM_REAL
+ b_veloc(:,iglob) = 0._CUSTOM_REAL
+ b_accel(:,iglob) = 0._CUSTOM_REAL
+ enddo; enddo; enddo
+ enddo
+
+ call read_abs(0,b_PML_field,b_reclen_PML_field,NSTEP-it+1)
+
+ do iglob = 1, nglob_interface_PML_elastic
+ b_displ(1,iglob) = b_PML_field(1,iglob)
+ b_displ(2,iglob) = b_PML_field(2,iglob)
+ b_displ(3,iglob) = b_PML_field(3,iglob)
+
+ b_veloc(1,iglob) = b_PML_field(4,iglob)
+ b_veloc(2,iglob) = b_PML_field(5,iglob)
+ b_veloc(3,iglob) = b_PML_field(6,iglob)
+
+ b_accel(1,iglob) = b_PML_field(7,iglob)
+ b_accel(2,iglob) = b_PML_field(8,iglob)
+ b_accel(3,iglob) = b_PML_field(9,iglob)
+ enddo
+
+end subroutine read_field_on_pml_interface
+!
+!=====================================================================
+!
+subroutine save_potential_on_pml_interface(potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic,&
+ nglob_interface_PML_acoustic,b_PML_potential,b_reclen_PML_potential)
+ use specfem_par, only: NGLOB_AB,it
+ use constants, only: CUSTOM_REAL
+ integer, intent(in) :: nglob_interface_PML_acoustic,b_reclen_PML_potential
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(3,nglob_interface_PML_acoustic) :: b_PML_potential
+
+ integer :: iglob
+
+ do iglob = 1, nglob_interface_PML_acoustic
+ b_PML_potential(1,iglob) = potential_acoustic(iglob)
+ b_PML_potential(2,iglob) = potential_dot_acoustic(iglob)
+ b_PML_potential(3,iglob) = potential_dot_dot_acoustic(iglob)
+ enddo
+
+ call write_abs(1,b_PML_potential,b_reclen_PML_potential,it)
+
+end subroutine save_potential_on_pml_interface
+!
+!=====================================================================
+!
+subroutine 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)
+ use specfem_par, only: NGLOB_AB,NSPEC_AB,ibool,NSTEP,it
+ use pml_par, only: NSPEC_CPML,CPML_to_spec
+ use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ
+ integer, intent(in) :: nglob_interface_PML_acoustic,b_reclen_PML_potential
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(3,nglob_interface_PML_acoustic) :: b_PML_potential
+
+ integer :: iglob,ispec,ispec_pml,i,j,k
+
+ do ispec_pml = 1, NSPEC_CPML
+ ispec = CPML_to_spec(ispec_pml)
+ do i = 1, NGLLX; do j = 1, NGLLY; do k = 1, NGLLZ
+ iglob = ibool(i,j,k,ispec)
+ b_potential_acoustic(iglob) = 0._CUSTOM_REAL
+ b_potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ b_potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ enddo; enddo; enddo
+ enddo
+
+ call read_abs(1,b_PML_potential,b_reclen_PML_potential,NSTEP-it+1)
+
+ do iglob = 1, nglob_interface_PML_acoustic
+ b_potential_acoustic(iglob) = b_PML_potential(1,iglob)
+ b_potential_dot_acoustic(iglob) = b_PML_potential(2,iglob)
+ b_potential_dot_dot_acoustic(iglob) = b_PML_potential(3,iglob)
+ enddo
+
+end subroutine read_potential_on_pml_interface
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-05-20 23:35:24 UTC (rev 22122)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-05-21 19:43:56 UTC (rev 22123)
@@ -113,6 +113,9 @@
! need the array stored the points on interface between PML and interior computational domain
! --------------------------------------------------------------------------------------------
integer :: nglob_interface_PML_acoustic,nglob_interface_PML_elastic
- integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
+ integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
+ integer :: b_reclen_PML_field,b_reclen_PML_potential !ZN
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_PML_field,b_PML_potential
+
end module pml_par
More information about the CIG-COMMITS
mailing list