BUG: fix hard-coded 2x2x2 zone-count assumptions in octree neighbour-finding (issue 5402) - #5508
BUG: fix hard-coded 2x2x2 zone-count assumptions in octree neighbour-finding (issue 5402)#5508jurb33 wants to merge 5 commits into
Conversation
|
Hi! Welcome, and thanks for opening this pull request. We have some guidelines for new pull requests, and soon you'll hear back about the results of our tests and continuous integration checks. Thank you for your contribution! |
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
cphyc
left a comment
cphyc
left a comment
There was a problem hiding this comment.
Thanks @jurb33 for fixing this!
This looks reasonable to me, but I'm wondering whether it'd be possible to simplify the test? As far as I could understand, you're creating an octree from a collection of particles using an arbitrary nz.
Could we use the dedicated yt.testing.fake_octree_ds function instead? From what I could gather, I think the answer is not a straight yes, but if you could give it a try that'd be awesome!
Otherwise, the fixes look great.
NOTE: I edited your PR for cosmetic changes, so please git pull before adding anything new.
| morton = get_morton_indices(np.floor((centers - DLE) / dx).astype("uint64")) | ||
| morton.sort() | ||
|
|
||
| for nz in ((2, 2, 2), (2, 3, 4)): |
There was a problem hiding this comment.
I would suggest also testing that (1, 1, 1) doesn't work (it doesn't on main).
| for nz in ((2, 2, 2), (2, 3, 4)): | |
| for nz in ((1, 1, 1), (2, 2, 2), (2, 3, 4)): |
| octree = ParticleOctreeContainer((1, 1, 1), DLE, DRE, num_zones=nz) | ||
| octree.n_ref = 1 | ||
| octree.add(morton) | ||
| octree.finalize() | ||
|
|
||
| nzones = nz[0] * nz[1] * nz[2] | ||
| n_per_oct = (nz[0] + 2) * (nz[1] + 2) * (nz[2] + 2) # +1 ghost zone each side |
There was a problem hiding this comment.
This is fairly convoluted - do you think it'd be possible to use instead yt.load_octree(..., num_zones=(1, 2, 3))?
EDIT: I just tested, and the short answer is no... because load_octree expects num_zones to be a scalar.
| i1 = 2 + self.n_ghost_zones | ||
| # Loop over cells in and directly around oct | ||
| for i in range(i0, i1): | ||
| for i in range(i0, self.nz[0] + self.n_ghost_zones): |
There was a problem hiding this comment.
How confident are you that this shouldn't be the following?
| for i in range(i0, self.nz[0] + self.n_ghost_zones): | |
| for i in range(i0, self.nz[2] + self.n_ghost_zones): |
(and self.nz[0] for the innermost loop)
There was a problem hiding this comment.
This is correct - looking at these lines:
yt/geometry/oct_visitors.pyx:391
local_oct &= (0 <= ishift[i] < self.nz[i])
yt/geometry/oct_visitors.pyx:415
self.neigh_ind[i] = <np.uint32_t>(ishift[i]) % self.nz[i]
is using the same ishift[i] self.nz[i] pairing convention everywhere else in oct_visitors.pyx
Using self.nz[2] would break the test case which is what this was intended to fix. (ran this locally)
As far as the test implementation, yes, I can make these changes. May be a bit since school is starting up for me and I am pursuing different projects.
Thank you for reviewing my changes!
3 participants
Addresses issue #5402.
PR Summary
Neighbour finding in
oct_container.pyx/oct_visitors.pyx(used for ghost zone construction) assumed every oct is 2x2x2 cells, failing when non cubic zones encountered. The fix is changing how the index of neighbour cells are calculated in several places (issue list was not exhaustive)Changes:
BaseNeighbourVisitor.set_neighbour_info: checking "am I still in this oct" and neighbour cell index used%2 and <=1, now usesself.nz[I]for the actual width.NeighbourCellIndexVisitor.visit/NeighbourCellVisitor.visit: Used a shared bound2 + n_ghost_zonesto check each axis, and hardcoded8as the no cell condition. now usesself.nz[axis]andself.nz[0,1,2]**3(not real python, but you get the idea)oct_container.pyx: the output forfill_octcellindex_neighbours/file_index_octs_with_ghost_zoneswere sized as4**3regardless of actual cell count.Two of the cases above were not present in the original ticket
set_neighbour_info.Testing:
Added
test_octcellindex_neighbours_num_zones:yt/data_objects/tests/test_octree.pywhich builds the error case withnum_zones=(2,2,2) and (2,3,4)and asserts the expected cell count. Old code fails this, where this change fixes situations like these. All previous tests pass unchanged.While tracing function calls, I noticed in a higher level function that
yt/frontends/ramses/data_structures.py:fcoordscalls
oct_container.pyx:fill_octcellindex_neighbourslike
oh.fill_octcellindex_neighbours(self.selector, self._num_ghost_zones)which is actually
num_octsindef fill_octcellindex_neighbours(self, selector, num_octs=-1, domain_id=-1)passing ghost_zones incorrectly. This could cause overflow of a single oct and give incorrect results. Maybe a Issue should be opened up for this?
PR Checklist