[cig-commits] r22060 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon May 13 23:37:06 PDT 2013
Author: xie.zhinan
Date: 2013-05-13 23:37:06 -0700 (Mon, 13 May 2013)
New Revision: 22060
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
Log:
commit part of code aiming for implementing PML for adjoint inverse simulation
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -300,7 +300,7 @@
endif
if( ACOUSTIC_SIMULATION) then
deallocate(rmass_acoustic)
- deallocate(rmass_acoustic_interface) !ZN
+ deallocate(rmass_acoustic_interface)
deallocate(rmassz_acoustic)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -389,10 +389,14 @@
endif
write(IMAIN,*)
- if(STACEY_ABSORBING_CONDITIONS) then
- write(IMAIN,*) 'incorporating absorbing conditions'
+ if(STACEY_ABSORBING_CONDITIONS) then
+ write(IMAIN,*) 'incorporating Stacey absorbing conditions'
else
- write(IMAIN,*) 'no absorbing condition'
+ if(PML_CONDITIONS) then
+ write(IMAIN,*) 'incorporating absorbing conditions of perfectly matched layer'
+ else
+ write(IMAIN,*) 'no absorbing condition'
+ endif
endif
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -148,6 +148,11 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: K_store_x, K_store_y, K_store_z
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: alpha_store
+ ! array recording the points on interface shared by PML and interior computational domain
+ logical, dimension(:), allocatable :: mask_ibool_interior_domain
+ integer :: nglob_interface_PML_acoustic,nglob_interface_PML_elastic
+ integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
+
! moho (optional)
integer :: nspec2D_moho_ext
integer, dimension(:), allocatable :: ibelm_moho
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -35,7 +35,10 @@
CPML_width_x,CPML_width_y,CPML_width_z,NPOWER,K_MAX_PML, &
CUSTOM_REAL,SIZE_REAL,NGLLX,NGLLY,NGLLZ,nspec_cpml,PML_INSTEAD_OF_FREE_SURFACE, &
IMAIN,FOUR_THIRDS,CPML_REGIONS,f0_FOR_PML,PI, &
- CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
+ CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ,&
+ SIMULATION_TYPE,SAVE_FORWARD,nspec => NSPEC_AB,is_CPML,&
+ mask_ibool_interior_domain,nglob_interface_PML_acoustic,points_interface_PML_acoustic,&
+ nglob_interface_PML_elastic,points_interface_PML_elastic
use create_regions_mesh_ext_par, only: rhostore,rho_vp,ispec_is_acoustic,ispec_is_elastic, &
ELASTIC_SIMULATION, ACOUSTIC_SIMULATION
@@ -255,7 +258,6 @@
call max_all_all_cr(vp_max,vp_max_all)
-
! user output
if( myrank == 0 ) then
write(IMAIN,*)
@@ -322,17 +324,11 @@
alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- ! avoid d_x to be less than zero due to
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ ! avoid d_x to be less than zero
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -345,16 +341,10 @@
alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer x"
endif
@@ -389,16 +379,10 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else if( ystore(iglob) - y_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -411,16 +395,10 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer y"
@@ -457,15 +435,10 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
endif
else if( zstore(iglob) - z_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
@@ -479,16 +452,10 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer z"
endif
@@ -534,26 +501,14 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. ystore(iglob) - y_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -577,28 +532,16 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. ystore(iglob) - y_origin>0._CUSTOM_REAL ) then
- ! gets abscissa of current grid point along the damping profile
+ ! gets abscissa of current grid point along the damping profilenspec
abscissa_in_PML_x = xoriginleft - xstore(iglob)
! determines distance to C-PML/mesh interface
@@ -620,26 +563,14 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin <0._CUSTOM_REAL .and. ystore(iglob) - y_origin <0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -663,26 +594,14 @@
alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer xy"
endif
@@ -730,26 +649,13 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
@@ -775,26 +681,14 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
@@ -819,26 +713,13 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
@@ -864,26 +745,13 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer xz"
endif
@@ -930,25 +798,13 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
endif
else if( ystore(iglob) - y_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
@@ -974,26 +830,14 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else if( ystore(iglob) - y_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
@@ -1018,26 +862,13 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( ystore(iglob) - y_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
@@ -1063,26 +894,14 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer yz"
endif
@@ -1141,35 +960,17 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
endif
else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
@@ -1208,37 +1009,18 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
-
else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
@@ -1277,36 +1059,17 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
@@ -1347,36 +1110,18 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
ystore(iglob) - y_origin>0._CUSTOM_REAL .and. &
zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
@@ -1414,36 +1159,17 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
@@ -1482,36 +1208,18 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
@@ -1549,36 +1257,17 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
-
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
endif
else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
@@ -1617,36 +1306,18 @@
alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_x < 0._CUSTOM_REAL .or. K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL; K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1._CUSTOM_REAL ) then
- d_x = 0._CUSTOM_REAL
- K_x = 1._CUSTOM_REAL
+ if( d_y < 0._CUSTOM_REAL .or. K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL; K_y = 1._CUSTOM_REAL
endif
- if( d_y < 0._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
+ if( d_z < 0._CUSTOM_REAL .or. K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL; K_z = 1._CUSTOM_REAL
endif
- if( K_y < 1._CUSTOM_REAL ) then
- d_y = 0._CUSTOM_REAL
- K_y = 1._CUSTOM_REAL
- endif
-
- if( d_z < 0._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
- if( K_z < 1._CUSTOM_REAL ) then
- d_z = 0._CUSTOM_REAL
- K_z = 1._CUSTOM_REAL
- endif
-
else
stop "there is error in mesh of CPML-layer xyz"
endif
@@ -1669,6 +1340,103 @@
enddo
enddo !ispec_CPML
+! --------------------------------------------------------------------------------------------
+! for adjoint tomography
+! create the array store the points on interface between PML and interior computational domain
+! --------------------------------------------------------------------------------------------
+
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+
+ !mask all points belong interior computational domain
+ allocate(mask_ibool_interior_domain(NGLOB_AB),stat=ier)
+ if(ier /= 0) stop 'error allocating array mask_ibool_interior_domain'
+ mask_ibool_interior_domain = .false.
+ do ispec=1,nspec
+ if(.not. is_CPML(ispec)) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ mask_ibool_interior_domain(ibool(i,j,k,ispec))=.true.
+ enddo; enddo; enddo
+ endif
+ enddo
+
+ !------------------------------------------------------
+ ! begin of acoustic domain
+ !------------------------------------------------------
+ nglob_interface_PML_acoustic = 0
+
+ if(ACOUSTIC_SIMULATION) then
+
+ do ispec=1,nspec
+ if( ispec_is_acoustic(ispec) .and. is_CPML(ispec)) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
+ nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
+ endif
+ enddo; enddo; enddo
+ endif
+ enddo
+
+ if(nglob_interface_PML_acoustic > 0)then
+ allocate(points_interface_PML_acoustic(nglob_interface_PML_acoustic),stat=ier)
+ if(ier /= 0) stop 'error allocating array points_interface_PML_acoustic'
+ points_interface_PML_acoustic = 0
+ nglob_interface_PML_acoustic = 0
+ do ispec=1,nspec
+ if( ispec_is_acoustic(ispec) .and. is_CPML(ispec)) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
+ nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
+ points_interface_PML_acoustic(nglob_interface_PML_acoustic) = ibool(i,j,k,ispec)
+ endif
+ enddo; enddo; enddo
+ endif
+ enddo
+ endif
+
+ endif
+ !------------------------------------------------------
+ ! end of acoustic domains
+ !------------------------------------------------------
+
+ !------------------------------------------------------
+ ! begin of elastic domains
+ !------------------------------------------------------
+ nglob_interface_PML_elastic = 0
+
+ if(ELASTIC_SIMULATION) then
+
+ do ispec=1,nspec
+ if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
+ nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
+ endif
+ enddo; enddo; enddo
+ endif
+ enddo
+
+ if(nglob_interface_PML_elastic > 0)then
+ allocate(points_interface_PML_elastic(nglob_interface_PML_elastic),stat=ier)
+ if(ier /= 0) stop 'error allocating array points_interface_PML_elastic'
+ points_interface_PML_elastic = 0
+ nglob_interface_PML_elastic = 0
+ do ispec=1,nspec
+ if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
+ do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
+ if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
+ nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
+ points_interface_PML_elastic(nglob_interface_PML_elastic) = ibool(i,j,k,ispec)
+ endif
+ enddo; enddo; enddo
+ endif
+ enddo
+ endif
+
+ endif
+ !------------------------------------------------------
+ ! end of elastic domains
+ !------------------------------------------------------
+ endif
end subroutine pml_set_local_dampingcoeff
!
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -36,7 +36,11 @@
CPML_regions,is_CPML,nspec_cpml_tot, &
d_store_x,d_store_y,d_store_z,k_store_x,k_store_y,k_store_z,alpha_store, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top,PML_CONDITIONS
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top,PML_CONDITIONS,&
+ !for adjoint tomography
+ SIMULATION_TYPE,SAVE_FORWARD,mask_ibool_interior_domain, &
+ nglob_interface_PML_acoustic,points_interface_PML_acoustic,&
+ nglob_interface_PML_elastic,points_interface_PML_elastic
use create_regions_mesh_ext_par
implicit none
@@ -98,7 +102,7 @@
! acoustic
if( ACOUSTIC_SIMULATION ) then
write(IOUT) rmass_acoustic
- write(IOUT) rmass_acoustic_interface !ZN
+ write(IOUT) rmass_acoustic_interface
endif
! this array is needed for acoustic simulations but also for elastic simulations with CPML,
@@ -147,6 +151,16 @@
write(IOUT) k_store_y
write(IOUT) k_store_z
write(IOUT) alpha_store
+ ! --------------------------------------------------------------------------------------------
+ ! for adjoint tomography
+ ! save the array stored the points on interface between PML and interior computational domain
+ ! --------------------------------------------------------------------------------------------
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+ write(IOUT) nglob_interface_PML_acoustic
+ write(IOUT) nglob_interface_PML_elastic
+ if(nglob_interface_PML_acoustic > 0) write(IOUT) points_interface_PML_acoustic
+ if(nglob_interface_PML_elastic > 0) write(IOUT) points_interface_PML_elastic
+ endif
endif
! absorbing boundary surface
@@ -337,6 +351,20 @@
deallocate(k_store_y,stat=ier); if( ier /= 0 ) stop 'error deallocating array d_store_y'
deallocate(k_store_z,stat=ier); if( ier /= 0 ) stop 'error deallocating array d_store_z'
deallocate(alpha_store,stat=ier); if( ier /= 0 ) stop 'error deallocating array alpha_store'
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+ deallocate(mask_ibool_interior_domain,stat=ier)
+ if(ier /= 0) stop 'error deallocating array mask_ibool_interior_domain'
+
+ if(nglob_interface_PML_acoustic > 0) then
+ deallocate(points_interface_PML_acoustic,stat=ier)
+ if( ier /= 0 ) stop 'error deallocating array points_interface_PML_acoustic'
+ endif
+
+ if(nglob_interface_PML_elastic > 0) then
+ deallocate(points_interface_PML_elastic,stat=ier)
+ if( ier /= 0 ) stop 'error deallocating array points_interface_PML_elastic'
+ endif
+ endif
endif
end subroutine save_arrays_solver_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -108,4 +108,11 @@
! C-PML contribution to update displacement on elastic/acoustic interface
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: rmemory_coupling_ac_el_displ
+ ! --------------------------------------------------------------------------------------------
+ ! for adjoint tomography
+ ! 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
+
end module pml_par
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-05-14 01:16:51 UTC (rev 22059)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-05-14 06:37:06 UTC (rev 22060)
@@ -374,6 +374,22 @@
read(27) k_store_y
read(27) k_store_z
read(27) alpha_store
+
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+ read(27) nglob_interface_PML_acoustic
+ read(27) nglob_interface_PML_elastic
+ if(nglob_interface_PML_acoustic > 0) then
+ allocate(points_interface_PML_acoustic(nglob_interface_PML_acoustic),stat=ier)
+ if(ier /= 0) stop 'error allocating array points_interface_PML_acoustic'
+ read(27) points_interface_PML_acoustic
+ endif
+
+ if(nglob_interface_PML_elastic > 0) then
+ allocate(points_interface_PML_elastic(nglob_interface_PML_elastic),stat=ier)
+ if(ier /= 0) stop 'error allocating array points_interface_PML_elastic'
+ read(27) points_interface_PML_elastic
+ endif
+ endif
endif
! absorbing boundary surface
More information about the CIG-COMMITS
mailing list