[cig-commits] [commit] devel: updates interpolation tool by adding separate flags for separating elements on 410/650 and moho discontinuities (74780ef)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Nov 25 06:56:20 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/0e5b55c6f30be94583639fd325373eecd6facc6d...8be3e0b0267c8d4cf5af3bc26e8903da17bc4fd1
>---------------------------------------------------------------
commit 74780ef3fdf1a13fc083a7feb2c1a5e64d7b903a
Author: daniel peter <peterda at ethz.ch>
Date: Wed Nov 19 15:59:47 2014 +0100
updates interpolation tool by adding separate flags for separating elements on 410/650 and moho discontinuities
>---------------------------------------------------------------
74780ef3fdf1a13fc083a7feb2c1a5e64d7b903a
src/tomography/interpolate_model.F90 | 178 +++++++++++++++++-----------
src/tomography/interpolate_model_kdtree.f90 | 32 ++++-
2 files changed, 134 insertions(+), 76 deletions(-)
diff --git a/src/tomography/interpolate_model.F90 b/src/tomography/interpolate_model.F90
index a0d3279..699f5fe 100644
--- a/src/tomography/interpolate_model.F90
+++ b/src/tomography/interpolate_model.F90
@@ -98,6 +98,15 @@
logical,parameter :: TREE_INTERNAL_GLL_POINTS = .true.
logical,parameter :: TREE_MID_POINTS = .false.
+ ! special case for elements
+ ! around 410-km discontinuity where internal topography distorts meshes
+ logical,parameter :: DO_SEPARATION_410_650 = .true.
+ ! around surface (due to moho-stretching)
+ logical,parameter :: DO_SEPARATION_TOPO = .true.
+
+ ! use closest point value in case of large differences
+ logical,parameter :: USE_FALLBACK = .false.
+
!-------------------------------------------------------------------
! Gauss-Lobatto-Legendre points of integration and weights
@@ -375,13 +384,24 @@
print*,' model2 = ',NGLLX*NGLLY*NGLLZ*NSPEC_CRUST_MANTLE*nparams*dble(CUSTOM_REAL)/1024./1024.,'MB'
print*
print*,'total mpi processes: ',sizeprocs
+ print*
if (DO_BRUTE_FORCE_SEARCH) then
print*,'location search by : brute-force approach'
else
+ print*,'location search by : kd-tree search'
if (USE_MIDPOINT_SEARCH) then
- print*,'location search by : midpoint element search'
+ print*,'location search by : uses midpoints of elements only'
else
- print*,'location search by : gll element search'
+ print*,' uses internal gll points'
+ endif
+ if (DO_SEPARATION_410_650) then
+ print*,' uses element separation for 410-km/650-km discontinuity'
+ endif
+ if (DO_SEPARATION_TOPO) then
+ print*,' uses element separation for surface (moho) discontinuity'
+ endif
+ if (USE_FALLBACK) then
+ print*,' uses fall-back to model value of closest point in case of large differences'
endif
endif
print*
@@ -420,9 +440,9 @@
! model files
allocate( model1(NGLLX,NGLLY,NGLLZ,nspec_max_old,nparams,0:nproc_chunk1-1),stat=ier )
- if (ier /= 0) stop 'Error allocating model1'
+ if (ier /= 0) stop 'Error allocating initial model1'
allocate( model2(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE,nparams),stat=ier )
- if (ier /= 0) stop 'Error allocating model2'
+ if (ier /= 0) stop 'Error allocating target model2'
! statistics
allocate( model_maxdiff(nparams),stat=ier)
if (ier /= 0) stop 'Error allocating model_maxdiff'
@@ -551,6 +571,10 @@
! user output
if (myrank == 0) then
print*, 'loading source model ... '
+ do iker = 1,nparams
+ print*, ' for parameter: ',trim(fname(iker))
+ enddo
+ print*
endif
! reads in old model files
@@ -571,10 +595,8 @@
! reads in model slices
do iker = 1,nparams
- ! user output
- if (myrank == 0) then
- print *, ' for parameter: ',trim(fname(iker))
- endif
+ ! debug user output
+ !if (myrank == 0) print *, ' for parameter: ',trim(fname(iker))
! opens model file
write(m_file,'(a,i6.6,a)') trim(input_model_dir)//'proc',rank,'_reg1_'//trim(fname(iker))//'.bin'
open(IIN,file=trim(m_file),status='old',form='unformatted',iostat=ier)
@@ -833,7 +855,7 @@
call get_model_values_kdtree(ispec,nspec,nglob,ibool2,x2,y2,z2,nparams,model2, &
nspec_max_old,nglob_max_old,nproc_chunk1,ibool1,x1,y1,z1,model1, &
iaddx,iaddy,iaddr,xigll,yigll,zigll,typical_size,myrank,model_maxdiff, &
- USE_MIDPOINT_SEARCH)
+ USE_MIDPOINT_SEARCH,DO_SEPARATION_410_650,DO_SEPARATION_TOPO,USE_FALLBACK)
endif
enddo ! ispec
if (myrank == 0) print*
@@ -1028,7 +1050,7 @@
subroutine get_model_values_kdtree(ispec,nspec,nglob,ibool2,x2,y2,z2,nparams,model2, &
nspec_max_old,nglob_max_old,nproc_chunk1,ibool1,x1,y1,z1,model1, &
iaddx,iaddy,iaddr,xigll,yigll,zigll,typical_size,myrank,model_maxdiff, &
- USE_MIDPOINT_SEARCH)
+ USE_MIDPOINT_SEARCH,DO_SEPARATION_410_650,DO_SEPARATION_TOPO,USE_FALLBACK)
use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGNOD,MIDX,MIDY,MIDZ,R_EARTH_KM,R_EARTH
@@ -1066,6 +1088,8 @@
real(kind=CUSTOM_REAL),dimension(nparams),intent(inout) :: model_maxdiff
logical,intent(in) :: USE_MIDPOINT_SEARCH
+ logical,intent(in) :: DO_SEPARATION_410_650,DO_SEPARATION_TOPO
+ logical,intent(in) :: USE_FALLBACK
! local parameters
integer :: i,j,k,iglob,iker
@@ -1103,19 +1127,9 @@
! (s362ani: 650 topography perturbations have min/max ~ -14/+19 km)
double precision,parameter :: R650_MARGIN = 50000.d0
- !------------------------------------------------------
-
- ! warn about large model value differences
+ ! debug warning about large model value differences
logical,parameter :: DO_WARNING = .false.
- ! special case for elements
- ! e.g. around 410-km discontinuity and surface (due to moho-stretching) where internal topography distorts meshes
- logical,parameter :: DO_SPECIAL_SEPARATION = .true.
-
- ! use closest point value in case of large differences
- logical,parameter :: USE_FALLBACK = .false.
- !------------------------------------------------------
-
! checks given ispec
if (ispec < 1 .or. ispec > nspec) then
print*,'Error: rank ',myrank,' has invalid ispec'
@@ -1128,7 +1142,7 @@
is_critical = .false.
! searches for element using mid-point
- if (USE_MIDPOINT_SEARCH .or. DO_SPECIAL_SEPARATION) then
+ if (USE_MIDPOINT_SEARCH .or. DO_SEPARATION_410_650 .or. DO_SEPARATION_TOPO) then
iglob = ibool2(MIDX,MIDY,MIDZ,ispec)
xyz_target(1) = x2(iglob)
@@ -1150,8 +1164,8 @@
!if (myrank == 0 .and. iglob < 100) &
! print*,'dist_min kdtree midpoint: ',dist_min * R_EARTH_KM,'(km)',typical_size * R_EARTH_KM
- ! special case for 410-km discontinuity
- if (DO_SPECIAL_SEPARATION) then
+ ! special case for 410-km/650-km discontinuity
+ if (DO_SEPARATION_410_650) then
! point radius
r = sqrt(xyz_target(1)*xyz_target(1) + xyz_target(2)*xyz_target(2) + xyz_target(3)*xyz_target(3))
@@ -1172,23 +1186,35 @@
! elements within around 650 km depth
is_critical = .true.
endif
+ endif
- if (is_critical) then
- ! stores mid-point radius
- mid_radius = r
-
- ! element height: size along a vertical edge
- ! top point
- iglob = ibool2(1,1,NGLLZ,ispec)
- r = sqrt(x2(iglob)*x2(iglob) + y2(iglob)*y2(iglob) + z2(iglob)*z2(iglob))
- ! bottom point
- iglob = ibool2(1,1,1,ispec)
- elem_height = r - sqrt(x2(iglob)*x2(iglob) + y2(iglob)*y2(iglob) + z2(iglob)*z2(iglob))
- ! debug
- !if (myrank == 0) print*,'element height: ',elem_height * R_EARTH_KM,'(km)','radius: ',mid_radius*R_EARTH_KM
+ ! special case for surface (moho) discontinuity
+ if (DO_SEPARATION_TOPO) then
+ ! point radius
+ r = sqrt(xyz_target(1)*xyz_target(1) + xyz_target(2)*xyz_target(2) + xyz_target(3)*xyz_target(3))
+
+ ! surface
+ if (r >= R220/R_EARTH) then
+ ! elements close to surface
+ is_critical = .true.
endif
endif
+ if (is_critical) then
+ ! stores mid-point radius
+ mid_radius = r
+
+ ! element height: size along a vertical edge
+ ! top point
+ iglob = ibool2(1,1,NGLLZ,ispec)
+ r = sqrt(x2(iglob)*x2(iglob) + y2(iglob)*y2(iglob) + z2(iglob)*z2(iglob))
+ ! bottom point
+ iglob = ibool2(1,1,1,ispec)
+ elem_height = r - sqrt(x2(iglob)*x2(iglob) + y2(iglob)*y2(iglob) + z2(iglob)*z2(iglob))
+ ! debug
+ !if (myrank == 0) print*,'element height: ',elem_height * R_EARTH_KM,'(km)','radius: ',mid_radius*R_EARTH_KM
+ endif
+
endif
! loops over all element gll points
@@ -1205,54 +1231,64 @@
! kd-search for this single GLL point
if (.not. USE_MIDPOINT_SEARCH) then
+ search_internal = .false.
+
! avoids getting values from "wrong" side on 410-km discontinuity,etc.
- if (DO_SPECIAL_SEPARATION) then
+ if (DO_SEPARATION_410_650) then
if (is_critical) then
! gll point radius
r = sqrt(x_target*x_target + y_target*y_target + z_target*z_target)
! takes corresponding internal gll point for element search
- search_internal = .false.
- ! surface elements
- if (r >= (RTOP - RTOP_MARGIN)/R_EARTH) search_internal = .true.
! 410-km discontinuity
if (r >= (R410 - R410_MARGIN)/R_EARTH .and. r <= (R410 + R410_MARGIN)/R_EARTH) search_internal = .true.
! 650-km discontinuity
if (r >= (R650 - R650_MARGIN)/R_EARTH .and. r <= (R650 + R650_MARGIN)/R_EARTH) search_internal = .true.
+ endif
+ endif
- ! avoid using nodes at upper/lower/..outer surfaces
- if (search_internal) then
- ! new search point indices
- ii = i
- jj = j
- kk = k
-
- ! takes corresponding internal one for element search
- if (i == 1) then
- ii = 2
- else if (i == NGLLX) then
- ii = NGLLX - 1
- endif
- if (j == 1) then
- jj = 2
- else if (j == NGLLY) then
- jj = NGLLY - 1
- endif
- if (k == 1) then
- kk = 2
- else if (k == NGLLZ) then
- kk = NGLLZ - 1
- endif
-
- ! target point location
- iglob = ibool2(ii,jj,kk,ispec)
- x_target = x2(iglob)
- y_target = y2(iglob)
- z_target = z2(iglob)
- endif
+ if (DO_SEPARATION_TOPO) then
+ if (is_critical) then
+ ! gll point radius
+ r = sqrt(x_target*x_target + y_target*y_target + z_target*z_target)
+
+ ! takes corresponding internal gll point for element search
+ ! surface elements
+ if (r >= (RTOP - RTOP_MARGIN)/R_EARTH) search_internal = .true.
endif
endif
+ ! avoid using nodes at upper/lower/..outer surfaces
+ if (search_internal) then
+ ! new search point indices
+ ii = i
+ jj = j
+ kk = k
+
+ ! takes corresponding internal one for element search
+ if (i == 1) then
+ ii = 2
+ else if (i == NGLLX) then
+ ii = NGLLX - 1
+ endif
+ if (j == 1) then
+ jj = 2
+ else if (j == NGLLY) then
+ jj = NGLLY - 1
+ endif
+ if (k == 1) then
+ kk = 2
+ else if (k == NGLLZ) then
+ kk = NGLLZ - 1
+ endif
+
+ ! target point location
+ iglob = ibool2(ii,jj,kk,ispec)
+ x_target = x2(iglob)
+ y_target = y2(iglob)
+ z_target = z2(iglob)
+ endif
+
! kdtree search for each single GLL point
xyz_target(1) = x_target
xyz_target(2) = y_target
diff --git a/src/tomography/interpolate_model_kdtree.f90 b/src/tomography/interpolate_model_kdtree.f90
index 613a246..7d5a98e 100644
--- a/src/tomography/interpolate_model_kdtree.f90
+++ b/src/tomography/interpolate_model_kdtree.f90
@@ -578,7 +578,6 @@ contains
! local parameters
double precision :: dist
- double precision,dimension(3) :: xyz
! debug
!if (node%idim == 0) then
@@ -596,8 +595,7 @@ contains
if (node%ipoint < 1 ) stop 'Error searched node has wrong point index'
! squared distance to associated data point
- xyz(:) = xyz_target(:) - points_data(:,node%ipoint)
- dist = xyz(1) * xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+ dist = get_distance_squared(xyz_target(:),points_data(:,node%ipoint))
if (dist < dist_min) then
! debug
!if (ipoint_min < 1 ) then
@@ -642,7 +640,7 @@ contains
if (associated(node%right) ) then
! checks right node as a final node
if (node%right%idim == 0 ) then
- dist = sum((xyz_target(:) - points_data(:,node%right%ipoint))**2)
+ dist = get_distance_squared(xyz_target(:),points_data(:,node%right%ipoint))
if (dist <= dist_min) then
! stores minimum point
dist_min = dist
@@ -660,7 +658,7 @@ contains
if (associated(node%left) ) then
! checks left node as a final node
if (node%left%idim == 0 ) then
- dist = sum((xyz_target(:) - points_data(:,node%left%ipoint))**2)
+ dist = get_distance_squared(xyz_target(:),points_data(:,node%left%ipoint))
if (dist <= dist_min) then
! stores minimum point
dist_min = dist
@@ -690,4 +688,28 @@ contains
end subroutine kdtree_set_verbose
+
+!===================================================================================================
+
+ double precision function get_distance_squared(xyz0,xyz1)
+
+ implicit none
+ double precision,dimension(3),intent(in) :: xyz0,xyz1
+
+ ! local parameters
+ double precision :: dist
+ !double precision,dimension(3) :: xyz
+
+ ! calculates distance (squared) between 2 points
+ dist = sum((xyz0(:) - xyz1(:))**2)
+
+ ! explicit alternative calculation
+ !xyz(:) = xyz0(:) - xyz1(:)
+ !dist = xyz(1) * xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+
+ ! returns distance squared
+ get_distance_squared = dist
+
+ end function get_distance_squared
+
end module
More information about the CIG-COMMITS
mailing list