Problem plotting PLUTO Sedov 1D Spherical

dutta-alankar commented 1 year ago

On trying to load and then plot 1D spherical using the following code,

import yt
ds = yt.load("./data.%04d.vtk"%1)
p = yt.ProfilePlot(
    ds, "radius", [("gas", "density")], weight_field=None, accumulation=True
)

I get following error and I h'm attaching the data. pluto_sedov.zip

/home/alankar/Documents/myrepos/yt_idefix/yt_idefix/_io/commons.py:72: RuntimeWarning: overflow encountered in square
  r = np.sqrt(xcart[:, 0, 0] ** 2 + ycart[:, 0, 0] ** 2)
/home/alankar/Documents/myrepos/yt_idefix/yt_idefix/_io/commons.py:77: RuntimeWarning: overflow encountered in square
  ycart[0, :, 0] / np.sqrt(xcart[0, :, 0] ** 2 + ycart[0, :, 0] ** 2)
yt : [INFO     ] 2023-02-27 19:11:01,119 Parameters: current_time              = 0.5
yt : [INFO     ] 2023-02-27 19:11:01,119 Parameters: domain_dimensions         = [400   1   1]
yt : [INFO     ] 2023-02-27 19:11:01,119 Parameters: domain_left_edge          = [       nan 1.57079637 0.78539819]
yt : [INFO     ] 2023-02-27 19:11:01,120 Parameters: domain_right_edge         = [       nan 2.57079637 1.78539819]
yt : [INFO     ] 2023-02-27 19:11:01,120 Parameters: cosmological_simulation   = 0
test <module 'yt_idefix._io.vtk_io' from '/home/alankar/Documents/myrepos/yt_idefix/yt_idefix/_io/vtk_io.py'>
/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/numpy/lib/arraysetops.py:89: RuntimeWarning: invalid value encountered in subtract
  return ary[1:] - ary[:-1]
Traceback (most recent call last):
  File "/home/alankar/Documents/myrepos/data-pluto/tests/pluto_sedov/test.py", line 6, in <module>
    p = yt.ProfilePlot(
  File "/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/yt/visualization/profile_plotter.py", line 211, in __init__
    create_profile(
  File "/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/yt/data_objects/profiles.py", line 1271, in create_profile
    bin_fields = data_source._determine_fields(bin_fields)
  File "/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/yt/data_objects/data_containers.py", line 1473, in _determine_fields
    finfo = self.ds._get_field_info(ftype, fname)
  File "/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/yt/data_objects/static_output.py", line 858, in _get_field_info
    field_info, candidates = self._get_field_info_helper(ftype, fname)
  File "/home/alankar/Documents/myrepos/yt_idefix/venv/lib/python3.9/site-packages/yt/data_objects/static_output.py", line 978, in _get_field_info_helper
    raise YTFieldNotFound(field=INPUT, ds=self)
yt.utilities.exceptions.YTFieldNotFound: Could not find field ('unknown', 'radius') in data.0001.vtk.

dutta-alankar commented 1 year ago

Perhaps this is because 1d spherical is not treated in the same way as 1d cartesian.

dutta-alankar commented 1 year ago

This might be due to bug in PLUTO itself. See my commit in h5_io.py (something similar is needed for vtk files as well) which somewhat fixes this, but this is not ideal. Perhaps parsing grid.out in this special case of 1d spherical/polar is the better idea.

I'm also putting in the relevant problematic lines in PLUTO from hdf5_io.c (also present in vtk io codes):

#if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
      DIM_EXPAND(node_coords[IDIR][kk][jj][ii] = (float)x1;  ,
               node_coords[JDIR][kk][jj][ii] = (float)x2;  ,
               node_coords[KDIR][kk][jj][ii] = (float)x3;)
      #elif GEOMETRY == POLAR
      DIM_EXPAND(node_coords[IDIR][kk][jj][ii] = (float)(x1*cos(x2)); ,
               node_coords[JDIR][kk][jj][ii] = (float)(x1*sin(x2)); ,
               node_coords[KDIR][kk][jj][ii] = (float)(x3);)
      #elif GEOMETRY == SPHERICAL
        #if DIMENSIONS <= 2
        DIM_EXPAND(node_coords[IDIR][kk][jj][ii] = (float)(x1*sin(x2)); ,
                 node_coords[JDIR][kk][jj][ii] = (float)(x1*cos(x2)); ,
                 node_coords[KDIR][kk][jj][ii] = 0.0;)
        #else
        DIM_EXPAND(node_coords[IDIR][kk][jj][ii] = (float)(x1*sin(x2)*cos(x3)); ,
                 node_coords[JDIR][kk][jj][ii] = (float)(x1*sin(x2)*sin(x3)); ,
                 node_coords[KDIR][kk][jj][ii] = (float)(x1*cos(x2));)
        #endif
      #else

Since for 1D left edge of x2=0 in spherical for the default range 0 to 1, all node coordinates become 0 as sin(x2)=0.

neutrinoceros commented 1 year ago

I'll look into this tonight if I can

neutrinoceros commented 1 year ago

So I see two issues to be discussed here: 1) I am actually not sure "radius" is a canonically valid field name in spherical coordinates in yt ("r" works) 2) In the log you provided we see that domain edges are nan on both sides in the first direction (which I assume is "r"). Inspection shows that cartesian coordinates that are read from file already contain invalid data (nan as well as what looks like uninitialised memory). Since the vtk format is very strict about precision and endianness, I'd be enclined to think the data is actually corrupted and that this isn't a bug in yt/yt_idefix but indeed a bug in Pluto itself, as you already hinted to.

neutrinoceros commented 1 year ago

@dutta-alankar are you able to visualise this dataset with any other method than yt ? I'm not convinced this is actionable here.

xshaokun commented 1 year ago

I checked the vtk io part in write_vtk.c ( version 4.4-patch2), and I'm surprised that the coordinate information was ignored to update before being written into vtk files in the case of 1D spherical geometry.

  if (node_coord == NULL) node_coord = ARRAY_2D(nx1 + INCLUDE_IDIR, 3, float);

  sprintf(header,"POINTS %d float\n", (nx1+INCLUDE_IDIR)*(nx2+INCLUDE_JDIR)*(nx3+INCLUDE_KDIR));
  VTK_HEADER_WRITE_STRING(header);

/* -- Write structured grid information -- */

  x1 = x2 = x3 = 0.0;
  for (k = 0; k < nx3 + INCLUDE_KDIR; k++){
  for (j = 0; j < nx2 + INCLUDE_JDIR; j++){ 
    for (i = 0; i < nx1 + INCLUDE_IDIR; i++){
      DIM_EXPAND(x1 = grid->xl_glob[IDIR][IBEG + i];  ,
                 x2 = grid->xl_glob[JDIR][JBEG + j];  ,
                 x3 = grid->xl_glob[KDIR][KBEG + k];)

      #if (GEOMETRY == CARTESIAN) || (GEOMETRY == CYLINDRICAL)
      node_coord[i][0] = x1;
      node_coord[i][1] = x2;
      node_coord[i][2] = x3;
      #elif GEOMETRY == POLAR
      node_coord[i][0] = x1*cos(x2);
      node_coord[i][1] = x1*sin(x2);
      node_coord[i][2] = x3;
      #elif GEOMETRY == SPHERICAL
      #if DIMENSIONS == 2
      node_coord[i][0] = x1*sin(x2);
      node_coord[i][1] = x1*cos(x2);
      node_coord[i][2] = 0.0;
      #elif DIMENSIONS == 3
      node_coord[i][0] = x1*sin(x2)*cos(x3);
      node_coord[i][1] = x1*sin(x2)*sin(x3);
      node_coord[i][2] = x1*cos(x2);
      #endif
      #endif

      if (IsLittleEndian()){
        SWAP_VAR(node_coord[i][0]);
        SWAP_VAR(node_coord[i][1]);
        SWAP_VAR(node_coord[i][2]);
      }
    }
    VTK_HEADER_WRITE_FLTARR(node_coord[0],3*(nx1+INCLUDE_IDIR));
  }}

Referring to the source code above, node_coord was directly written into vtk files with default values (i.e. zero in my output, based on Test_Problem/HD/Sedov/defitions_03.h)

So in this case, we have to read the grid coordinates from grid.out. @neutrinoceros @dutta-alankar

dutta-alankar commented 1 year ago

I checked the vtk io part in write_vtk.c ( version 4.4-patch2), and I'm surprised that the coordinate information was ignored to update before being written into vtk files in the case of 1D spherical geometry.

  if (node_coord == NULL) node_coord = ARRAY_2D(nx1 + INCLUDE_IDIR, 3, float);

  sprintf(header,"POINTS %d float\n", (nx1+INCLUDE_IDIR)*(nx2+INCLUDE_JDIR)*(nx3+INCLUDE_KDIR));
  VTK_HEADER_WRITE_STRING(header);

/* -- Write structured grid information -- */

  x1 = x2 = x3 = 0.0;
  for (k = 0; k < nx3 + INCLUDE_KDIR; k++){
  for (j = 0; j < nx2 + INCLUDE_JDIR; j++){ 
    for (i = 0; i < nx1 + INCLUDE_IDIR; i++){
      DIM_EXPAND(x1 = grid->xl_glob[IDIR][IBEG + i];  ,
                 x2 = grid->xl_glob[JDIR][JBEG + j];  ,
                 x3 = grid->xl_glob[KDIR][KBEG + k];)

      #if (GEOMETRY == CARTESIAN) || (GEOMETRY == CYLINDRICAL)
      node_coord[i][0] = x1;
      node_coord[i][1] = x2;
      node_coord[i][2] = x3;
      #elif GEOMETRY == POLAR
      node_coord[i][0] = x1*cos(x2);
      node_coord[i][1] = x1*sin(x2);
      node_coord[i][2] = x3;
      #elif GEOMETRY == SPHERICAL
      #if DIMENSIONS == 2
      node_coord[i][0] = x1*sin(x2);
      node_coord[i][1] = x1*cos(x2);
      node_coord[i][2] = 0.0;
      #elif DIMENSIONS == 3
      node_coord[i][0] = x1*sin(x2)*cos(x3);
      node_coord[i][1] = x1*sin(x2)*sin(x3);
      node_coord[i][2] = x1*cos(x2);
      #endif
      #endif

      if (IsLittleEndian()){
        SWAP_VAR(node_coord[i][0]);
        SWAP_VAR(node_coord[i][1]);
        SWAP_VAR(node_coord[i][2]);
      }
    }
    VTK_HEADER_WRITE_FLTARR(node_coord[0],3*(nx1+INCLUDE_IDIR));
  }}

Referring to the source code above, node_coord was directly written into vtk files with default values (i.e. zero in my output, based on Test_Problem/HD/Sedov/defitions_03.h)

So in this case, we have to read the grid coordinates from grid.out. @neutrinoceros @dutta-alankar

@xshaokun Yes, this is why parsing grid.out is what I was suggesting for 1d with any geometry that is non-cartesian. This issue is present both in vtk and hdf5 dumps.

dutta-alankar commented 1 year ago

@dutta-alankar are you able to visualise this dataset with any other method than yt ? I'm not convinced this is actionable here.

Sorry I missed relying on this. The field values are fine if I read it using just h5py in python and then read in the grid values from grid.out manually. But since 1d is so easy to deal with the manual coding isn't much of an issue.

neutrinoceros commented 1 year ago

I'm surprised that the coordinate information was ignored to update before being written into vtk files in the case of 1D spherical geometry.

ok then this is clearly a bug in Pluto. I'm okay with to reading grid.out but I would prefer it be kept as a fallback strategy , and we don't do it unless we absolutely must. I think we should be able to run a couple validation steps to decide wether the grid data from the .vtk file is sound or not, for instance : are all numbers valid floats (no infs, no NaNs) and are they sorted ?

neutrinoceros commented 1 year ago

Oh and can someone please report this bug to Andrea Mignone ?

neutrinoceros commented 1 year ago

(FTR the bug has been fixed in Idefix a year and a half ago, but it seems that we missed the opportunity to report upstream)

xshaokun commented 1 year ago

It is strange that the values of grid data in the .vtk file provided by @dutta-alankar are not all zero, and even include nan. I checked the values by the following script:

import numpy as np
with open("data.0001.vtk", "rb") as f:
    f.seek(0)
    for _ in range(6):
        line = next(f).decode()
        print(line)
    coords = np.fromfile(f, dtype=">f", count=3*401)
    print(coords[2::3])

I guess it was supposed to be all zero like mine since the geometry setups are the same. @dutta-alankar Do you have any idea about this?

xshaokun commented 1 year ago

Oh and can someone please report this bug to Andrea Mignone ?

I can report this to the google group of Pluto. But I'm not sure if it'll be addressed since the group seems not active. Or should I just email it to Andrea Mignone directly?

neutrinoceros commented 1 year ago

t is strange that the values of grid data in the .vtk file provided by @dutta-alankar are not all zero, and even include nan.

If it's uninitialised memory it can be anything that was just sitting in RAM, not necessarily zeros.

Or should I just email it to Andrea Mignone directly?

From what I could gather yes, emailing Andrea directly is the way to go.

dutta-alankar commented 1 year ago

Oh and can someone please report this bug to Andrea Mignone ?

I can report this to the google group of Pluto. But I'm not sure if it'll be addressed since the group seems not active. Or should I just email it to Andrea Mignone directly?

Yes the PLUTO google group is extremely inactive. E-mailing Andrea is a good option.

dutta-alankar commented 1 year ago

It is strange that the values of grid data in the .vtk file provided by @dutta-alankar are not all zero, and even include nan. I checked the values by the following script:
import numpy as np
with open("data.0001.vtk", "rb") as f:
    f.seek(0)
    for _ in range(6):
        line = next(f).decode()
        print(line)
    coords = np.fromfile(f, dtype=">f", count=3*401)
    print(coords[2::3])
I guess it was supposed to be all zero like mine since the geometry setups are the same. @dutta-alankar Do you have any idea about this?

This has probably something to do with the conversion from non-cartesian to cartesian. Can you tell me what was geometry and problem for which this is hapenning?

xshaokun commented 1 year ago

Can you tell me what was geometry and problem for which this is hapenning?

It's actually your output provided in this issue, but never mind, I got it now.

BTW, I just sent an email to Andrea and cc you two. Did you receive it? @neutrinoceros @dutta-alankar I'd like to check it because I'm not sure I've sent it successfully. The email doesn't appear in my sent box.

neutrinoceros commented 1 year ago

I did not receive it.

dutta-alankar commented 1 year ago

Can you tell me what was geometry and problem for which this is hapenning?

It's actually your output provided in this issue, but never mind, I got it now.

BTW, I just sent an email to Andrea and cc you two. Did you receive it? @neutrinoceros @dutta-alankar I'd like to check it because I'm not sure I've sent it successfully. The email doesn't appear in my sent box.

I didn't receive any email. My email id is: alankardutta@iisc.ac.in

dutta-alankar commented 1 year ago

Can you tell me what was geometry and problem for which this is hapenning?

It's actually your output provided in this issue, but never mind, I got it now.

BTW, I just sent an email to Andrea and cc you two. Did you receive it? @neutrinoceros @dutta-alankar I'd like to check it because I'm not sure I've sent it successfully. The email doesn't appear in my sent box.

Got your mail right now.

xshaokun commented 1 year ago

Well, I just resent the email. It should succeed this time. @neutrinoceros @dutta-alankar

dutta-alankar commented 1 year ago

Well, I just resent the email. It should succeed this time. @neutrinoceros @dutta-alankar

@xshaokun Just wanted explicitly have this mentioned here so that we don't forget testing it later on. This issue is present for not only spherical 1d but any non-cartesian 1d setups. Also not only is vtk io affected but also hdf5 io.

xshaokun commented 1 year ago

This issue is present for not only spherical 1d but any non-cartesian 1d setups. Also not only is vtk io affected but also hdf5 io.

Oh sorry, I only focused on the uninitialized issue and overlooked other parts, since it is well-treated in your xdmf loader.

I still haven't figured out how this goes wrong. @dutta-alankar Could you please send another email quoting mine to describe and report this?

dutta-alankar commented 1 year ago

Looking at the code I can understand the occurrence of zeros but cannot figure out the NaNs. This needs a bit more investigation.

neutrinoceros commented 1 year ago

To summarize

we've identified that the root issue was with PLUTO itself
it only affects 1D simulations, which I expect to be pretty cheap to rerun if needed.
working around it on our side is possible in principle but would require what feels like a disproportionate amount of work: we'd need to validate coordinates read from vtks, and reconstruct them from the parameter file in case they are found to be invalid, which is complicated in the general case.

I'm inclined to think that adding minimal validation and a warning linking to this discussion would suffice to close this issue.

dutta-alankar commented 1 year ago

To summarize

we've identified that the root issue was with PLUTO itself

it only affects 1D simulations, which I expect to be pretty cheap to rerun if needed.

working around it on our side is possible in principle but would require what feels like a disproportionate amount of work: we'd need to validate coordinates read from vtks, and reconstruct them from the parameter file in case they are found to be invalid, which is complicated in the general case.

I'm inclined to think that adding minimal validation and a warning linking to this discussion would suffice to close this issue.

I think just a warning is a good idea.

neutrinoceros / yt_idefix

Problem plotting PLUTO Sedov 1D Spherical #210