[cig-commits] r21755 - in seismo/3D/SPECFEM3D/trunk/src: shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Apr 7 18:05:40 PDT 2013
Author: dkomati1
Date: 2013-04-07 18:05:40 -0700 (Sun, 07 Apr 2013)
New Revision: 21755
Modified:
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
Log:
fixed element size calculation to avoid including the diagonals (which is non standard)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-04-08 00:27:31 UTC (rev 21754)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-04-08 01:05:40 UTC (rev 21755)
@@ -927,8 +927,8 @@
subroutine get_GLL_minmaxdistance(distance_min,distance_max,ispec, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
-! calculates the min/max distances between neighboring GLL points within
-! the specified element (ispec)
+! calculates the min/max distances between neighboring GLL points within the specified element (ispec);
+! we purposely do not include the distance along the diagonals of the element, only along its three coordinate axes.
implicit none
@@ -943,41 +943,54 @@
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
! local parameters
- real(kind=CUSTOM_REAL) :: dx,x0,y0,z0
- integer :: i,j,k,ii,jj,kk,iglob_a,iglob_b
+ real(kind=CUSTOM_REAL) :: dist
+ integer :: i,j,k,iglob_a,iglob_b
! initializes
distance_min = HUGEVAL
distance_max = -HUGEVAL
- ! loops over all GLL points
- do k=1,NGLLZ-1
- do j=1,NGLLY-1
+ ! loops over all GLL points along X
+ do k=1,NGLLZ
+ do j=1,NGLLY
do i=1,NGLLX-1
- iglob_a = ibool(i,j,k,ispec)
- x0 = xstore(iglob_a)
- y0 = ystore(iglob_a)
- z0 = zstore(iglob_a)
+ iglob_a = ibool(i ,j,k,ispec)
+ iglob_b = ibool(i+1,j,k,ispec)
+ dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ + ( ystore(iglob_a) - ystore(iglob_b) )**2 &
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ if( dist < distance_min) distance_min = dist
+ if( dist > distance_max) distance_max = dist
+ enddo
+ enddo
+ enddo
- ! loops over nearest neighbor points
- ! maybe a faster method could be found...
- do kk=k-1,k+1
- do jj=j-1,j+1
- do ii=i-1,i+1
- if( ii < 1 .or. jj < 1 .or. kk < 1 ) cycle
- ! distances between points
- iglob_b = ibool(ii,jj,kk,ispec)
- if( iglob_a /= iglob_b) then
- dx = sqrt( ( x0 - xstore(iglob_b) )**2 &
- + ( y0 - ystore(iglob_b) )**2 &
- + ( z0 - zstore(iglob_b) )**2 )
- if( dx < distance_min) distance_min = dx
- if( dx > distance_max) distance_max = dx
- endif
- enddo
- enddo
- enddo
+ ! loops over all GLL points along Y
+ do k=1,NGLLZ
+ do i=1,NGLLX
+ do j=1,NGLLY-1
+ iglob_a = ibool(i,j ,k,ispec)
+ iglob_b = ibool(i,j+1,k,ispec)
+ dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ + ( ystore(iglob_a) - ystore(iglob_b) )**2 &
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ if( dist < distance_min) distance_min = dist
+ if( dist > distance_max) distance_max = dist
+ enddo
+ enddo
+ enddo
+ ! loops over all GLL points along Z
+ do i=1,NGLLX
+ do j=1,NGLLY
+ do k=1,NGLLZ-1
+ iglob_a = ibool(i,j,k ,ispec)
+ iglob_b = ibool(i,j,k+1,ispec)
+ dist = sqrt( ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ + ( ystore(iglob_a) - ystore(iglob_b) )**2 &
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2 )
+ if( dist < distance_min) distance_min = dist
+ if( dist > distance_max) distance_max = dist
enddo
enddo
enddo
@@ -991,7 +1004,8 @@
subroutine get_elem_minmaxsize(elemsize_min,elemsize_max,ispec, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
-! calculates the min/max size of the specified element (ispec)
+! calculates the min/max size of an edge of the specified element (ispec);
+! we purposely do not include the distance along the diagonals of the element, only the size of its edges.
implicit none
@@ -1006,49 +1020,56 @@
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
! local parameters
- real(kind=CUSTOM_REAL) :: dx,x0,y0,z0
- integer :: i,j,k,icorner,jcorner,iglob_a,iglob_b
+ real(kind=CUSTOM_REAL) :: dist
+ integer :: i,j,k
+ integer :: i1,i2,j1,j2,k1,k2,ibool1,ibool2
- ! corners indices of reference cube faces
- ! shapes of arrays below
- integer,dimension(2),parameter :: corner_shape = (/3,NGNOD_EIGHT_CORNERS/)
- integer,dimension(3,NGNOD_EIGHT_CORNERS),parameter :: corner_ijk = &
- reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ, &
- NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),corner_shape)
-
! initializes
elemsize_min = HUGEVAL
elemsize_max = -HUGEVAL
- ! loops over corners
- do icorner=1,NGNOD_EIGHT_CORNERS
- i = corner_ijk(1,icorner)
- j = corner_ijk(2,icorner)
- k = corner_ijk(3,icorner)
+ ! loops over the four edges that are along X
+ i1 = 1
+ i2 = NGLLX
+ do j = 1, NGLLY, NGLLY-1
+ do k = 1, NGLLZ, NGLLZ-1
+ ibool1 = ibool(i1,j,k,ispec)
+ ibool2 = ibool(i2,j,k,ispec)
+ dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ + (ystore(ibool1) - ystore(ibool2))**2 &
+ + (zstore(ibool1) - zstore(ibool2))**2)
+ if(dist < elemsize_min) elemsize_min = dist
+ if(dist > elemsize_max) elemsize_max = dist
+ enddo
+ enddo
- iglob_a = ibool(i,j,k,ispec)
- x0 = xstore(iglob_a)
- y0 = ystore(iglob_a)
- z0 = zstore(iglob_a)
+ ! loops over the four edges that are along Y
+ j1 = 1
+ j2 = NGLLY
+ do i = 1, NGLLX, NGLLX-1
+ do k = 1, NGLLZ, NGLLZ-1
+ ibool1 = ibool(i,j1,k,ispec)
+ ibool2 = ibool(i,j2,k,ispec)
+ dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ + (ystore(ibool1) - ystore(ibool2))**2 &
+ + (zstore(ibool1) - zstore(ibool2))**2)
+ if(dist < elemsize_min) elemsize_min = dist
+ if(dist > elemsize_max) elemsize_max = dist
+ enddo
+ enddo
- ! loops over all other corners
- do jcorner = icorner+1,NGNOD_EIGHT_CORNERS
- i = corner_ijk(1,jcorner)
- j = corner_ijk(2,jcorner)
- k = corner_ijk(3,jcorner)
-
- ! coordinates
- iglob_b = ibool(i,j,k,ispec)
-
- ! distances between points
- if( iglob_a /= iglob_b) then
- dx = sqrt( ( x0 - xstore(iglob_b) )**2 &
- + ( y0 - ystore(iglob_b) )**2 &
- + ( z0 - zstore(iglob_b) )**2 )
- if( dx < elemsize_min) elemsize_min = dx
- if( dx > elemsize_max) elemsize_max = dx
- endif
-
+ ! loops over the four edges that are along Z
+ k1 = 1
+ k2 = NGLLZ
+ do i = 1, NGLLX, NGLLX-1
+ do j = 1, NGLLY, NGLLY-1
+ ibool1 = ibool(i,j,k1,ispec)
+ ibool2 = ibool(i,j,k2,ispec)
+ dist = sqrt((xstore(ibool1) - xstore(ibool2))**2 &
+ + (ystore(ibool1) - ystore(ibool2))**2 &
+ + (zstore(ibool1) - zstore(ibool2))**2)
+ if(dist < elemsize_min) elemsize_min = dist
+ if(dist > elemsize_max) elemsize_max = dist
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-08 00:27:31 UTC (rev 21754)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-08 01:05:40 UTC (rev 21755)
@@ -52,7 +52,7 @@
! mask of C-PML elements for the global mesh
logical, dimension(:), allocatable :: CPML_mask_ibool
- ! thickness of C-PML layers read in "absorbing_cpml_file.asc"
+ ! thickness of C-PML layers
real(CUSTOM_REAL) :: CPML_width,CPML_width_x,CPML_width_y,CPML_width_z
! C-PML damping profile arrays
More information about the CIG-COMMITS
mailing list