BUG: Fix intersections between SPH particles and regions when the boundary is periodic

[jzuhone]

PR Summary

This PR fixes the issue raised by #3916 and seen in other contexts where sometimes particles which properly belong to a rectangular region straddling a periodic boundary would get missed.

I think something was not quite right in the overall logic, but I have to be honest I am not certain. What I do know for sure is that instead of checking the particle position itself for particles that lay outside the boundary of a region, we should be checking the "left_edge" and "right_edge" of the particle in question (defined by the +/- the radius away from the position), and then computing the radius of the particle to the edge. That was inspired by considering the region/cell intersection implementation:

yt/yt/geometry/_selection_routines/region_selector.pxi

Lines 81 to 92 in 6c0780f

	@cython.boundscheck(False)
	@cython.wraparound(False)
	@cython.cdivision(True)
	cdef int select_bbox(self, np.float64_t left_edge[3],
	np.float64_t right_edge[3]) nogil:
	cdef int i
	for i in range(3):
	if (right_edge[i] < self.left_edge[i] and \
	left_edge[i] >= self.right_edge_shift[i]) or \
	left_edge[i] >= self.right_edge[i]:
	return 0
	return 1

This PR implements that logic, and makes the code logic clearer. It also uses the already defined periodic_difference in the superclass to compute the distance between the particle position and the region edge. However, I noticed during this exercise that we have this same logic implemented in two different places (both of which are used):

Here:

yt/yt/geometry/_selection_routines/selector_object.pxi

Lines 266 to 278 in 6c0780f

	@cython.boundscheck(False)
	@cython.wraparound(False)
	@cython.cdivision(True)
	cdef np.float64_t periodic_difference(self, np.float64_t x1, np.float64_t x2, int d) nogil:
	# domain_width is already in code units, and we assume what is fed in
	# is too.
	cdef np.float64_t rel = x1 - x2
	if self.periodicity[d]:
	if rel > self.domain_width[d] * 0.5:
	rel -= self.domain_width[d]
	elif rel < -self.domain_width[d] * 0.5:
	rel += self.domain_width[d]
	return rel

and here:

yt/yt/geometry/selection_routines.pxd

Lines 81 to 89 in 6c0780f

	cdef inline np.float64_t _periodic_dist(np.float64_t x1, np.float64_t x2,
	np.float64_t dw, bint periodic) nogil:
	cdef np.float64_t rel = x1 - x2
	if not periodic: return rel
	if rel > dw * 0.5:
	rel -= dw
	elif rel < -dw * 0.5:
	rel += dw
	return rel

Should I presume the latter would be faster?

Also, it's possible this change may break some answer tests.

PR Checklist

• New features are documented, with docstrings and narrative docs
• Adds a test for any bugs fixed. Adds tests for new features.

[matthewturk]

THANK YOU, John.

My guess is that the latter probably is faster, but there may be cases where we can't call it, or where type conversion may end up being more expensive than we'd expect.

[jzuhone]

By the way, this is my modified script from @neutrinoceros:

import numpy as np
import yt


def fake_ds(hsml_factor):
    npart = 3**3
    x = np.empty(npart)
    y = np.empty(npart)
    z = np.empty(npart)

    tot = 0
    for i in range(0, 3):
        for j in range(0, 3):
            for k in range(0, 3):
                x[tot] = i + 0.5
                y[tot] = j + 0.5
                z[tot] = k + 0.5
                tot += 1

    data = {
        "particle_position_x": (x, "cm"),
        "particle_position_y": (y, "cm"),
        "particle_position_z": (z, "cm"),
        "particle_mass": (np.ones(npart), "g"),
        "particle_velocity_x": (np.zeros(npart), "cm/s"),
        "particle_velocity_y": (np.zeros(npart), "cm/s"),
        "particle_velocity_z": (np.zeros(npart), "cm/s"),
        "smoothing_length": (0.05 * np.ones(npart) * hsml_factor, "cm"),
        "density": (np.ones(npart), "g/cm**3"),
        "temperature": (np.ones(npart), "K"),
    }

    bbox = np.array([[0, 3], [0, 3], [0, 3]])

    return yt.load_particles(data=data, length_unit=1.0, bbox=bbox)


ds = fake_ds(10)

for w in (3, 2, 1.5, 1.0, 0.5):
    center = np.array([2.8]*3)
    proj = yt.ProjectionPlot(
        ds, "z", ("gas", "density"), center=center, width=w,
    )
    proj.save(f"/tmp/test_{w=}.png")

    le = center-0.5*w
    re = center+0.5*w
    box = ds.box(le, re)
    print(box["io","smoothing_length"].size)
    proj = yt.ProjectionPlot(
        ds, "z", ("gas", "density"), center=center, data_source=box)
    proj.save(f"/tmp/test_box_{w=}.png")

which adds more widths and it also adds a test where we keep the image fixed at the domain width but make the region itself smaller. As expected, the number of particles plotted decreases.

And the results:

[neutrinoceros]

Very nice, thanks for solving this. I will not be able to review this properly until tomorrow, but if Matt wants to push the button, I don't want to make you guys wait for me when it looks "obviously right".

[jzuhone]

Ping @benopp99 let's test on that dataset of yours

[matthewturk]

I think this looks right.

[benopp99]

I'm going to have to refigure out how to access CAMELS data now that I don't have my account.

…

-Ben

On Thu, Aug 4, 2022 at 11:37 AM Matthew Turk ***@***.***> wrote: ***@***.**** approved this pull request. I think this looks right. — Reply to this email directly, view it on GitHub <#4050 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMEF72NKQYF5HQHIB4MRWOTVXP5XDANCNFSM55TFEXWQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[benopp99]

Okay, this link now should work: https://wormhole.app/MjqOK#npxkeFwTqrXr9MTxT3yTFg On Thu, Aug 4, 2022 at 11:42 AM Benjamin Oppenheimer ***@***.***> wrote:

…

I'm going to have to refigure out how to access CAMELS data now that I don't have my account. -Ben On Thu, Aug 4, 2022 at 11:37 AM Matthew Turk ***@***.***> wrote: > ***@***.**** approved this pull request. > > I think this looks right. > > — > Reply to this email directly, view it on GitHub > <#4050 (review)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AMEF72NKQYF5HQHIB4MRWOTVXP5XDANCNFSM55TFEXWQ> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

[jzuhone]

@benopp99 I can confirm this code fixes your plots.

[benopp99]

Outstanding! Glad I could find a bug (an important one) stress testing YT on simulations. This will be quite useful. Let me know when the bug fix is in an updated version. Best, Ben

…

On Thu, Aug 4, 2022 at 5:17 PM John ZuHone ***@***.***> wrote: @benopp99 <https://github.com/benopp99> I can confirm this code fixes your plots. — Reply to this email directly, view it on GitHub <#4050 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMEF72LIK2EIHVQMELLD4ZLVXRFSJANCNFSM55TFEXWQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[jzuhone]

@benopp99 you can always check out this branch until then :) but it will be merged into main no later than tomorrow I assume.

[jzuhone]

Ok, I added a test for this. Provided it passes, this is ready to go.

[neutrinoceros]

LGTM, nice craft !