Athena++ Problem Loading SMR in Spherical-Polar Coordinates

[smressle]

Bug report

Bug summary

yt returns an error when attempting to construct the mesh for athena++ data files when using static mesh refinement in spherical-polar coordinates.

Code for reproduction

ds = yt.load('Implode.out2.00010.athdf')
print(ds.field_list)

data set: http://use.yt/upload/3f780fa3

Actual outcome

/global/scratch/smressle/star_cluster/restart_grmhd/vis/python/athena_script.py in <module>()
----> 1 print(ds.field_list)

~/yt/yt/data_objects/static_output.py in field_list(self)
    543     @property
    544     def field_list(self):
--> 545         return self.index.field_list
    546 
    547     def create_field_info(self):

~/yt/yt/data_objects/static_output.py in index(self)
    501                 raise RuntimeError("You should not instantiate Dataset.")
    502             self._instantiated_index = self._index_class(
--> 503                 self, dataset_type=self.dataset_type)
    504             # Now we do things that we need an instantiated index for
    505             # ...first off, we create our field_info now.

~/yt/yt/frontends/athena_pp/data_structures.py in __init__(self, ds, dataset_type)
     66     def __init__(self, ds, dataset_type = 'athena_pp'):
     67         self._handle = ds._handle
---> 68         super(AthenaPPLogarithmicIndex, self).__init__(ds, dataset_type)
     69         self.index_filename = self.dataset.filename
     70         self.directory = os.path.dirname(self.dataset.filename)

~/yt/yt/geometry/unstructured_mesh_handler.py in __init__(self, ds, dataset_type)
     34         self.directory = os.path.dirname(self.index_filename)
     35         self.float_type = np.float64
---> 36         super(UnstructuredIndex, self).__init__(ds, dataset_type)
     37 
     38     def _setup_geometry(self):

~/yt/yt/geometry/geometry_handler.py in __init__(self, ds, dataset_type)
     48 
     49         mylog.debug("Setting up domain geometry.")
---> 50         self._setup_geometry()
     51 
     52         mylog.debug("Initializing data grid data IO")

~/yt/yt/geometry/unstructured_mesh_handler.py in _setup_geometry(self)
     38     def _setup_geometry(self):
     39         mylog.debug("Initializing Unstructured Mesh Geometry Handler.")
---> 40         self._initialize_mesh()
     41 
     42     def get_smallest_dx(self):

~/yt/yt/frontends/athena_pp/data_structures.py in _initialize_mesh(self)
    107         pbar = get_pbar("Constructing meshes", num_meshes)
    108         for i in range(num_meshes):
--> 109             ob = bc[i][0]
    110             x = x1f[ob,:]
    111             y = x2f[ob,:]

IndexError: too many indices for array
#

Expected outcome

Something like this

Constructing meshes: 100%|###################| 420/420 [00:04<00:00, 101.56it/s]
 [('athena_pp', 'Bcc1'), ('athena_pp', 'Bcc2'), ('athena_pp', 'Bcc3'), ('athena_pp', 'press'), ('athena_pp', 'rho'), ('athena_pp', 'vel1'), ('athena_pp', 'vel2'), ('athena_pp', 'vel3')] 

Version Information

• Operating System:
• Python Version:
• yt version:
• Other Libraries (if applicable): Python 2.7.5 yt 3.4-dev Scientific Linux 7.8

Installed yt and python from source.

[neutrinoceros]

Hi @smressle, thank you for reporting this. I see that you installed a pretty old version of yt (and Python), but I report that, with your dataset, I'm able to reproduce this even with Python 3.9.5 + the current dev version of yt. I will now dig into this.

[neutrinoceros]

@jzuhone it looks like you're the person who last touched the part of the code that fail (5yrs ago). It's a bit hard for me to tell wether the code is buggy (despite it being apparently very stable), or if something about the format of this particular datafile has grown out of compatibility with your frontend. From what I'm seeing, the conditional where the crash happens should not be hit in the first place, though I'm not 100% sure.

[smressle]

Thanks for looking into this! In case it is relevant, I am using Athena++ v 1.1.1 https://github.com/PrincetonUniversity/athena-public-version/releases/tag/v1.1.1

[jzuhone]

This is a known issue because this dataset is using a non-uniform cell-spacing. I tried to hack something up which worked for this using a hexahedral mesh, but there were always corner cases which didn't work and it is also very slow.

The real fix for this needs yt to support cell sizes which are constant in logspace, which I believe @neutrinoceros and @matthewturk started working on.

We can try to get this version working for @smressle's dataset, but I won't be able to get to it very soon.

[neutrinoceros]

Thank you for this insight. Indeed logspacing is something Matt has been working on and I'm most interested in it, but I also believe that it's a relatively long way there, sorry about that @smressle

[smressle]

That's OK, thanks both for looking into it for me. Hopefully I can find another way to deal with the data.

Cheers! Sean

On Tue, Jun 29, 2021 at 2:05 PM Clément Robert @.***> wrote:

Thank you for this insight. Indeed logspacing is something Matt has been working on and I'm most interested in it, but I also believe that it's a relatively long way there, sorry about that @smressle https://github.com/smressle

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/yt-project/yt/issues/3389#issuecomment-870805150, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6S4RW5RKBJI2B4WMFWZUTTVIDNNANCNFSM47QEA47A .

[matthewturk]

@smressle I'm sorry that this hasn't been fixed yet, but we have been working on it. I hope to provide more information soon!

[smressle]

No worries, Matthew. Thankfully I was able to find a workaround for my specific application. It's a difficult problem!

On Tue, Jan 11, 2022 at 11:36 AM Matthew Turk @.***> wrote:

@smressle https://github.com/smressle I'm sorry that this hasn't been fixed yet, but we have been working on it. I hope to provide more information soon!

— Reply to this email directly, view it on GitHub https://github.com/yt-project/yt/issues/3389#issuecomment-1010147278, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6S4RUCJLJASNNPXGWUOVTUVRL7ZANCNFSM47QEA47A . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you were mentioned.Message ID: @.***>

[jzuhone]

This has likely been resolved by https://github.com/yt-project/yt/pull/4562. Given how long it's been, I don't know if @smressle wants to try again, but I'm closing this unless it needs reopening.

[jzuhone]

This has likely been resolved by https://github.com/yt-project/yt/pull/4562. Given how long it's been, I don't know if @smressle wants to try again, but I'm closing this unless it needs reopening.

[neutrinoceros]

Confirmed. This now works correctly !

[smressle]

Awesome, thank you all! Yes I no longer am in immediate need of this feature but I have heard of others that are and I am sure it will come in handy in the future.

On Mon, Jul 10, 2023 at 9:45 AM Clément Robert @.***> wrote:

Confirmed. This now works correctly !

— Reply to this email directly, view it on GitHub https://github.com/yt-project/yt/issues/3389#issuecomment-1629001364, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6S4RQSODRQ623RH7L3APTXPQBN3ANCNFSM47QEA47A . You are receiving this because you were mentioned.Message ID: @.***>