BUG: Order of Cylindrical Coordinates for Athena++ Incorrect

[forrestglines]

Bug report

Bug summary

The order of coordinates for cylindrical datasets from Athena++ seem to be incorrect. I ran across this issue writing #4323

I have a hot-fix for just my purposes to reorder the coordinates. However, changing the geometry to Geometry.POLAR might be a better fix once #4790 is done. I can provide my hot-fix on request.

Code for reproduction

Here's the athdf file for a 2D Keplerian disk with uniform r-coordinates modified from the input deck athinput.disk_cyl included with Athena++:

I try to plot the disk with

import yt
ds = yt.load("disk.out1.00000.athdf")
slc = yt.SlicePlot(ds, "z", ("gas", "density"))
slc.show()

Actual outcome

Expected outcome

Whereas if I force the correct coordinate ordering by changing the Athena++ frontend to respect $(r,\theta,z)$ I get the expected disk

Version Information

Using current main commit: d5b3335

[neutrinoceros]

Hi @forrestglines, thank you for reporting. I think this is reasonably decoupled with #4790 so feel free to open a PR with your fix at any point !

[matthewturk]

It might be worth checking that we're using the ds.coordinates.axis_name dict wherever we need to; that's a potential source of problems in the selection and pixelization routines. (i.e., assuming that 0 == r, 1 == theta, 2 == phi, rather than using the axis_names to map)

[forrestglines]

Hi @forrestglines, thank you for reporting. I think this is reasonably decoupled with #4790 so feel free to open a PR with your fix at any point !

Thanks for the input! This issue with coordinates is the same issue I encountered developing the Parthenon backend. The two data formats (and the AthenaPK code) are extremely similar so I'd like the frontends to behave similarly. I see two options preferable options to fix both:

1. We go with Geometry.CYLINDRICAL and change axis_order for both Athena++, as I do in #4801, and for Parthenon as I originally did in #4323

2. We make both Athena++ and Parthenon use Geometry.POLAR and we fix #4790.

I'm a bit inclined to go with the first option to keep Athena++ as Geometry.CYLINDRICAL and have Parthenon match that solution. I'll still help out to debug #4790.

[forrestglines]

It might be worth checking that we're using the ds.coordinates.axis_name dict wherever we need to; that's a potential source of problems in the selection and pixelization routines. (i.e., assuming that 0 == r, 1 == theta, 2 == phi, rather than using the axis_names to map)

Thanks @matthewturk! This is a suggestion for implementation for the selection and pixelization routines and what's likely causing issues in #4790, right? In this issue I'm having the trouble that yt by default uses r == 0, z == 1, theta == 2 for cylindrical dataset when Athena++ cylindrical datasets are r == 0, theta == 1, z == 2.

Does the fix in #4801 look congruent with yt standards/expectations? It only changes the assumptions about Athena++ cylindrical dataset but hopefully no other assumptions in the code.

[neutrinoceros]

I see two options to fix both:

I don't think those options should be mutually exclusive: if using Geometry.CYLINDRICAL works for you and avoids the current bugs with Geometry.POLAR, that's a reasonable thing to do as a workaround, and we can always switch both frontends to use the more correct solution (POLAR) once known issues are addressed. What do you think ?

[forrestglines]

I don't think those options should be mutually exclusive: if using Geometry.CYLINDRICAL works for you and avoids the current bugs with Geometry.POLAR, that's a reasonable thing to do as a workaround, and we can always switch both frontends to use the more correct solution (POLAR) once known issues are addressed. What do you think ?

I think this is the right solution, this is good enough for both frontends. When the issues with POLAR are addressed we can switch both frontends.