[cig-commits] r21813 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Apr 11 03:25:12 PDT 2013
Author: xie.zhinan
Date: 2013-04-11 03:25:12 -0700 (Thu, 11 Apr 2013)
New Revision: 21813
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/part_decompose_mesh.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/read_partition_files.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/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
Log:
add one simple but not general way to compute the width of pml layer
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -100,7 +100,6 @@
integer :: nspec_cpml
integer, dimension(:), allocatable :: CPML_to_spec, CPML_regions
logical, dimension(:), allocatable :: CPML_mask_ibool
- real(kind=CUSTOM_REAL) :: CPML_width
! moho surface (optional)
integer :: nspec2D_moho
@@ -658,7 +657,7 @@
if( ier /= 0 .or. .not. PML_CONDITIONS) then
nspec_cpml = 0
else
- read(98,*) nspec_cpml, CPML_width
+ read(98,*) nspec_cpml
endif
! sanity check
@@ -1083,7 +1082,7 @@
glob2loc_nodes_parts, glob2loc_nodes, part, NGNOD2D)
! writes out C-PML elements indices, CPML-regions and thickness of C-PML layer
- call write_cpml_database(IIN_database, ipart, nspec, nspec_cpml, CPML_width, CPML_to_spec, &
+ call write_cpml_database(IIN_database, ipart, nspec, nspec_cpml, CPML_to_spec, &
CPML_regions, CPML_mask_ibool, glob2loc_elmnts, part)
! gets number of MPI interfaces
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/part_decompose_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/part_decompose_mesh.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/part_decompose_mesh.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -974,8 +974,8 @@
!--------------------------------------------------
! Write C-PML elements indices, CPML-regions and thickness of C-PML layer
! pertaining to iproc partition in the corresponding Database
- !--------------------------------------------------
- subroutine write_cpml_database(IIN_database, iproc, nspec, nspec_cpml, CPML_width, CPML_to_spec, &
+ !--------------------------------------------------
+ subroutine write_cpml_database(IIN_database, iproc, nspec, nspec_cpml, CPML_to_spec, &
CPML_regions, CPML_mask_ibool, glob2loc_elmnts, part)
integer, intent(in) :: IIN_database
@@ -988,8 +988,6 @@
logical, dimension(nspec), intent(in) :: CPML_mask_ibool
- real(kind=CUSTOM_REAL), intent(in) :: CPML_width
-
integer, dimension(:), pointer :: glob2loc_elmnts
integer, dimension(1:nspec), intent(in) :: part
@@ -1011,9 +1009,6 @@
write(IIN_database) nspec_cpml_local
- ! writes thickness of C-PML layers for the global mesh
- write(IIN_database) CPML_width
-
! writes C-PML regions and C-PML spectral elements global indexing
do i=1,nspec_cpml
! #id_cpml_regions = 1 : X_surface C-PML
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -138,8 +138,8 @@
! mask of C-PML elements for the global mesh
logical, dimension(:), allocatable :: CPML_mask_ibool
- ! thickness of C-PML layers
- real(kind=CUSTOM_REAL) :: CPML_width,CPML_width_x,CPML_width_y,CPML_width_z
+ ! thickness of C-PML layers in each direction
+ real(kind=CUSTOM_REAL) :: CPML_width_x,CPML_width_y,CPML_width_z
! C-PML damping profile arrays
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: d_store_x, d_store_y, d_store_z
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-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -32,7 +32,7 @@
use generate_databases_par, only: ibool,NGLOB_AB,d_store_x,d_store_y,d_store_z, &
K_store_x,K_store_y,K_store_z,alpha_store,CPML_to_spec, &
- CPML_width,CPML_width_x,CPML_width_y,CPML_width_z,NPOWER,K_MAX_PML, &
+ CPML_width_x,CPML_width_y,CPML_width_z,NPOWER,K_MAX_PML, &
CUSTOM_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
@@ -59,6 +59,15 @@
real(kind=CUSTOM_REAL) :: x_min,x_min_all,y_min,y_min_all,z_min,z_min_all,&
x_max,x_max_all,y_max,y_max_all,z_max,z_max_all,&
x_origin,y_origin,z_origin
+ real(kind=CUSTOM_REAL) :: CPML_width_x_left, CPML_width_x_right,&
+ CPML_width_y_front,CPML_width_y_back,&
+ CPML_width_z_top,CPML_width_z_bottom,&
+ CPML_x_left, CPML_x_right,&
+ CPML_y_front,CPML_y_back,&
+ CPML_z_top,CPML_z_bottom,&
+ CPML_width_x_left_max_all, CPML_width_x_right_max_all,&
+ CPML_width_y_front_max_all,CPML_width_y_back_max_all,&
+ CPML_width_z_top_max_all,CPML_width_z_bottom_max_all
! stores damping profiles
allocate(d_store_x(NGLLX,NGLLY,NGLLZ,nspec_cpml),stat=ier)
@@ -91,12 +100,6 @@
ALPHA_MAX_PML = PI*f0_FOR_PML ! ELASTIC from Festa and Vilotte (2005)
ALPHA_MAX_PML = PI*f0_FOR_PML*2.0 ! ACOUSTIC from Festa and Vilotte (2005)
- !the idea of fix PML width is bad when element sizes are not the same in x,y,z direction
-
- CPML_width_x = CPML_width
- CPML_width_y = CPML_width
- CPML_width_z = CPML_width
-
x_min = minval(xstore(:))
x_max = maxval(xstore(:))
y_min = minval(ystore(:))
@@ -123,15 +126,96 @@
y_origin = (y_min_all + y_max_all)/2.0
z_origin = (z_max_all + z_min_all)/2.0
- ! determines equations of C-PML/mesh interface planes
- xoriginleft = x_min_all + CPML_width_x
- xoriginright = x_max_all - CPML_width_x
- yoriginback = y_min_all + CPML_width_y
- yoriginfront = y_max_all - CPML_width_y
- zoriginbottom = z_min_all + CPML_width_z
+!Calculation of CPML_width_x,CPML_width_y,CPML_width_Z
+!we assume CPML_width_x,CPML_width_y,CPML_width_Z are constants inside PML layer
+ CPML_width_x_left = 0.0
+ CPML_width_x_right = 0.0
+ CPML_width_y_front = 0.0
+ CPML_width_y_back = 0.0
+ CPML_width_z_top = 0.0
+ CPML_width_z_bottom = 0.0
+
+ CPML_x_right = x_max_all
+ CPML_x_left = x_min_all
+ CPML_y_front = y_max_all
+ CPML_y_back = y_min_all
+ CPML_z_top = z_max_all
+ CPML_z_bottom = z_min_all
+
+ do ispec_CPML=1,nspec_cpml
+ ispec = CPML_to_spec(ispec_CPML)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if( CPML_regions(ispec_CPML) == CPML_X_ONLY ) then
+ if(xstore(iglob) - x_origin > 0.d0)then
+ if(xstore(iglob) - x_origin <= CPML_x_right - x_origin )then
+ CPML_x_right = xstore(iglob)
+ endif
+ else
+ if(abs(xstore(iglob) - x_origin) <= abs(CPML_x_left-x_origin))then
+ CPML_x_left = xstore(iglob)
+ endif
+ endif
+ endif
+
+ if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
+ if(ystore(iglob) - y_origin > 0.d0)then
+ if(ystore(iglob) - y_origin <= CPML_y_front - y_origin )then
+ CPML_y_front = ystore(iglob)
+ endif
+ else
+ if(abs(ystore(iglob) - y_origin) <= abs(CPML_y_back-y_origin))then
+ CPML_y_back = ystore(iglob)
+ endif
+ endif
+ endif
+
+ if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
+ if(zstore(iglob) - z_origin > 0.d0)then
+ if(zstore(iglob) - z_origin <= CPML_z_top - z_origin )then
+ CPML_z_top = zstore(iglob)
+ endif
+ else
+ if(abs(zstore(iglob) - z_origin) <= abs(CPML_z_bottom-z_origin))then
+ CPML_z_bottom = zstore(iglob)
+ endif
+ endif
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ CPML_width_x_right = x_max_all - CPML_x_right
+ CPML_width_x_left = CPML_x_left - x_min_all
+ CPML_width_y_front = y_max_all - CPML_y_front
+ CPML_width_y_back = CPML_y_back - y_min_all
+ CPML_width_z_top = z_max_all - CPML_z_top
+ CPML_width_z_bottom = CPML_z_bottom - z_min_all
+
+ call max_all_all_cr(CPML_width_x_left,CPML_width_x_left_max_all)
+ call max_all_all_cr(CPML_width_x_right,CPML_width_x_right_max_all)
+ call max_all_all_cr(CPML_width_y_front,CPML_width_y_front_max_all)
+ call max_all_all_cr(CPML_width_y_back,CPML_width_y_back_max_all)
+ call max_all_all_cr(CPML_width_z_top,CPML_width_z_top_max_all)
+ call max_all_all_cr(CPML_width_z_bottom,CPML_width_z_bottom_max_all)
+
+ xoriginleft = x_min_all + CPML_width_x_left_max_all
+ xoriginright = x_max_all - CPML_width_x_right_max_all
+ yoriginback = y_min_all + CPML_width_y_front_max_all
+ yoriginfront = y_max_all - CPML_width_y_back_max_all
+ zoriginbottom = z_min_all + CPML_width_z_bottom_max_all
+
+ CPML_width_x = max(CPML_width_x_left_max_all,CPML_width_x_right_max_all)
+ CPML_width_y = max(CPML_width_y_front_max_all,CPML_width_y_back_max_all)
+ CPML_width_z = max(CPML_width_z_bottom_max_all,CPML_width_z_top_max_all )
+
if( PML_INSTEAD_OF_FREE_SURFACE ) then
- zorigintop = z_max_all - CPML_width_z
+ zorigintop = z_max_all - CPML_width_z_top_max_all
endif
! user output
@@ -186,7 +270,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -204,7 +288,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -243,7 +327,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -261,7 +345,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -300,7 +384,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -318,7 +402,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -356,7 +440,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -373,7 +457,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -391,7 +475,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -408,7 +492,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -426,7 +510,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -443,7 +527,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -461,7 +545,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -478,7 +562,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -518,7 +602,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -535,7 +619,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -554,7 +638,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -571,7 +655,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -590,7 +674,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -607,7 +691,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -626,7 +710,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -643,7 +727,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -682,7 +766,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -699,7 +783,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -718,7 +802,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -735,7 +819,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -754,7 +838,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -771,7 +855,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -790,7 +874,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -807,7 +891,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -845,7 +929,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -862,7 +946,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -879,7 +963,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -898,7 +982,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -915,7 +999,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -932,7 +1016,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -951,7 +1035,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -968,7 +1052,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -985,7 +1069,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -1004,7 +1088,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -1021,7 +1105,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -1038,7 +1122,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -1057,7 +1141,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -1074,7 +1158,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -1091,7 +1175,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -1110,7 +1194,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -1127,7 +1211,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -1144,7 +1228,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -1163,7 +1247,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -1180,7 +1264,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -1197,7 +1281,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
@@ -1216,7 +1300,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
- if( abscissa_in_PML_x >= 0.d0 ) then
+ if( abscissa_in_PML_x >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_x / CPML_width_x
@@ -1233,7 +1317,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
- if( abscissa_in_PML_y >= 0.d0 ) then
+ if( abscissa_in_PML_y >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_y / CPML_width_y
@@ -1250,7 +1334,7 @@
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
- if( abscissa_in_PML_z >= 0.d0 ) then
+ if( abscissa_in_PML_z >= -0.1d0 ) then
! determines distance to C-PML/mesh interface
dist = abscissa_in_PML_z / CPML_width_z
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -228,9 +228,6 @@
! checks that the sum of C-PML elements over all partitions is correct
if( myrank == 0 .and. nspec_cpml_tot /= num_cpml ) stop 'error while summing C-PML elements over all partitions'
- ! reads thickness of C-PML layers for the global mesh
- read(IIN) CPML_width
-
! reads C-PML regions and C-PML spectral elements global indexing
allocate(CPML_to_spec(nspec_cpml),stat=ier)
if(ier /= 0) stop 'error allocating array CPML_to_spec'
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -32,7 +32,8 @@
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
SAVE_MESH_FILES,ANISOTROPY)
- use generate_databases_par, only: nspec_cpml,CPML_width,CPML_to_spec,CPML_regions,CPML_mask_ibool,nspec_cpml_tot, &
+ use generate_databases_par, only: nspec_cpml,CPML_width_x,CPML_width_y,CPML_width_z,CPML_to_spec,&
+ CPML_regions,CPML_mask_ibool,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
@@ -131,7 +132,9 @@
! C-PML absorbing boundary conditions
write(IOUT) nspec_cpml
- write(IOUT) CPML_width
+ write(IOUT) CPML_width_x
+ write(IOUT) CPML_width_y
+ write(IOUT) CPML_width_z
if( nspec_cpml > 0 ) then
write(IOUT) CPML_regions
write(IOUT) CPML_to_spec
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -53,7 +53,8 @@
logical, dimension(:), allocatable :: CPML_mask_ibool
! thickness of C-PML layers
- real(CUSTOM_REAL) :: CPML_width,CPML_width_x,CPML_width_y,CPML_width_z
+!ZN real(CUSTOM_REAL) :: CPML_width,CPML_width_x,CPML_width_y,CPML_width_z
+ real(CUSTOM_REAL) :: CPML_width_x,CPML_width_y,CPML_width_z
! C-PML damping profile arrays
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: d_store_x, d_store_y, d_store_z
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -703,11 +703,6 @@
! local parameters
integer :: ispec,ispec_CPML,NSPEC_CPML_GLOBAL
- ! sets thickness of C-PML layers
- CPML_width_x = CPML_width
- CPML_width_y = CPML_width
- CPML_width_z = CPML_width
-
call sum_all_i(NSPEC_CPML,NSPEC_CPML_GLOBAL)
! user output
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-11 03:39:49 UTC (rev 21812)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-11 10:25:12 UTC (rev 21813)
@@ -334,7 +334,9 @@
! C-PML absorbing boundary conditions
read(27) NSPEC_CPML
- read(27) CPML_width
+ read(27) CPML_width_x
+ read(27) CPML_width_y
+ read(27) CPML_width_z
if( PML_CONDITIONS .and. NSPEC_CPML > 0 ) then
allocate(CPML_regions(NSPEC_CPML),stat=ier)
if(ier /= 0) stop 'error allocating array CPML_regions'
More information about the CIG-COMMITS
mailing list