[Yt-svn] yt: Adding support for arbitrary periods for all tools that use ...
hg at spacepope.org
hg at spacepope.org
Mon Dec 13 14:43:34 PST 2010
hg Repository: yt
details: yt/rev/9d0d3dbb249e
changeset: 3602:9d0d3dbb249e
user: Stephen Skory <stephenskory at yahoo.com>
date:
Mon Dec 13 15:43:13 2010 -0700
description:
Adding support for arbitrary periods for all tools that use the fortran
kdtree: parallelHF, halo merger tree, and the two point functions toolkit.
diffstat:
yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py | 1 +
yt/analysis_modules/halo_merger_tree/merger_tree.py | 5 +
yt/analysis_modules/two_point_functions/two_point_functions.py | 3 +-
yt/utilities/kdtree/fKD.f90 | 6 +-
yt/utilities/kdtree/fKD.v | 5 +-
yt/utilities/kdtree/fKD_source.f90 | 47 ++++++---
6 files changed, 45 insertions(+), 22 deletions(-)
diffs (truncated from 331 to 300 lines):
diff -r 309549706702 -r 9d0d3dbb249e yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py
--- a/yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/analysis_modules/halo_finding/parallel_hop/parallel_hop_interface.py Mon Dec 13 15:43:13 2010 -0700
@@ -358,6 +358,7 @@
fKD.sort = True # Slower, but needed in _connect_chains
fKD.rearrange = self.rearrange # True is faster, but uses more memory
# Now call the fortran.
+ fKD.period = self.period
create_tree(0)
self.__max_memory()
yt_counters("init kd tree")
diff -r 309549706702 -r 9d0d3dbb249e yt/analysis_modules/halo_merger_tree/merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/merger_tree.py Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/analysis_modules/halo_merger_tree/merger_tree.py Mon Dec 13 15:43:13 2010 -0700
@@ -164,6 +164,7 @@
self.dm_only = dm_only
self.refresh = refresh
self.sleep = sleep # How long to wait between db sync checks.
+ self.period = -1
if self.sleep <= 0.:
self.sleep = 5
# MPI stuff
@@ -215,6 +216,9 @@
for cycle, file in enumerate(self.restart_files):
gc.collect()
pf = load(file)
+ # get the period, only once is sufficient
+ if self.period == -1:
+ self.period = pf.domain_right_edge - pf.domain_left_edge
# If the halos are already found, skip this data step, unless
# refresh is True.
dir = os.path.dirname(file)
@@ -365,6 +369,7 @@
fKD.nn = 5
fKD.sort = True
fKD.rearrange = True
+ fKD.period = self.period
create_tree(0)
# Find the parent points from the database.
diff -r 309549706702 -r 9d0d3dbb249e yt/analysis_modules/two_point_functions/two_point_functions.py
--- a/yt/analysis_modules/two_point_functions/two_point_functions.py Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/analysis_modules/two_point_functions/two_point_functions.py Mon Dec 13 15:43:13 2010 -0700
@@ -29,7 +29,7 @@
from yt.utilities.parallel_tools.parallel_analysis_interface import ParallelAnalysisInterface, parallel_blocking_call, parallel_root_only
try:
- from yt.extensions.kdtree import *
+ from yt.utilities.kdtree import *
except ImportError:
mylog.debug("The Fortran kD-Tree did not import correctly.")
@@ -306,6 +306,7 @@
fKD.nn = 1
fKD.sort = False
fKD.rearrange = True
+ fKD.period = self.period
create_tree(0)
def _build_sort_array(self):
diff -r 309549706702 -r 9d0d3dbb249e yt/utilities/kdtree/fKD.f90
--- a/yt/utilities/kdtree/fKD.f90 Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/utilities/kdtree/fKD.f90 Mon Dec 13 15:43:13 2010 -0700
@@ -19,11 +19,11 @@
! create a kd tree object
if (ID == 1) then
- t1%tree2 => kdtree2_create(t1%pos,sort=t1%sort,rearrange=t1%rearrange)
+ t1%tree2 => kdtree2_create(t1%pos,sort=t1%sort,rearrange=t1%rearrange,period=t1%period)
elseif (ID == 2) then
- t2%tree2 => kdtree2_create(t2%pos,sort=t2%sort,rearrange=t2%rearrange)
+ t2%tree2 => kdtree2_create(t2%pos,sort=t2%sort,rearrange=t2%rearrange,period=t2%period)
else
- tree2 => kdtree2_create(pos,sort=sort,rearrange=rearrange)
+ tree2 => kdtree2_create(pos,sort=sort,rearrange=rearrange,period=period)
end if
return
diff -r 309549706702 -r 9d0d3dbb249e yt/utilities/kdtree/fKD.v
--- a/yt/utilities/kdtree/fKD.v Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/utilities/kdtree/fKD.v Mon Dec 13 15:43:13 2010 -0700
@@ -24,6 +24,7 @@
tree2 _kdtree2
sort logical /.false./
rearrange logical /.true./
+period(:) _real
radius real # the unsquared radius for radius searches
radius_n integer # the number of results to return
nfound integer # number of neighbors within radius
@@ -97,7 +98,8 @@
root _tree_node
#type(tree_node), pointer :: root => null()
# root pointer of the tree
-
+period(:) _real
+# the an array giving the periodicity of the data.
%%%% tree_set:
# This contains the objects that any given kD tree might wish to use.
@@ -121,6 +123,7 @@
tree2 _kdtree2
sort logical /.false./
rearrange logical /.true./
+period(:) _real
radius real # the unsquared radius for radius searches
radius_n integer # the number of results to return
nfound integer # number of neighbors within radius
diff -r 309549706702 -r 9d0d3dbb249e yt/utilities/kdtree/fKD_source.f90
--- a/yt/utilities/kdtree/fKD_source.f90 Mon Dec 13 13:58:37 2010 -0800
+++ b/yt/utilities/kdtree/fKD_source.f90 Mon Dec 13 15:43:13 2010 -0700
@@ -593,6 +593,7 @@
! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0
integer :: nalloc ! how much allocated for results(:)?
logical :: rearrange ! are the data rearranged or original?
+ real, pointer :: period(:)
! did the # of points found overflow the storage provided?
logical :: overflow
real(kdkind), pointer :: qv(:) ! query vector
@@ -609,7 +610,7 @@
contains
- function kdtree2_create(input_data,dim,sort,rearrange) result (mr)
+ function kdtree2_create(input_data,dim,sort,rearrange,period) result (mr)
!
! create the actual tree structure, given an input array of data.
!
@@ -640,6 +641,7 @@
! ..
! .. Array Arguments ..
real(kdkind), target :: input_data(:,:)
+ real, target :: period(:)
!
integer :: i
! ..
@@ -666,7 +668,7 @@
end if
call build_tree(mr)
-
+
if (present(sort)) then
mr%sort = sort
else
@@ -678,6 +680,9 @@
else
mr%rearrange = .true.
endif
+
+ mr%period => period
+
if (mr%rearrange) then
allocate(mr%rearranged_data(mr%dimen,mr%n))
@@ -1052,6 +1057,7 @@
sr%ind => tp%ind
sr%rearrange = tp%rearrange
+ sr%period => tp%period
if (tp%rearrange) then
sr%Data => tp%rearranged_data
else
@@ -1095,6 +1101,7 @@
sr%ind => tp%ind
sr%rearrange = tp%rearrange
+ sr%period => tp%period
if (sr%rearrange) then
sr%Data => tp%rearranged_data
@@ -1147,6 +1154,7 @@
sr%overflow = .false.
sr%ind => tp%ind
sr%rearrange= tp%rearrange
+ sr%period => tp%period
if (tp%rearrange) then
sr%Data => tp%rearranged_data
@@ -1212,6 +1220,7 @@
sr%ind => tp%ind
sr%rearrange = tp%rearrange
+ sr%period => tp%period
if (tp%rearrange) then
sr%Data => tp%rearranged_data
@@ -1266,6 +1275,7 @@
! for counting.
sr%ind => tp%ind
sr%rearrange = tp%rearrange
+ sr%period => tp%period
if (tp%rearrange) then
sr%Data => tp%rearranged_data
else
@@ -1315,6 +1325,7 @@
sr%ind => tp%ind
sr%rearrange = tp%rearrange
+ sr%period => tp%period
if (sr%rearrange) then
sr%Data => tp%rearranged_data
@@ -1352,7 +1363,7 @@
return
end subroutine validate_query_storage
- function square_distance(d, iv,qv) result (res)
+ function square_distance(d, iv,qv,period) result (res)
! distance between iv[1:n] and qv[1:n]
! .. Function Return Value ..
! re-implemented to improve vectorization.
@@ -1365,11 +1376,12 @@
! ..
! .. Array Arguments ..
real(kdkind) :: iv(:),qv(:)
+ real, intent(in) :: period(:)
! ..
! ..
! .. Periodicity added by S Skory
! res = sum( (iv(1:d)-qv(1:d))**2 )
- d_min = min( abs(iv(1:d) - qv(1:d)) , 1. - abs(iv(1:d) - qv(1:d)) )
+ d_min = min( abs(iv(1:d) - qv(1:d)) , period - abs(iv(1:d) - qv(1:d)) )
res = sum(d_min(1:d)**2)
end function square_distance
@@ -1410,9 +1422,9 @@
! node edge is closer, rather than just doing a simple less than or
! greater than to the cut in this dimension.
dis_left = min( (node%box(cut_dim)%lower - qval)**2, &
- (1 - abs(node%box(cut_dim)%lower - qval))**2)
+ (sr%period(cut_dim) - abs(node%box(cut_dim)%lower - qval))**2)
dis_right = min( (node%box(cut_dim)%upper - qval)**2, &
- (1 - abs(node%box(cut_dim)%upper - qval))**2)
+ (sr%period(cut_dim) - abs(node%box(cut_dim)%upper - qval))**2)
if (qval < node%cut_val) then
@@ -1447,7 +1459,7 @@
box => node%box(1:)
do i=1,sr%dimen
if (i .ne. cut_dim) then
- dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper)
+ dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper,sr%period(i))
if (dis > ballsize) then
return
endif
@@ -1482,17 +1494,18 @@
! return
! end function dis2_from_bnd
- real(kdkind) function dis2_from_bnd(x,amin,amax) result (res)
+ real(kdkind) function dis2_from_bnd(x,amin,amax,period) result (res)
! Periodicity added by S Skory
real(kdkind), intent(in) :: x, amin,amax
+ real, intent(in) :: period
real :: dxmax, dxmin
if ((x < amax) .and. (x > amin)) then
res = 0.0
return
else
- dxmax = (min( abs(x - amax), 1. - abs(x - amax)))**2
- dxmin = (min( abs(x - amin), 1. - abs(x - amin)))**2
+ dxmax = (min( abs(x - amax), period - abs(x - amax)))**2
+ dxmin = (min( abs(x - amin), period - abs(x - amin)))**2
res = min(dxmax, dxmin)
return
endif
@@ -1520,7 +1533,7 @@
do i=1,dimen
l = node%box(i)%lower
u = node%box(i)%upper
- dis = dis + (dis2_from_bnd(sr%qv(i),l,u))
+ dis = dis + (dis2_from_bnd(sr%qv(i),l,u,sr%period(i)))
if (dis > ballsize) then
res = .false.
return
@@ -1574,7 +1587,7 @@
do k = 1,dimen
!sd = sd + (data(k,i) - qv(k))**2
! Periodicity by S Skory
- sdtemp = min( abs(data(k,i) - qv(k)) , 1. - abs(data(k,i) - qv(k)))
+ sdtemp = min( abs(data(k,i) - qv(k)) , sr%period(i) - abs(data(k,i) - qv(k)))
sd = sd + sdtemp**2
if (sd>ballsize) cycle mainloop
end do
@@ -1584,7 +1597,7 @@
sd = 0.0
do k = 1,dimen
!sd = sd + (data(k,indexofi) - qv(k))**2
- sdtemp = min( abs(data(k,indexofi) - qv(k)), 1. - abs(data(k,indexofi) - qv(k)))
+ sdtemp = min( abs(data(k,indexofi) - qv(k)), sr%period(i) - abs(data(k,indexofi) - qv(k)))
sd = sd + sdtemp**2
if (sd>ballsize) cycle mainloop
end do
@@ -1702,7 +1715,7 @@
do k = 1,dimen
!sd = sd + (data(k,i) - qv(k))**2
! Periodicity S Skory
- sdtemp = min( abs(data(k,i) - qv(k)) , 1. - abs(data(k,i) - qv(k)))
More information about the yt-svn
mailing list