BUG: fix incorrect default limits for slice plots along curvilinear directions (phi, theta)

[neutrinoceros]

PR Summary

addresses #3529 for the phi and theta cases
TODO:

• fix "phi" slices
• fix "theta" slices
• cleanup

Here's the result from #3529's script with this branch:

[neutrinoceros]

My current solution works for the exact dataset that I wrote the test for, however it is correct in the general case. What's missing is a way to calculate the min and max for elevation (z) and algebraic cylindrical radius (R), composing only with left and right edges of the domain.

[neutrinoceros]

Got a much more general solution for the phi-oriented case:

Now there's the issue of duplicated code, but I'll deal with that when I get every case to work (including theta-oriented plots)

[neutrinoceros]

Sigh... this solution isn't general enough either, because it will always set Rmax=rmax
Will come back to this later

[neutrinoceros]

I think I got the the right approach now

[neutrinoceros]

While I was at it, I also fixed boxing in the "theta"-oriented case

import yt

ds = yt.load_sample("bw_spherical_2d", unit_system="code")
p = yt.SlicePlot(ds, "theta", "density")
p.save("/tmp/raw.png")

main branch:

this branch:

Now I need to come back one more time to clean up this messy code before I can open this to review.
edit: the axis ticks are also inverted in the "fixed" case, this is also something I'll need to address before marking this as ready.

[neutrinoceros]

I'm having a rough time figuring out how to fix the "theta-normal" case. In my previous comment, I noted that axis ticks were x/y inverted, despite the width and centre being correctly computed and the overall result being correct (maybe not the aspect, but well, that would be another sign of a single problem).
I think it way have something to do with :

• the fact that SphercialCoordinateHandler._sanitize_center returns a tuple (display_center) with 3 elements instead of 2, the second one being completely unimportant except for the fact that something else downstream looks for a third element
• the definitions of SphericalCoordinateHandler._x_pairs (.y_pairs respectively) which don't look easy to retro-engineer. I think they are supposed to be "right-handed" and "left-handed" direction pairs but that assumption doesn't hold, and I really don't get how they are currently written.
  @matthewturk , you seem to be the most recent contributor here, any chance you would spot something terribly wrong here ? (or maybe provide some insight on how those pairs are defined ?)
  

  yt/yt/geometry/coordinates/spherical_coordinates.py

  Line 273 in c531b73

  _x_pairs = (("r", "theta"), ("theta", "r"), ("phi", "r"))

[neutrinoceros]

I found a pretty ugly hack that fixes the issue:

I'm really not found of the current solution as a maintainer, but it seems valuable (maybe not viable) if it works.
On the user side I don't have any complaints left so I'll open this to review and hope for the best, maybe we can sort this out.
edit : while this appears to be working fine with AMRVAC data, it completely breaks with another spherical dataset I'm using here : https://github.com/neutrinoceros/yt_idefix/tree/main/tests/data/FargoMHDSpherical
and I don't even get a proper failure: it's just halting my code forever or creating an infinite loop... anyway. More work needed again, I guess.

[neutrinoceros]

Alright, figured out the last (fatal) issue I had with my Idefix dataset came from some floating point truncation in a case where the domain right edge in the azimuthal direction was supposed to be 2 PI, so I added a small tolerance margin, large enough for single precision

Here's the plot for this dataset with yt's main branch

and this branch

[neutrinoceros]

So I've simplified this as much as I could but I wasn't able to fix the r-case. Still this PR fixes 2/3 of the problematic cases reported in #3529 so I'm going to open it for review.

[matthewturk]

Thank you for doing this. Some comments...

[matthewturk]

Your comment makes me think that we need to have a discussion about moving some of these functions around a bit. It seems that you've likely already come up to speed on what everything does and would be a good person to speculate about potential refactorings that touch the coordinate handlers. Could we set up a time to chat?

[neutrinoceros]

Actually I went with a very dumb approach of touching everything in a sensible way until something worked. Now I wouldn't say I was able to understand most of it, and I don't think I would be in a position to champion such a refactor, but yes, let's set up a time to discuss that over zoom :)

[neutrinoceros]

@matthewturk I've backported functools.cached_property from Python 3.8 to help mitigate your concerns about running the exact same computations multiple times.
I also added a test but I'm not 100% convinced it will check what I want.

[neutrinoceros]

@yt-fido test this please

[matthewturk]

Just to be clear, my concerns were not articulated as well as I'd like. I am not really concerned about the performance impact as much as just the conceptual nature of it being a singly-calculated item that could be done just once.

[neutrinoceros]

yup, agreed that it didn't make much sense. Now, as a bonus, we'll be able to start using cached_property before we even drop python 3.7, yay.

[Xarthisius]

@yt-fido test this please

[matthewturk]

This looks good to me. Any reason not to merge?

[neutrinoceros]

Thank you for asking. I needed to unlink the linked issue since it always resolves 2/3 cases, but this can go in as is AFAIC, I don't think I'll be able to manage the remaining (theta-phi) soon

[matthewturk]

@meeseeksdev backport to yt-4.0.x