[cig-commits] [commit] devel: uses vectorized loops in smoothing tool; renames smooth_sem.f90 to smooth_sem.F90 (58e8efb)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 07:35:48 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/47b844e0941e7031b9b9330a3ed6e3a5e6b6407a...58e8efb773c4727b4a02f43d1ddc2cde8680f383
>---------------------------------------------------------------
commit 58e8efb773c4727b4a02f43d1ddc2cde8680f383
Author: daniel peter <peterda at ethz.ch>
Date: Thu Dec 11 13:01:17 2014 +0100
uses vectorized loops in smoothing tool; renames smooth_sem.f90 to smooth_sem.F90
>---------------------------------------------------------------
58e8efb773c4727b4a02f43d1ddc2cde8680f383
src/meshfem3D/model_ppm.f90 | 3 +-
src/shared/smooth_weights_vec.F90 | 10 +-
src/tomography/{smooth_sem.f90 => smooth_sem.F90} | 331 +++++++++++-----------
3 files changed, 180 insertions(+), 164 deletions(-)
diff --git a/src/meshfem3D/model_ppm.f90 b/src/meshfem3D/model_ppm.f90
index 477865b..88c3e51 100644
--- a/src/meshfem3D/model_ppm.f90
+++ b/src/meshfem3D/model_ppm.f90
@@ -908,7 +908,8 @@
! checks distance between centers of elements
if (dist_h > sigma_h3 .or. dist_v > sigma_v3 ) cycle
- factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:) ! integration factors
+ ! integration factors
+ factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:)
! loop over GLL points of the elements in current slice (ispec)
do k = 1, NGLLZ
diff --git a/src/shared/smooth_weights_vec.F90 b/src/shared/smooth_weights_vec.F90
index 78f002d..24ad272 100644
--- a/src/shared/smooth_weights_vec.F90
+++ b/src/shared/smooth_weights_vec.F90
@@ -33,14 +33,13 @@
subroutine smoothing_weights_vec(x0,y0,z0,sigma_h2,sigma_v2,exp_val,&
xx_elem,yy_elem,zz_elem)
-
use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGLLCUBE
implicit none
+ real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,sigma_h2,sigma_v2
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: exp_val
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(in) :: xx_elem, yy_elem, zz_elem
- real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,sigma_h2,sigma_v2
! local parameters
real(kind=CUSTOM_REAL) :: dist_h,dist_v
@@ -117,19 +116,19 @@
! Gaussian function
val = - dist_h*sigma_h2_inv - dist_v*sigma_v2_inv
- !exp_val(i,j,k) = exp(val) ! * factor(i,j,k)
+ !exp_val(INDEX_IJK) = exp(val) ! * factor(INDEX_IJK)
! limits to single precision
if (val < - 86.0) then
! smaller than numerical precision: exp(-86) < 1.e-37
exp_val(INDEX_IJK) = 0.0_CUSTOM_REAL
else
- exp_val(INDEX_IJK) = exp(val) ! * factor(i,j,k)
+ exp_val(INDEX_IJK) = exp(val) ! * factor(INDEX_IJK)
endif
! debug
!if (debug) then
- ! print*,i,j,k,'smoothing:',dist_v,dist_h,sigma_h2,sigma_v2,ratio,theta,'val',- dist_h/sigma_h2 - dist_v/sigma_v2
+ ! print*,INDEX_IJK,'smoothing:',dist_v,dist_h,sigma_h2,sigma_v2,ratio,theta,'val',- dist_h/sigma_h2 - dist_v/sigma_v2
!endif
ENDDO_LOOP_IJK
@@ -200,3 +199,4 @@
!dist_h = sqrt( vx*vx + vy*vy + vz*vz )
end subroutine get_distance_vec
+
diff --git a/src/tomography/smooth_sem.f90 b/src/tomography/smooth_sem.F90
similarity index 79%
rename from src/tomography/smooth_sem.f90
rename to src/tomography/smooth_sem.F90
index 17b02b4..25988d4 100644
--- a/src/tomography/smooth_sem.f90
+++ b/src/tomography/smooth_sem.F90
@@ -25,6 +25,10 @@
!
!=====================================================================
+! we switch between vectorized and non-vectorized version by using pre-processor flag FORCE_VECTORIZATION
+! and macros INDEX_IJK, DO_LOOP_IJK, ENDDO_LOOP_IJK defined in config.fh
+#include "config.fh"
+
! smooth_sem_globe
!
! this program can be used for smoothing a (summed) event kernel,
@@ -62,7 +66,7 @@ program smooth_sem_globe
! algorithm uses vector components in radial/horizontal direction
use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,NX_BATHY,NY_BATHY,IIN,IOUT, &
- GAUSSALPHA,GAUSSBETA,PI,TWO_PI,R_EARTH_KM,MAX_STRING_LEN,DEGREES_TO_RADIANS,SIZE_INTEGER
+ GAUSSALPHA,GAUSSBETA,PI,TWO_PI,R_EARTH_KM,MAX_STRING_LEN,DEGREES_TO_RADIANS,SIZE_INTEGER,NGLLCUBE
use kdtree_search
@@ -99,17 +103,22 @@ program smooth_sem_globe
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore, ystore, zstore
real(kind=CUSTOM_REAL), dimension(NSPEC_AB) :: cx0, cy0, cz0, cx, cy, cz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: exp_val
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor
+
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) :: sigma_h, sigma_h2, sigma_h3, sigma_v, sigma_v2, sigma_v3
- real(kind=CUSTOM_REAL) :: x0, y0, z0, norm, norm_h, norm_v, element_size_m
+ real(kind=CUSTOM_REAL) :: norm, norm_h, norm_v, element_size_m
+ real(kind=CUSTOM_REAL) :: x0, y0, z0
+ real(kind=CUSTOM_REAL) :: center_x0, center_y0, center_z0
+ real(kind=CUSTOM_REAL) :: center_x, center_y, center_z
+
real(kind=CUSTOM_REAL) :: max_old, max_new, min_old, min_new
real(kind=CUSTOM_REAL) :: max_old_all, max_new_all, min_old_all, min_new_all
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: exp_val
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor
-
integer :: sizeprocs,ier,myrank,ichunk, ixi, ieta, iglob,nums,ival
- integer :: i,j,k,ispec,iproc,ispec2,inum
+ integer :: ispec,iproc,ispec2,inum
+ integer :: i,j,k
character(len=MAX_STRING_LEN) :: arg(5)
character(len=MAX_STRING_LEN) :: kernel_filename, topo_dir, scratch_file_dir
@@ -148,11 +157,18 @@ program smooth_sem_globe
double precision :: tCPU
double precision, external :: wtime
+#ifdef FORCE_VECTORIZATION
+! in this vectorized version we have to assume that N_SLS == 3 in order to be able to unroll and thus suppress
+! an inner loop that would otherwise prevent vectorization; this is safe in practice in all cases because N_SLS == 3
+! in all known applications, and in the main program we check that N_SLS == 3 if FORCE_VECTORIZATION is used and we stop
+ integer :: ijk
+#endif
+
!------------------------------------------------------------------------------------------
! USER PARAMETERS
! number of steps to reach 100 percent, i.e. 10 outputs info for every 10 percent
- integer,parameter :: NSTEP_PERCENT_INFO = 100
+ integer,parameter :: NSTEP_PERCENT_INFO = 10
! brute-force search for neighboring elements, if set to .false. a kd-tree search will be used
logical,parameter :: DO_BRUTE_FORCE_SEARCH = .false.
@@ -379,34 +395,34 @@ program smooth_sem_globe
! get the location of the center of the elements
do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- xx0(i,j,k,ispec) = xstore(iglob)
- yy0(i,j,k,ispec) = ystore(iglob)
- zz0(i,j,k,ispec) = zstore(iglob)
- ! build jacobian
- ! get derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- ! compute the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- - xiyl*(etaxl*gammazl-etazl*gammaxl) &
- + xizl*(etaxl*gammayl-etayl*gammaxl))
+ DO_LOOP_IJK
+
+ iglob = ibool(INDEX_IJK,ispec)
+ xx0(INDEX_IJK,ispec) = xstore(iglob)
+ yy0(INDEX_IJK,ispec) = ystore(iglob)
+ zz0(INDEX_IJK,ispec) = zstore(iglob)
+
+ ! build jacobian
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(INDEX_IJK,ispec)
+ xiyl = xiy(INDEX_IJK,ispec)
+ xizl = xiz(INDEX_IJK,ispec)
+ etaxl = etax(INDEX_IJK,ispec)
+ etayl = etay(INDEX_IJK,ispec)
+ etazl = etaz(INDEX_IJK,ispec)
+ gammaxl = gammax(INDEX_IJK,ispec)
+ gammayl = gammay(INDEX_IJK,ispec)
+ gammazl = gammaz(INDEX_IJK,ispec)
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ jacobian(INDEX_IJK,ispec) = jacobianl
+
+ ENDDO_LOOP_IJK
- jacobian(i,j,k,ispec) = jacobianl
- enddo
- enddo
- enddo
cx0(ispec) = (xx0(1,1,1,ispec) + xx0(NGLLX,NGLLY,NGLLZ,ispec)) * 0.5_CUSTOM_REAL
cy0(ispec) = (yy0(1,1,1,ispec) + yy0(NGLLX,NGLLY,NGLLZ,ispec)) * 0.5_CUSTOM_REAL
cz0(ispec) = (zz0(1,1,1,ispec) + zz0(NGLLX,NGLLY,NGLLZ,ispec)) * 0.5_CUSTOM_REAL
@@ -528,33 +544,33 @@ program smooth_sem_globe
! get the location of the center of the elements
do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- xx(i,j,k,ispec) = xstore(iglob)
- yy(i,j,k,ispec) = ystore(iglob)
- zz(i,j,k,ispec) = zstore(iglob)
-
- ! build jacobian
- ! get derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- ! compute the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- - xiyl*(etaxl*gammazl-etazl*gammaxl) &
- + xizl*(etaxl*gammayl-etayl*gammaxl))
- jacobian(i,j,k,ispec) = jacobianl
- enddo
- enddo
- enddo
+
+ DO_LOOP_IJK
+
+ iglob = ibool(INDEX_IJK,ispec)
+ xx(INDEX_IJK,ispec) = xstore(iglob)
+ yy(INDEX_IJK,ispec) = ystore(iglob)
+ zz(INDEX_IJK,ispec) = zstore(iglob)
+
+ ! build jacobian
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(INDEX_IJK,ispec)
+ xiyl = xiy(INDEX_IJK,ispec)
+ xizl = xiz(INDEX_IJK,ispec)
+ etaxl = etax(INDEX_IJK,ispec)
+ etayl = etay(INDEX_IJK,ispec)
+ etazl = etaz(INDEX_IJK,ispec)
+ gammaxl = gammax(INDEX_IJK,ispec)
+ gammayl = gammay(INDEX_IJK,ispec)
+ gammazl = gammaz(INDEX_IJK,ispec)
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(INDEX_IJK,ispec) = jacobianl
+
+ ENDDO_LOOP_IJK
+
! calculate element center location
cx(ispec) = (xx(1,1,1,ispec) + xx(NGLLX,NGLLY,NGLLZ,ispec)) * 0.5_CUSTOM_REAL
cy(ispec) = (yy(1,1,1,ispec) + yy(NGLLX,NGLLY,NGLLZ,ispec)) * 0.5_CUSTOM_REAL
@@ -627,11 +643,16 @@ program smooth_sem_globe
! initializes
num_elem_local = 0
+ ! element center position
+ center_x0 = cx0(ispec)
+ center_y0 = cy0(ispec)
+ center_z0 = cz0(ispec)
+
! sets number of elements to loop over
if (.not. DO_BRUTE_FORCE_SEARCH) then
- xyz_target(1) = cx0(ispec)
- xyz_target(2) = cy0(ispec)
- xyz_target(3) = cz0(ispec)
+ xyz_target(1) = center_x0
+ xyz_target(2) = center_y0
+ xyz_target(3) = center_z0
! counts nearest neighbors
if (DO_SEARCH_ELLIP) then
@@ -695,34 +716,27 @@ program smooth_sem_globe
! --- only double loop over the elements in the search radius ---
do ielem = 1, num_elem_local
-
+ ! search element
if (.not. DO_BRUTE_FORCE_SEARCH) then
! kd-tree search elements
ispec2 = kdtree_search_index(ielem)
-
- ! (spherical search radius)
- ! calculates horizontal and vertical distance between two element centers
- ! vector approximation
- call get_distance_vec(dist_h,dist_v,cx0(ispec),cy0(ispec),cz0(ispec),&
- cx(ispec2),cy(ispec2),cz(ispec2))
-
- ! note: distances and sigmah, sigmav are normalized by R_EARTH_KM
-
- ! checks distance between centers of elements
- if ( dist_h > sigma_h3 .or. dist_v > sigma_v3 ) cycle
else
! brute-force search
- ! takes only elements in search radius
- ! calculates horizontal and vertical distance between two element centers
ispec2 = ielem
- call get_distance_vec(dist_h,dist_v,cx0(ispec),cy0(ispec),cz0(ispec),&
- cx(ispec2),cy(ispec2),cz(ispec2))
+ endif
- ! note: distances and sigmah, sigmav are normalized by R_EARTH_KM
+ ! search element center position
+ center_x = cx(ispec2)
+ center_y = cy(ispec2)
+ center_z = cz(ispec2)
- ! checks distance between centers of elements
- if ( dist_h > sigma_h3 .or. dist_v > sigma_v3 ) cycle
- endif
+ ! takes only elements in search radius
+ ! calculates horizontal and vertical distance between two element centers
+ call get_distance_vec(dist_h,dist_v,center_x0,center_y0,center_z0,center_x,center_y,center_z)
+
+ ! checks distance between centers of elements
+ ! note: distances and sigmah, sigmav are normalized by R_EARTH_KM
+ if ( dist_h > sigma_h3 .or. dist_v > sigma_v3 ) cycle
! integration factors:
! uses volume assigned to GLL points
@@ -730,54 +744,57 @@ program smooth_sem_globe
! no volume
!factor(:,:,:) = 1.0_CUSTOM_REAL
- ! loop over GLL points of the elements in current slice (ispec)
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
-
- ! reference location
- ! current point (i,j,k,ispec) location, cartesian coordinates
- x0 = xx0(i,j,k,ispec)
- y0 = yy0(i,j,k,ispec)
- z0 = zz0(i,j,k,ispec)
-
- ! calculate weights based on gaussian smoothing
- call smoothing_weights_vec(x0,y0,z0,sigma_h2,sigma_v2,exp_val,&
- xx(:,:,:,ispec2),yy(:,:,:,ispec2),zz(:,:,:,ispec2))
-
- ! debug
- if (DEBUG .and. myrank == 0) then
- if (i == 2 .and. j == 3 .and. k == 4 .and. ispec == tmp_ispec_dbg) then
- tmp_bk(:,:,:,ispec2) = exp_val(:,:,:)
- print*,'debug',myrank,'ispec',ispec,'ispec2',ispec2,'dist:',dist_h,dist_v
- print*,'debug exp',minval(exp_val),maxval(exp_val)
- print*,'debug kernel',minval(kernel(:,:,:,ispec2)),maxval(kernel(:,:,:,ispec2))
- endif
+ ! loops over all GLL points of reference element in current slice (ispec) and search elements (ispec2)
+ DO_LOOP_IJK
+
+ ! reference location
+ ! current point (i,j,k,ispec) location, cartesian coordinates
+ x0 = xx0(INDEX_IJK,ispec)
+ y0 = yy0(INDEX_IJK,ispec)
+ z0 = zz0(INDEX_IJK,ispec)
+
+ ! calculate weights based on gaussian smoothing
+ call smoothing_weights_vec(x0,y0,z0,sigma_h2,sigma_v2,exp_val,&
+ xx(:,:,:,ispec2),yy(:,:,:,ispec2),zz(:,:,:,ispec2))
+
+ ! debug
+ if (DEBUG .and. myrank == 0) then
+ if (ispec == tmp_ispec_dbg) then
+#ifdef FORCE_VECTORIZATION
+ if (ijk == (2 + 3*NGLLX + 4*NGLLY*NGLLX)) then
+#else
+ if (i == 2 .and. j == 3 .and. k == 4) then
+#endif
+ tmp_bk(:,:,:,ispec2) = exp_val(:,:,:)
+ print*,'debug',myrank,'ispec',ispec,'ispec2',ispec2,'dist:',dist_h,dist_v
+ print*,'debug exp',minval(exp_val),maxval(exp_val)
+ print*,'debug kernel',minval(kernel(:,:,:,ispec2)),maxval(kernel(:,:,:,ispec2))
endif
+ endif
+ endif
- ! adds GLL integration weights
- exp_val(:,:,:) = exp_val(:,:,:) * factor(:,:,:)
+ ! adds GLL integration weights
+ exp_val(:,:,:) = exp_val(:,:,:) * factor(:,:,:)
- ! adds contribution of element ispec2 to smoothed kernel values
- tk(i,j,k,ispec) = tk(i,j,k,ispec) + sum(exp_val(:,:,:) * kernel(:,:,:,ispec2))
+ ! adds contribution of element ispec2 to smoothed kernel values
+ tk(INDEX_IJK,ispec) = tk(INDEX_IJK,ispec) + sum(exp_val(:,:,:) * kernel(:,:,:,ispec2))
- ! normalization, integrated values of gaussian smoothing function
- bk(i,j,k,ispec) = bk(i,j,k,ispec) + sum(exp_val(:,:,:))
+ ! normalization, integrated values of gaussian smoothing function
+ bk(INDEX_IJK,ispec) = bk(INDEX_IJK,ispec) + sum(exp_val(:,:,:))
- ! checks number
- !if (isNaN(tk(i,j,k,ispec))) then
- ! print*,'Error tk NaN: ',tk(i,j,k,ispec)
- ! print*,'rank:',myrank
- ! print*,'i,j,k,ispec:',i,j,k,ispec
- ! print*,'tk: ',tk(i,j,k,ispec),'bk:',bk(i,j,k,ispec)
- ! print*,'sum exp_val: ',sum(exp_val(:,:,:)),'sum factor:',sum(factor(:,:,:))
- ! print*,'sum kernel:',sum(kernel(:,:,:,ispec2))
- ! call exit_mpi(myrank, 'Error NaN')
- !endif
+ ! checks number
+ !if (isNaN(tk(INDEX_IJK,ispec))) then
+ ! print*,'Error tk NaN: ',tk(INDEX_IJK,ispec)
+ ! print*,'rank:',myrank
+ ! print*,'INDEX_IJK,ispec:',INDEX_IJK,ispec
+ ! print*,'tk: ',tk(INDEX_IJK,ispec),'bk:',bk(INDEX_IJK,ispec)
+ ! print*,'sum exp_val: ',sum(exp_val(:,:,:)),'sum factor:',sum(factor(:,:,:))
+ ! print*,'sum kernel:',sum(kernel(:,:,:,ispec2))
+ ! call exit_mpi(myrank, 'Error NaN')
+ !endif
+
+ ENDDO_LOOP_IJK
- enddo
- enddo
- enddo ! (i,j,k)
enddo ! (ielem)
! frees search results
@@ -794,8 +811,7 @@ program smooth_sem_globe
if (ispec == tmp_ispec_dbg) then
! outputs gaussian weighting function
write(filename,'(a,i4.4,a,i6.6)') trim(scratch_file_dir)//'/search_elem',tmp_ispec_dbg,'_gaussian_proc',iproc
- call write_VTK_data_elem_cr(NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
- ibool,tmp_bk,filename)
+ call write_VTK_data_elem_cr(NSPEC_AB,NGLOB_AB,xstore,ystore,zstore,ibool,tmp_bk,filename)
print*,'file written: ',trim(filename)//'.vtk'
endif
endif
@@ -838,42 +854,41 @@ program smooth_sem_globe
! compute the smoothed kernel values
kernel_smooth(:,:,:,:) = 0.0_CUSTOM_REAL
do ispec = 1, NSPEC_AB
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
-
- ! checks the normalization criterion
- ! e.g. sigma_h 160km, sigma_v 40km:
- ! norm (not squared sigma_h ) ~ 0.001
- ! norm ( squared sigma_h) ~ 6.23 * e-5
- !if (abs(bk(i,j,k,ispec) - norm) > 1.e-4 ) then
- ! print *, 'Problem norm here --- ', myrank, ispec, i, j, k, bk(i,j,k,ispec), norm
- ! !call exit_mpi(myrank, 'Error computing Gaussian function on the grid')
- !endif
- !debug
- !if (tk(i,j,k,ispec) == 0.0_CUSTOM_REAL .and. myrank == 0) then
- ! print*,myrank,'zero tk: ',i,j,k,ispec,'tk:',tk(i,j,k,ispec),'bk:',bk(i,j,k,ispec)
- !endif
+ DO_LOOP_IJK
+ ! checks the normalization criterion
+ ! e.g. sigma_h 160km, sigma_v 40km:
+ ! norm (not squared sigma_h ) ~ 0.001
+ ! norm ( squared sigma_h) ~ 6.23 * e-5
+ !if (abs(bk(INDEX_IJK,ispec) - norm) > 1.e-4 ) then
+ ! print *, 'Problem norm here --- ', myrank, ispec, i, j, k, bk(INDEX_IJK,ispec), norm
+ ! !call exit_mpi(myrank, 'Error computing Gaussian function on the grid')
+ !endif
- ! avoids division by zero
- if (bk(i,j,k,ispec) == 0.0_CUSTOM_REAL) bk(i,j,k,ispec) = 1.0_CUSTOM_REAL
+ !debug
+ !if (tk(INDEX_IJK,ispec) == 0.0_CUSTOM_REAL .and. myrank == 0) then
+ ! print*,myrank,'zero tk: ',INDEX_IJK,ispec,'tk:',tk(INDEX_IJK,ispec),'bk:',bk(INDEX_IJK,ispec)
+ !endif
- ! normalizes smoothed kernel values by integral value of gaussian weighting
- kernel_smooth(i,j,k,ispec) = tk(i,j,k,ispec) / bk(i,j,k,ispec)
- ! checks number (isNaN check)
- if (kernel_smooth(i,j,k,ispec) /= kernel_smooth(i,j,k,ispec)) then
- print*,'Error kernel_smooth value not a number: ',kernel_smooth(i,j,k,ispec)
- print*,'rank:',myrank
- print*,'i,j,k,ispec:',i,j,k,ispec
- print*,'tk: ',tk(i,j,k,ispec),'bk:',bk(i,j,k,ispec),'norm:',norm
- call exit_mpi(myrank, 'Error kernel value is NaN')
- endif
- enddo
- enddo
- enddo
+ ! avoids division by zero
+ if (bk(INDEX_IJK,ispec) == 0.0_CUSTOM_REAL) bk(INDEX_IJK,ispec) = 1.0_CUSTOM_REAL
+
+ ! normalizes smoothed kernel values by integral value of gaussian weighting
+ kernel_smooth(INDEX_IJK,ispec) = tk(INDEX_IJK,ispec) / bk(INDEX_IJK,ispec)
+
+ ! checks number (isNaN check)
+ if (kernel_smooth(INDEX_IJK,ispec) /= kernel_smooth(INDEX_IJK,ispec)) then
+ print*,'Error kernel_smooth value not a number: ',kernel_smooth(INDEX_IJK,ispec)
+ print*,'rank:',myrank
+ print*,'INDEX_IJK,ispec:',INDEX_IJK,ispec
+ print*,'tk: ',tk(INDEX_IJK,ispec),'bk:',bk(INDEX_IJK,ispec),'norm:',norm
+ call exit_mpi(myrank, 'Error kernel value is NaN')
+ endif
+
+ ENDDO_LOOP_IJK
+
enddo
! statistics
More information about the CIG-COMMITS
mailing list