From 333436c9605d7be6615d0ea3e386e206c494c9a9 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 21 Jan 2020 15:49:59 +0000 Subject: [PATCH 01/61] Add support for fast neighbour search --- setup.py | 5 + yt/frontends/ramses/data_structures.py | 2 +- yt/frontends/ramses/io_utils.pyx | 2 +- yt/geometry/fake_octree.pyx | 6 +- yt/geometry/oct_container.pxd | 2 - yt/geometry/oct_container.pyx | 5 +- yt/geometry/ramses_oct_container.pxd | 15 +++ yt/geometry/ramses_oct_container.pyx | 139 +++++++++++++++++++++++++ 8 files changed, 165 insertions(+), 11 deletions(-) create mode 100644 yt/geometry/ramses_oct_container.pxd create mode 100644 yt/geometry/ramses_oct_container.pyx diff --git a/setup.py b/setup.py index 0c3f8fea9d0..9123917a441 100644 --- a/setup.py +++ b/setup.py @@ -97,6 +97,11 @@ def _compile( "yt/utilities/lib/tsearch.c"], include_dirs=["yt/utilities/lib"], libraries=std_libs), + Extension("yt.geometry.ramses_oct_container", + ["yt/geometry/ramses_oct_container.pyx", + "yt/utilities/lib/tsearch.c"], + include_dirs=["yt/utilities/lib"], + libraries=std_libs), Extension("yt.geometry.oct_visitors", ["yt/geometry/oct_visitors.pyx"], include_dirs=["yt/utilities/lib/"], diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index a65262209ea..2335987c80e 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -27,7 +27,7 @@ from .particle_handlers import get_particle_handlers from .field_handlers import get_field_handlers from yt.utilities.cython_fortran_utils import FortranFile as fpu -from yt.geometry.oct_container import \ +from yt.geometry.ramses_oct_container import \ RAMSESOctreeContainer from yt.arraytypes import blankRecordArray diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 160cdb62b39..12fc0918f33 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -2,7 +2,7 @@ cimport cython cimport numpy as np import numpy as np from yt.utilities.cython_fortran_utils cimport FortranFile -from yt.geometry.oct_container cimport RAMSESOctreeContainer +from yt.geometry.ramses_oct_container cimport RAMSESOctreeContainer from yt.utilities.exceptions import YTIllDefinedAMRData ctypedef np.int32_t INT32_t diff --git a/yt/geometry/fake_octree.pyx b/yt/geometry/fake_octree.pyx index e69dc92165f..375f6a5856f 100644 --- a/yt/geometry/fake_octree.pyx +++ b/yt/geometry/fake_octree.pyx @@ -13,11 +13,11 @@ from oct_visitors cimport cind import numpy as np cimport cython -from oct_container cimport Oct, RAMSESOctreeContainer +from oct_container cimport Oct, SparseOctreeContainer # Create a balanced octree by a random walk that recursively # subdivides -def create_fake_octree(RAMSESOctreeContainer oct_handler, +def create_fake_octree(SparseOctreeContainer oct_handler, long max_noct, long max_level, np.ndarray[np.int32_t, ndim=1] ndd, @@ -44,7 +44,7 @@ def create_fake_octree(RAMSESOctreeContainer oct_handler, return cur_leaf -cdef long subdivide(RAMSESOctreeContainer oct_handler, +cdef long subdivide(SparseOctreeContainer oct_handler, Oct *parent, int ind[3], int dd[3], long cur_leaf, long cur_level, diff --git a/yt/geometry/oct_container.pxd b/yt/geometry/oct_container.pxd index 4f7c482f075..3a3571464c4 100644 --- a/yt/geometry/oct_container.pxd +++ b/yt/geometry/oct_container.pxd @@ -94,8 +94,6 @@ cdef class SparseOctreeContainer(OctreeContainer): cdef void key_to_ipos(self, np.int64_t key, np.int64_t pos[3]) cdef np.int64_t ipos_to_key(self, int pos[3]) nogil -cdef class RAMSESOctreeContainer(SparseOctreeContainer): - pass cdef extern from "tsearch.h" nogil: void *tsearch(const void *key, void **rootp, diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 4f077a13712..e4695918370 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -210,7 +210,7 @@ cdef class OctreeContainer: ind32[i] = ind[i] self.get_root(ind32, &next) # We want to stop recursing when there's nowhere else to go - while next != NULL and level <= max_level: + while next != NULL and level < max_level: level += 1 for i in range(3): ipos[i] = (ipos[i] << 1) + ind[i] @@ -888,9 +888,6 @@ cdef class SparseOctreeContainer(OctreeContainer): # called. if self.root_nodes != NULL: free(self.root_nodes) -cdef class RAMSESOctreeContainer(SparseOctreeContainer): - pass - cdef class ARTOctreeContainer(OctreeContainer): def __init__(self, oct_domain_dimensions, domain_left_edge, domain_right_edge, partial_coverage = 0, diff --git a/yt/geometry/ramses_oct_container.pxd b/yt/geometry/ramses_oct_container.pxd new file mode 100644 index 00000000000..f399e87e75e --- /dev/null +++ b/yt/geometry/ramses_oct_container.pxd @@ -0,0 +1,15 @@ +""" +RAMSES Oct definitions file + + + + +""" +from oct_container cimport SparseOctreeContainer, OctInfo +from .oct_visitors cimport OctVisitor, Oct, cind +from yt.utilities.lib.fp_utils cimport * +cimport numpy as np + +cdef class RAMSESOctreeContainer(SparseOctreeContainer): + cdef Oct neighbor_in_direction(self, OctInfo *oinfo, np.int64_t *nneighbors, + Oct *o, bint periodicity[3]) \ No newline at end of file diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx new file mode 100644 index 00000000000..ac7aec06495 --- /dev/null +++ b/yt/geometry/ramses_oct_container.pyx @@ -0,0 +1,139 @@ +cimport cython +cimport oct_visitors +from selection_routines cimport SelectorObject, AlwaysSelector +cimport numpy as np +import numpy as np + +cdef class FillFileIndices(oct_visitors.OctVisitor): + cdef np.int64_t[:,:,:,:] cell_inds + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + if selected == 0: return + self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index + self.index += 1 + +cdef class NeighborVisitor(oct_visitors.OctVisitor): + cdef np.int64_t[:,:,:,:] cell_inds + cdef np.int64_t[:,:,:,:] neigh_cell_inds + cdef int idim # 0,1,2 for x,y,z + cdef int direction # +1 for +x, -1 for -x + cdef np.uint8_t neigh_ind[3] + cdef RAMSESOctreeContainer octree + cdef OctInfo oi + cdef Oct *neighbour + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void set_neighbour_oct(self): + cdef int i + cdef np.float64_t c, dx + cdef np.int64_t ipos + cdef np.float64_t fcoords[3] + cdef Oct *neighbour + dx = 1.0 / ((1 << self.oref) << self.level) + # Compute position of neighbouring cell + for i in range(3): + c = ((self.pos[i] << self.oref) + self.ind[i]) + if i == self.idim: + fcoords[i] = (c + 0.5 + self.direction) * dx / self.octree.nn[i] + else: + fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] + + # Use octree to find neighbour + neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) + + # Extra step - compute cell position in neighbouring oct (and store in oi.ipos) + if self.oi.level == self.level - 1: + for i in range(3): + ipos = (((self.pos[i] << self.oref) + self.ind[i])) >> 1 + if i == self.idim: + ipos += self.direction + #print('oi.level=%s level=%s oi.ipos[%s]<<2=%s ipos=%s' % ( + # self.oi.level, self.level, i, self.oi.ipos[i]<<2, ipos)) + if (self.oi.ipos[i] << 1) == ipos: + self.oi.ipos[i] = 0 + else: + self.oi.ipos[i] = 1 + self.neighbour = neighbour + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef int i + cdef np.int64_t cell_ind + cdef bint other_oct # True if the neighbouring cell lies in another oct + + # Note that we provide an index even if the cell is not selected. + if selected == 0: return + # Index of neighbouring cell within its oct + for i in range(3): + if i == self.idim: + self.neigh_ind[i] = (self.ind[i] + self.direction) + other_oct = self.neigh_ind[i] != 0 and self.neigh_ind[i] != 1 + if other_oct: + self.neigh_ind[i] %= 2 + else: + self.neigh_ind[i] = self.ind[i] + + if not other_oct: + # Simple case: the neighbouring cell is within the oct + cell_ind = self.cell_inds[o.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] + else: + # Complicated case: the cell is in a neighbouring oct + if self.last != o.domain_ind: + self.set_neighbour_oct() + self.last = o.domain_ind + + if self.neighbour != NULL: + if self.oi.level == self.level -1: + # Need to find cell position in neighbouring oct + for i in range(3): + self.neigh_ind[i] = self.oi.ipos[i] + elif self.oi.level != self.level: + print('FUUUUUCK %s %s' % (self.oi.level, self.level)) + cell_ind = self.cell_inds[self.neighbour.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] + else: + cell_ind = -1 + self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind + + + +cdef class RAMSESOctreeContainer(SparseOctreeContainer): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + cdef Oct neighbor_in_direction(self, OctInfo *oi, np.int64_t *nneighbors, Oct *o, + bint periodicity[3]): + pass + + def neighbors_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): + """Return index on file of all neighbors in a given direction""" + cdef SelectorObject always_selector = AlwaysSelector(None) + cdef FillFileIndices visitor + + cdef int num_cells = selector.count_oct_cells(self, -1) + + # Get the on-file index of each cell + cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((num_cells//8, 2, 2, 2), dtype="int64") + cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) + visitor = FillFileIndices(self, -1) + visitor.cell_inds = cell_inds + + self.visit_all_octs(selector, visitor) + + # Revisit the tree, now querying the neighbour in a given direction + cdef NeighborVisitor n_visitor + n_visitor = NeighborVisitor(self, -1) + n_visitor.idim = idim + n_visitor.direction = -direction + n_visitor.cell_inds = cell_inds + n_visitor.neigh_cell_inds = neigh_cell_inds + n_visitor.octree = self + n_visitor.last = -1 + self.visit_all_octs(always_selector, n_visitor) + + return np.asarray(cell_inds), np.asarray(neigh_cell_inds) From 6ce041687566ea9e09df93a8f7120825052769c7 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 22 Jan 2020 09:34:57 +0000 Subject: [PATCH 02/61] More generic code --- yt/geometry/ramses_oct_container.pxd | 4 +- yt/geometry/ramses_oct_container.pyx | 75 ++++++++++++++++++++-------- 2 files changed, 56 insertions(+), 23 deletions(-) diff --git a/yt/geometry/ramses_oct_container.pxd b/yt/geometry/ramses_oct_container.pxd index f399e87e75e..6d5dddd8997 100644 --- a/yt/geometry/ramses_oct_container.pxd +++ b/yt/geometry/ramses_oct_container.pxd @@ -11,5 +11,5 @@ from yt.utilities.lib.fp_utils cimport * cimport numpy as np cdef class RAMSESOctreeContainer(SparseOctreeContainer): - cdef Oct neighbor_in_direction(self, OctInfo *oinfo, np.int64_t *nneighbors, - Oct *o, bint periodicity[3]) \ No newline at end of file + cdef Oct neighbour_in_direction(self, OctInfo *oinfo, np.int64_t *nneighbors, + Oct *o, bint periodicity[3]) diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index ac7aec06495..2cf71d7cbbc 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -14,12 +14,11 @@ cdef class FillFileIndices(oct_visitors.OctVisitor): self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index self.index += 1 -cdef class NeighborVisitor(oct_visitors.OctVisitor): - cdef np.int64_t[:,:,:,:] cell_inds - cdef np.int64_t[:,:,:,:] neigh_cell_inds +cdef class BaseNeighbourVisitor(oct_visitors.OctVisitor): cdef int idim # 0,1,2 for x,y,z cdef int direction # +1 for +x, -1 for -x cdef np.uint8_t neigh_ind[3] + cdef np.int64_t[:,:,:,:] cell_inds cdef RAMSESOctreeContainer octree cdef OctInfo oi cdef Oct *neighbour @@ -51,8 +50,6 @@ cdef class NeighborVisitor(oct_visitors.OctVisitor): ipos = (((self.pos[i] << self.oref) + self.ind[i])) >> 1 if i == self.idim: ipos += self.direction - #print('oi.level=%s level=%s oi.ipos[%s]<<2=%s ipos=%s' % ( - # self.oi.level, self.level, i, self.oi.ipos[i]<<2, ipos)) if (self.oi.ipos[i] << 1) == ipos: self.oi.ipos[i] = 0 else: @@ -62,13 +59,13 @@ cdef class NeighborVisitor(oct_visitors.OctVisitor): @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): + cdef np.int64_t get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): cdef int i cdef np.int64_t cell_ind cdef bint other_oct # True if the neighbouring cell lies in another oct # Note that we provide an index even if the cell is not selected. - if selected == 0: return + # if selected == 0: return -1 # Index of neighbouring cell within its oct for i in range(3): if i == self.idim: @@ -89,47 +86,83 @@ cdef class NeighborVisitor(oct_visitors.OctVisitor): self.last = o.domain_ind if self.neighbour != NULL: - if self.oi.level == self.level -1: - # Need to find cell position in neighbouring oct + if self.oi.level == self.level - 1: + # Position within neighbouring oct is stored in oi.ipos for i in range(3): self.neigh_ind[i] = self.oi.ipos[i] elif self.oi.level != self.level: - print('FUUUUUCK %s %s' % (self.oi.level, self.level)) + print('This should not happen! %s %s' % (self.oi.level, self.level)) + return -1 cell_ind = self.cell_inds[self.neighbour.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] else: cell_ind = -1 + return cell_ind + + cdef inline np.uint8_t neighbour_rind(self): + cdef int d = (1 << self.oref) + return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) + +cdef class NeighbourVisitor(BaseNeighbourVisitor): + cdef np.int64_t[:,:,:,:] neigh_cell_inds + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef np.int64_t cell_ind + cell_ind = self.get_neighbour_cell_index(o, selected) self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind - +cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): + cdef np.uint8_t[:] shifted_levels + cdef np.int64_t[:] shifted_file_inds + cdef np.uint8_t[:] shifted_cell_inds + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef np.int64_t neigbour_cell_index + if selected == 0: return + # Note: only selected items have an index + neighbour_index = self.get_neighbour_cell_index(o, selected) + self.shifted_levels[self.index] = self.oi.level + self.shifted_file_inds[self.index] = self.neighbour.file_ind + self.shifted_cell_inds[self.index] = self.neighbour_rind() + self.index += 1 + # if neighbour_index > -1: + # self.shifted_levels[neighbour_index] = self.level + # self.shifted_file_inds[neighbour_index] = o.file_ind + # self.shifted_cell_inds[neighbour_index] = self.rind() cdef class RAMSESOctreeContainer(SparseOctreeContainer): @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) - cdef Oct neighbor_in_direction(self, OctInfo *oi, np.int64_t *nneighbors, Oct *o, + cdef Oct neighbour_in_direction(self, OctInfo *oi, np.int64_t *nneighbours, Oct *o, bint periodicity[3]): pass - def neighbors_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): - """Return index on file of all neighbors in a given direction""" + def neighbours_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): + """Return index on file of all neighbours in a given direction""" cdef SelectorObject always_selector = AlwaysSelector(None) - cdef FillFileIndices visitor cdef int num_cells = selector.count_oct_cells(self, -1) # Get the on-file index of each cell - cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((num_cells//8, 2, 2, 2), dtype="int64") - cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) + cdef FillFileIndices visitor + cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") + visitor = FillFileIndices(self, -1) visitor.cell_inds = cell_inds self.visit_all_octs(selector, visitor) - # Revisit the tree, now querying the neighbour in a given direction - cdef NeighborVisitor n_visitor - n_visitor = NeighborVisitor(self, -1) + # Store the index of the neighbour + cdef NeighbourVisitor n_visitor + cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) + n_visitor = NeighbourVisitor(self, -1) n_visitor.idim = idim - n_visitor.direction = -direction + n_visitor.direction = direction n_visitor.cell_inds = cell_inds n_visitor.neigh_cell_inds = neigh_cell_inds n_visitor.octree = self From 2b0baa360e7c5838d5ad36c8cccc119190b5dc19 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 22 Jan 2020 14:47:50 +0000 Subject: [PATCH 03/61] Remove useless variable on oct container --- yt/geometry/oct_container.pyx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index e4695918370..544c88af7c3 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -728,14 +728,13 @@ cdef class OctreeContainer: cdef np.ndarray[np.float64_t, ndim=2] source cdef np.ndarray[np.float64_t, ndim=1] dest cdef int i - cdef np.int64_t local_filled = 0 + cdef str key for key in dest_fields: dest = dest_fields[key] source = source_fields[key] for i in range(levels.shape[0]): if levels[i] != level: continue dest[i + offset] = source[file_inds[i], cell_inds[i]] - local_filled += 1 def finalize(self): cdef SelectorObject selector = selection_routines.AlwaysSelector(None) From 8cf8a565564a1ac31c1c12e6a82101922b25e37d Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 22 Jan 2020 14:48:59 +0000 Subject: [PATCH 04/61] First working (not crashing version) --- yt/geometry/ramses_oct_container.pyx | 85 +++++++++++++++++++++++++--- 1 file changed, 77 insertions(+), 8 deletions(-) diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index 2cf71d7cbbc..f6b6e98ef55 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -13,7 +13,7 @@ cdef class FillFileIndices(oct_visitors.OctVisitor): if selected == 0: return self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index self.index += 1 - + cdef class BaseNeighbourVisitor(oct_visitors.OctVisitor): cdef int idim # 0,1,2 for x,y,z cdef int direction # +1 for +x, -1 for -x @@ -125,14 +125,15 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): if selected == 0: return # Note: only selected items have an index neighbour_index = self.get_neighbour_cell_index(o, selected) - self.shifted_levels[self.index] = self.oi.level - self.shifted_file_inds[self.index] = self.neighbour.file_ind - self.shifted_cell_inds[self.index] = self.neighbour_rind() + self.shifted_levels[self.index] = self.level + if self.neighbour != NULL: + # Note: we store the local level, not the remote one + self.shifted_file_inds[self.index] = self.neighbour.file_ind + self.shifted_cell_inds[self.index] = self.neighbour_rind() + else: + self.shifted_file_inds[self.index] = -1 + self.shifted_cell_inds[self.index] = 255 # -1 on uint8 self.index += 1 - # if neighbour_index > -1: - # self.shifted_levels[neighbour_index] = self.level - # self.shifted_file_inds[neighbour_index] = o.file_ind - # self.shifted_cell_inds[neighbour_index] = self.rind() cdef class RAMSESOctreeContainer(SparseOctreeContainer): @cython.boundscheck(False) @@ -170,3 +171,71 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): self.visit_all_octs(always_selector, n_visitor) return np.asarray(cell_inds), np.asarray(neigh_cell_inds) + + def file_index_octs_with_shift(self, SelectorObject selector, int domain_id, + int idim, int direction, int num_cells = -1): + """Return index on file of all neighbours in a given direction""" + # We create oct arrays of the correct size + cdef np.int64_t i + if num_cells < 0: + num_cells = selector.count_oct_cells(self, domain_id) + + # TODO remove this requirement + cdef np.ndarray[np.int64_t, ndim=4] cell_inds + cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") - 1 + + # Fill value of each cell with its neighbouring value + cdef FillFileIndicesRNeighbour neigh_visitor + cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels + cdef np.ndarray[np.uint8_t, ndim=1] shifted_cell_inds + cdef np.ndarray[np.int64_t, ndim=1] shifted_file_inds + shifted_levels = np.zeros(num_cells, dtype="uint8") + shifted_file_inds = np.zeros(num_cells, dtype="int64") + shifted_cell_inds = np.zeros(num_cells, dtype="uint8") + + if self.fill_style == "r": + neigh_visitor = FillFileIndicesRNeighbour(self, domain_id) + # input: index of neighbouring cells + neigh_visitor.cell_inds = cell_inds + # output: level, file_ind and cell_ind of the neighbouring cells + neigh_visitor.shifted_levels = shifted_levels + neigh_visitor.shifted_file_inds = shifted_file_inds + neigh_visitor.shifted_cell_inds = shifted_cell_inds + # direction to explore and extra parameters of the visitor + neigh_visitor.idim = idim + neigh_visitor.direction = direction + neigh_visitor.octree = self + neigh_visitor.last = -1 + elif self.fill_style == "o": + raise NotImplementedError('C-style filling with spatial offset has not been implemented.') + else: + raise RuntimeError + self.visit_all_octs(selector, neigh_visitor) + return shifted_levels, shifted_cell_inds, shifted_file_inds + + def file_index_octs(self, SelectorObject selector, int domain_id, + num_cells = -1, spatial_offset=(0, 0, 0)): + + + cdef int i, idim, direction + cdef bint do_spatial_offset + cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds + cdef np.ndarray source_shifted + + do_spatial_offset = False + for i in range(3): + if spatial_offset[i] == 1 or spatial_offset[i] == -1: + idim = i + direction = spatial_offset[i] + if do_spatial_offset: + raise Exception( + 'ERROR: You can only specify one spatial offset direction, got [%s, %s, %s]!' % + (spatial_offset[0], spatial_offset[1], spatial_offset[2])) + do_spatial_offset = True + + if not do_spatial_offset: + return super(RAMSESOctreeContainer, self).file_index_octs( + selector, domain_id, num_cells) + else: + return self.file_index_octs_with_shift( + selector, domain_id, idim, direction, num_cells) From ab4410eac4c787dfd91d92cf9a4404c334167991 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 09:16:35 +0100 Subject: [PATCH 05/61] Does not work yet, but hey it compiles --- yt/geometry/oct_container.pxd | 8 +- yt/geometry/oct_visitors.pxd | 32 +++++ yt/geometry/oct_visitors.pyx | 119 +++++++++++++++- yt/geometry/ramses_oct_container.pyx | 195 +++++---------------------- 4 files changed, 186 insertions(+), 168 deletions(-) diff --git a/yt/geometry/oct_container.pxd b/yt/geometry/oct_container.pxd index 3a3571464c4..eb500d7fd56 100644 --- a/yt/geometry/oct_container.pxd +++ b/yt/geometry/oct_container.pxd @@ -12,7 +12,7 @@ cimport numpy as np from yt.utilities.lib.fp_utils cimport * cimport oct_visitors cimport selection_routines -from .oct_visitors cimport OctVisitor, Oct, cind +from .oct_visitors cimport OctVisitor, Oct, cind, OctInfo from libc.stdlib cimport bsearch, qsort, realloc, malloc, free from libc.math cimport floor from yt.utilities.lib.allocation_container cimport \ @@ -27,12 +27,6 @@ cdef struct OctKey: np.int64_t *indices np.int64_t pcount -cdef struct OctInfo: - np.float64_t left_edge[3] - np.float64_t dds[3] - np.int64_t ipos[3] - np.int32_t level - cdef struct OctList cdef struct OctList: diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 125c69a523e..89dec648108 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -17,6 +17,12 @@ cdef struct Oct: np.int64_t domain # (opt) addl int index Oct **children # Up to 8 long +cdef struct OctInfo: + np.float64_t left_edge[3] + np.float64_t dds[3] + np.int64_t ipos[3] + np.int32_t level + cdef struct OctPadded: np.int64_t file_ind np.int64_t domain_ind @@ -136,3 +142,29 @@ cdef inline int cind(int i, int j, int k) nogil: # THIS ONLY WORKS FOR CHILDREN. It is not general for zones. return (((i*2)+j)*2+k) +from oct_container cimport OctreeContainer + +# cimport oct_container +cdef class BaseNeighbourVisitor(OctVisitor): + cdef int idim # 0,1,2 for x,y,z + cdef int direction # +1 for +x, -1 for -x + cdef np.int8_t[:] neigh_ind + cdef Oct *neighbour + cdef OctreeContainer octree + cdef OctInfo oi + + cdef void set_neighbour_oct(self) + cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected) + + cdef inline np.uint8_t neighbour_rind(self): + cdef int d = (1 << self.oref) + return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) + +# # # cdef class NeighbourVisitor(BaseNeighbourVisitor): +# # # cdef np.int64_t[:,:,:,:] neigh_cell_inds + +cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): + cdef np.uint8_t[:] shifted_levels + cdef np.int64_t[:] shifted_file_inds + cdef np.uint8_t[:] shifted_cell_inds + diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index f8212b6dff8..ce754f40f80 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -12,7 +12,7 @@ cimport numpy import numpy from yt.utilities.lib.fp_utils cimport * from libc.stdlib cimport malloc, free -from yt.geometry.oct_container cimport OctreeContainer +from yt.geometry.oct_container cimport OctreeContainer, OctInfo from yt.utilities.lib.geometry_utils cimport encode_morton_64bit # Now some visitor functions @@ -336,3 +336,120 @@ cdef class MortonIndexOcts(OctVisitor): np.uint64(coord[1]), np.uint64(coord[2])) self.index += 1 + + +cdef class BaseNeighbourVisitor(OctVisitor): + # cdef OctInfo oi + # cdef OctreeContainer octree + + def __init__(self, OctreeContainer octree, int domain_id = -1): + self.octree = octree + super(BaseNeighbourVisitor, self).__init__(octree, domain_id) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void set_neighbour_oct(self): + cdef int i + cdef np.float64_t c, dx + cdef np.int64_t ipos + cdef np.float64_t fcoords[3] + cdef Oct *neighbour + dx = 1.0 / ((1 << self.oref) << self.level) + # Compute position of neighbouring cell + for i in range(3): + c = ((self.pos[i] << self.oref) + self.ind[i]) + if i == self.idim: + fcoords[i] = (c + 0.5 + self.direction) * dx / self.octree.nn[i] + else: + fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] + + # Use octree to find neighbour + neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) + + # Extra step - compute cell position in neighbouring oct (and store in oi.ipos) + if self.oi.level == self.level - 1: + for i in range(3): + ipos = (((self.pos[i] << self.oref) + self.ind[i])) >> 1 + if i == self.idim: + ipos += self.direction + if (self.oi.ipos[i] << 1) == ipos: + self.oi.ipos[i] = 0 + else: + self.oi.ipos[i] = 1 + self.neighbour = neighbour + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): + cdef int i + # cdef np.int64_t cell_ind + cdef bint other_oct # True if the neighbouring cell lies in another oct + + # Note that we provide an index even if the cell is not selected. + # if selected == 0: return -1 + # Index of neighbouring cell within its oct + for i in range(3): + if i == self.idim: + self.neigh_ind[i] = (self.ind[i] + self.direction) + other_oct = self.neigh_ind[i] < 0 or self.neigh_ind[i] > 1 + if other_oct: + self.neigh_ind[i] %= 2 + else: + self.neigh_ind[i] = self.ind[i] + + if not other_oct: + # Simple case: the neighbouring cell is within the oct + # cell_ind = self.cell_inds[o.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] + self.neighbour = o + else: + # Complicated case: the cell is in a neighbouring oct + if self.last != o.domain_ind: + self.set_neighbour_oct() + self.last = o.domain_ind + + if self.neighbour != NULL: + if self.oi.level == self.level - 1: + # Position within neighbouring oct is stored in oi.ipos + for i in range(3): + self.neigh_ind[i] = self.oi.ipos[i] + if not ( 0<= self.oi.ipos[i] <= 1): + print('WHAAAAT?!', self.oi.ipos[i]) + elif self.oi.level != self.level: + print('This should not happen! %s %s' % (self.oi.level, self.level)) + return # -1 + # cell_ind = self.cell_inds[self.neighbour.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] + # else: + # cell_ind = -1 + # return cell_ind + +# cdef class NeighbourVisitor(BaseNeighbourVisitor): +# @cython.boundscheck(False) +# @cython.wraparound(False) +# @cython.initializedcheck(False) +# cdef void visit(self, Oct* o, np.uint8_t selected): +# cdef np.int64_t cell_ind +# # cell_ind = self.get_neighbour_cell_index(o, selected) +# cell_ind = 0 +# self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind + +cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef np.int64_t neigbour_cell_index + if selected == 0: return + # Note: only selected items have an index + neighbour_index = self.get_neighbour_cell_index(o, selected) + self.shifted_levels[self.index] = self.level + if self.neighbour != NULL: + # Note: we store the local level, not the remote one + self.shifted_file_inds[self.index] = self.neighbour.file_ind + self.shifted_cell_inds[self.index] = self.neighbour_rind() + else: + self.shifted_file_inds[self.index] = -1 + self.shifted_cell_inds[self.index] = 255 # -1 on uint8 + self.index += 1 \ No newline at end of file diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index f6b6e98ef55..edf6e60b8d0 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -1,139 +1,18 @@ cimport cython -cimport oct_visitors +from oct_visitors cimport FillFileIndicesRNeighbour from selection_routines cimport SelectorObject, AlwaysSelector cimport numpy as np import numpy as np -cdef class FillFileIndices(oct_visitors.OctVisitor): - cdef np.int64_t[:,:,:,:] cell_inds - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - if selected == 0: return - self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index - self.index += 1 - -cdef class BaseNeighbourVisitor(oct_visitors.OctVisitor): - cdef int idim # 0,1,2 for x,y,z - cdef int direction # +1 for +x, -1 for -x - cdef np.uint8_t neigh_ind[3] - cdef np.int64_t[:,:,:,:] cell_inds - cdef RAMSESOctreeContainer octree - cdef OctInfo oi - cdef Oct *neighbour - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void set_neighbour_oct(self): - cdef int i - cdef np.float64_t c, dx - cdef np.int64_t ipos - cdef np.float64_t fcoords[3] - cdef Oct *neighbour - dx = 1.0 / ((1 << self.oref) << self.level) - # Compute position of neighbouring cell - for i in range(3): - c = ((self.pos[i] << self.oref) + self.ind[i]) - if i == self.idim: - fcoords[i] = (c + 0.5 + self.direction) * dx / self.octree.nn[i] - else: - fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] - - # Use octree to find neighbour - neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) - - # Extra step - compute cell position in neighbouring oct (and store in oi.ipos) - if self.oi.level == self.level - 1: - for i in range(3): - ipos = (((self.pos[i] << self.oref) + self.ind[i])) >> 1 - if i == self.idim: - ipos += self.direction - if (self.oi.ipos[i] << 1) == ipos: - self.oi.ipos[i] = 0 - else: - self.oi.ipos[i] = 1 - self.neighbour = neighbour - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef np.int64_t get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): - cdef int i - cdef np.int64_t cell_ind - cdef bint other_oct # True if the neighbouring cell lies in another oct - - # Note that we provide an index even if the cell is not selected. - # if selected == 0: return -1 - # Index of neighbouring cell within its oct - for i in range(3): - if i == self.idim: - self.neigh_ind[i] = (self.ind[i] + self.direction) - other_oct = self.neigh_ind[i] != 0 and self.neigh_ind[i] != 1 - if other_oct: - self.neigh_ind[i] %= 2 - else: - self.neigh_ind[i] = self.ind[i] - - if not other_oct: - # Simple case: the neighbouring cell is within the oct - cell_ind = self.cell_inds[o.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] - else: - # Complicated case: the cell is in a neighbouring oct - if self.last != o.domain_ind: - self.set_neighbour_oct() - self.last = o.domain_ind - - if self.neighbour != NULL: - if self.oi.level == self.level - 1: - # Position within neighbouring oct is stored in oi.ipos - for i in range(3): - self.neigh_ind[i] = self.oi.ipos[i] - elif self.oi.level != self.level: - print('This should not happen! %s %s' % (self.oi.level, self.level)) - return -1 - cell_ind = self.cell_inds[self.neighbour.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] - else: - cell_ind = -1 - return cell_ind - - cdef inline np.uint8_t neighbour_rind(self): - cdef int d = (1 << self.oref) - return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) - -cdef class NeighbourVisitor(BaseNeighbourVisitor): - cdef np.int64_t[:,:,:,:] neigh_cell_inds - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - cdef np.int64_t cell_ind - cell_ind = self.get_neighbour_cell_index(o, selected) - self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind - -cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): - cdef np.uint8_t[:] shifted_levels - cdef np.int64_t[:] shifted_file_inds - cdef np.uint8_t[:] shifted_cell_inds - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - cdef np.int64_t neigbour_cell_index - if selected == 0: return - # Note: only selected items have an index - neighbour_index = self.get_neighbour_cell_index(o, selected) - self.shifted_levels[self.index] = self.level - if self.neighbour != NULL: - # Note: we store the local level, not the remote one - self.shifted_file_inds[self.index] = self.neighbour.file_ind - self.shifted_cell_inds[self.index] = self.neighbour_rind() - else: - self.shifted_file_inds[self.index] = -1 - self.shifted_cell_inds[self.index] = 255 # -1 on uint8 - self.index += 1 +# cdef class FillFileIndices(oct_visitors.OctVisitor): +# cdef np.int64_t[:,:,:,:] cell_inds +# @cython.boundscheck(False) +# @cython.wraparound(False) +# @cython.initializedcheck(False) +# cdef void visit(self, Oct* o, np.uint8_t selected): +# if selected == 0: return +# self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index +# self.index += 1 cdef class RAMSESOctreeContainer(SparseOctreeContainer): @cython.boundscheck(False) @@ -143,34 +22,34 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): bint periodicity[3]): pass - def neighbours_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): - """Return index on file of all neighbours in a given direction""" - cdef SelectorObject always_selector = AlwaysSelector(None) + # def neighbours_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): + # """Return index on file of all neighbours in a given direction""" + # cdef SelectorObject always_selector = AlwaysSelector(None) - cdef int num_cells = selector.count_oct_cells(self, -1) + # cdef int num_cells = selector.count_oct_cells(self, -1) - # Get the on-file index of each cell - cdef FillFileIndices visitor - cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") + # # Get the on-file index of each cell + # cdef FillFileIndices visitor + # cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") - visitor = FillFileIndices(self, -1) - visitor.cell_inds = cell_inds + # visitor = FillFileIndices(self, -1) + # visitor.cell_inds = cell_inds - self.visit_all_octs(selector, visitor) + # self.visit_all_octs(selector, visitor) + + # # Store the index of the neighbour + # cdef NeighbourVisitor n_visitor + # cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) + # n_visitor = NeighbourVisitor(self, -1) + # n_visitor.idim = idim + # n_visitor.direction = direction + # n_visitor.cell_inds = cell_inds + # n_visitor.neigh_cell_inds = neigh_cell_inds + # n_visitor.octree = self + # n_visitor.last = -1 + # self.visit_all_octs(always_selector, n_visitor) - # Store the index of the neighbour - cdef NeighbourVisitor n_visitor - cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) - n_visitor = NeighbourVisitor(self, -1) - n_visitor.idim = idim - n_visitor.direction = direction - n_visitor.cell_inds = cell_inds - n_visitor.neigh_cell_inds = neigh_cell_inds - n_visitor.octree = self - n_visitor.last = -1 - self.visit_all_octs(always_selector, n_visitor) - - return np.asarray(cell_inds), np.asarray(neigh_cell_inds) + # return np.asarray(cell_inds), np.asarray(neigh_cell_inds) def file_index_octs_with_shift(self, SelectorObject selector, int domain_id, int idim, int direction, int num_cells = -1): @@ -180,10 +59,6 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): if num_cells < 0: num_cells = selector.count_oct_cells(self, domain_id) - # TODO remove this requirement - cdef np.ndarray[np.int64_t, ndim=4] cell_inds - cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") - 1 - # Fill value of each cell with its neighbouring value cdef FillFileIndicesRNeighbour neigh_visitor cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels @@ -195,8 +70,6 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): if self.fill_style == "r": neigh_visitor = FillFileIndicesRNeighbour(self, domain_id) - # input: index of neighbouring cells - neigh_visitor.cell_inds = cell_inds # output: level, file_ind and cell_ind of the neighbouring cells neigh_visitor.shifted_levels = shifted_levels neigh_visitor.shifted_file_inds = shifted_file_inds @@ -210,7 +83,9 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): raise NotImplementedError('C-style filling with spatial offset has not been implemented.') else: raise RuntimeError + print('visiting all octs') self.visit_all_octs(selector, neigh_visitor) + print('visited') return shifted_levels, shifted_cell_inds, shifted_file_inds def file_index_octs(self, SelectorObject selector, int domain_id, From b7c0c2f362c58c1342be1cca22f04d33cf34acde Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 22 Jan 2020 16:38:10 +0000 Subject: [PATCH 06/61] Wiring cython code with ghost zones --- yt/data_objects/data_containers.py | 12 ++--- yt/data_objects/octree_subset.py | 3 +- yt/frontends/ramses/data_structures.py | 69 ++++++++++++++++++++++++-- 3 files changed, 73 insertions(+), 11 deletions(-) diff --git a/yt/data_objects/data_containers.py b/yt/data_objects/data_containers.py index 8ebfe3988fd..4b9270c60ac 100644 --- a/yt/data_objects/data_containers.py +++ b/yt/data_objects/data_containers.py @@ -352,13 +352,11 @@ def _generate_spatial_fluid(self, field, ngz): rv = self.ds.arr(np.empty(wogz.ires.size, dtype="float64"), units) outputs.append(rv) - if gz._type_name == 'octree_subset': - raise NotImplementedError - else: - ind += wogz.select( - self.selector, - gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz], - rv, ind) + data = gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz] + ind += wogz.select( + self.selector, + data, + rv, ind) if accumulate: rv = uconcatenate(outputs) return rv diff --git a/yt/data_objects/octree_subset.py b/yt/data_objects/octree_subset.py index dbe297e5d0c..f98876c5e34 100644 --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -37,9 +37,10 @@ class OctreeSubset(YTSelectionContainer): _cell_count = -1 _block_reorder = None - def __init__(self, base_region, domain, ds, over_refine_factor = 1): + def __init__(self, base_region, domain, ds, over_refine_factor = 1, num_ghost_zones = 0): super(OctreeSubset, self).__init__(ds, None) self._num_zones = 1 << (over_refine_factor) + self._num_ghost_zones = num_ghost_zones self._oref = over_refine_factor self.domain = domain self.domain_id = domain.domain_id diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 2335987c80e..c3209b42948 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -182,7 +182,9 @@ class RAMSESDomainSubset(OctreeSubset): _domain_offset = 1 _block_reorder = "F" - def fill(self, fd, fields, selector, file_handler): + _base_grid = None + + def _fill_no_ghostzones(self, fd, fields, selector, file_handler): ndim = self.ds.dimensionality # Here we get a copy of the file, which we skip through and read the # bits we want. @@ -198,13 +200,73 @@ def fill(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, 'float64') - fill_hydro(fd, file_handler.offset, file_handler.level_count, levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler) return tr + def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zones): + ndim = self.ds.dimensionality + # Here we get a copy of the file, which we skip through and read the + # bits we want. + oct_handler = self.oct_handler + all_fields = [f for ft, f in file_handler.field_list] + fields = [f for ft, f in fields] + + oct_count = selector.count_octs(self.oct_handler, self.domain_id) + cell_count = oct_count * self._num_zones**ndim + iwidth = self._num_zones + num_ghost_zones * 2 + + tr = {} + tr_all = {} + for field in fields: + tr_all[field] = np.zeros((oct_count, iwidth, iwidth, iwidth), 'float64') + tr[field] = np.zeros(cell_count, 'float64') + + # Compute the index to read with a positive and negative shift in all dimensions + for idim in range(ndim): + ishift_all = [num_ghost_zones]*ndim + for shift in range(-num_ghost_zones, num_ghost_zones+1, 2): + ishift_all[idim] = num_ghost_zones + shift + import yt + new_selector = yt.geometry.selection_routines.OctreeSubsetSelector(self) + if shift == 0: + continue + levels, cell_inds, file_inds = self.oct_handler.file_index_octs_with_shift( + new_selector, self.domain_id, cell_count, idim, shift) + + # Initializing data container + for field in fields: + tr[field][:] = 0 + + fill_hydro(fd, file_handler.offset, + file_handler.level_count, levels, cell_inds, + file_inds, ndim, all_fields, fields, tr, + oct_handler) + _slice = tuple( + [slice(None)] + + [slice(i, i+self._num_zones) for i in ishift_all]) + for field in fields: + tr_all[field][_slice] = \ + tr[field].reshape(oct_count, 2, 2, 2) + for field in fields: + tr_all[field] = tr_all[field].reshape(-1) + return tr_all + + def fill(self, fd, fields, selector, file_handler): + if self._num_ghost_zones == 0: + return self._fill_no_ghostzones(fd, fields, selector, file_handler) + else: + return self._fill_with_ghostzones(fd, fields, selector, file_handler, self._num_ghost_zones) + + def retrieve_ghost_zones(self, ngz, fields, smoothed=False): + new_subset = RAMSESDomainSubset(self.base_region, self.domain, self.ds, num_ghost_zones=ngz) + new_subset._base_grid = self + + return new_subset + + class RAMSESIndex(OctreeIndex): def __init__(self, ds, dataset_type='ramses'): @@ -255,13 +317,14 @@ def _detect_output_fields(self): self.field_list = self.particle_field_list + self.fluid_field_list def _identify_base_chunk(self, dobj): + ngz = dobj._num_ghost_zones if getattr(dobj, "_chunk_info", None) is None: domains = [dom for dom in self.domains if dom.included(dobj.selector)] base_region = getattr(dobj, "base_region", dobj) if len(domains) > 1: mylog.debug("Identified %s intersecting domains", len(domains)) - subsets = [RAMSESDomainSubset(base_region, domain, self.dataset) + subsets = [RAMSESDomainSubset(base_region, domain, self.dataset, num_ghost_zones=ngz) for domain in domains] dobj._chunk_info = subsets dobj._current_chunk = list(self._chunk_all(dobj))[0] From 8da0a17318b713548f8ecfb3f79bbbb772bdf22a Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 22 Jan 2020 20:43:13 +0000 Subject: [PATCH 07/61] Cleanup --- yt/geometry/oct_visitors.pxd | 3 ++- yt/geometry/oct_visitors.pyx | 30 ++++++------------------------ 2 files changed, 8 insertions(+), 25 deletions(-) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 89dec648108..18493856878 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -16,6 +16,7 @@ cdef struct Oct: np.int64_t domain_ind # index within the global set of domains np.int64_t domain # (opt) addl int index Oct **children # Up to 8 long + Oct *parent cdef struct OctInfo: np.float64_t left_edge[3] @@ -153,7 +154,7 @@ cdef class BaseNeighbourVisitor(OctVisitor): cdef OctreeContainer octree cdef OctInfo oi - cdef void set_neighbour_oct(self) + cdef void set_neighbour_oct(self, Oct* o) cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected) cdef inline np.uint8_t neighbour_rind(self): diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index ce754f40f80..3d4d639a977 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -8,8 +8,8 @@ Oct visitor functions cimport cython -cimport numpy -import numpy +cimport numpy as np +import numpy as np from yt.utilities.lib.fp_utils cimport * from libc.stdlib cimport malloc, free from yt.geometry.oct_container cimport OctreeContainer, OctInfo @@ -339,17 +339,15 @@ cdef class MortonIndexOcts(OctVisitor): cdef class BaseNeighbourVisitor(OctVisitor): - # cdef OctInfo oi - # cdef OctreeContainer octree - def __init__(self, OctreeContainer octree, int domain_id = -1): self.octree = octree + self.neigh_ind = np.zeros(3, np.int8) super(BaseNeighbourVisitor, self).__init__(octree, domain_id) @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) - cdef void set_neighbour_oct(self): + cdef void set_neighbour_oct(self, Oct *o): cdef int i cdef np.float64_t c, dx cdef np.int64_t ipos @@ -385,7 +383,6 @@ cdef class BaseNeighbourVisitor(OctVisitor): @cython.cdivision(True) cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): cdef int i - # cdef np.int64_t cell_ind cdef bint other_oct # True if the neighbouring cell lies in another oct # Note that we provide an index even if the cell is not selected. @@ -402,12 +399,11 @@ cdef class BaseNeighbourVisitor(OctVisitor): if not other_oct: # Simple case: the neighbouring cell is within the oct - # cell_ind = self.cell_inds[o.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] self.neighbour = o else: # Complicated case: the cell is in a neighbouring oct if self.last != o.domain_ind: - self.set_neighbour_oct() + self.set_neighbour_oct(o) self.last = o.domain_ind if self.neighbour != NULL: @@ -419,21 +415,7 @@ cdef class BaseNeighbourVisitor(OctVisitor): print('WHAAAAT?!', self.oi.ipos[i]) elif self.oi.level != self.level: print('This should not happen! %s %s' % (self.oi.level, self.level)) - return # -1 - # cell_ind = self.cell_inds[self.neighbour.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] - # else: - # cell_ind = -1 - # return cell_ind - -# cdef class NeighbourVisitor(BaseNeighbourVisitor): -# @cython.boundscheck(False) -# @cython.wraparound(False) -# @cython.initializedcheck(False) -# cdef void visit(self, Oct* o, np.uint8_t selected): -# cdef np.int64_t cell_ind -# # cell_ind = self.get_neighbour_cell_index(o, selected) -# cell_ind = 0 -# self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind + return cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): @cython.boundscheck(False) From adb21bab0433481ca9a28c76a29a7cff5047890e Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 23 Jan 2020 14:48:14 +0000 Subject: [PATCH 08/61] First working (not crashing) version --- yt/frontends/ramses/data_structures.py | 73 +++++++++++++++++++++++--- yt/geometry/oct_container.pyx | 9 +++- yt/geometry/oct_visitors.pxd | 9 +++- yt/geometry/oct_visitors.pyx | 66 +++++++++++++++++------ yt/geometry/ramses_oct_container.pyx | 67 ++++++++++++++--------- 5 files changed, 172 insertions(+), 52 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index c3209b42948..b50a3d76f28 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -10,6 +10,8 @@ setdefaultattr from yt.geometry.oct_geometry_handler import \ OctreeIndex +from yt.geometry.selection_routines import OctreeSubsetSelector + from yt.geometry.geometry_handler import \ YTDataChunk from yt.data_objects.static_output import \ @@ -182,7 +184,18 @@ class RAMSESDomainSubset(OctreeSubset): _domain_offset = 1 _block_reorder = "F" - _base_grid = None + _base_domain = None + + def __init__(self, base_region, domain, ds, over_refine_factor=1, num_ghost_zones=0, + base_grid=None): + super(RAMSESDomainSubset, self).__init__(base_region, domain, ds, over_refine_factor, num_ghost_zones) + + self._base_grid = base_grid + + if num_ghost_zones > 0: + # Create a base domain *with no self._base_domain.fwidth + base_domain = RAMSESDomainSubset(ds.all_data(), domain, ds, over_refine_factor) + self._base_domain = base_domain def _fill_no_ghostzones(self, fd, fields, selector, file_handler): ndim = self.ds.dimensionality @@ -221,7 +234,7 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo tr = {} tr_all = {} for field in fields: - tr_all[field] = np.zeros((oct_count, iwidth, iwidth, iwidth), 'float64') + tr_all[field] = np.full((oct_count, iwidth, iwidth, iwidth), np.nan, 'float64') tr[field] = np.zeros(cell_count, 'float64') # Compute the index to read with a positive and negative shift in all dimensions @@ -229,12 +242,10 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo ishift_all = [num_ghost_zones]*ndim for shift in range(-num_ghost_zones, num_ghost_zones+1, 2): ishift_all[idim] = num_ghost_zones + shift - import yt - new_selector = yt.geometry.selection_routines.OctreeSubsetSelector(self) if shift == 0: continue levels, cell_inds, file_inds = self.oct_handler.file_index_octs_with_shift( - new_selector, self.domain_id, cell_count, idim, shift) + selector, self.domain_id, idim, shift, cell_count) # Initializing data container for field in fields: @@ -254,6 +265,55 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo tr_all[field] = tr_all[field].reshape(-1) return tr_all + @property + def fwidth(self): + fwidth = super(RAMSESDomainSubset, self).fwidth + if self._num_ghost_zones > 0: + fwidth = super(RAMSESDomainSubset, self).fwidth.reshape(-1, 8, 3) + n_oct = fwidth.shape[0] + new_fwidth = np.zeros((n_oct, self.nz**3, 3), dtype=fwidth.dtype) + new_fwidth[:, :, :] = fwidth[:, 0:1, :] + fwidth = new_fwidth.reshape(-1, 3) + return fwidth + + @property + def fcoords(self): + fcoords = super(RAMSESDomainSubset, self).fcoords + num_ghost_zones = self._num_ghost_zones + if num_ghost_zones == 0: + return fcoords + + fcoords_base = self._base_domain.fcoords + oct_selector = OctreeSubsetSelector(self) + oh = self.oct_handler + + n_oct = fcoords_base.size // 3 // 8 + new_fcoords = np.full((n_oct, 4, 4, 4, 3), np.nan) + + icell = oh.fill_index(oct_selector) + Ncell = icell.size + + for idim in range(3): + for idir in (-1, 1): + ishift_all = [1, 1, 1] + ishift_all[idim] += idir + + nicell = oh.neighbours_in_direction(idim, idir, icell).reshape(-1) + + tmp = np.full((n_oct * 2 * 2 * 2, 3), np.nan) + + oh.copy_neighbour_data( + icell.reshape(-1), nicell, + fcoords_base, tmp, Ncell) + + _slice = tuple([slice(None)] + + [slice(i, i+2) for i in ishift_all] + + [slice(None)]) + new_fcoords[_slice] = tmp.reshape(n_oct, 2, 2, 2, -1) + new_fcoords = self.ds.arr(new_fcoords, fcoords_base.units) + return new_fcoords + + def fill(self, fd, fields, selector, file_handler): if self._num_ghost_zones == 0: return self._fill_no_ghostzones(fd, fields, selector, file_handler) @@ -261,8 +321,7 @@ def fill(self, fd, fields, selector, file_handler): return self._fill_with_ghostzones(fd, fields, selector, file_handler, self._num_ghost_zones) def retrieve_ghost_zones(self, ngz, fields, smoothed=False): - new_subset = RAMSESDomainSubset(self.base_region, self.domain, self.ds, num_ghost_zones=ngz) - new_subset._base_grid = self + new_subset = RAMSESDomainSubset(self.base_region, self.domain, self.ds, num_ghost_zones=ngz, base_grid=self) return new_subset diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 544c88af7c3..73d349547f7 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -566,6 +566,7 @@ cdef class OctreeContainer: cdef int ind[3] cdef int nb = 0 cdef Oct *cur + cdef Oct *parent cdef np.float64_t pp[3] cdef np.float64_t cp[3] cdef np.float64_t dds[3] @@ -574,6 +575,7 @@ cdef class OctreeContainer: cdef OctAllocationContainer *cont = self.domains.get_cont(curdom - 1) cdef int initial = cont.n_assigned cdef int in_boundary = 0 + parent = NULL # How do we bootstrap ourselves? for p in range(no): #for every oct we're trying to add find the @@ -606,10 +608,12 @@ cdef class OctreeContainer: ind[i] = 1 cp[i] += dds[i]/2.0 # Check if it has not been allocated + parent = cur cur = self.next_child(curdom, ind, cur) # Now we should be at the right level cur.domain = curdom cur.file_ind = p + cur.parent = parent return cont.n_assigned - initial + nb def allocate_domains(self, domain_counts): @@ -734,7 +738,10 @@ cdef class OctreeContainer: source = source_fields[key] for i in range(levels.shape[0]): if levels[i] != level: continue - dest[i + offset] = source[file_inds[i], cell_inds[i]] + if file_inds[i] < 0: + dest[i + offset] = np.nan + else: + dest[i + offset] = source[file_inds[i], cell_inds[i]] def finalize(self): cdef SelectorObject selector = selection_routines.AlwaysSelector(None) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 18493856878..675bacdf8c5 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -145,11 +145,15 @@ cdef inline int cind(int i, int j, int k) nogil: from oct_container cimport OctreeContainer +cdef class StoreIndex(OctVisitor): + cdef np.int64_t[:,:,:,:] cell_inds + # cimport oct_container cdef class BaseNeighbourVisitor(OctVisitor): cdef int idim # 0,1,2 for x,y,z cdef int direction # +1 for +x, -1 for -x cdef np.int8_t[:] neigh_ind + cdef bint other_oct cdef Oct *neighbour cdef OctreeContainer octree cdef OctInfo oi @@ -161,8 +165,9 @@ cdef class BaseNeighbourVisitor(OctVisitor): cdef int d = (1 << self.oref) return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) -# # # cdef class NeighbourVisitor(BaseNeighbourVisitor): -# # # cdef np.int64_t[:,:,:,:] neigh_cell_inds +cdef class NeighbourVisitor(BaseNeighbourVisitor): + cdef np.int64_t[:,:,:,:] cell_inds + cdef np.int64_t[:,:,:,:] neigh_cell_inds cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.uint8_t[:] shifted_levels diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 3d4d639a977..18d5dc64042 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -337,6 +337,15 @@ cdef class MortonIndexOcts(OctVisitor): np.uint64(coord[2])) self.index += 1 +cdef class StoreIndex(OctVisitor): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + if not selected: return + self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index + + self.index += 1 cdef class BaseNeighbourVisitor(OctVisitor): def __init__(self, OctreeContainer octree, int domain_id = -1): @@ -358,7 +367,7 @@ cdef class BaseNeighbourVisitor(OctVisitor): for i in range(3): c = ((self.pos[i] << self.oref) + self.ind[i]) if i == self.idim: - fcoords[i] = (c + 0.5 + self.direction) * dx / self.octree.nn[i] + fcoords[i] = (c + 0.5 + 2*self.direction) * dx / self.octree.nn[i] else: fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] @@ -384,6 +393,12 @@ cdef class BaseNeighbourVisitor(OctVisitor): cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): cdef int i cdef bint other_oct # True if the neighbouring cell lies in another oct + cdef np.int64_t cell_ind + + # Compute information about neighbour once per oct + if self.last != o.domain_ind: + self.set_neighbour_oct(o) + self.last = o.domain_ind # Note that we provide an index even if the cell is not selected. # if selected == 0: return -1 @@ -393,29 +408,48 @@ cdef class BaseNeighbourVisitor(OctVisitor): self.neigh_ind[i] = (self.ind[i] + self.direction) other_oct = self.neigh_ind[i] < 0 or self.neigh_ind[i] > 1 if other_oct: - self.neigh_ind[i] %= 2 + # trick here: we want modulo with positive remainder, but neigh_ind may be negative so cast + # it to unsigned int *before* applying modulo. + self.neigh_ind[i] = (self.neigh_ind[i]) % 2 else: self.neigh_ind[i] = self.ind[i] - if not other_oct: - # Simple case: the neighbouring cell is within the oct - self.neighbour = o - else: - # Complicated case: the cell is in a neighbouring oct - if self.last != o.domain_ind: - self.set_neighbour_oct(o) - self.last = o.domain_ind - + self.other_oct = other_oct + if other_oct: if self.neighbour != NULL: if self.oi.level == self.level - 1: # Position within neighbouring oct is stored in oi.ipos for i in range(3): self.neigh_ind[i] = self.oi.ipos[i] - if not ( 0<= self.oi.ipos[i] <= 1): - print('WHAAAAT?!', self.oi.ipos[i]) elif self.oi.level != self.level: print('This should not happen! %s %s' % (self.oi.level, self.level)) - return + self.neighbour = NULL + +cdef class NeighbourVisitor(BaseNeighbourVisitor): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef np.int64_t cell_ind + cdef Oct *neighbour_oct + cdef bint ok + + self.get_neighbour_cell_index(o, selected) + if not self.other_oct: + neighbour_oct = o + ok = True + elif self.neighbour != NULL: + neighbour_oct = self.neighbour + ok = True + else: + ok = False + + if ok: + cell_ind = self.cell_inds[neighbour_oct.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] + else: + cell_ind = -1 + + self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): @cython.boundscheck(False) @@ -425,7 +459,7 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.int64_t neigbour_cell_index if selected == 0: return # Note: only selected items have an index - neighbour_index = self.get_neighbour_cell_index(o, selected) + self.get_neighbour_cell_index(o, selected) self.shifted_levels[self.index] = self.level if self.neighbour != NULL: # Note: we store the local level, not the remote one @@ -434,4 +468,4 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): else: self.shifted_file_inds[self.index] = -1 self.shifted_cell_inds[self.index] = 255 # -1 on uint8 - self.index += 1 \ No newline at end of file + self.index += 1 diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index edf6e60b8d0..8f471775f11 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -1,5 +1,5 @@ cimport cython -from oct_visitors cimport FillFileIndicesRNeighbour +from oct_visitors cimport FillFileIndicesRNeighbour, StoreIndex, NeighbourVisitor from selection_routines cimport SelectorObject, AlwaysSelector cimport numpy as np import numpy as np @@ -22,34 +22,51 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): bint periodicity[3]): pass - # def neighbours_in_direction(self, int idim, int direction, SelectorObject selector = AlwaysSelector(None)): - # """Return index on file of all neighbours in a given direction""" - # cdef SelectorObject always_selector = AlwaysSelector(None) + def fill_index(self, SelectorObject selector = AlwaysSelector(None)): + # Get the on-file index of each cell + cdef StoreIndex visitor - # cdef int num_cells = selector.count_oct_cells(self, -1) + cdef np.int64_t[:, :, :, :] cell_inds, - # # Get the on-file index of each cell - # cdef FillFileIndices visitor - # cdef np.ndarray[np.int64_t, ndim=4] cell_inds = np.zeros((self.nocts, 2, 2, 2), dtype="int64") + cell_inds = np.full((self.nocts, 2, 2, 2), -1, dtype=np.int64) - # visitor = FillFileIndices(self, -1) - # visitor.cell_inds = cell_inds + visitor = StoreIndex(self, -1) + visitor.cell_inds = cell_inds - # self.visit_all_octs(selector, visitor) - - # # Store the index of the neighbour - # cdef NeighbourVisitor n_visitor - # cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.empty_like(cell_inds) - # n_visitor = NeighbourVisitor(self, -1) - # n_visitor.idim = idim - # n_visitor.direction = direction - # n_visitor.cell_inds = cell_inds - # n_visitor.neigh_cell_inds = neigh_cell_inds - # n_visitor.octree = self - # n_visitor.last = -1 - # self.visit_all_octs(always_selector, n_visitor) + self.visit_all_octs(selector, visitor) - # return np.asarray(cell_inds), np.asarray(neigh_cell_inds) + return np.asarray(cell_inds) + + def neighbours_in_direction(self, int idim, int direction, + np.int64_t[:, :, :, :] cell_inds): + """Return index on file of all neighbours in a given direction""" + cdef SelectorObject always_selector = AlwaysSelector(None) + + # Store the index of the neighbour + cdef NeighbourVisitor n_visitor + cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.full_like(cell_inds, -1) + n_visitor = NeighbourVisitor(self, -1) + n_visitor.idim = idim + n_visitor.direction = direction + n_visitor.cell_inds = cell_inds + n_visitor.neigh_cell_inds = neigh_cell_inds + n_visitor.octree = self + n_visitor.last = -1 + self.visit_all_octs(always_selector, n_visitor) + + return np.asarray(neigh_cell_inds) + + #@cython.boundscheck(False) + @cython.wraparound(False) + def copy_neighbour_data(self, + np.int64_t[:] icell, np.int64_t[:] nicell, + np.float64_t[:, :] input, np.float64_t[:, :] output, + int N,): + cdef int i + + for i in range(N): + if nicell[i] > -1 and icell[i] > -1: + output[icell[i], :] = input[nicell[i], :] def file_index_octs_with_shift(self, SelectorObject selector, int domain_id, int idim, int direction, int num_cells = -1): @@ -83,9 +100,7 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): raise NotImplementedError('C-style filling with spatial offset has not been implemented.') else: raise RuntimeError - print('visiting all octs') self.visit_all_octs(selector, neigh_visitor) - print('visited') return shifted_levels, shifted_cell_inds, shifted_file_inds def file_index_octs(self, SelectorObject selector, int domain_id, From bc077f4f23cd1c6d742ec5119b344daf7158bdec Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 23 Jan 2020 16:53:40 +0000 Subject: [PATCH 09/61] Fix octree accessor --- yt/geometry/oct_visitors.pyx | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 18d5dc64042..d294575528b 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -456,16 +456,25 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): @cython.wraparound(False) @cython.initializedcheck(False) cdef void visit(self, Oct* o, np.uint8_t selected): - cdef np.int64_t neigbour_cell_index + cdef np.int64_t neigh_file_ind + cdef np.uint8_t neigh_cell_ind + if selected == 0: return # Note: only selected items have an index self.get_neighbour_cell_index(o, selected) - self.shifted_levels[self.index] = self.level - if self.neighbour != NULL: - # Note: we store the local level, not the remote one - self.shifted_file_inds[self.index] = self.neighbour.file_ind - self.shifted_cell_inds[self.index] = self.neighbour_rind() + if not self.other_oct: + neigh_file_ind = o.file_ind + neigh_cell_ind = self.neighbour_rind() + elif self.neighbour != NULL: + neigh_file_ind = self.neighbour.file_ind + neigh_cell_ind = self.neighbour_rind() else: - self.shifted_file_inds[self.index] = -1 - self.shifted_cell_inds[self.index] = 255 # -1 on uint8 + neigh_file_ind = -1 + neigh_cell_ind = 255 + + self.shifted_levels[self.index] = self.level + # Note: we store the local level, not the remote one + self.shifted_file_inds[self.index] = neigh_file_ind + self.shifted_cell_inds[self.index] = neigh_cell_ind + self.index += 1 From 266ac1fb9dab1dd8644e3bcb41bf10c6cc994345 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 09:17:30 +0100 Subject: [PATCH 10/61] Forget about nan when doing gradients --- yt/fields/fluid_fields.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 0a213e43409..f62881a0747 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -65,7 +65,7 @@ def _cell_mass(field, data): def _sound_speed(field, data): tr = data.ds.gamma * data[ftype, "pressure"] / data[ftype, "density"] return np.sqrt(tr) - + registry.add_field((ftype, "sound_speed"), sampling_type="local", function=_sound_speed, @@ -75,7 +75,7 @@ def _radial_mach_number(field, data): """ Radial component of M{|v|/c_sound} """ tr = data[ftype, "radial_velocity"] / data[ftype, "sound_speed"] return np.abs(tr) - + registry.add_field((ftype, "radial_mach_number"), sampling_type="local", function=_radial_mach_number, @@ -93,7 +93,7 @@ def _kin_energy(field, data): def _mach_number(field, data): """ M{|v|/c_sound} """ return data[ftype, "velocity_magnitude"] / data[ftype, "sound_speed"] - + registry.add_field((ftype, "mach_number"), sampling_type="local", function=_mach_number, @@ -213,7 +213,7 @@ def func(field, data): f = data[grad_field][slice_3dr]/ds[slice_3d] f -= data[grad_field][slice_3dl]/ds[slice_3d] new_field = np.zeros_like(data[grad_field], dtype=np.float64) - new_field = data.ds.arr(new_field, f.units) + new_field = data.ds.arr(new_field, vr.units / ds.units) new_field[slice_3d] = f return new_field return func From 0d8bf3acd37a5614661e668b10af351eefc6dcc9 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 23 Jan 2020 16:54:02 +0000 Subject: [PATCH 11/61] Accept unyt_array as center --- yt/geometry/coordinates/coordinate_handler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/geometry/coordinates/coordinate_handler.py b/yt/geometry/coordinates/coordinate_handler.py index 6e7f2f00bd0..b674cada4ce 100644 --- a/yt/geometry/coordinates/coordinate_handler.py +++ b/yt/geometry/coordinates/coordinate_handler.py @@ -7,7 +7,7 @@ fix_unitary, \ iterable from yt.units.yt_array import \ - YTArray, YTQuantity + YTArray, YTQuantity, unyt_array from yt.utilities.exceptions import \ YTCoordinateNotImplemented, \ YTInvalidWidthError @@ -213,7 +213,7 @@ def sanitize_width(self, axis, width, depth): width = (w[0], w[1]) elif iterable(width): width = validate_iterable_width(width, self.ds) - elif isinstance(width, YTQuantity): + elif isinstance(width, (YTQuantity, unyt_array)): width = (width, width) elif isinstance(width, Number): width = (self.ds.quan(width, 'code_length'), From 59923e22e2648ecc895527979ba9332c8dab5388 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 09:18:29 +0100 Subject: [PATCH 12/61] Gradient should be block order-aware --- yt/fields/fluid_fields.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index f62881a0747..5966f8612ff 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -205,16 +205,25 @@ def grad_func(axi, ax): slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:] slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:] def func(field, data): + block_reorder = getattr(data, '_block_reorder', 'C') + if block_reorder == 'F': + field_data = data[grad_field].swapaxes(0, 2) + else: + field_data = data[grad_field] ds = div_fac * data[ftype, "d%s" % ax] if ax == "theta": ds *= data[ftype, "r"] if ax == "phi": ds *= data[ftype, "r"] * np.sin(data[ftype, "theta"]) - f = data[grad_field][slice_3dr]/ds[slice_3d] - f -= data[grad_field][slice_3dl]/ds[slice_3d] + f = field_data[slice_3dr]/ds[slice_3d] + f -= field_data[slice_3dl]/ds[slice_3d] new_field = np.zeros_like(data[grad_field], dtype=np.float64) new_field = data.ds.arr(new_field, vr.units / ds.units) new_field[slice_3d] = f + + if block_reorder: + new_field = new_field.swapaxes(0, 2) + return new_field return func From 6d091d5277c8318563e693fff76e7ef5a16eef90 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 23 Jan 2020 18:43:46 +0000 Subject: [PATCH 13/61] cAdd comments --- yt/geometry/oct_visitors.pyx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index d294575528b..00c183cc199 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -337,6 +337,7 @@ cdef class MortonIndexOcts(OctVisitor): np.uint64(coord[2])) self.index += 1 +# Store cell index cdef class StoreIndex(OctVisitor): @cython.boundscheck(False) @cython.wraparound(False) @@ -425,6 +426,7 @@ cdef class BaseNeighbourVisitor(OctVisitor): print('This should not happen! %s %s' % (self.oi.level, self.level)) self.neighbour = NULL +# Store neighbouring cell index in current cell cdef class NeighbourVisitor(BaseNeighbourVisitor): @cython.boundscheck(False) @cython.wraparound(False) @@ -451,6 +453,7 @@ cdef class NeighbourVisitor(BaseNeighbourVisitor): self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind +# Store file position + cell of neighbouring cell in current cell cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): @cython.boundscheck(False) @cython.wraparound(False) From 08fe5a91e90ac58ba45b6efcacda0c9fef0b82e4 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 23 Jan 2020 18:44:04 +0000 Subject: [PATCH 14/61] Support reading boundaries between CPUs --- yt/frontends/ramses/data_structures.py | 17 ++++--- yt/frontends/ramses/io_utils.pyx | 67 +++++++++++++------------- 2 files changed, 44 insertions(+), 40 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index b50a3d76f28..6c729231636 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -213,8 +213,9 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, 'float64') - fill_hydro(fd, file_handler.offset, - file_handler.level_count, levels, cell_inds, + fill_hydro(fd, file_handler.offset[self.domain_id-1:self.domain_id], + file_handler.level_count[self.domain_id-1:self.domain_id], + levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler) return tr @@ -251,10 +252,13 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo for field in fields: tr[field][:] = 0 - fill_hydro(fd, file_handler.offset, - file_handler.level_count, levels, cell_inds, - file_inds, ndim, all_fields, fields, tr, - oct_handler) + fill_hydro( + fd, + file_handler.offset, file_handler.level_count, + levels, cell_inds, + file_inds, ndim, all_fields, fields, tr, + oct_handler + ) _slice = tuple( [slice(None)] + [slice(i, i+self._num_zones) for i in ishift_all]) @@ -325,7 +329,6 @@ def retrieve_ghost_zones(self, ngz, fields, smoothed=False): return new_subset - class RAMSESIndex(OctreeIndex): def __init__(self, ds, dataset_type='ramses'): diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 12fc0918f33..9414c2bf352 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -83,7 +83,7 @@ def read_amr(FortranFile f, dict headers, @cython.nonecheck(False) cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers): - cdef np.ndarray[np.int64_t, ndim=1] offset, level_count + cdef np.ndarray[np.int64_t, ndim=2] offset, level_count cdef INT64_t ndim, twotondim, nlevelmax, n_levels, nboundary, ncpu, ncpu_and_bound cdef INT64_t ilevel, icpu, skip_len cdef INT32_t file_ilevel, file_ncache @@ -101,12 +101,11 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n skip_len = twotondim * nvar # It goes: level, CPU, 8-variable (1 oct) - offset = np.zeros(n_levels, dtype=np.int64) - offset -= 1 - level_count = np.zeros(n_levels, dtype=np.int64) + offset = np.full((ncpu, n_levels), -1, dtype=np.int64) + level_count = np.zeros((ncpu, n_levels), dtype=np.int64) - cdef np.int64_t[:] level_count_view = level_count - cdef np.int64_t[:] offset_view = offset + cdef np.int64_t[:,:] level_count_view = level_count + cdef np.int64_t[:,:] offset_view = offset for ilevel in range(nlevelmax): for icpu in range(ncpu_and_bound): @@ -121,9 +120,9 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n 'from data (%s) is not coherent with the expected (%s)', f.name, file_ilevel, ilevel) - if icpu + 1 == domain_id and ilevel >= min_level: - offset[ilevel - min_level] = f.tell() - level_count[ilevel - min_level] = file_ncache + if ilevel >= min_level: + offset_view[icpu, ilevel - min_level] = f.tell() + level_count_view[icpu, ilevel - min_level] = file_ncache f.skip(skip_len) return offset, level_count @@ -133,45 +132,47 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n @cython.cdivision(True) @cython.nonecheck(False) def fill_hydro(FortranFile f, - np.ndarray[np.int64_t, ndim=1] offsets, - np.ndarray[np.int64_t, ndim=1] level_count, + np.ndarray[np.int64_t, ndim=2] offsets, + np.ndarray[np.int64_t, ndim=2] level_count, np.ndarray[np.uint8_t, ndim=1] levels, np.ndarray[np.uint8_t, ndim=1] cell_inds, np.ndarray[np.int64_t, ndim=1] file_inds, INT64_t ndim, list all_fields, list fields, dict tr, RAMSESOctreeContainer oct_handler): - cdef INT64_t ilevel, ifield, nfields, noffset cdef INT64_t offset cdef dict tmp cdef str field - cdef INT64_t twotondim, i + cdef INT64_t twotondim + cdef int ilevel, icpu, ifield, nfields, noffset, nc cdef np.ndarray[np.uint8_t, ndim=1] mask twotondim = 2**ndim nfields = len(all_fields) - noffset = len(offsets) + ncpu = offsets.shape[0] + noffset = offsets.shape[1] mask = np.array([(field in fields) for field in all_fields], dtype=np.uint8) # Loop over levels for ilevel in range(noffset): - offset = offsets[ilevel] - if offset == -1: - continue - f.seek(offset) - nc = level_count[ilevel] - tmp = {} - # Initalize temporary data container for io - for field in all_fields: - tmp[field] = np.empty((nc, twotondim), dtype="float64") - - for i in range(twotondim): - # Read the selected fields - for ifield in range(nfields): - if not mask[ifield]: - f.skip() - else: - tmp[all_fields[ifield]][:, i] = f.read_vector('d') # i-th cell - - oct_handler.fill_level(ilevel, levels, cell_inds, file_inds, tr, tmp) + for icpu in range(ncpu): + offset = offsets[icpu, ilevel] + if offset == -1: + continue + f.seek(offset) + nc = level_count[icpu, ilevel] + tmp = {} + # Initalize temporary data container for io + for field in all_fields: + tmp[field] = np.empty((nc, twotondim), dtype="float64") + + for i in range(twotondim): + # Read the selected fields + for ifield in range(nfields): + if not mask[ifield]: + f.skip() + else: + tmp[all_fields[ifield]][:, i] = f.read_vector('d') # i-th cell + + oct_handler.fill_level(ilevel, levels, cell_inds, file_inds, tr, tmp) From d6b508b4df10315caf1c0d4e4ba5ef8166708f4e Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 28 Jan 2020 09:02:22 +0000 Subject: [PATCH 15/61] Filling hydro with boundary information --- yt/frontends/ramses/data_structures.py | 31 +++++++++++++++++---- yt/frontends/ramses/io_utils.pyx | 19 +++++++++---- yt/geometry/oct_container.pyx | 4 +-- yt/geometry/oct_visitors.pxd | 1 + yt/geometry/oct_visitors.pyx | 9 ++++-- yt/geometry/ramses_oct_container.pyx | 38 +++++++++++++++++++++++--- 6 files changed, 83 insertions(+), 19 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 6c729231636..55776805fe6 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -213,8 +213,9 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, 'float64') - fill_hydro(fd, file_handler.offset[self.domain_id-1:self.domain_id], - file_handler.level_count[self.domain_id-1:self.domain_id], + fill_hydro(fd, file_handler.offset, + file_handler.level_count, + [self.domain_id-1], levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler) @@ -228,10 +229,19 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo all_fields = [f for ft, f in file_handler.field_list] fields = [f for ft, f in fields] - oct_count = selector.count_octs(self.oct_handler, self.domain_id) + # Select all cells, including those in *other domain* + selector = OctreeSubsetSelector(self) + selector_with_edge = OctreeSubsetSelector(self) + selector_with_edge.domain_id = -1 + + oct_count = selector_with_edge.count_octs(self.oct_handler, -1) cell_count = oct_count * self._num_zones**ndim iwidth = self._num_zones + num_ghost_zones * 2 + oct_count0 = selector.count_octs(self.oct_handler, self.domain_id) + cell_count0 = oct_count0 * self._num_zones**ndim + ncpus = file_handler.offset.shape[0] + tr = {} tr_all = {} for field in fields: @@ -245,9 +255,15 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo ishift_all[idim] = num_ghost_zones + shift if shift == 0: continue - levels, cell_inds, file_inds = self.oct_handler.file_index_octs_with_shift( - selector, self.domain_id, idim, shift, cell_count) + levels0, cell_inds0, file_inds0, domain0 = \ + self.oct_handler.file_index_octs_with_shift( + selector, self.domain_id, idim, shift, cell_count0) + + levels, cell_inds, file_inds, domain = \ + self.oct_handler.file_index_octs_with_shift( + selector_with_edge, -1, idim, shift, cell_count) + # import ipdb; ipdb.set_trace() # Initializing data container for field in fields: tr[field][:] = 0 @@ -255,9 +271,12 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo fill_hydro( fd, file_handler.offset, file_handler.level_count, + #file_handler.offset, file_handler.level_count, + list(range(ncpus)), levels, cell_inds, file_inds, ndim, all_fields, fields, tr, - oct_handler + oct_handler, + domains=domain.astype(np.int64) ) _slice = tuple( [slice(None)] + diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 9414c2bf352..0993668b0b2 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -134,29 +134,34 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n def fill_hydro(FortranFile f, np.ndarray[np.int64_t, ndim=2] offsets, np.ndarray[np.int64_t, ndim=2] level_count, + list cpu_enumerator, np.ndarray[np.uint8_t, ndim=1] levels, np.ndarray[np.uint8_t, ndim=1] cell_inds, np.ndarray[np.int64_t, ndim=1] file_inds, INT64_t ndim, list all_fields, list fields, dict tr, - RAMSESOctreeContainer oct_handler): + RAMSESOctreeContainer oct_handler, + np.ndarray[np.int64_t, ndim=1] domains=np.array([], dtype='int64'), + int domain=-1): cdef INT64_t offset cdef dict tmp cdef str field cdef INT64_t twotondim - cdef int ilevel, icpu, ifield, nfields, noffset, nc + cdef int ilevel, icpu, ifield, nfields, noffset, nc, ncpu_selected cdef np.ndarray[np.uint8_t, ndim=1] mask twotondim = 2**ndim nfields = len(all_fields) ncpu = offsets.shape[0] noffset = offsets.shape[1] + ncpu_selected = len(cpu_enumerator) mask = np.array([(field in fields) for field in all_fields], dtype=np.uint8) # Loop over levels for ilevel in range(noffset): - for icpu in range(ncpu): + # Loop over cpu domains + for icpu in cpu_enumerator: offset = offsets[icpu, ilevel] if offset == -1: continue @@ -174,5 +179,9 @@ def fill_hydro(FortranFile f, f.skip() else: tmp[all_fields[ifield]][:, i] = f.read_vector('d') # i-th cell - - oct_handler.fill_level(ilevel, levels, cell_inds, file_inds, tr, tmp) + if ncpu_selected > 1: + oct_handler.fill_level_with_domain( + ilevel, levels, cell_inds, file_inds, domains, tr, tmp, domain=icpu) + else: + oct_handler.fill_level( + ilevel, levels, cell_inds, file_inds, tr, tmp) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 73d349547f7..02f00f19437 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -669,9 +669,9 @@ cdef class OctreeContainer: file_inds = np.zeros(num_cells, dtype="int64") cell_inds = np.zeros(num_cells, dtype="uint8") for i in range(num_cells): - levels[i] = 100 + levels[i] = 255 file_inds[i] = -1 - cell_inds[i] = 9 + cell_inds[i] = 8 cdef oct_visitors.FillFileIndicesO visitor_o cdef oct_visitors.FillFileIndicesR visitor_r if self.fill_style == "r": diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 675bacdf8c5..5e0f8a395c3 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -173,4 +173,5 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.uint8_t[:] shifted_levels cdef np.int64_t[:] shifted_file_inds cdef np.uint8_t[:] shifted_cell_inds + cdef np.int32_t[:] neigh_domain diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 00c183cc199..a2cc19f6e2a 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -461,23 +461,28 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef void visit(self, Oct* o, np.uint8_t selected): cdef np.int64_t neigh_file_ind cdef np.uint8_t neigh_cell_ind + cdef np.int32_t neigh_domain if selected == 0: return # Note: only selected items have an index self.get_neighbour_cell_index(o, selected) if not self.other_oct: - neigh_file_ind = o.file_ind + neigh_domain = 0 + neigh_file_ind = o.domain neigh_cell_ind = self.neighbour_rind() elif self.neighbour != NULL: + neigh_domain = self.neighbour.domain neigh_file_ind = self.neighbour.file_ind neigh_cell_ind = self.neighbour_rind() else: + neigh_domain = -1 neigh_file_ind = -1 - neigh_cell_ind = 255 + neigh_cell_ind = 8 self.shifted_levels[self.index] = self.level # Note: we store the local level, not the remote one self.shifted_file_inds[self.index] = neigh_file_ind self.shifted_cell_inds[self.index] = neigh_cell_ind + self.neigh_domain[self.index] = neigh_domain self.index += 1 diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index 8f471775f11..bc001b2e836 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -81,9 +81,11 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels cdef np.ndarray[np.uint8_t, ndim=1] shifted_cell_inds cdef np.ndarray[np.int64_t, ndim=1] shifted_file_inds - shifted_levels = np.zeros(num_cells, dtype="uint8") - shifted_file_inds = np.zeros(num_cells, dtype="int64") - shifted_cell_inds = np.zeros(num_cells, dtype="uint8") + cdef np.ndarray[np.int32_t, ndim=1] neigh_domain + shifted_levels = np.full(num_cells, 255, dtype="uint8") + shifted_file_inds = np.full(num_cells, -1, dtype="int64") + shifted_cell_inds = np.full(num_cells, 8, dtype="uint8") + neigh_domain = np.full(num_cells, -1, dtype="int32") if self.fill_style == "r": neigh_visitor = FillFileIndicesRNeighbour(self, domain_id) @@ -91,6 +93,7 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): neigh_visitor.shifted_levels = shifted_levels neigh_visitor.shifted_file_inds = shifted_file_inds neigh_visitor.shifted_cell_inds = shifted_cell_inds + neigh_visitor.neigh_domain = neigh_domain # direction to explore and extra parameters of the visitor neigh_visitor.idim = idim neigh_visitor.direction = direction @@ -101,7 +104,7 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): else: raise RuntimeError self.visit_all_octs(selector, neigh_visitor) - return shifted_levels, shifted_cell_inds, shifted_file_inds + return shifted_levels, shifted_cell_inds, shifted_file_inds, neigh_domain def file_index_octs(self, SelectorObject selector, int domain_id, num_cells = -1, spatial_offset=(0, 0, 0)): @@ -129,3 +132,30 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): else: return self.file_index_octs_with_shift( selector, domain_id, idim, direction, num_cells) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def fill_level_with_domain( + self, int level, + np.uint8_t[:] levels, + np.uint8_t[:] cell_inds, + np.int64_t[:] file_inds, + np.int64_t[:] domains, + dest_fields, source_fields, + np.int64_t domain, + np.int64_t offset = 0 + ): + cdef np.ndarray[np.float64_t, ndim=2] source + cdef np.ndarray[np.float64_t, ndim=1] dest + cdef int i + cdef str key + for key in dest_fields: + dest = dest_fields[key] + source = source_fields[key] + for i in range(levels.shape[0]): + if levels[i] != level or domains[i] != domain: continue + if file_inds[i] < 0: + dest[i + offset] = np.nan + else: + dest[i + offset] = source[file_inds[i], cell_inds[i]] From fb6462aa2309cb0e0f373de9ee8ac4a73acd3582 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 09:10:45 +0000 Subject: [PATCH 16/61] Use Fortran ordering to read in for faster results --- yt/frontends/ramses/io_utils.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 0993668b0b2..0cb1398ea72 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -170,7 +170,7 @@ def fill_hydro(FortranFile f, tmp = {} # Initalize temporary data container for io for field in all_fields: - tmp[field] = np.empty((nc, twotondim), dtype="float64") + tmp[field] = np.empty((nc, twotondim), dtype="float64", order='F') for i in range(twotondim): # Read the selected fields From d23b9a1649ec9ef2413d70803d25d4a5babe8fc1 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 12:57:04 +0000 Subject: [PATCH 17/61] Working version? This is not crashing and there are very few glitches. --- yt/frontends/ramses/data_structures.py | 67 +++---------- yt/frontends/ramses/io_utils.pyx | 13 +-- yt/geometry/oct_container.pyx | 2 +- yt/geometry/oct_visitors.pxd | 5 + yt/geometry/oct_visitors.pyx | 124 ++++++++++++++++++++++++- yt/geometry/ramses_oct_container.pyx | 102 ++++++++++---------- 6 files changed, 198 insertions(+), 115 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 55776805fe6..7f2dd60c51a 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -228,65 +228,24 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo oct_handler = self.oct_handler all_fields = [f for ft, f in file_handler.field_list] fields = [f for ft, f in fields] + tr = {} - # Select all cells, including those in *other domain* - selector = OctreeSubsetSelector(self) - selector_with_edge = OctreeSubsetSelector(self) - selector_with_edge.domain_id = -1 - - oct_count = selector_with_edge.count_octs(self.oct_handler, -1) - cell_count = oct_count * self._num_zones**ndim - iwidth = self._num_zones + num_ghost_zones * 2 + cell_count = selector.count_octs(self.oct_handler, self.domain_id) * self.nz**ndim - oct_count0 = selector.count_octs(self.oct_handler, self.domain_id) - cell_count0 = oct_count0 * self._num_zones**ndim - ncpus = file_handler.offset.shape[0] + levels, cell_inds, file_inds, domains = self.oct_handler.compute_domain_mapper( + selector, self.domain_id, cell_count) - tr = {} - tr_all = {} + # Initializing data container for field in fields: - tr_all[field] = np.full((oct_count, iwidth, iwidth, iwidth), np.nan, 'float64') tr[field] = np.zeros(cell_count, 'float64') - - # Compute the index to read with a positive and negative shift in all dimensions - for idim in range(ndim): - ishift_all = [num_ghost_zones]*ndim - for shift in range(-num_ghost_zones, num_ghost_zones+1, 2): - ishift_all[idim] = num_ghost_zones + shift - if shift == 0: - continue - levels0, cell_inds0, file_inds0, domain0 = \ - self.oct_handler.file_index_octs_with_shift( - selector, self.domain_id, idim, shift, cell_count0) - - levels, cell_inds, file_inds, domain = \ - self.oct_handler.file_index_octs_with_shift( - selector_with_edge, -1, idim, shift, cell_count) - - # import ipdb; ipdb.set_trace() - # Initializing data container - for field in fields: - tr[field][:] = 0 - - fill_hydro( - fd, - file_handler.offset, file_handler.level_count, - #file_handler.offset, file_handler.level_count, - list(range(ncpus)), - levels, cell_inds, - file_inds, ndim, all_fields, fields, tr, - oct_handler, - domains=domain.astype(np.int64) - ) - _slice = tuple( - [slice(None)] + - [slice(i, i+self._num_zones) for i in ishift_all]) - for field in fields: - tr_all[field][_slice] = \ - tr[field].reshape(oct_count, 2, 2, 2) - for field in fields: - tr_all[field] = tr_all[field].reshape(-1) - return tr_all + fill_hydro(fd, file_handler.offset, + file_handler.level_count, + [self.domain_id-1], + levels, cell_inds, + file_inds, ndim, all_fields, fields, tr, + oct_handler, + domains=domains) + return tr @property def fwidth(self): diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 0cb1398ea72..fb74ed48216 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -141,32 +141,33 @@ def fill_hydro(FortranFile f, INT64_t ndim, list all_fields, list fields, dict tr, RAMSESOctreeContainer oct_handler, - np.ndarray[np.int64_t, ndim=1] domains=np.array([], dtype='int64'), - int domain=-1): + np.ndarray[np.int32_t, ndim=1] domains=np.array([], dtype='int32')): cdef INT64_t offset cdef dict tmp cdef str field cdef INT64_t twotondim - cdef int ilevel, icpu, ifield, nfields, noffset, nc, ncpu_selected + cdef int ilevel, icpu, ifield, nfields, nlevels, nc, ncpu_selected cdef np.ndarray[np.uint8_t, ndim=1] mask twotondim = 2**ndim nfields = len(all_fields) ncpu = offsets.shape[0] - noffset = offsets.shape[1] + nlevels = offsets.shape[1] ncpu_selected = len(cpu_enumerator) mask = np.array([(field in fields) for field in all_fields], dtype=np.uint8) # Loop over levels - for ilevel in range(noffset): + for ilevel in range(nlevels): # Loop over cpu domains for icpu in cpu_enumerator: + nc = level_count[icpu, ilevel] + if nc == 0: + continue offset = offsets[icpu, ilevel] if offset == -1: continue f.seek(offset) - nc = level_count[icpu, ilevel] tmp = {} # Initalize temporary data container for io for field in all_fields: diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 02f00f19437..437462c333f 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -548,7 +548,7 @@ cdef class OctreeContainer: def domain_ind(self, selector, int domain_id = -1): cdef np.ndarray[np.int64_t, ndim=1] ind # Here's where we grab the masked items. - ind = np.zeros(self.nocts, 'int64') - 1 + ind = np.full(self.nocts, -1, 'int64') cdef oct_visitors.IndexOcts visitor visitor = oct_visitors.IndexOcts(self, domain_id) visitor.oct_index = ind diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 5e0f8a395c3..1616b36b943 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -175,3 +175,8 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.uint8_t[:] shifted_cell_inds cdef np.int32_t[:] neigh_domain +cdef class NeighbourCellVisitor(BaseNeighbourVisitor): + cdef np.uint8_t[:] shifted_levels + cdef np.int64_t[:] shifted_file_inds + cdef np.uint8_t[:] shifted_cell_inds + cdef np.int32_t[:] neigh_domain \ No newline at end of file diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index a2cc19f6e2a..3e5c708907d 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -471,7 +471,7 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): neigh_file_ind = o.domain neigh_cell_ind = self.neighbour_rind() elif self.neighbour != NULL: - neigh_domain = self.neighbour.domain + neigh_domain = self.neighbour.domain neigh_file_ind = self.neighbour.file_ind neigh_cell_ind = self.neighbour_rind() else: @@ -486,3 +486,125 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): self.neigh_domain[self.index] = neigh_domain self.index += 1 + +# Store file position + cell of neighbouring cell in current cell +cdef class NeighbourCellNeighbourCellVisitor(BaseNeighbourVisitor): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void set_neighbour_info(self, Oct *o, int ishift[3]): + cdef int i + cdef np.float64_t c, dx + cdef np.int64_t ipos + cdef np.float64_t fcoords[3] + cdef Oct *neighbour + cdef bint local_oct + cdef bint other_oct + dx = 1.0 / ((1 << self.oref) << self.level) + local_oct = True + + # Compute position of neighbouring cell + for i in range(3): + c = (self.pos[i] << self.oref) + fcoords[i] = (c + 0.5 + ishift[i]) * dx / self.octree.nn[i] + local_oct &= (0 <= ishift[i] <= 1) + other_oct = not local_oct + + # Use octree to find neighbour + if other_oct: + neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) + else: + neighbour = o + self.oi.level = self.level + for i in range(3): + self.oi.ipos[i] = (self.pos[i] << self.oref) + ishift[i] + + # Extra step - compute cell position in neighbouring oct (and store in oi.ipos) + if self.oi.level == self.level - 1: + for i in range(3): + ipos = (((self.pos[i] << self.oref) + ishift[i])) >> 1 + if (self.oi.ipos[i] << 1) == ipos: + self.oi.ipos[i] = 0 + else: + self.oi.ipos[i] = 1 + self.neighbour = neighbour + + # Index of neighbouring cell within its oct + for i in range(3): + self.neigh_ind[i] = (ishift[i]) % 2 + + self.other_oct = other_oct + if other_oct: + if self.neighbour != NULL: + if self.oi.level == self.level - 1: + # Position within neighbouring oct is stored in oi.ipos + for i in range(3): + self.neigh_ind[i] = self.oi.ipos[i] + elif self.oi.level != self.level: + print('This should not happen! %s %s' % (self.oi.level, self.level)) + self.neighbour = NULL + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef int i, j, k + cdef int ishift[3] + cdef np.int64_t neigh_file_ind + cdef np.uint8_t neigh_cell_ind + cdef np.int32_t neigh_domain + cdef np.uint8_t neigh_level + if selected == 0: return + # Work at oct level + if self.last == o.domain_ind: return + + self.last = o.domain_ind + + # Loop over cells in and directly around oct + for i in range(-1, 3): + ishift[0] = i + for j in range(-1, 3): + ishift[1] = j + for k in range(-1, 3): + ishift[2] = k + self.set_neighbour_info(o, ishift) + + if not self.other_oct: + neigh_level = self.level + neigh_domain = o.domain + neigh_file_ind = o.file_ind + neigh_cell_ind = self.neighbour_rind() + elif self.neighbour != NULL: + neigh_level = self.oi.level + neigh_domain = self.neighbour.domain + neigh_file_ind = self.neighbour.file_ind + neigh_cell_ind = self.neighbour_rind() + else: + neigh_level = 255 + neigh_domain = -1 + neigh_file_ind = -1 + neigh_cell_ind = 8 + + self.shifted_levels[self.index] = neigh_level + self.shifted_file_inds[self.index] = neigh_file_ind + self.shifted_cell_inds[self.index] = neigh_cell_ind + self.neigh_domain[self.index] = neigh_domain + + self.index += 1# Compute mapping from index in domain to index in domain with buffer zones +cdef class DomainMapper(BaseNeighbourVisitor): + # Intended to map all octs ids to only octs in local domain + # Should be used in conjunction with OctreeSubsetSelector + def __init__(self, OctreeContainer octree, int domain_id): + super(DomainMapper, self).__init__(octree, domain_id) + self.index_in = 0 + + #@cython.boundscheck(False) + #@cython.wraparound(False) + #@cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + if self.last != o.domain_ind: + self.last = o.domain_ind + if self.marked[self.index]: + self.mapping[self.index] = self.index_in + self.index_in += 1 + self.index += 1 diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index bc001b2e836..40a7d7d94e2 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -1,18 +1,9 @@ cimport cython -from oct_visitors cimport FillFileIndicesRNeighbour, StoreIndex, NeighbourVisitor -from selection_routines cimport SelectorObject, AlwaysSelector +from oct_visitors cimport StoreIndex, NeighbourVisitor, NeighbourCellVisitor +from selection_routines cimport SelectorObject, AlwaysSelector, OctreeSubsetSelector cimport numpy as np import numpy as np -# cdef class FillFileIndices(oct_visitors.OctVisitor): -# cdef np.int64_t[:,:,:,:] cell_inds -# @cython.boundscheck(False) -# @cython.wraparound(False) -# @cython.initializedcheck(False) -# cdef void visit(self, Oct* o, np.uint8_t selected): -# if selected == 0: return -# self.cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = self.index -# self.index += 1 cdef class RAMSESOctreeContainer(SparseOctreeContainer): @cython.boundscheck(False) @@ -68,44 +59,6 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): if nicell[i] > -1 and icell[i] > -1: output[icell[i], :] = input[nicell[i], :] - def file_index_octs_with_shift(self, SelectorObject selector, int domain_id, - int idim, int direction, int num_cells = -1): - """Return index on file of all neighbours in a given direction""" - # We create oct arrays of the correct size - cdef np.int64_t i - if num_cells < 0: - num_cells = selector.count_oct_cells(self, domain_id) - - # Fill value of each cell with its neighbouring value - cdef FillFileIndicesRNeighbour neigh_visitor - cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels - cdef np.ndarray[np.uint8_t, ndim=1] shifted_cell_inds - cdef np.ndarray[np.int64_t, ndim=1] shifted_file_inds - cdef np.ndarray[np.int32_t, ndim=1] neigh_domain - shifted_levels = np.full(num_cells, 255, dtype="uint8") - shifted_file_inds = np.full(num_cells, -1, dtype="int64") - shifted_cell_inds = np.full(num_cells, 8, dtype="uint8") - neigh_domain = np.full(num_cells, -1, dtype="int32") - - if self.fill_style == "r": - neigh_visitor = FillFileIndicesRNeighbour(self, domain_id) - # output: level, file_ind and cell_ind of the neighbouring cells - neigh_visitor.shifted_levels = shifted_levels - neigh_visitor.shifted_file_inds = shifted_file_inds - neigh_visitor.shifted_cell_inds = shifted_cell_inds - neigh_visitor.neigh_domain = neigh_domain - # direction to explore and extra parameters of the visitor - neigh_visitor.idim = idim - neigh_visitor.direction = direction - neigh_visitor.octree = self - neigh_visitor.last = -1 - elif self.fill_style == "o": - raise NotImplementedError('C-style filling with spatial offset has not been implemented.') - else: - raise RuntimeError - self.visit_all_octs(selector, neigh_visitor) - return shifted_levels, shifted_cell_inds, shifted_file_inds, neigh_domain - def file_index_octs(self, SelectorObject selector, int domain_id, num_cells = -1, spatial_offset=(0, 0, 0)): @@ -142,20 +95,63 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): np.uint8_t[:] cell_inds, np.int64_t[:] file_inds, np.int64_t[:] domains, - dest_fields, source_fields, - np.int64_t domain, + dict dest_fields, dict source_fields, + np.int32_t domain, np.int64_t offset = 0 ): cdef np.ndarray[np.float64_t, ndim=2] source cdef np.ndarray[np.float64_t, ndim=1] dest - cdef int i + cdef np.float64_t tmp + cdef int i, count cdef str key for key in dest_fields: dest = dest_fields[key] source = source_fields[key] + count = 0 for i in range(levels.shape[0]): if levels[i] != level or domains[i] != domain: continue + count += 1 if file_inds[i] < 0: dest[i + offset] = np.nan else: - dest[i + offset] = source[file_inds[i], cell_inds[i]] + # print(f'\t{i}: Accessing source {file_inds[i]}:{cell_inds[i]} source.shape=({source.shape[0]},{source.shape[1]})') + tmp =source[file_inds[i], cell_inds[i]] + dest[i + offset] = tmp # source[file_inds[i], cell_inds[i]] + return count + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def file_index_octs_with_shift( + self, SelectorObject selector, int domain_id, + int num_cells = -1): + cdef np.int64_t i + cdef int num_octs + if num_cells < 0: + num_octs = selector.count_octs(self, domain_id) + num_cells = num_octs * 4**3 + cdef NeighbourCellVisitor visitor + + cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels + cdef np.ndarray[np.uint8_t, ndim=1] shifted_cell_inds + cdef np.ndarray[np.int64_t, ndim=1] shifted_file_inds + cdef np.ndarray[np.int32_t, ndim=1] neigh_domain + shifted_levels = np.full(num_cells, 255, dtype="uint8") + shifted_file_inds = np.full(num_cells, -1, dtype="int64") + shifted_cell_inds = np.full(num_cells, 8, dtype="uint8") + neigh_domain = np.full(num_cells, -1, dtype="int32") + + visitor = NeighbourCellVisitor(self, -1) + # output: level, file_ind and cell_ind of the neighbouring cells + visitor.shifted_levels = shifted_levels + visitor.shifted_file_inds = shifted_file_inds + visitor.shifted_cell_inds = shifted_cell_inds + visitor.neigh_domain = neigh_domain + # direction to explore and extra parameters of the visitor + visitor.octree = self + visitor.last = -1 + + # Compute indices + self.visit_all_octs(selector, visitor) + + return shifted_levels, shifted_cell_inds, shifted_file_inds, neigh_domain \ No newline at end of file From 5859895159e44dc857f473fe81a755e630413dfe Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 12:57:38 +0000 Subject: [PATCH 18/61] Correct erroneous reference --- yt/geometry/oct_visitors.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 3e5c708907d..36474c5ac1e 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -467,8 +467,8 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): # Note: only selected items have an index self.get_neighbour_cell_index(o, selected) if not self.other_oct: - neigh_domain = 0 - neigh_file_ind = o.domain + neigh_domain = o.domain + neigh_file_ind = o.file_ind neigh_cell_ind = self.neighbour_rind() elif self.neighbour != NULL: neigh_domain = self.neighbour.domain From cac9d1968bb1a3543442bc8244cb9fe1ccc15e6f Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 12:57:56 +0000 Subject: [PATCH 19/61] Catch early bug --- yt/frontends/ramses/io_utils.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index fb74ed48216..7ca7d106b9d 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -74,6 +74,8 @@ def read_amr(FortranFile f, dict headers, count_boundary = 1) if n > 0: max_level = max(ilevel - min_level, max_level) + if n != ng: + raise Exception('Expected %s octs, got %s' % (ng, n)) return max_level From dbeccb83d53209f959293e0e41a70630bba1e7a2 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:03:46 +0000 Subject: [PATCH 20/61] More explicit names --- yt/frontends/ramses/data_structures.py | 2 +- yt/geometry/ramses_oct_container.pyx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 7f2dd60c51a..c6258ba38ef 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -232,7 +232,7 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo cell_count = selector.count_octs(self.oct_handler, self.domain_id) * self.nz**ndim - levels, cell_inds, file_inds, domains = self.oct_handler.compute_domain_mapper( + levels, cell_inds, file_inds, domains = self.oct_handler.file_index_octs_with_ghost_zones( selector, self.domain_id, cell_count) # Initializing data container diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index 40a7d7d94e2..af390d68b82 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -122,7 +122,7 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) - def file_index_octs_with_shift( + def file_index_octs_with_ghost_zones( self, SelectorObject selector, int domain_id, int num_cells = -1): cdef np.int64_t i From 2ddaebb0f5a02828104149ce0d9c7b0b377df512 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:11:16 +0000 Subject: [PATCH 21/61] Correct discrepancy between .pxd and .pyx file --- yt/geometry/oct_visitors.pxd | 4 +++- yt/geometry/oct_visitors.pyx | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 1616b36b943..84199f27230 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -179,4 +179,6 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): cdef np.uint8_t[:] shifted_levels cdef np.int64_t[:] shifted_file_inds cdef np.uint8_t[:] shifted_cell_inds - cdef np.int32_t[:] neigh_domain \ No newline at end of file + cdef np.int32_t[:] neigh_domain + + cdef void set_neighbour_info(self, Oct *o, int ishift[3]) \ No newline at end of file diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 36474c5ac1e..ebb14f7cee7 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -488,7 +488,7 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): self.index += 1 # Store file position + cell of neighbouring cell in current cell -cdef class NeighbourCellNeighbourCellVisitor(BaseNeighbourVisitor): +cdef class NeighbourCellVisitor(BaseNeighbourVisitor): @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) From 3c20d534d2c6d249dc6d63856fae3e3d599e423b Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:19:23 +0000 Subject: [PATCH 22/61] UNSURE ABOUT THIS ONE, REMOVE ME IF YOU HAVE WEIRD RESULTS --- yt/frontends/ramses/io_utils.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 7ca7d106b9d..e666b68deb3 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -184,7 +184,7 @@ def fill_hydro(FortranFile f, tmp[all_fields[ifield]][:, i] = f.read_vector('d') # i-th cell if ncpu_selected > 1: oct_handler.fill_level_with_domain( - ilevel, levels, cell_inds, file_inds, domains, tr, tmp, domain=icpu) + ilevel, levels, cell_inds, file_inds, domains, tr, tmp, domain=icpu+1) else: oct_handler.fill_level( ilevel, levels, cell_inds, file_inds, tr, tmp) From 12657695638595058733407f13e20de8d38520e0 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:33:20 +0000 Subject: [PATCH 23/61] Read boundary regions as well --- yt/frontends/ramses/data_structures.py | 3 ++- yt/geometry/ramses_oct_container.pyx | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index c6258ba38ef..62b894541b8 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -223,6 +223,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zones): ndim = self.ds.dimensionality + ncpu = self.ds.parameters['ncpu'] # Here we get a copy of the file, which we skip through and read the # bits we want. oct_handler = self.oct_handler @@ -240,7 +241,7 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo tr[field] = np.zeros(cell_count, 'float64') fill_hydro(fd, file_handler.offset, file_handler.level_count, - [self.domain_id-1], + list(range(ncpu)), levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler, diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx index af390d68b82..3d80bf4ffc5 100644 --- a/yt/geometry/ramses_oct_container.pyx +++ b/yt/geometry/ramses_oct_container.pyx @@ -94,7 +94,7 @@ cdef class RAMSESOctreeContainer(SparseOctreeContainer): np.uint8_t[:] levels, np.uint8_t[:] cell_inds, np.int64_t[:] file_inds, - np.int64_t[:] domains, + np.int32_t[:] domains, dict dest_fields, dict source_fields, np.int32_t domain, np.int64_t offset = 0 From 41f05017a9cf769a3896bc5abd05eb3a9d014e50 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:46:52 +0000 Subject: [PATCH 24/61] Cannot assume that keys are string by default (can be tuple) --- yt/geometry/oct_container.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 437462c333f..63d80994bb6 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -732,7 +732,7 @@ cdef class OctreeContainer: cdef np.ndarray[np.float64_t, ndim=2] source cdef np.ndarray[np.float64_t, ndim=1] dest cdef int i - cdef str key + for key in dest_fields: dest = dest_fields[key] source = source_fields[key] From a39b4e7518446da4289d23cff14820627457ddfd Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 13:58:45 +0000 Subject: [PATCH 25/61] Inform gradient about ordering --- yt/data_objects/octree_subset.py | 1 + yt/fields/fluid_fields.py | 6 +++--- yt/frontends/ramses/data_structures.py | 1 + 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/yt/data_objects/octree_subset.py b/yt/data_objects/octree_subset.py index f98876c5e34..bf02807cf05 100644 --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -36,6 +36,7 @@ class OctreeSubset(YTSelectionContainer): _domain_offset = 0 _cell_count = -1 _block_reorder = None + _gradient_swap_axes = False # Set to True if one should swap the axes when computing gradient def __init__(self, base_region, domain, ds, over_refine_factor = 1, num_ghost_zones = 0): super(OctreeSubset, self).__init__(ds, None) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 5966f8612ff..89cdf112f53 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -205,8 +205,8 @@ def grad_func(axi, ax): slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:] slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:] def func(field, data): - block_reorder = getattr(data, '_block_reorder', 'C') - if block_reorder == 'F': + reorder = getattr(data, '_gradient_swap_axes', False) + if reorder: field_data = data[grad_field].swapaxes(0, 2) else: field_data = data[grad_field] @@ -221,7 +221,7 @@ def func(field, data): new_field = data.ds.arr(new_field, vr.units / ds.units) new_field[slice_3d] = f - if block_reorder: + if reorder: new_field = new_field.swapaxes(0, 2) return new_field diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 62b894541b8..2c829e4f1b0 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -183,6 +183,7 @@ class RAMSESDomainSubset(OctreeSubset): _domain_offset = 1 _block_reorder = "F" + _gradient_swap_axes = True _base_domain = None From 463cd9a0d57185508a1b2f8c92063a6413b9ee4d Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 09:19:12 +0100 Subject: [PATCH 26/61] Use less ambiguous name ds usually refers to a dataset. --- yt/fields/fluid_fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 89cdf112f53..9a060259ee2 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -218,7 +218,7 @@ def func(field, data): f = field_data[slice_3dr]/ds[slice_3d] f -= field_data[slice_3dl]/ds[slice_3d] new_field = np.zeros_like(data[grad_field], dtype=np.float64) - new_field = data.ds.arr(new_field, vr.units / ds.units) + new_field = data.ds.arr(new_field, vr.units / dt.units) new_field[slice_3d] = f if reorder: From 295695ebeb265ed74b141ab3e17872557f23d42e Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 16:28:34 +0000 Subject: [PATCH 27/61] Use block_reorder to decide when transposing --- yt/data_objects/octree_subset.py | 1 - yt/fields/fluid_fields.py | 6 +++--- yt/frontends/ramses/data_structures.py | 1 - 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/yt/data_objects/octree_subset.py b/yt/data_objects/octree_subset.py index bf02807cf05..f98876c5e34 100644 --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -36,7 +36,6 @@ class OctreeSubset(YTSelectionContainer): _domain_offset = 0 _cell_count = -1 _block_reorder = None - _gradient_swap_axes = False # Set to True if one should swap the axes when computing gradient def __init__(self, base_region, domain, ds, over_refine_factor = 1, num_ghost_zones = 0): super(OctreeSubset, self).__init__(ds, None) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 9a060259ee2..ce927240a2a 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -205,8 +205,8 @@ def grad_func(axi, ax): slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:] slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:] def func(field, data): - reorder = getattr(data, '_gradient_swap_axes', False) - if reorder: + block_reorder = getattr(data, '_block_reorder', None) + if block_reorder == 'F': field_data = data[grad_field].swapaxes(0, 2) else: field_data = data[grad_field] @@ -221,7 +221,7 @@ def func(field, data): new_field = data.ds.arr(new_field, vr.units / dt.units) new_field[slice_3d] = f - if reorder: + if block_reorder == 'F': new_field = new_field.swapaxes(0, 2) return new_field diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 2c829e4f1b0..62b894541b8 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -183,7 +183,6 @@ class RAMSESDomainSubset(OctreeSubset): _domain_offset = 1 _block_reorder = "F" - _gradient_swap_axes = True _base_domain = None From 9c04156e07b38a5397057ecb1afac4b43a4c1751 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 16:28:48 +0000 Subject: [PATCH 28/61] ("index, "ones") should return data with units --- yt/fields/geometric_fields.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/yt/fields/geometric_fields.py b/yt/fields/geometric_fields.py index 3b8fa041fa3..d6c73b3b3b7 100644 --- a/yt/fields/geometric_fields.py +++ b/yt/fields/geometric_fields.py @@ -94,9 +94,10 @@ def _zeros(field, data): def _ones(field, data): """Returns one for all cells""" arr = np.ones(data.ires.shape, dtype="float64") + tmp = data.apply_units(arr, field.units) if data._spatial: - return data._reshape_vals(arr) - return data.apply_units(arr, field.units) + return data._reshape_vals(tmp) + return tmp registry.add_field(("index", "ones"), sampling_type="cell", From 7c8fe2fd86e789d74451a5a1e10001d827b79e98 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 09:19:42 +0100 Subject: [PATCH 29/61] Add unit test --- yt/frontends/ramses/tests/test_outputs.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index 109ded6608c..7b55031ffaa 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -461,3 +461,17 @@ def test_magnetic_field_aliasing(): 'magnetic_field_divergence']: assert ('gas',field) in ds.derived_field_list ad[('gas',field)] + +output_00080 = "output_00080/info_00080.txt" +@requires_file(output_00080) +def test_field_accession(): + ds = yt.load(output_00080) + fields = [ + ('gas', 'density'), # basic ones + ('gas' ,'pressure'), + ('gas', 'pressure_gradient_magnitude'), # requires ghost zones + ] + # Check accessing gradient works for a variety of spatial domains + for reg in (ds.all_data(), ds.sphere([.1]*3, .01), ds.sphere([.5]*3, 0.05), ds.box([.1]*3, [.2]*3)): + for field in fields: + reg[field] From c42a2fb0bac87fb1795e607afe71b7e13153e6ca Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 16:52:20 +0000 Subject: [PATCH 30/61] Add new answer test --- yt/frontends/ramses/tests/test_outputs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index 7b55031ffaa..1538344e694 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -17,7 +17,7 @@ import yt import numpy as np -_fields = ("temperature", "density", "velocity_magnitude") +_fields = ("temperature", "density", "velocity_magnitude", "pressure_gradient_magnitude") output_00080 = "output_00080/info_00080.txt" @requires_ds(output_00080) From 7a7807dd69a1ca4d9ddf19a64087ce7d4547c53c Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 17:15:47 +0000 Subject: [PATCH 31/61] Move non-specific code into oct_container.pyx --- setup.py | 5 - yt/frontends/ramses/data_structures.py | 2 +- yt/frontends/ramses/io_utils.pyx | 2 +- yt/geometry/oct_container.pxd | 2 + yt/geometry/oct_container.pyx | 159 ++++++++++++++++++++++++- yt/geometry/oct_visitors.pxd | 8 +- yt/geometry/oct_visitors.pyx | 29 +---- yt/geometry/ramses_oct_container.pxd | 15 --- yt/geometry/ramses_oct_container.pyx | 157 ------------------------ 9 files changed, 167 insertions(+), 212 deletions(-) delete mode 100644 yt/geometry/ramses_oct_container.pxd delete mode 100644 yt/geometry/ramses_oct_container.pyx diff --git a/setup.py b/setup.py index 9123917a441..0c3f8fea9d0 100644 --- a/setup.py +++ b/setup.py @@ -97,11 +97,6 @@ def _compile( "yt/utilities/lib/tsearch.c"], include_dirs=["yt/utilities/lib"], libraries=std_libs), - Extension("yt.geometry.ramses_oct_container", - ["yt/geometry/ramses_oct_container.pyx", - "yt/utilities/lib/tsearch.c"], - include_dirs=["yt/utilities/lib"], - libraries=std_libs), Extension("yt.geometry.oct_visitors", ["yt/geometry/oct_visitors.pyx"], include_dirs=["yt/utilities/lib/"], diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 62b894541b8..3bf55d9e344 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -29,7 +29,7 @@ from .particle_handlers import get_particle_handlers from .field_handlers import get_field_handlers from yt.utilities.cython_fortran_utils import FortranFile as fpu -from yt.geometry.ramses_oct_container import \ +from yt.geometry.oct_container import \ RAMSESOctreeContainer from yt.arraytypes import blankRecordArray diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index e666b68deb3..1975a744eee 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -2,7 +2,7 @@ cimport cython cimport numpy as np import numpy as np from yt.utilities.cython_fortran_utils cimport FortranFile -from yt.geometry.ramses_oct_container cimport RAMSESOctreeContainer +from yt.geometry.oct_container cimport RAMSESOctreeContainer from yt.utilities.exceptions import YTIllDefinedAMRData ctypedef np.int32_t INT32_t diff --git a/yt/geometry/oct_container.pxd b/yt/geometry/oct_container.pxd index eb500d7fd56..5b6ece4c1b5 100644 --- a/yt/geometry/oct_container.pxd +++ b/yt/geometry/oct_container.pxd @@ -88,6 +88,8 @@ cdef class SparseOctreeContainer(OctreeContainer): cdef void key_to_ipos(self, np.int64_t key, np.int64_t pos[3]) cdef np.int64_t ipos_to_key(self, int pos[3]) nogil +cdef class RAMSESOctreeContainer(OctreeContainer): + pass cdef extern from "tsearch.h" nogil: void *tsearch(const void *key, void **rootp, diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 63d80994bb6..3c5e450cd3f 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -10,10 +10,10 @@ Oct container cimport cython cimport numpy as np import numpy as np -from selection_routines cimport SelectorObject +from selection_routines cimport SelectorObject, AlwaysSelector from libc.math cimport floor, ceil -cimport selection_routines -from yt.geometry.oct_visitors cimport OctPadded +from yt.geometry.oct_visitors cimport OctPadded, NeighbourCellVisitor, StoreIndex, NeighbourVisitor + ORDER_MAX = 20 _ORDER_MAX = ORDER_MAX @@ -76,7 +76,7 @@ cdef class OctreeContainer: header['right_edge'], over_refine = header['over_refine'], partial_coverage = header['partial_coverage']) # NOTE: We do not allow domain/file indices to be specified. - cdef SelectorObject selector = selection_routines.AlwaysSelector(None) + cdef SelectorObject selector = AlwaysSelector(None) cdef oct_visitors.LoadOctree visitor visitor = oct_visitors.LoadOctree(obj, -1) cdef int i, j, k, n @@ -471,7 +471,7 @@ cdef class OctreeContainer: right_edge = (self.DRE[0], self.DRE[1], self.DRE[2]), over_refine = self.oref, partial_coverage = self.partial_coverage) - cdef SelectorObject selector = selection_routines.AlwaysSelector(None) + cdef SelectorObject selector = AlwaysSelector(None) # domain_id = -1 here, because we want *every* oct cdef oct_visitors.StoreOctree visitor visitor = oct_visitors.StoreOctree(self, -1) @@ -743,8 +743,155 @@ cdef class OctreeContainer: else: dest[i + offset] = source[file_inds[i], cell_inds[i]] + def fill_index(self, SelectorObject selector = AlwaysSelector(None)): + """Get the on-file index of each cell""" + cdef StoreIndex visitor + + cdef np.int64_t[:, :, :, :] cell_inds, + + cell_inds = np.full((self.nocts, 2, 2, 2), -1, dtype=np.int64) + + visitor = StoreIndex(self, -1) + visitor.cell_inds = cell_inds + + self.visit_all_octs(selector, visitor) + + return np.asarray(cell_inds) + + def neighbours_in_direction(self, int idim, int direction, + np.int64_t[:, :, :, :] cell_inds): + """Return index on file of all neighbours in a given direction""" + cdef SelectorObject always_selector = AlwaysSelector(None) + + # Store the index of the neighbour + cdef NeighbourVisitor n_visitor + cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.full_like(cell_inds, -1) + n_visitor = NeighbourVisitor(self, -1) + n_visitor.idim = idim + n_visitor.direction = direction + n_visitor.cell_inds = cell_inds + n_visitor.neigh_cell_inds = neigh_cell_inds + n_visitor.octree = self + n_visitor.last = -1 + self.visit_all_octs(always_selector, n_visitor) + + return np.asarray(neigh_cell_inds) + + @cython.boundscheck(False) + @cython.wraparound(False) + def copy_neighbour_data(self, + np.int64_t[:] icell, np.int64_t[:] nicell, + np.float64_t[:, :] input, np.float64_t[:, :] output, + int N): + """Copy data from neighbouring cell into current one""" + cdef int i + + for i in range(N): + if nicell[i] > -1 and icell[i] > -1: + output[icell[i], :] = input[nicell[i], :] + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def fill_level_with_domain( + self, int level, + np.uint8_t[:] levels, + np.uint8_t[:] cell_inds, + np.int64_t[:] file_inds, + np.int32_t[:] domains, + dict dest_fields, dict source_fields, + np.int32_t domain, + np.int64_t offset = 0 + ): + """Similar to fill_level but accepts a domain argument. + + This is particularly useful for frontends that have buffer zones at CPU boundaries. + These buffer oct cells have a different domain than the local one and + are usually not read, but one has to read them e.g. to compute ghost zones. + """ + cdef np.ndarray[np.float64_t, ndim=2] source + cdef np.ndarray[np.float64_t, ndim=1] dest + cdef int i, count + + for key in dest_fields: + dest = dest_fields[key] + source = source_fields[key] + count = 0 + for i in range(levels.shape[0]): + if levels[i] != level or domains[i] != domain: continue + count += 1 + if file_inds[i] < 0: + dest[i + offset] = np.nan + else: + dest[i + offset] = source[file_inds[i], cell_inds[i]] + return count + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + def file_index_octs_with_ghost_zones( + self, SelectorObject selector, int domain_id, + int num_cells = -1): + """Similar as file_index_octs, but return as well the level, cell index, + file index and domain of the neighbouring cells. + + Arguments + --------- + selector : SelectorObject + The selector object. It is expected to select all cells for a selected oct. + domain_id : int + The domain to select. Set to -1 to select all domains. + num_cells : int, optional + The total number of cells (accounting for the ghost zones) + + Returns + ------- + shifted + + +---+---+---+---+ + | | | | | + |---+---+---+---| + | | x | x | | + |---+---+---+---| + | | x | x | | + |---+---+---+---| + | | | | | + +---+---+---+---+ + + """ + cdef np.int64_t i + cdef int num_octs + if num_cells < 0: + num_octs = selector.count_octs(self, domain_id) + num_cells = num_octs * 4**3 + cdef NeighbourCellVisitor visitor + + cdef np.ndarray[np.uint8_t, ndim=1] levels + cdef np.ndarray[np.uint8_t, ndim=1] cell_inds + cdef np.ndarray[np.int64_t, ndim=1] file_inds + cdef np.ndarray[np.int32_t, ndim=1] domains + levels = np.full(num_cells, 255, dtype="uint8") + file_inds = np.full(num_cells, -1, dtype="int64") + cell_inds = np.full(num_cells, 8, dtype="uint8") + domains = np.full(num_cells, -1, dtype="int32") + + visitor = NeighbourCellVisitor(self, -1) + # output: level, file_ind and cell_ind of the neighbouring cells + visitor.levels = levels + visitor.file_inds = file_inds + visitor.cell_inds = cell_inds + visitor.domains = domains + # direction to explore and extra parameters of the visitor + visitor.octree = self + visitor.last = -1 + + # Compute indices + self.visit_all_octs(selector, visitor) + + return levels, cell_inds, file_inds, domains + def finalize(self): - cdef SelectorObject selector = selection_routines.AlwaysSelector(None) + cdef SelectorObject selector = AlwaysSelector(None) cdef oct_visitors.AssignDomainInd visitor visitor = oct_visitors.AssignDomainInd(self, 1) self.visit_all_octs(selector, visitor) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 84199f27230..7caa7fef47e 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -176,9 +176,9 @@ cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.int32_t[:] neigh_domain cdef class NeighbourCellVisitor(BaseNeighbourVisitor): - cdef np.uint8_t[:] shifted_levels - cdef np.int64_t[:] shifted_file_inds - cdef np.uint8_t[:] shifted_cell_inds - cdef np.int32_t[:] neigh_domain + cdef np.uint8_t[:] levels + cdef np.int64_t[:] file_inds + cdef np.uint8_t[:] cell_inds + cdef np.int32_t[:] domains cdef void set_neighbour_info(self, Oct *o, int ishift[3]) \ No newline at end of file diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index ebb14f7cee7..6d47d4d363d 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -585,26 +585,9 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): neigh_file_ind = -1 neigh_cell_ind = 8 - self.shifted_levels[self.index] = neigh_level - self.shifted_file_inds[self.index] = neigh_file_ind - self.shifted_cell_inds[self.index] = neigh_cell_ind - self.neigh_domain[self.index] = neigh_domain - - self.index += 1# Compute mapping from index in domain to index in domain with buffer zones -cdef class DomainMapper(BaseNeighbourVisitor): - # Intended to map all octs ids to only octs in local domain - # Should be used in conjunction with OctreeSubsetSelector - def __init__(self, OctreeContainer octree, int domain_id): - super(DomainMapper, self).__init__(octree, domain_id) - self.index_in = 0 - - #@cython.boundscheck(False) - #@cython.wraparound(False) - #@cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - if self.last != o.domain_ind: - self.last = o.domain_ind - if self.marked[self.index]: - self.mapping[self.index] = self.index_in - self.index_in += 1 - self.index += 1 + self.levels[self.index] = neigh_level + self.file_inds[self.index] = neigh_file_ind + self.cell_inds[self.index] = neigh_cell_ind + self.domains[self.index] = neigh_domain + + self.index += 1 \ No newline at end of file diff --git a/yt/geometry/ramses_oct_container.pxd b/yt/geometry/ramses_oct_container.pxd deleted file mode 100644 index 6d5dddd8997..00000000000 --- a/yt/geometry/ramses_oct_container.pxd +++ /dev/null @@ -1,15 +0,0 @@ -""" -RAMSES Oct definitions file - - - - -""" -from oct_container cimport SparseOctreeContainer, OctInfo -from .oct_visitors cimport OctVisitor, Oct, cind -from yt.utilities.lib.fp_utils cimport * -cimport numpy as np - -cdef class RAMSESOctreeContainer(SparseOctreeContainer): - cdef Oct neighbour_in_direction(self, OctInfo *oinfo, np.int64_t *nneighbors, - Oct *o, bint periodicity[3]) diff --git a/yt/geometry/ramses_oct_container.pyx b/yt/geometry/ramses_oct_container.pyx deleted file mode 100644 index 3d80bf4ffc5..00000000000 --- a/yt/geometry/ramses_oct_container.pyx +++ /dev/null @@ -1,157 +0,0 @@ -cimport cython -from oct_visitors cimport StoreIndex, NeighbourVisitor, NeighbourCellVisitor -from selection_routines cimport SelectorObject, AlwaysSelector, OctreeSubsetSelector -cimport numpy as np -import numpy as np - - -cdef class RAMSESOctreeContainer(SparseOctreeContainer): - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef Oct neighbour_in_direction(self, OctInfo *oi, np.int64_t *nneighbours, Oct *o, - bint periodicity[3]): - pass - - def fill_index(self, SelectorObject selector = AlwaysSelector(None)): - # Get the on-file index of each cell - cdef StoreIndex visitor - - cdef np.int64_t[:, :, :, :] cell_inds, - - cell_inds = np.full((self.nocts, 2, 2, 2), -1, dtype=np.int64) - - visitor = StoreIndex(self, -1) - visitor.cell_inds = cell_inds - - self.visit_all_octs(selector, visitor) - - return np.asarray(cell_inds) - - def neighbours_in_direction(self, int idim, int direction, - np.int64_t[:, :, :, :] cell_inds): - """Return index on file of all neighbours in a given direction""" - cdef SelectorObject always_selector = AlwaysSelector(None) - - # Store the index of the neighbour - cdef NeighbourVisitor n_visitor - cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.full_like(cell_inds, -1) - n_visitor = NeighbourVisitor(self, -1) - n_visitor.idim = idim - n_visitor.direction = direction - n_visitor.cell_inds = cell_inds - n_visitor.neigh_cell_inds = neigh_cell_inds - n_visitor.octree = self - n_visitor.last = -1 - self.visit_all_octs(always_selector, n_visitor) - - return np.asarray(neigh_cell_inds) - - #@cython.boundscheck(False) - @cython.wraparound(False) - def copy_neighbour_data(self, - np.int64_t[:] icell, np.int64_t[:] nicell, - np.float64_t[:, :] input, np.float64_t[:, :] output, - int N,): - cdef int i - - for i in range(N): - if nicell[i] > -1 and icell[i] > -1: - output[icell[i], :] = input[nicell[i], :] - - def file_index_octs(self, SelectorObject selector, int domain_id, - num_cells = -1, spatial_offset=(0, 0, 0)): - - - cdef int i, idim, direction - cdef bint do_spatial_offset - cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds - cdef np.ndarray source_shifted - - do_spatial_offset = False - for i in range(3): - if spatial_offset[i] == 1 or spatial_offset[i] == -1: - idim = i - direction = spatial_offset[i] - if do_spatial_offset: - raise Exception( - 'ERROR: You can only specify one spatial offset direction, got [%s, %s, %s]!' % - (spatial_offset[0], spatial_offset[1], spatial_offset[2])) - do_spatial_offset = True - - if not do_spatial_offset: - return super(RAMSESOctreeContainer, self).file_index_octs( - selector, domain_id, num_cells) - else: - return self.file_index_octs_with_shift( - selector, domain_id, idim, direction, num_cells) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - def fill_level_with_domain( - self, int level, - np.uint8_t[:] levels, - np.uint8_t[:] cell_inds, - np.int64_t[:] file_inds, - np.int32_t[:] domains, - dict dest_fields, dict source_fields, - np.int32_t domain, - np.int64_t offset = 0 - ): - cdef np.ndarray[np.float64_t, ndim=2] source - cdef np.ndarray[np.float64_t, ndim=1] dest - cdef np.float64_t tmp - cdef int i, count - cdef str key - for key in dest_fields: - dest = dest_fields[key] - source = source_fields[key] - count = 0 - for i in range(levels.shape[0]): - if levels[i] != level or domains[i] != domain: continue - count += 1 - if file_inds[i] < 0: - dest[i + offset] = np.nan - else: - # print(f'\t{i}: Accessing source {file_inds[i]}:{cell_inds[i]} source.shape=({source.shape[0]},{source.shape[1]})') - tmp =source[file_inds[i], cell_inds[i]] - dest[i + offset] = tmp # source[file_inds[i], cell_inds[i]] - return count - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - def file_index_octs_with_ghost_zones( - self, SelectorObject selector, int domain_id, - int num_cells = -1): - cdef np.int64_t i - cdef int num_octs - if num_cells < 0: - num_octs = selector.count_octs(self, domain_id) - num_cells = num_octs * 4**3 - cdef NeighbourCellVisitor visitor - - cdef np.ndarray[np.uint8_t, ndim=1] shifted_levels - cdef np.ndarray[np.uint8_t, ndim=1] shifted_cell_inds - cdef np.ndarray[np.int64_t, ndim=1] shifted_file_inds - cdef np.ndarray[np.int32_t, ndim=1] neigh_domain - shifted_levels = np.full(num_cells, 255, dtype="uint8") - shifted_file_inds = np.full(num_cells, -1, dtype="int64") - shifted_cell_inds = np.full(num_cells, 8, dtype="uint8") - neigh_domain = np.full(num_cells, -1, dtype="int32") - - visitor = NeighbourCellVisitor(self, -1) - # output: level, file_ind and cell_ind of the neighbouring cells - visitor.shifted_levels = shifted_levels - visitor.shifted_file_inds = shifted_file_inds - visitor.shifted_cell_inds = shifted_cell_inds - visitor.neigh_domain = neigh_domain - # direction to explore and extra parameters of the visitor - visitor.octree = self - visitor.last = -1 - - # Compute indices - self.visit_all_octs(selector, visitor) - - return shifted_levels, shifted_cell_inds, shifted_file_inds, neigh_domain \ No newline at end of file From 4fa0a57a657f5d390ccb5e76214a1c0f6e3b529b Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 17:21:49 +0000 Subject: [PATCH 32/61] Remove unused parent accessor --- yt/geometry/oct_container.pyx | 4 ---- yt/geometry/oct_visitors.pxd | 1 - 2 files changed, 5 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 3c5e450cd3f..bdf9eab8798 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -566,7 +566,6 @@ cdef class OctreeContainer: cdef int ind[3] cdef int nb = 0 cdef Oct *cur - cdef Oct *parent cdef np.float64_t pp[3] cdef np.float64_t cp[3] cdef np.float64_t dds[3] @@ -575,7 +574,6 @@ cdef class OctreeContainer: cdef OctAllocationContainer *cont = self.domains.get_cont(curdom - 1) cdef int initial = cont.n_assigned cdef int in_boundary = 0 - parent = NULL # How do we bootstrap ourselves? for p in range(no): #for every oct we're trying to add find the @@ -608,12 +606,10 @@ cdef class OctreeContainer: ind[i] = 1 cp[i] += dds[i]/2.0 # Check if it has not been allocated - parent = cur cur = self.next_child(curdom, ind, cur) # Now we should be at the right level cur.domain = curdom cur.file_ind = p - cur.parent = parent return cont.n_assigned - initial + nb def allocate_domains(self, domain_counts): diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 7caa7fef47e..177412b6132 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -16,7 +16,6 @@ cdef struct Oct: np.int64_t domain_ind # index within the global set of domains np.int64_t domain # (opt) addl int index Oct **children # Up to 8 long - Oct *parent cdef struct OctInfo: np.float64_t left_edge[3] From 9b4e5c7aadbbe70eac2abcd05e1f87b2577222d6 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 17:22:00 +0000 Subject: [PATCH 33/61] RAMSES octree is sparse, restore this --- yt/geometry/oct_container.pxd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/geometry/oct_container.pxd b/yt/geometry/oct_container.pxd index 5b6ece4c1b5..75ada2b97e7 100644 --- a/yt/geometry/oct_container.pxd +++ b/yt/geometry/oct_container.pxd @@ -88,7 +88,7 @@ cdef class SparseOctreeContainer(OctreeContainer): cdef void key_to_ipos(self, np.int64_t key, np.int64_t pos[3]) cdef np.int64_t ipos_to_key(self, int pos[3]) nogil -cdef class RAMSESOctreeContainer(OctreeContainer): +cdef class RAMSESOctreeContainer(SparseOctreeContainer): pass cdef extern from "tsearch.h" nogil: From b5b2bfa134dfee2c4c5a4be4aee25f899094ab09 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 30 Jan 2020 17:48:50 +0000 Subject: [PATCH 34/61] Update docstring --- yt/geometry/oct_container.pyx | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index bdf9eab8798..2592573ab43 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -800,7 +800,7 @@ cdef class OctreeContainer: np.int64_t offset = 0 ): """Similar to fill_level but accepts a domain argument. - + This is particularly useful for frontends that have buffer zones at CPU boundaries. These buffer oct cells have a different domain than the local one and are usually not read, but one has to read them e.g. to compute ghost zones. @@ -828,7 +828,7 @@ cdef class OctreeContainer: def file_index_octs_with_ghost_zones( self, SelectorObject selector, int domain_id, int num_cells = -1): - """Similar as file_index_octs, but return as well the level, cell index, + """Similar as file_index_octs, but return as well the level, cell index, file index and domain of the neighbouring cells. Arguments @@ -838,11 +838,30 @@ cdef class OctreeContainer: domain_id : int The domain to select. Set to -1 to select all domains. num_cells : int, optional - The total number of cells (accounting for the ghost zones) + The total number of cells (accounting for a 1-cell thick ghost zone layer). Returns ------- - shifted + levels : uint8, shape (num_cells,) + The level of each cell of the super oct + cell_inds : uint8, shape (num_cells, ) + The index of each cell of the super oct within its own oct + file_inds : int64, shape (num_cells, ) + The on-file position of the cell. See notes below. + domains : int32, shape (num_cells) + The domain to which the cells belongs. See notes below. + + Notes + ----- + + The algorithm constructs a "super-oct" around each oct (see sketch below, + where the original oct cells are marked with an x). + + Note that for sparse octrees (such as RAMSES'), the neighbouring cells + may belong to another domain (this is stored in `domains`). If the dataset + provides buffer zones between domains (such as RAMSES), this may be stored + locally and can be accessed directly. + +---+---+---+---+ | | | | | From c3262c74a8ca7a2cacb77a4a85035bcc2f303d88 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Sat, 8 Feb 2020 14:14:02 +0000 Subject: [PATCH 35/61] Only query fcoords w/o ghost zones --- yt/frontends/ramses/data_structures.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 3bf55d9e344..1473d5e1a22 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -261,10 +261,9 @@ def fwidth(self): @property def fcoords(self): - fcoords = super(RAMSESDomainSubset, self).fcoords num_ghost_zones = self._num_ghost_zones if num_ghost_zones == 0: - return fcoords + return super(RAMSESDomainSubset, self).fcoords fcoords_base = self._base_domain.fcoords oct_selector = OctreeSubsetSelector(self) From 5e5af6fb700fa379ee62316b61f1b08e2a71e358 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 11 Feb 2020 10:57:04 +0000 Subject: [PATCH 36/61] Assume periodicity --- yt/geometry/oct_visitors.pyx | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 6d47d4d363d..454b18c7ed6 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -372,6 +372,12 @@ cdef class BaseNeighbourVisitor(OctVisitor): else: fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] + # Assuming periodicity + if fcoords[i] < 0: + fcoords[i] += 1 + elif fcoords[i] > 1: + fcoords[i] -= 1 + # Use octree to find neighbour neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) @@ -507,6 +513,11 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): for i in range(3): c = (self.pos[i] << self.oref) fcoords[i] = (c + 0.5 + ishift[i]) * dx / self.octree.nn[i] + # Assuming periodicity + if fcoords[i] < 0: + fcoords[i] += 1 + elif fcoords[i] > 1: + fcoords[i] -= 1 local_oct &= (0 <= ishift[i] <= 1) other_oct = not local_oct From 7266517957cb9aa0c75ec8f263a541b7759e7e5a Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 11 Feb 2020 16:47:04 +0000 Subject: [PATCH 37/61] Add warning --- yt/frontends/ramses/data_structures.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 1473d5e1a22..4a65d43acfe 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -2,6 +2,7 @@ import numpy as np import stat import weakref +import warnings from collections import defaultdict from glob import glob @@ -193,6 +194,8 @@ def __init__(self, base_region, domain, ds, over_refine_factor=1, num_ghost_zone self._base_grid = base_grid if num_ghost_zones > 0: + if not all(ds.periodicity): + warnings.warn('Ghost zones will wrongly assume the domain to be periodic.') # Create a base domain *with no self._base_domain.fwidth base_domain = RAMSESDomainSubset(ds.all_data(), domain, ds, over_refine_factor) self._base_domain = base_domain From 9fc52dfd595864f1bfb0f469f08c92bbc060db51 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 12 Feb 2020 16:16:09 +0000 Subject: [PATCH 38/61] Vastly simplify code --- yt/frontends/ramses/data_structures.py | 33 ++--- yt/geometry/oct_container.pyx | 64 ++++++--- yt/geometry/oct_visitors.pxd | 10 +- yt/geometry/oct_visitors.pyx | 187 ++++++------------------- 4 files changed, 102 insertions(+), 192 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 4a65d43acfe..e14918c8d1e 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -268,36 +268,21 @@ def fcoords(self): if num_ghost_zones == 0: return super(RAMSESDomainSubset, self).fcoords - fcoords_base = self._base_domain.fcoords - oct_selector = OctreeSubsetSelector(self) oh = self.oct_handler - n_oct = fcoords_base.size // 3 // 8 - new_fcoords = np.full((n_oct, 4, 4, 4, 3), np.nan) + indices = oh.fill_index(self.selector).reshape(-1, 8) + oinds, cinds = oh.fill_octcellindex_neighbours(self.selector) - icell = oh.fill_index(oct_selector) - Ncell = icell.size + oinds = oinds.reshape(-1, 64) + cinds = cinds.reshape(-1, 64) - for idim in range(3): - for idir in (-1, 1): - ishift_all = [1, 1, 1] - ishift_all[idim] += idir + inds = indices[oinds, cinds] - nicell = oh.neighbours_in_direction(idim, idir, icell).reshape(-1) - - tmp = np.full((n_oct * 2 * 2 * 2, 3), np.nan) - - oh.copy_neighbour_data( - icell.reshape(-1), nicell, - fcoords_base, tmp, Ncell) - - _slice = tuple([slice(None)] + - [slice(i, i+2) for i in ishift_all] + - [slice(None)]) - new_fcoords[_slice] = tmp.reshape(n_oct, 2, 2, 2, -1) - new_fcoords = self.ds.arr(new_fcoords, fcoords_base.units) - return new_fcoords + fcoords = self.ds.arr( + oh.fcoords(self.selector)[inds].reshape(-1, 3), + 'unitary') + return fcoords def fill(self, fd, fields, selector, file_handler): if self._num_ghost_zones == 0: diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 2592573ab43..ae7421f73f4 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -12,7 +12,7 @@ cimport numpy as np import numpy as np from selection_routines cimport SelectorObject, AlwaysSelector from libc.math cimport floor, ceil -from yt.geometry.oct_visitors cimport OctPadded, NeighbourCellVisitor, StoreIndex, NeighbourVisitor +from yt.geometry.oct_visitors cimport OctPadded, NeighbourCellVisitor, StoreIndex, NeighbourCellIndexVisitor ORDER_MAX = 20 @@ -743,7 +743,7 @@ cdef class OctreeContainer: """Get the on-file index of each cell""" cdef StoreIndex visitor - cdef np.int64_t[:, :, :, :] cell_inds, + cdef np.int64_t[:, :, :, :] cell_inds cell_inds = np.full((self.nocts, 2, 2, 2), -1, dtype=np.int64) @@ -754,24 +754,48 @@ cdef class OctreeContainer: return np.asarray(cell_inds) - def neighbours_in_direction(self, int idim, int direction, - np.int64_t[:, :, :, :] cell_inds): - """Return index on file of all neighbours in a given direction""" - cdef SelectorObject always_selector = AlwaysSelector(None) - - # Store the index of the neighbour - cdef NeighbourVisitor n_visitor - cdef np.ndarray[np.int64_t, ndim=4] neigh_cell_inds = np.full_like(cell_inds, -1) - n_visitor = NeighbourVisitor(self, -1) - n_visitor.idim = idim - n_visitor.direction = direction - n_visitor.cell_inds = cell_inds - n_visitor.neigh_cell_inds = neigh_cell_inds - n_visitor.octree = self - n_visitor.last = -1 - self.visit_all_octs(always_selector, n_visitor) - - return np.asarray(neigh_cell_inds) + def fill_octcellindex_neighbours(self, SelectorObject selector, int num_octs = -1, domain_id = -1): + """Compute the oct and cell indices of all the cells within all selected octs, extended + by one cell in all directions (for ghost zones computations). + + Parameters + ---------- + selector : SelectorObject + Selector for the octs to compute neighbour of + num_octs : int, optional + The number of octs to read in + domain_id : int, optional + The domain to perform the selection over + + Returns + ------- + oct_inds : int64 ndarray (nocts*8, ) + The on-domain index of the octs containing each cell + cell_inds : uint8 array (nocts*8, ) + The index of the cell in its parent oct + + Note + ---- + oct_inds/cell_inds + """ + if num_octs == -1: + num_octs = selector.count_octs(self, domain_id) + + cdef NeighbourCellIndexVisitor visitor + + cdef np.uint8_t[:] cell_inds + cdef np.int64_t[:] oct_inds + + cell_inds = np.full(num_octs*4**3, 8, dtype=np.uint8) + oct_inds = np.full(num_octs*4**3, -1, dtype=np.int64) + + visitor = NeighbourCellIndexVisitor(self, -1) + visitor.cell_inds = cell_inds + visitor.domain_inds = oct_inds + + self.visit_all_octs(selector, visitor) + + return np.asarray(oct_inds), np.asarray(cell_inds) @cython.boundscheck(False) @cython.wraparound(False) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 177412b6132..791f2c82fe1 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -157,16 +157,15 @@ cdef class BaseNeighbourVisitor(OctVisitor): cdef OctreeContainer octree cdef OctInfo oi - cdef void set_neighbour_oct(self, Oct* o) - cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected) + cdef void set_neighbour_info(self, Oct *o, int ishift[3]) cdef inline np.uint8_t neighbour_rind(self): cdef int d = (1 << self.oref) return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) -cdef class NeighbourVisitor(BaseNeighbourVisitor): - cdef np.int64_t[:,:,:,:] cell_inds - cdef np.int64_t[:,:,:,:] neigh_cell_inds +cdef class NeighbourCellIndexVisitor(BaseNeighbourVisitor): + cdef np.uint8_t[:] cell_inds + cdef np.int64_t[:] domain_inds cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): cdef np.uint8_t[:] shifted_levels @@ -180,4 +179,3 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): cdef np.uint8_t[:] cell_inds cdef np.int32_t[:] domains - cdef void set_neighbour_info(self, Oct *o, int ishift[3]) \ No newline at end of file diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 454b18c7ed6..e35c8c14c8d 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -354,147 +354,6 @@ cdef class BaseNeighbourVisitor(OctVisitor): self.neigh_ind = np.zeros(3, np.int8) super(BaseNeighbourVisitor, self).__init__(octree, domain_id) - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void set_neighbour_oct(self, Oct *o): - cdef int i - cdef np.float64_t c, dx - cdef np.int64_t ipos - cdef np.float64_t fcoords[3] - cdef Oct *neighbour - dx = 1.0 / ((1 << self.oref) << self.level) - # Compute position of neighbouring cell - for i in range(3): - c = ((self.pos[i] << self.oref) + self.ind[i]) - if i == self.idim: - fcoords[i] = (c + 0.5 + 2*self.direction) * dx / self.octree.nn[i] - else: - fcoords[i] = (c + 0.5) * dx / self.octree.nn[i] - - # Assuming periodicity - if fcoords[i] < 0: - fcoords[i] += 1 - elif fcoords[i] > 1: - fcoords[i] -= 1 - - # Use octree to find neighbour - neighbour = self.octree.get(fcoords, &self.oi, max_level=self.level) - - # Extra step - compute cell position in neighbouring oct (and store in oi.ipos) - if self.oi.level == self.level - 1: - for i in range(3): - ipos = (((self.pos[i] << self.oref) + self.ind[i])) >> 1 - if i == self.idim: - ipos += self.direction - if (self.oi.ipos[i] << 1) == ipos: - self.oi.ipos[i] = 0 - else: - self.oi.ipos[i] = 1 - self.neighbour = neighbour - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - @cython.cdivision(True) - cdef void get_neighbour_cell_index(self, Oct* o, np.uint8_t selected): - cdef int i - cdef bint other_oct # True if the neighbouring cell lies in another oct - cdef np.int64_t cell_ind - - # Compute information about neighbour once per oct - if self.last != o.domain_ind: - self.set_neighbour_oct(o) - self.last = o.domain_ind - - # Note that we provide an index even if the cell is not selected. - # if selected == 0: return -1 - # Index of neighbouring cell within its oct - for i in range(3): - if i == self.idim: - self.neigh_ind[i] = (self.ind[i] + self.direction) - other_oct = self.neigh_ind[i] < 0 or self.neigh_ind[i] > 1 - if other_oct: - # trick here: we want modulo with positive remainder, but neigh_ind may be negative so cast - # it to unsigned int *before* applying modulo. - self.neigh_ind[i] = (self.neigh_ind[i]) % 2 - else: - self.neigh_ind[i] = self.ind[i] - - self.other_oct = other_oct - if other_oct: - if self.neighbour != NULL: - if self.oi.level == self.level - 1: - # Position within neighbouring oct is stored in oi.ipos - for i in range(3): - self.neigh_ind[i] = self.oi.ipos[i] - elif self.oi.level != self.level: - print('This should not happen! %s %s' % (self.oi.level, self.level)) - self.neighbour = NULL - -# Store neighbouring cell index in current cell -cdef class NeighbourVisitor(BaseNeighbourVisitor): - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - cdef np.int64_t cell_ind - cdef Oct *neighbour_oct - cdef bint ok - - self.get_neighbour_cell_index(o, selected) - if not self.other_oct: - neighbour_oct = o - ok = True - elif self.neighbour != NULL: - neighbour_oct = self.neighbour - ok = True - else: - ok = False - - if ok: - cell_ind = self.cell_inds[neighbour_oct.domain_ind, self.neigh_ind[2], self.neigh_ind[1], self.neigh_ind[0]] - else: - cell_ind = -1 - - self.neigh_cell_inds[o.domain_ind, self.ind[2], self.ind[1], self.ind[0]] = cell_ind - -# Store file position + cell of neighbouring cell in current cell -cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.initializedcheck(False) - cdef void visit(self, Oct* o, np.uint8_t selected): - cdef np.int64_t neigh_file_ind - cdef np.uint8_t neigh_cell_ind - cdef np.int32_t neigh_domain - - if selected == 0: return - # Note: only selected items have an index - self.get_neighbour_cell_index(o, selected) - if not self.other_oct: - neigh_domain = o.domain - neigh_file_ind = o.file_ind - neigh_cell_ind = self.neighbour_rind() - elif self.neighbour != NULL: - neigh_domain = self.neighbour.domain - neigh_file_ind = self.neighbour.file_ind - neigh_cell_ind = self.neighbour_rind() - else: - neigh_domain = -1 - neigh_file_ind = -1 - neigh_cell_ind = 8 - - self.shifted_levels[self.index] = self.level - # Note: we store the local level, not the remote one - self.shifted_file_inds[self.index] = neigh_file_ind - self.shifted_cell_inds[self.index] = neigh_cell_ind - self.neigh_domain[self.index] = neigh_domain - - self.index += 1 - -# Store file position + cell of neighbouring cell in current cell -cdef class NeighbourCellVisitor(BaseNeighbourVisitor): @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) @@ -555,6 +414,49 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): print('This should not happen! %s %s' % (self.oi.level, self.level)) self.neighbour = NULL +# Store neighbouring cell index in current cell +cdef class NeighbourCellIndexVisitor(BaseNeighbourVisitor): + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef void visit(self, Oct* o, np.uint8_t selected): + cdef int i, j, k + cdef int ishift[3] + cdef np.uint8_t neigh_cell_ind + cdef np.int64_t neigh_domain_ind + if selected == 0: return + # Work at oct level + if self.last == o.domain_ind: return + + self.last = o.domain_ind + + # Loop over cells in and directly around oct + for i in range(-1, 3): + ishift[0] = i + for j in range(-1, 3): + ishift[1] = j + for k in range(-1, 3): + ishift[2] = k + self.set_neighbour_info(o, ishift) + + if not self.other_oct: + neigh_domain_ind = o.domain_ind + neigh_cell_ind = self.neighbour_rind() + elif self.neighbour != NULL: + neigh_domain_ind = self.neighbour.domain_ind + neigh_cell_ind = self.neighbour_rind() + else: + neigh_domain_ind = -1 + neigh_cell_ind = 8 + + self.cell_inds[self.index] = neigh_cell_ind + self.domain_inds[self.index] = neigh_domain_ind + + self.index += 1 + +# Store file position + cell of neighbouring cell in current cell +cdef class NeighbourCellVisitor(BaseNeighbourVisitor): + @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) @@ -601,4 +503,5 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): self.cell_inds[self.index] = neigh_cell_ind self.domains[self.index] = neigh_domain - self.index += 1 \ No newline at end of file + self.index += 1 + From 8ec03540d3db06f8cb094e98a980108126ab7a0a Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 12 Feb 2020 16:16:21 +0000 Subject: [PATCH 39/61] Remove useless import --- yt/frontends/ramses/data_structures.py | 1 - 1 file changed, 1 deletion(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index e14918c8d1e..b71ccef109f 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -11,7 +11,6 @@ setdefaultattr from yt.geometry.oct_geometry_handler import \ OctreeIndex -from yt.geometry.selection_routines import OctreeSubsetSelector from yt.geometry.geometry_handler import \ YTDataChunk From aa2605dfc406f3f6da4d083d55fa45ec2935c27d Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 12 Feb 2020 16:25:22 +0000 Subject: [PATCH 40/61] Remove unused visitor --- yt/geometry/oct_visitors.pxd | 6 ------ 1 file changed, 6 deletions(-) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 791f2c82fe1..7a77a653299 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -167,12 +167,6 @@ cdef class NeighbourCellIndexVisitor(BaseNeighbourVisitor): cdef np.uint8_t[:] cell_inds cdef np.int64_t[:] domain_inds -cdef class FillFileIndicesRNeighbour(BaseNeighbourVisitor): - cdef np.uint8_t[:] shifted_levels - cdef np.int64_t[:] shifted_file_inds - cdef np.uint8_t[:] shifted_cell_inds - cdef np.int32_t[:] neigh_domain - cdef class NeighbourCellVisitor(BaseNeighbourVisitor): cdef np.uint8_t[:] levels cdef np.int64_t[:] file_inds From 11285a006ae8ce726ef49ff365681af832c19e7a Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 19 Feb 2020 10:54:22 +0000 Subject: [PATCH 41/61] Remove useless method --- yt/geometry/oct_container.pyx | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index ae7421f73f4..36d4b8dfb51 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -797,19 +797,6 @@ cdef class OctreeContainer: return np.asarray(oct_inds), np.asarray(cell_inds) - @cython.boundscheck(False) - @cython.wraparound(False) - def copy_neighbour_data(self, - np.int64_t[:] icell, np.int64_t[:] nicell, - np.float64_t[:, :] input, np.float64_t[:, :] output, - int N): - """Copy data from neighbouring cell into current one""" - cdef int i - - for i in range(N): - if nicell[i] > -1 and icell[i] > -1: - output[icell[i], :] = input[nicell[i], :] - @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) From 6fde2e5823ba9d76065667c7c764a0bf45564403 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 29 May 2020 09:42:01 +0100 Subject: [PATCH 42/61] Address comments --- yt/frontends/ramses/data_structures.py | 7 ++++--- yt/frontends/ramses/io_utils.pyx | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index b71ccef109f..5dbcf6e0319 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -254,7 +254,7 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo def fwidth(self): fwidth = super(RAMSESDomainSubset, self).fwidth if self._num_ghost_zones > 0: - fwidth = super(RAMSESDomainSubset, self).fwidth.reshape(-1, 8, 3) + fwidth = fwidth.reshape(-1, 8, 3) n_oct = fwidth.shape[0] new_fwidth = np.zeros((n_oct, self.nz**3, 3), dtype=fwidth.dtype) new_fwidth[:, :, :] = fwidth[:, 0:1, :] @@ -344,14 +344,15 @@ def _detect_output_fields(self): self.field_list = self.particle_field_list + self.fluid_field_list def _identify_base_chunk(self, dobj): - ngz = dobj._num_ghost_zones if getattr(dobj, "_chunk_info", None) is None: domains = [dom for dom in self.domains if dom.included(dobj.selector)] base_region = getattr(dobj, "base_region", dobj) if len(domains) > 1: mylog.debug("Identified %s intersecting domains", len(domains)) - subsets = [RAMSESDomainSubset(base_region, domain, self.dataset, num_ghost_zones=ngz) + subsets = [RAMSESDomainSubset( + base_region, domain, self.dataset, + num_ghost_zones=dobj._num_ghost_zones) for domain in domains] dobj._chunk_info = subsets dobj._current_chunk = list(self._chunk_all(dobj))[0] diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 1975a744eee..d745655cbfe 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -75,7 +75,7 @@ def read_amr(FortranFile f, dict headers, if n > 0: max_level = max(ilevel - min_level, max_level) if n != ng: - raise Exception('Expected %s octs, got %s' % (ng, n)) + raise ValueError('Expected %s octs, got %s' % (ng, n)) return max_level From b91838990279b03466f5c0caa5eab576eaae256e Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 14:23:38 +0100 Subject: [PATCH 43/61] Fix mistake --- yt/fields/fluid_fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index ce927240a2a..f9abc0785d7 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -218,7 +218,7 @@ def func(field, data): f = field_data[slice_3dr]/ds[slice_3d] f -= field_data[slice_3dl]/ds[slice_3d] new_field = np.zeros_like(data[grad_field], dtype=np.float64) - new_field = data.ds.arr(new_field, vr.units / dt.units) + new_field = data.ds.arr(new_field, field_data.units / ds.units) new_field[slice_3d] = f if block_reorder == 'F': From 7978ddcaeb3844b0cddc0cf94e990275fd479a71 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Wed, 17 Jun 2020 14:23:50 +0100 Subject: [PATCH 44/61] No need to format warning --- yt/fields/fluid_fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index f9abc0785d7..25134b65f40 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -189,7 +189,7 @@ def setup_gradient_fields(registry, grad_field, field_units, slice_info = None): geom = registry.ds.geometry if is_curvilinear(geom): - mylog.warning("In %s geometry, gradient fields may contain artifacts near cartesian axes." % geom) + mylog.warning("In %s geometry, gradient fields may contain artifacts near cartesian axes.", geom) assert(isinstance(grad_field, tuple)) ftype, fname = grad_field From 86a314e51d077146068ee9f6d34f9f0b1be1ee43 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 23 Jun 2020 10:01:30 +0100 Subject: [PATCH 45/61] regenerate answers --- tests/tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tests.yaml b/tests/tests.yaml index 72f35e51cea..8df95daf9ec 100644 --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -105,7 +105,7 @@ answer_tests: - yt/frontends/boxlib/tests/test_outputs.py:test_NyxDataset - yt/frontends/boxlib/tests/test_outputs.py:test_WarpXDataset - local_ramses_002: + local_ramses_003: - yt/frontends/ramses/tests/test_outputs.py local_ytdata_006: From 6d4adf8f9d6fea44c8dfb224247c20d7dc113ebe Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Tue, 23 Jun 2020 15:12:20 +0100 Subject: [PATCH 46/61] Support non periodic boundaries? --- yt/frontends/ramses/io_utils.pyx | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index d745655cbfe..edc17da082e 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -74,8 +74,6 @@ def read_amr(FortranFile f, dict headers, count_boundary = 1) if n > 0: max_level = max(ilevel - min_level, max_level) - if n != ng: - raise ValueError('Expected %s octs, got %s' % (ng, n)) return max_level @@ -103,8 +101,8 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n skip_len = twotondim * nvar # It goes: level, CPU, 8-variable (1 oct) - offset = np.full((ncpu, n_levels), -1, dtype=np.int64) - level_count = np.zeros((ncpu, n_levels), dtype=np.int64) + offset = np.full((ncpu_and_bound, n_levels), -1, dtype=np.int64) + level_count = np.zeros((ncpu_and_bound, n_levels), dtype=np.int64) cdef np.int64_t[:,:] level_count_view = level_count cdef np.int64_t[:,:] offset_view = offset From 6c318dcc263ea9ef5edc8c8881f57bc6f78468df Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:38:37 +0100 Subject: [PATCH 47/61] Add kwa to select calls --- yt/data_objects/data_containers.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/yt/data_objects/data_containers.py b/yt/data_objects/data_containers.py index 4b9270c60ac..f2f98892bfa 100644 --- a/yt/data_objects/data_containers.py +++ b/yt/data_objects/data_containers.py @@ -340,7 +340,11 @@ def _generate_spatial_fluid(self, field, ngz): outputs.append(rv) ind = 0 # Does this work with mesh? with o._activate_cache(): - ind += o.select(self.selector, self[field], rv, ind) + ind += o.select( + self.selector, + source=self[field], + dest=rv, + offset=ind) else: chunks = self.index._chunk(self, "spatial", ngz = ngz) for i, chunk in enumerate(chunks): @@ -355,8 +359,9 @@ def _generate_spatial_fluid(self, field, ngz): data = gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz] ind += wogz.select( self.selector, - data, - rv, ind) + source=data, + dest=rv, + offset=ind) if accumulate: rv = uconcatenate(outputs) return rv From 236134b5d90bb5cd17bb50fbde8e57a33f742aae Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:42:16 +0100 Subject: [PATCH 48/61] =?UTF-8?q?ds=20=E2=86=92=20dx?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- yt/fields/fluid_fields.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 25134b65f40..b421670c7ab 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -210,15 +210,15 @@ def func(field, data): field_data = data[grad_field].swapaxes(0, 2) else: field_data = data[grad_field] - ds = div_fac * data[ftype, "d%s" % ax] + dx = div_fac * data[ftype, f"d{ax}"] if ax == "theta": - ds *= data[ftype, "r"] + dx *= data[ftype, "r"] if ax == "phi": - ds *= data[ftype, "r"] * np.sin(data[ftype, "theta"]) - f = field_data[slice_3dr]/ds[slice_3d] - f -= field_data[slice_3dl]/ds[slice_3d] + dx *= data[ftype, "r"] * np.sin(data[ftype, "theta"]) + f = field_data[slice_3dr]/dx[slice_3d] + f -= field_data[slice_3dl]/dx[slice_3d] new_field = np.zeros_like(data[grad_field], dtype=np.float64) - new_field = data.ds.arr(new_field, field_data.units / ds.units) + new_field = data.ds.arr(new_field, field_data.units / dx.units) new_field[slice_3d] = f if block_reorder == 'F': From 550e887213150398aba68530831289ad8f729352 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:45:36 +0100 Subject: [PATCH 49/61] Return "arr" not "tmp" variable --- yt/fields/geometric_fields.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yt/fields/geometric_fields.py b/yt/fields/geometric_fields.py index d6c73b3b3b7..3a8ec9d55d3 100644 --- a/yt/fields/geometric_fields.py +++ b/yt/fields/geometric_fields.py @@ -93,11 +93,11 @@ def _zeros(field, data): def _ones(field, data): """Returns one for all cells""" - arr = np.ones(data.ires.shape, dtype="float64") - tmp = data.apply_units(arr, field.units) + tmp = np.ones(data.ires.shape, dtype="float64") + arr = data.apply_units(tmp, field.units) if data._spatial: - return data._reshape_vals(tmp) - return tmp + return data._reshape_vals(arr) + return arr registry.add_field(("index", "ones"), sampling_type="cell", From 0fc513160d6a7fdd63cadc34790e50149f54f049 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:46:11 +0100 Subject: [PATCH 50/61] Each line for its argument --- yt/geometry/oct_container.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 36d4b8dfb51..8b9e88e055a 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -806,7 +806,8 @@ cdef class OctreeContainer: np.uint8_t[:] cell_inds, np.int64_t[:] file_inds, np.int32_t[:] domains, - dict dest_fields, dict source_fields, + dict dest_fields, + dict source_fields, np.int32_t domain, np.int64_t offset = 0 ): From f7d9e4ae5ac57f064ccc511ed077df52b7df8dca Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:50:33 +0100 Subject: [PATCH 51/61] Minor nits in docstrings --- yt/geometry/oct_container.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 8b9e88e055a..4ea1e38da4c 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -840,13 +840,13 @@ cdef class OctreeContainer: def file_index_octs_with_ghost_zones( self, SelectorObject selector, int domain_id, int num_cells = -1): - """Similar as file_index_octs, but return as well the level, cell index, - file index and domain of the neighbouring cells. + """Similar as file_index_octs, but returns the level, cell index, + file index and domain of the neighbouring cells as well. Arguments --------- selector : SelectorObject - The selector object. It is expected to select all cells for a selected oct. + The selector object. It is expected to select all cells for a given oct. domain_id : int The domain to select. Set to -1 to select all domains. num_cells : int, optional From 9ca0d0e13af6a4348660fa74b9ba17154680d363 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:50:59 +0100 Subject: [PATCH 52/61] Using python-friendly syntax --- yt/geometry/oct_container.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index 4ea1e38da4c..ea8c64058f8 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -825,8 +825,8 @@ cdef class OctreeContainer: dest = dest_fields[key] source = source_fields[key] count = 0 - for i in range(levels.shape[0]): - if levels[i] != level or domains[i] != domain: continue + for i, (lev, dom) in enumerate(zip(levels, domains)): + if lev != level or dom != domain: continue count += 1 if file_inds[i] < 0: dest[i + offset] = np.nan From a92d4778323beea20ec798cb43b6157173208dbc Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 10:55:40 +0100 Subject: [PATCH 53/61] Remove duplicates in testing file --- yt/frontends/ramses/tests/test_outputs.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py index 1538344e694..38d491cab71 100644 --- a/yt/frontends/ramses/tests/test_outputs.py +++ b/yt/frontends/ramses/tests/test_outputs.py @@ -399,7 +399,6 @@ def check_unit(array, unit): check_unit(ds.r[('gas','Electron_number_density')],'cm**(-3)') -ramses_rt = "ramses_rt_00088/output_00088/info_00088.txt" @requires_file(ramses_rt) def test_ramses_mixed_files(): # Test that one can use derived fields that depend on different @@ -462,7 +461,6 @@ def test_magnetic_field_aliasing(): assert ('gas',field) in ds.derived_field_list ad[('gas',field)] -output_00080 = "output_00080/info_00080.txt" @requires_file(output_00080) def test_field_accession(): ds = yt.load(output_00080) From 0bf76b2808ab00471a1129c96ca34ec782f65d6f Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 11:00:02 +0100 Subject: [PATCH 54/61] _block_reorder -> _block_order --- yt/data_objects/octree_subset.py | 6 +++--- yt/fields/fluid_fields.py | 8 +++++--- yt/frontends/ramses/data_structures.py | 2 +- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/yt/data_objects/octree_subset.py b/yt/data_objects/octree_subset.py index dbce321e657..c66359bf943 100644 --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -35,7 +35,7 @@ class OctreeSubset(YTSelectionContainer): _con_args = ('base_region', 'domain', 'ds') _domain_offset = 0 _cell_count = -1 - _block_reorder = None + _block_order = 'C' def __init__(self, base_region, domain, ds, over_refine_factor = 1, num_ghost_zones = 0): super(OctreeSubset, self).__init__(ds, None) @@ -456,8 +456,8 @@ def __init__(self, ind, block_slice): def __getitem__(self, key): bs = self.block_slice rv = bs.octree_subset[key][:,:,:,self.ind].T - if bs.octree_subset._block_reorder: - rv = rv.copy(order=bs.octree_subset._block_reorder) + if bs.octree_subset._block_order: + rv = rv.copy(order=bs.octree_subset._block_order) return rv @property diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index b421670c7ab..44af8c961ad 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -205,8 +205,10 @@ def grad_func(axi, ax): slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:] slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:] def func(field, data): - block_reorder = getattr(data, '_block_reorder', None) - if block_reorder == 'F': + block_order = data._block_order + if block_order == 'F': + # Fortran-ordering: we need to swap axes here and + # reswap below field_data = data[grad_field].swapaxes(0, 2) else: field_data = data[grad_field] @@ -221,7 +223,7 @@ def func(field, data): new_field = data.ds.arr(new_field, field_data.units / dx.units) new_field[slice_3d] = f - if block_reorder == 'F': + if block_order == 'F': new_field = new_field.swapaxes(0, 2) return new_field diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 1e24bdfbc7f..217f5192621 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -181,7 +181,7 @@ def included(self, selector): class RAMSESDomainSubset(OctreeSubset): _domain_offset = 1 - _block_reorder = "F" + _block_order = "F" _base_domain = None From 19058d8c88fbde284dae230232acb24043e66929 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Thu, 2 Jul 2020 11:03:43 +0100 Subject: [PATCH 55/61] Field detector does not have a block_order argument --- yt/fields/fluid_fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index 44af8c961ad..ddb2d0fe257 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -205,7 +205,7 @@ def grad_func(axi, ax): slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:] slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:] def func(field, data): - block_order = data._block_order + block_order = getattr(data, '_block_order', 'C') if block_order == 'F': # Fortran-ordering: we need to swap axes here and # reswap below From 6d2b43978f53e8425953270f20418ef5c0229483 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 3 Jul 2020 09:30:56 +0100 Subject: [PATCH 56/61] Reorder if necessary --- yt/data_objects/octree_subset.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/yt/data_objects/octree_subset.py b/yt/data_objects/octree_subset.py index c66359bf943..cd006c1be5b 100644 --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -455,9 +455,9 @@ def __init__(self, ind, block_slice): def __getitem__(self, key): bs = self.block_slice - rv = bs.octree_subset[key][:,:,:,self.ind].T - if bs.octree_subset._block_order: - rv = rv.copy(order=bs.octree_subset._block_order) + rv = np.require( + bs.octree_subset[key][:,:,:,self.ind].T, + requirements=bs.octree_subset._block_order) return rv @property From 3c8b75ccf6d2a84588c487ccdbdd68fda6cfc631 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 17 Jul 2020 10:01:31 +0100 Subject: [PATCH 57/61] Take comments into account --- yt/frontends/ramses/data_structures.py | 25 +++++++++++++------ yt/frontends/ramses/io_utils.pyx | 1 + yt/geometry/coordinates/coordinate_handler.py | 4 +-- yt/geometry/oct_container.pyx | 21 +++++++--------- yt/geometry/oct_visitors.pxd | 14 +++++------ yt/geometry/oct_visitors.pyx | 24 ++++++++++++++++-- 6 files changed, 58 insertions(+), 31 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index a0f1babab14..47ff2992790 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -1,7 +1,6 @@ import os import numpy as np import weakref -import warnings from collections import defaultdict from glob import glob @@ -193,10 +192,15 @@ def __init__(self, base_region, domain, ds, over_refine_factor=1, num_ghost_zone if num_ghost_zones > 0: if not all(ds.periodicity): - warnings.warn('Ghost zones will wrongly assume the domain to be periodic.') + mylog.warn('Ghost zones will wrongly assume the domain to be periodic.') # Create a base domain *with no self._base_domain.fwidth base_domain = RAMSESDomainSubset(ds.all_data(), domain, ds, over_refine_factor) self._base_domain = base_domain + elif num_ghost_zones < 0: + raise RuntimeError( + 'Cannot initialize a domain subset with a negative number of ghost zones,' + ' was called with num_ghost_zones=%s' % num_ghost_zones + ) def _fill_no_ghostzones(self, fd, fields, selector, file_handler): ndim = self.ds.dimensionality @@ -214,9 +218,10 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, 'float64') + cpu_list = [self.domain_id - 1] fill_hydro(fd, file_handler.offset, file_handler.level_count, - [self.domain_id-1], + cpu_list, levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler) @@ -240,9 +245,10 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, 'float64') + cpu_list = list(range(ncpu)) fill_hydro(fd, file_handler.offset, file_handler.level_count, - list(range(ncpu)), + cpu_list, levels, cell_inds, file_inds, ndim, all_fields, fields, tr, oct_handler, @@ -255,6 +261,9 @@ def fwidth(self): if self._num_ghost_zones > 0: fwidth = fwidth.reshape(-1, 8, 3) n_oct = fwidth.shape[0] + # new_fwidth contains the fwidth of the oct+ghost zones + # this is a constant array in each oct, so we simply copy + # the oct value using numpy fancy-indexing new_fwidth = np.zeros((n_oct, self.nz**3, 3), dtype=fwidth.dtype) new_fwidth[:, :, :] = fwidth[:, 0:1, :] fwidth = new_fwidth.reshape(-1, 3) @@ -269,12 +278,12 @@ def fcoords(self): oh = self.oct_handler indices = oh.fill_index(self.selector).reshape(-1, 8) - oinds, cinds = oh.fill_octcellindex_neighbours(self.selector) + oct_inds, cell_inds = oh.fill_octcellindex_neighbours(self.selector) - oinds = oinds.reshape(-1, 64) - cinds = cinds.reshape(-1, 64) + oct_inds = oct_inds.reshape(-1, 64) + cell_inds = cell_inds.reshape(-1, 64) - inds = indices[oinds, cinds] + inds = indices[oct_inds, cell_inds] fcoords = self.ds.arr( oh.fcoords(self.selector)[inds].reshape(-1, 3), diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 39d4d86d43d..21abdb1609e 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -172,6 +172,7 @@ def fill_hydro(FortranFile f, f.seek(offset) tmp = {} # Initalize temporary data container for io + # note: we use Fortran ordering to reflect the in-file ordering for field in all_fields: tmp[field] = np.empty((nc, twotondim), dtype="float64", order='F') diff --git a/yt/geometry/coordinates/coordinate_handler.py b/yt/geometry/coordinates/coordinate_handler.py index bb24accfad2..d5521e804ac 100644 --- a/yt/geometry/coordinates/coordinate_handler.py +++ b/yt/geometry/coordinates/coordinate_handler.py @@ -7,7 +7,7 @@ fix_unitary, \ iterable from yt.units.yt_array import \ - YTArray, YTQuantity, unyt_array + YTArray, YTQuantity from yt.utilities.exceptions import \ YTCoordinateNotImplemented, \ YTInvalidWidthError @@ -213,7 +213,7 @@ def sanitize_width(self, axis, width, depth): width = (w[0], w[1]) elif iterable(width): width = validate_iterable_width(width, self.ds) - elif isinstance(width, (YTQuantity, unyt_array)): + elif isinstance(width, YTQuantity): width = (width, width) elif isinstance(width, Number): width = (self.ds.quan(width, 'code_length'), diff --git a/yt/geometry/oct_container.pyx b/yt/geometry/oct_container.pyx index a6498556aff..db2e13490bf 100644 --- a/yt/geometry/oct_container.pyx +++ b/yt/geometry/oct_container.pyx @@ -664,13 +664,10 @@ cdef class OctreeContainer: cdef np.ndarray[np.int64_t, ndim=1] file_inds if num_cells < 0: num_cells = selector.count_oct_cells(self, domain_id) - levels = np.zeros(num_cells, dtype="uint8") - file_inds = np.zeros(num_cells, dtype="int64") - cell_inds = np.zeros(num_cells, dtype="uint8") - for i in range(num_cells): - levels[i] = 255 - file_inds[i] = -1 - cell_inds[i] = 8 + # Initialize variables with dummy values + levels = np.full(num_cells, 255, dtype="uint8") + file_inds = np.full(num_cells, -1, dtype="int64") + cell_inds = np.full(num_cells, 8, dtype="uint8") cdef oct_visitors.FillFileIndicesO visitor_o cdef oct_visitors.FillFileIndicesR visitor_r if self.fill_style == "r": @@ -735,8 +732,8 @@ cdef class OctreeContainer: for key in dest_fields: dest = dest_fields[key] source = source_fields[key] - for i in range(levels.shape[0]): - if levels[i] != level: continue + for i, lvl in enumerate(levels): + if lvl != level: continue if file_inds[i] < 0: dest[i + offset] = np.nan else: @@ -774,7 +771,7 @@ cdef class OctreeContainer: ------- oct_inds : int64 ndarray (nocts*8, ) The on-domain index of the octs containing each cell - cell_inds : uint8 array (nocts*8, ) + cell_inds : uint8 ndarray (nocts*8, ) The index of the cell in its parent oct Note @@ -786,8 +783,8 @@ cdef class OctreeContainer: cdef NeighbourCellIndexVisitor visitor - cdef np.uint8_t[:] cell_inds - cdef np.int64_t[:] oct_inds + cdef np.uint8_t[::1] cell_inds + cdef np.int64_t[::1] oct_inds cell_inds = np.full(num_octs*4**3, 8, dtype=np.uint8) oct_inds = np.full(num_octs*4**3, -1, dtype=np.int64) diff --git a/yt/geometry/oct_visitors.pxd b/yt/geometry/oct_visitors.pxd index 7a77a653299..63920480cf6 100644 --- a/yt/geometry/oct_visitors.pxd +++ b/yt/geometry/oct_visitors.pxd @@ -151,7 +151,7 @@ cdef class StoreIndex(OctVisitor): cdef class BaseNeighbourVisitor(OctVisitor): cdef int idim # 0,1,2 for x,y,z cdef int direction # +1 for +x, -1 for -x - cdef np.int8_t[:] neigh_ind + cdef np.uint8_t neigh_ind[3] cdef bint other_oct cdef Oct *neighbour cdef OctreeContainer octree @@ -164,12 +164,12 @@ cdef class BaseNeighbourVisitor(OctVisitor): return (((self.neigh_ind[2]*d)+self.neigh_ind[1])*d+self.neigh_ind[0]) cdef class NeighbourCellIndexVisitor(BaseNeighbourVisitor): - cdef np.uint8_t[:] cell_inds - cdef np.int64_t[:] domain_inds + cdef np.uint8_t[::1] cell_inds + cdef np.int64_t[::1] domain_inds cdef class NeighbourCellVisitor(BaseNeighbourVisitor): - cdef np.uint8_t[:] levels - cdef np.int64_t[:] file_inds - cdef np.uint8_t[:] cell_inds - cdef np.int32_t[:] domains + cdef np.uint8_t[::1] levels + cdef np.int64_t[::1] file_inds + cdef np.uint8_t[::1] cell_inds + cdef np.int32_t[::1] domains diff --git a/yt/geometry/oct_visitors.pyx b/yt/geometry/oct_visitors.pyx index 220da61f014..cd951d26271 100644 --- a/yt/geometry/oct_visitors.pyx +++ b/yt/geometry/oct_visitors.pyx @@ -353,7 +353,9 @@ cdef class StoreIndex(OctVisitor): cdef class BaseNeighbourVisitor(OctVisitor): def __init__(self, OctreeContainer octree, int domain_id = -1): self.octree = octree - self.neigh_ind = np.zeros(3, np.int8) + self.neigh_ind[0] = 0 + self.neigh_ind[1] = 0 + self.neigh_ind[2] = 0 super(BaseNeighbourVisitor, self).__init__(octree, domain_id) @cython.boundscheck(False) @@ -418,6 +420,25 @@ cdef class BaseNeighbourVisitor(OctVisitor): # Store neighbouring cell index in current cell cdef class NeighbourCellIndexVisitor(BaseNeighbourVisitor): + # This piece of code is very much optimizable. Here are possible routes to achieve + # much better performance: + + # - Work oct by oct, which would reduce the number of neighbor lookup + # from 4³=64 to 3³=27, + # - Use faster neighbor lookup method(s). For now, all searches are started from + # the root mesh down to leaf nodes, but we could instead go up the tree from the + # central oct then down to find all neighbors (see e.g. + # https://geidav.wordpress.com/2017/12/02/advanced-octrees-4-finding-neighbor-nodes/). + # - Pre-compute the face-neighbors of all octs. + + # Note that for the last point, algorithms exist that generate the neighbors of all + # octs in O(1) time (https://link.springer.com/article/10.1007/s13319-015-0060-9) + # during the octree construction. + # Another possible solution would be to keep a unordered hash map of all the octs + # indexed by their (3-integers) position. With such structure, finding a neighbor + # takes O(1) time. This could even come as a replacement of the current + # pointer-based octree structure. + @cython.boundscheck(False) @cython.wraparound(False) @cython.initializedcheck(False) @@ -506,4 +527,3 @@ cdef class NeighbourCellVisitor(BaseNeighbourVisitor): self.domains[self.index] = neigh_domain self.index += 1 - From 13e5dfad9fa9f73e5ef87dd0a901c45fcc50222b Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 17 Jul 2020 10:06:26 +0100 Subject: [PATCH 58/61] Inform Cython about strides --- yt/frontends/ramses/io_utils.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 21abdb1609e..1b29199e176 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -106,8 +106,8 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n offset = np.full((ncpu_and_bound, n_levels), -1, dtype=np.int64) level_count = np.zeros((ncpu_and_bound, n_levels), dtype=np.int64) - cdef np.int64_t[:,:] level_count_view = level_count - cdef np.int64_t[:,:] offset_view = offset + cdef np.int64_t[:, ::1] level_count_view = level_count + cdef np.int64_t[:, ::1] offset_view = offset for ilevel in range(nlevelmax): for icpu in range(ncpu_and_bound): From b1995407092950c600d40d6974f692acb8d8187e Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Sat, 18 Jul 2020 13:13:16 +0100 Subject: [PATCH 59/61] tr -> data --- yt/frontends/ramses/data_structures.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 5b827340cd2..8e75ea137aa 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -224,7 +224,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): oct_handler = self.oct_handler all_fields = [f for ft, f in file_handler.field_list] fields = [f for ft, f in fields] - tr = {} + data = {} cell_count = selector.count_oct_cells(self.oct_handler, self.domain_id) levels, cell_inds, file_inds = self.oct_handler.file_index_octs( @@ -233,7 +233,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: - tr[field] = np.zeros(cell_count, "float64") + data[field] = np.zeros(cell_count, "float64") cpu_list = [self.domain_id - 1] fill_hydro( @@ -247,10 +247,10 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): ndim, all_fields, fields, - tr, + data, oct_handler, ) - return tr + return data def _fill_with_ghostzones( self, fd, fields, selector, file_handler, num_ghost_zones From 8ecce954c0414911a3986b1d9b7268b8dbf27d86 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Sat, 18 Jul 2020 13:18:48 +0100 Subject: [PATCH 60/61] Make isort happy --- yt/frontends/ramses/data_structures.py | 1 + 1 file changed, 1 insertion(+) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 8e75ea137aa..a862424f114 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -2,6 +2,7 @@ import weakref from collections import defaultdict from glob import glob + import numpy as np from yt.arraytypes import blankRecordArray From c033f18b37ac634693ba1d4b7d9834c61176a50a Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Sat, 18 Jul 2020 20:48:47 +0100 Subject: [PATCH 61/61] Fix flaking --- yt/data_objects/data_containers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/yt/data_objects/data_containers.py b/yt/data_objects/data_containers.py index 1e3b543ff16..e9a0f5cf6ac 100644 --- a/yt/data_objects/data_containers.py +++ b/yt/data_objects/data_containers.py @@ -368,7 +368,6 @@ def _generate_spatial_fluid(self, field, ngz): np.empty(wogz.ires.size, dtype="float64"), units ) outputs.append(rv) - data = gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz] ind += wogz.select( self.selector, source=gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz],