[Yt-svn] yt: 4 new changesets
hg at spacepope.org
hg at spacepope.org
Tue Feb 22 14:09:38 PST 2011
hg Repository: yt
details: yt/rev/120365de1086
changeset: 3759:120365de1086
user: Matthew Turk <matthewturk at gmail.com>
date:
Sat Feb 19 21:07:19 2011 -0500
description:
This speeds things up a small amount, but I'm still avoiding adding ray packets
directly to grids.
hg Repository: yt
details: yt/rev/0579873ce4b6
changeset: 3760:0579873ce4b6
user: Matthew Turk <matthewturk at gmail.com>
date:
Mon Feb 21 00:05:12 2011 -0500
description:
Attempting to use a linked list for adaptive healpix renderings. There is
still an infinite loop problem with it, however.
hg Repository: yt
details: yt/rev/f96012b2adb9
changeset: 3761:f96012b2adb9
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Feb 22 17:09:08 2011 -0500
description:
I have identified the problem with the linked lists for the adaptive
ray-tracing, and I believe this fixes it.
The collision detection to identify the next grid into which rays pass is still
an issue, but I believe that the current method (of assuming a radially sorted
progression) is the best option at the moment. Future versions will improve
upon this; I think an obvious next step is to perform an initial pass inside
integrate_brick, set up a linked list of entries, and then to iterate over only
that linked list.
Anyway, there may still be some major bugs, and it's still somewhat slow
(although it only takes 4 minutes to generate an effective 1650^2 image on my
machine) and I have yet to figure out how to visualize, but at least all the
rays are traveling now.
hg Repository: yt
details: yt/rev/4e7845e84144
changeset: 3762:4e7845e84144
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Feb 22 17:09:16 2011 -0500
description:
Merging
diffstat:
yt/frontends/enzo/io.py | 3 +
yt/utilities/_amr_utils/VolumeIntegrator.pyx | 110 +++++++++---
yt/utilities/parallel_tools/parallel_analysis_interface.py | 7 +-
yt/visualization/volume_rendering/camera.py | 15 +-
4 files changed, 101 insertions(+), 34 deletions(-)
diffs (292 lines):
diff -r 2164bd47f1f9 -r 4e7845e84144 yt/frontends/enzo/io.py
--- a/yt/frontends/enzo/io.py Sat Feb 19 17:46:24 2011 -0500
+++ b/yt/frontends/enzo/io.py Tue Feb 22 17:09:16 2011 -0500
@@ -159,6 +159,9 @@
return field.swapaxes(0,2)
def preload(self, grids, sets):
+ if len(grids) == 0:
+ data = None
+ return
# We need to deal with files first
files_keys = defaultdict(lambda: [])
pf_field_list = grids[0].pf.h.field_list
diff -r 2164bd47f1f9 -r 4e7845e84144 yt/utilities/_amr_utils/VolumeIntegrator.pyx
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Sat Feb 19 17:46:24 2011 -0500
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Feb 22 17:09:16 2011 -0500
@@ -408,6 +408,8 @@
for k in range(nk):
tv[offset + k] = fv[k]
+cdef struct AdaptiveRayPacket
+
cdef class PartitionedGrid:
cdef public object my_data
cdef public object LeftEdge
@@ -828,6 +830,8 @@
np.float64_t pos[3]
AdaptiveRayPacket *next
AdaptiveRayPacket *prev
+ AdaptiveRayPacket *brick_next
+ #int cgi
cdef class AdaptiveRaySource:
cdef np.float64_t center[3]
@@ -835,20 +839,35 @@
cdef AdaptiveRayPacket *first
cdef public np.float64_t normalization
cdef public int nrays
+ cdef public int max_nside
+ cdef AdaptiveRayPacket **packet_pointers
+ cdef AdaptiveRayPacket **lpacket_pointers
def __cinit__(self, center, rays_per_cell, initial_nside,
- np.float64_t normalization):
+ np.float64_t normalization, brick_list, int max_nside = 8096):
cdef int i
+ self.max_nside = max_nside
self.center[0] = center[0]
self.center[1] = center[1]
self.center[2] = center[2]
self.rays_per_cell = rays_per_cell
cdef AdaptiveRayPacket *ray = self.first
cdef AdaptiveRayPacket *last = NULL
+ cdef PartitionedGrid pg
cdef double v_dir[3]
+ cdef int nbricks = len(brick_list)
+ # You see, killbots have a preset kill limit. Knowing their weakness, I
+ # sent wave after wave of my own men at them until they reached their
+ # limit and shut down.
+ self.lpacket_pointers = <AdaptiveRayPacket **> \
+ malloc(sizeof(AdaptiveRayPacket*)*nbricks)
+ self.packet_pointers = <AdaptiveRayPacket **> \
+ malloc(sizeof(AdaptiveRayPacket*)*nbricks)
+ for i in range(nbricks):
+ self.lpacket_pointers[i] = self.packet_pointers[i] = NULL
self.normalization = normalization
self.nrays = 12*initial_nside*initial_nside
- for i in range(12*initial_nside*initial_nside):
+ for i in range(self.nrays):
# Initialize rays here
ray = <AdaptiveRayPacket *> malloc(sizeof(AdaptiveRayPacket))
ray.prev = last
@@ -861,9 +880,19 @@
ray.v_dir[2] = v_dir[2] * normalization
ray.value[0] = ray.value[1] = ray.value[2] = ray.value[3] = 0.0
ray.next = NULL
- if i == 0: self.first = ray
- if last != NULL: last.next = ray
+ #ray.cgi = 0
+ ray.pos[0] = self.center[0]
+ ray.pos[1] = self.center[1]
+ ray.pos[2] = self.center[2]
+ ray.brick_next = NULL
+ if last != NULL:
+ last.next = ray
+ last.brick_next = ray
+ else:
+ self.first = ray
last = ray
+ self.packet_pointers[0] = self.first
+ self.lpacket_pointers[0] = last
def __dealloc__(self):
cdef AdaptiveRayPacket *next
@@ -872,6 +901,7 @@
next = ray.next
free(ray)
ray = next
+ free(self.packet_pointers)
def get_rays(self):
cdef AdaptiveRayPacket *ray = self.first
@@ -894,26 +924,51 @@
ray = ray.next
return info, values
- def integrate_brick(self, PartitionedGrid pg, TransferFunctionProxy tf):
- cdef AdaptiveRayPacket *ray = self.first
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def integrate_brick(self, PartitionedGrid pg, TransferFunctionProxy tf,
+ int pgi, np.ndarray[np.float64_t, ndim=2] ledges,
+ np.ndarray[np.float64_t, ndim=2] redges):
cdef np.float64_t domega
domega = self.get_domega(pg.left_edge, pg.right_edge)
#print "dOmega", domega, self.nrays
cdef int count = 0
cdef int i
+ cdef AdaptiveRayPacket *ray = self.packet_pointers[pgi]
+ cdef AdaptiveRayPacket *next
while ray != NULL:
# Note that we may end up splitting a ray such that it ends up
# outside the brick!
#print count
count +=1
- if self.intersects(ray, pg):
- ray = self.refine_ray(ray, domega, pg.dds[0],
- pg.left_edge, pg.right_edge)
- pg.integrate_ray(self.center, ray.v_dir, ray.value,
- tf, &ray.t)
- for i in range(3):
- ray.pos[i] = ray.v_dir[i] * (ray.t + 1e-8) + self.center[i]
- ray = ray.next
+ #if count > 10+self.nrays or ray.cgi != pgi:
+ # raise RuntimeError
+ # We don't need to check for intersection anymore, as we are the
+ # Ruler of the planet Omicron Persei 8
+ #if self.intersects(ray, pg):
+ ray = self.refine_ray(ray, domega, pg.dds[0],
+ pg.left_edge, pg.right_edge)
+ pg.integrate_ray(self.center, ray.v_dir, ray.value,
+ tf, &ray.t)
+ for i in range(3):
+ ray.pos[i] = ray.v_dir[i] * (ray.t + 1e-8) + self.center[i]
+ # We set 'next' after the refinement has occurred
+ next = ray.brick_next
+ for i in range(pgi+1, ledges.shape[0]):
+ if ((ledges[i, 0] <= ray.pos[0] <= redges[i, 0]) and
+ (ledges[i, 1] <= ray.pos[1] <= redges[i, 1]) and
+ (ledges[i, 2] <= ray.pos[2] <= redges[i, 2])):
+ if self.lpacket_pointers[i] == NULL:
+ self.packet_pointers[i] = \
+ self.lpacket_pointers[i] = ray
+ ray.brick_next = NULL
+ else:
+ self.lpacket_pointers[i].brick_next = ray
+ self.lpacket_pointers[i] = ray
+ ray.brick_next = NULL
+ #ray.cgi = i
+ break
+ ray = next
cdef int intersects(self, AdaptiveRayPacket *ray, PartitionedGrid pg):
cdef np.float64_t pos[3]
@@ -942,18 +997,11 @@
r2[2] = r2[1] + (edge[k][2] - self.center[2])**2.0
max_r2 = fmax(max_r2, r2[2])
domega = 4.0 * 3.1415926 * max_r2 # Used to be / Nrays
- #print domega, max_r2**0.5, self.nrays,
- #print self.center[0],
- #print self.center[1],
- #print self.center[2],
- #print edge[0][0],
- #print edge[0][1],
- #print edge[0][2],
- #print edge[1][0],
- #print edge[1][1],
- #print edge[1][2]
return domega
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
cdef AdaptiveRayPacket *refine_ray(self, AdaptiveRayPacket *ray,
np.float64_t domega,
np.float64_t dx,
@@ -962,32 +1010,42 @@
# This should be recursive; we are not correctly applying split
# criteria multiple times.
cdef long Nrays = 12 * ray.nside * ray.nside
+ cdef int i, j
if domega/Nrays < dx*dx/self.rays_per_cell:
#print ray.nside, domega/Nrays, dx, (domega/Nrays * self.rays_per_cell)/(dx*dx)
return ray
- if ray.nside == 8192: return ray
+ if ray.nside >= self.max_nside: return ray
#print "Refining %s from %s to %s" % (ray.ipix, ray.nside, ray.nside*2)
# Now we make four new ones
cdef double v_dir[3]
cdef AdaptiveRayPacket *prev = ray.prev
+ # It is important to note here that brick_prev is a local variable for
+ # the newly created rays, not the previous ray in this brick, as that
+ # has already been passed on to its next brick
+ cdef AdaptiveRayPacket *brick_prev = NULL
for i in range(4):
new_ray = <AdaptiveRayPacket *> malloc(
sizeof(AdaptiveRayPacket))
new_ray.nside = ray.nside * 2
new_ray.ipix = ray.ipix * 4 + i
new_ray.t = ray.t
+ #new_ray.cgi = ray.cgi
new_ray.prev = prev
if new_ray.prev != NULL:
new_ray.prev.next = new_ray
- prev = new_ray
+ if brick_prev != NULL:
+ brick_prev.brick_next = new_ray
+ prev = brick_prev = new_ray
healpix_interface.pix2vec_nest(
new_ray.nside, new_ray.ipix, v_dir)
for j in range(3):
new_ray.v_dir[j] = v_dir[j] * self.normalization
new_ray.value[j] = ray.value[j]
+ new_ray.pos[j] = self.center[j] + ray.t * new_ray.v_dir[j]
new_ray.value[3] = ray.value[3]
new_ray.next = ray.next
+ new_ray.brick_next = ray.brick_next
if new_ray.next != NULL:
new_ray.next.prev = new_ray
if self.first == ray:
diff -r 2164bd47f1f9 -r 4e7845e84144 yt/utilities/parallel_tools/parallel_analysis_interface.py
--- a/yt/utilities/parallel_tools/parallel_analysis_interface.py Sat Feb 19 17:46:24 2011 -0500
+++ b/yt/utilities/parallel_tools/parallel_analysis_interface.py Tue Feb 22 17:09:16 2011 -0500
@@ -1041,13 +1041,16 @@
ncols = -1
size = 0
else:
- if len(data.shape) == 1:
+ if len(data) == 0:
+ ncols = -1
+ size = 0
+ elif len(data.shape) == 1:
ncols = 1
size = data.shape[0]
else:
ncols, size = data.shape
ncols = MPI.COMM_WORLD.allreduce(ncols, op=MPI.MAX)
- if data is None:
+ if size == 0:
data = na.empty((ncols,0), dtype='float64') # This only works for
size = data.shape[-1]
sizes = na.zeros(MPI.COMM_WORLD.size, dtype='int64')
diff -r 2164bd47f1f9 -r 4e7845e84144 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py Sat Feb 19 17:46:24 2011 -0500
+++ b/yt/visualization/volume_rendering/camera.py Tue Feb 22 17:09:16 2011 -0500
@@ -604,7 +604,7 @@
transfer_function = None, fields = None,
sub_samples = 5, log_fields = None, volume = None,
pf = None, use_kd=True, no_ghost=False,
- rays_per_cell = 5.1):
+ rays_per_cell = 0.1):
if pf is not None: self.pf = pf
self.center = na.array(center, dtype='float64')
self.radius = radius
@@ -625,8 +625,6 @@
self.rays_per_cell = rays_per_cell
def snapshot(self, fn = None):
- ray_source = AdaptiveRaySource(self.center, self.rays_per_cell,
- self.initial_nside, self.radius)
tfp = TransferFunctionProxy(self.transfer_function)
tfp.ns = self.sub_samples
self.volume.initialize_source()
@@ -634,9 +632,14 @@
pbar = get_pbar("Ray casting",
(self.volume.brick_dimensions + 1).prod(axis=-1).sum())
total_cells = 0
- bricks = [b for b in self.volume.traverse(None, self.center, None)]
- for brick in reversed(bricks):
- ray_source.integrate_brick(brick, tfp)
+ bricks = [b for b in self.volume.traverse(None, self.center, None)][::-1]
+ left_edges = na.array([b.LeftEdge for b in bricks])
+ right_edges = na.array([b.RightEdge for b in bricks])
+ ray_source = AdaptiveRaySource(self.center, self.rays_per_cell,
+ self.initial_nside, self.radius,
+ bricks)
+ for i,brick in enumerate(bricks):
+ ray_source.integrate_brick(brick, tfp, i, left_edges, right_edges)
total_cells += na.prod(brick.my_data[0].shape)
pbar.update(total_cells)
pbar.finish()
More information about the yt-svn
mailing list