Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Question: expected behavior of center argument in slices of spherical datasets? #4841

Open
chrishavlin opened this issue Feb 29, 2024 · 3 comments

Comments

@chrishavlin
Copy link
Contributor

chrishavlin commented Feb 29, 2024

I was looking into a question on yt-users and got a bit confused myself so wanted to see if anyone has any thoughts. Here is a link to the original question, but the implicit question is a more general question about slices in spherical datasets: should the center argument have an effect?

Right now it seems that the center argument is only used to determine the value of the normal and is otherwise ignored:

from yt.testing import fake_amr_ds
import yt
ds = fake_amr_ds(geometry='spherical')

slc = yt.SlicePlot(ds, 'phi', 'Density', window_size=(3,3))
slc.save('phi_normal_0.png')

phival = ds.domain_center[-1].d
slc = yt.SlicePlot(ds, 'phi', 'Density',              
             center=(0.8, 3.14/2, phival),  
             window_size=(3,3)
             )
slc.save('phi_normal_1.png')

phival = ds.domain_center[-1].d
slc = yt.SlicePlot(ds, 'phi', 'Density',              
             center=(0.8, 3.14/2, phival),             
             width = (0.2, 0.2), 
             window_size=(3,3)
             )
slc.save('phi_normal_2.png')

phi_normal_0
phi_normal_1
phi_normal_2

I naively expected the image center to be calculated by transforming the r and theta values provided to the image x-y coordinates (in this case that'd be (0.8, 0.0) in the 2nd and 3rd image above), but looking into the code it seems that the center is fixed to always be the domain center (I think the relevant section of code is in _sanitize_center). Is this a case of the behavior simply not being implemented yet or is there a different reason for always using the full domain center?

@neutrinoceros
Copy link
Member

neutrinoceros commented Mar 1, 2024

Related to #4362

This has been a problem for as long as spherical coordinates have been supported. To me, the issue is primarily that the API, which was originally designed for cartesian coordinates only, is ill suited for non-cartesian geometries; in particular, specifying a {center, width} combo is ambiguous in spherical geometries: does it result in a spherical selection (possibly off-center), or does it really just select a slab in the data cube1 ? what coordinates is center in ? native, or cartesian2 ? is the width akin to a radius, a diameter, an azimuth, a latitude ? Combined with the fact that angular coordinates are, under the hood, represented as length dimensions, none of these questions have obvious answers to me.
So far, my approach has been very conservative, but I think we can't really beat this problem without a broader redesign of non-cartesian coordinates.

I think we'd need a couple fundamental changes:

  • stop representing all coordinates as lengths. This means in particular that coord triplets returned by, e.g., ds.find_max, need to be tuples of scalar unyt_quantity3 instead of a single unyt_quantity (which cannot have more than 1 dimensionality)
  • NormalPlot API:
    • elect a coordinate system (Cartesian or data-native) that the center. I don't think it matters much what we choose to support as long as it's explicit. I do have a slight preference for choosing data-native because it's a bit easier to convert to cartesian than the other way around, which might avoid some foot guns to users.
    • forbid passing a scalar width for non-cartesian geometries3: only triplets should be accepted (or pairs in 2D)
    • clarify, in docs, that center + width always selects a slab in data-native coordinates

Footnotes

  1. currently, the latter is true

  2. the current default, (0,0,0), is ambiguous but degenerate, which is why it "works" by default, but user-provided value are ignored. Also note that this ambiguity+degeneracy is intrinsic in cartesian coordinates, which is why it was never a problem there.

  3. That would be a major breaking API change. So, time for yt 5.0 ?? 2

@chrishavlin
Copy link
Contributor Author

Thanks for the detailed thoughts, @neutrinoceros ! And searching back through issues and PRs, I see now that there's been a lot of related discussion and attempts circling around this problem (should have searched a bit more before posting my question...). In any case, it does seem like a bigger refactor is needed... I do have more thoughts on this, will try to come back later today or tomorrow to type them up!

@matthewturk
Copy link
Member

I wanted to mention two things -- the first is that I think we could absolutely explore a new method of designing these slices and whatnot.

The second is that we should explore how we could revive the "IndexArray" type that @jzuhone worked on a few years ago in unyt. This would be a big change, and we may decide it's too much. But the idea there is twofold -- one is to make all of the arrays used for indexing read-only -- which could be a huge performance increase in some areas of the code that constantly do conversions or checks that arrays are in code-length. But the other part is to make some types of arrays non-broadcastable along all dimensions, and allow for non-homogeneous units. (Making this a general purpose solution in unyt would likely be a lot harder than making it work in yt, to be clear. So it's also possible that we could implement it in yt for that one specific use case.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants