Coordinate transformations and periodic boundaries
Hi all, I’m working on some updates to my yt-based pyXSIM package and I’m stuck and I thought I would email this list for some guidance. I’m thinking of a situation where we have a data object (sphere, box, whatever) and it straddles a periodic boundary. I want to convert the coordinates such that the coordinates are translated with respect to some arbitrary origin (say, the center of a sphere but in theory it could be anywhere), but are also continuous, i.e, they do not wrap at the boundary. I’ve looked at the functions get_radius and get_periodic_rvec in yt/fields/field_functions.py, and based on that I have come up with the code below, but it doesn’t quite work for any arbitrary value of the “ctr” argument (the origin). I think it’s because the functions I mentioned only need to calculate absolute value differences in order to compute a radius. The x, y, and z arguments to the function below are the input coordinate arrays. Could anyone who is more familiar with this sort of thing point out what I should be doing? Best, John Z def get_periodic_coords(ds, ctr, x, y, z): coords = ds.arr(np.zeros((3, x.size)), "kpc") coords[0, :] = (x - ctr[0]).to("kpc") coords[1, :] = (y - ctr[1]).to("kpc") coords[2, :] = (z - ctr[2]).to("kpc") if sum(ds.periodicity) == 0: return coords dw = ds.domain_width.to("kpc") for i in range(3): if ds.periodicity[i]: c = coords[i, :] cdw = c - dw[i] mins = np.argmin([np.abs(c), np.abs(cdw)], axis=0) ii = np.where(mins == 1)[0] coords[i, ii] = coords[i, ii] - dw[i] return coords
Hi John, On Thu, Sep 8, 2016 at 10:22 AM, John ZuHone <jzuhone@gmail.com> wrote:
Hi all,
I’m working on some updates to my yt-based pyXSIM package and I’m stuck and I thought I would email this list for some guidance.
I’m thinking of a situation where we have a data object (sphere, box, whatever) and it straddles a periodic boundary. I want to convert the coordinates such that the coordinates are translated with respect to some arbitrary origin (say, the center of a sphere but in theory it could be anywhere), but are also continuous, i.e, they do not wrap at the boundary.
I’ve looked at the functions get_radius and get_periodic_rvec in yt/fields/field_functions.py, and based on that I have come up with the code below, but it doesn’t quite work for any arbitrary value of the “ctr” argument (the origin). I think it’s because the functions I mentioned only need to calculate absolute value differences in order to compute a radius.
The x, y, and z arguments to the function below are the input coordinate arrays. Could anyone who is more familiar with this sort of thing point out what I should be doing?
I think that Doug implemented something in the region selector that does exactly this; it computes whether or not something needs a coordinate shift, and then if so, it applies that. If you look through selection_routines.pyx for anything with "shift" in the name, it should be there. If that doesn't do what you're looking for, we can try to figure something else out.
Best,
John Z
def get_periodic_coords(ds, ctr, x, y, z): coords = ds.arr(np.zeros((3, x.size)), "kpc") coords[0, :] = (x - ctr[0]).to("kpc") coords[1, :] = (y - ctr[1]).to("kpc") coords[2, :] = (z - ctr[2]).to("kpc") if sum(ds.periodicity) == 0: return coords dw = ds.domain_width.to("kpc") for i in range(3): if ds.periodicity[i]: c = coords[i, :] cdw = c - dw[i] mins = np.argmin([np.abs(c), np.abs(cdw)], axis=0) ii = np.where(mins == 1)[0] coords[i, ii] = coords[i, ii] - dw[i] return coords
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
Hi Matt, Thanks for the pointer. I think the real difficulty here is that in theory the data object used to generate the photons could be literally anything, a sphere, box, disk, or some goofy shape defined by a region cut. The code in selection_routines.pyx that handles such periodic shifting appears (to me at least) to apply only for rectangular regions for which such an operation is relatively straightforward—but I’d need something that doesn’t necessarily know about “left edge” and “right edge”. You could special-case that for things like rectangular regions and spheres, perhaps, but for anything more general determining a “left edge” becomes a problem. For example, you simply couldn’t just use the minimum value of the coordinate because of course your region could straddle the domain boundary and your minimum coordinate could actually be on the “rightward” side of the object. Best, John
On Sep 9, 2016, at 1:12 PM, Matthew Turk <matthewturk@gmail.com> wrote:
Hi John,
On Thu, Sep 8, 2016 at 10:22 AM, John ZuHone <jzuhone@gmail.com <mailto:jzuhone@gmail.com>> wrote: Hi all,
I’m working on some updates to my yt-based pyXSIM package and I’m stuck and I thought I would email this list for some guidance.
I’m thinking of a situation where we have a data object (sphere, box, whatever) and it straddles a periodic boundary. I want to convert the coordinates such that the coordinates are translated with respect to some arbitrary origin (say, the center of a sphere but in theory it could be anywhere), but are also continuous, i.e, they do not wrap at the boundary.
I’ve looked at the functions get_radius and get_periodic_rvec in yt/fields/field_functions.py, and based on that I have come up with the code below, but it doesn’t quite work for any arbitrary value of the “ctr” argument (the origin). I think it’s because the functions I mentioned only need to calculate absolute value differences in order to compute a radius.
The x, y, and z arguments to the function below are the input coordinate arrays. Could anyone who is more familiar with this sort of thing point out what I should be doing?
I think that Doug implemented something in the region selector that does exactly this; it computes whether or not something needs a coordinate shift, and then if so, it applies that. If you look through selection_routines.pyx for anything with "shift" in the name, it should be there. If that doesn't do what you're looking for, we can try to figure something else out.
Best,
John Z
def get_periodic_coords(ds, ctr, x, y, z): coords = ds.arr(np.zeros((3, x.size)), "kpc") coords[0, :] = (x - ctr[0]).to("kpc") coords[1, :] = (y - ctr[1]).to("kpc") coords[2, :] = (z - ctr[2]).to("kpc") if sum(ds.periodicity) == 0: return coords dw = ds.domain_width.to <http://ds.domain_width.to/>("kpc") for i in range(3): if ds.periodicity[i]: c = coords[i, :] cdw = c - dw[i] mins = np.argmin([np.abs(c), np.abs(cdw)], axis=0) ii = np.where(mins == 1)[0] coords[i, ii] = coords[i, ii] - dw[i] return coords
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org <mailto:yt-dev@lists.spacepope.org> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org>
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org <mailto:yt-dev@lists.spacepope.org> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org>
On Fri, Sep 9, 2016 at 12:57 PM, John Zuhone <jzuhone@gmail.com> wrote:
Hi Matt,
Thanks for the pointer.
I think the real difficulty here is that in theory the data object used to generate the photons could be literally anything, a sphere, box, disk, or some goofy shape defined by a region cut.
The code in selection_routines.pyx that handles such periodic shifting appears (to me at least) to apply only for rectangular regions for which such an operation is relatively straightforward—but I’d need something that doesn’t necessarily know about “left edge” and “right edge”. You could special-case that for things like rectangular regions and spheres, perhaps, but for anything more general determining a “left edge” becomes a problem. For example, you simply couldn’t just use the minimum value of the coordinate because of course your region could straddle the domain boundary and your minimum coordinate could actually be on the “rightward” side of the object.
You could do it with a minimum enclosing bounding box for non-rectilinear data objects. I also have thoughts about ways to improve the cut_region data object (the only one that this will be non-trivial for) by making use of Meagan Lang's fast, periodic C KDTree that she will be adding to yt soon. Perhaps you could punt on that for now, and then come back to this when we improve the cut_region code? -Nathan
Best,
John
On Sep 9, 2016, at 1:12 PM, Matthew Turk <matthewturk@gmail.com> wrote:
Hi John,
On Thu, Sep 8, 2016 at 10:22 AM, John ZuHone <jzuhone@gmail.com> wrote:
Hi all,
I’m working on some updates to my yt-based pyXSIM package and I’m stuck and I thought I would email this list for some guidance.
I’m thinking of a situation where we have a data object (sphere, box, whatever) and it straddles a periodic boundary. I want to convert the coordinates such that the coordinates are translated with respect to some arbitrary origin (say, the center of a sphere but in theory it could be anywhere), but are also continuous, i.e, they do not wrap at the boundary.
I’ve looked at the functions get_radius and get_periodic_rvec in yt/fields/field_functions.py, and based on that I have come up with the code below, but it doesn’t quite work for any arbitrary value of the “ctr” argument (the origin). I think it’s because the functions I mentioned only need to calculate absolute value differences in order to compute a radius.
The x, y, and z arguments to the function below are the input coordinate arrays. Could anyone who is more familiar with this sort of thing point out what I should be doing?
I think that Doug implemented something in the region selector that does exactly this; it computes whether or not something needs a coordinate shift, and then if so, it applies that. If you look through selection_routines.pyx for anything with "shift" in the name, it should be there. If that doesn't do what you're looking for, we can try to figure something else out.
Best,
John Z
def get_periodic_coords(ds, ctr, x, y, z): coords = ds.arr(np.zeros((3, x.size)), "kpc") coords[0, :] = (x - ctr[0]).to("kpc") coords[1, :] = (y - ctr[1]).to("kpc") coords[2, :] = (z - ctr[2]).to("kpc") if sum(ds.periodicity) == 0: return coords dw = ds.domain_width.to("kpc") for i in range(3): if ds.periodicity[i]: c = coords[i, :] cdw = c - dw[i] mins = np.argmin([np.abs(c), np.abs(cdw)], axis=0) ii = np.where(mins == 1)[0] coords[i, ii] = coords[i, ii] - dw[i] return coords
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
_______________________________________________ yt-dev mailing list yt-dev@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
participants (4)
-
John Zuhone -
John ZuHone
-
Matthew Turk -
Nathan Goldbaum