Shape error when reading idefix-vtk V2.2.00

[aborderies]

Last week I updated my idefix repository to version 2.2.00, which adds native coordinates to vtk files. Since then, I've had a problem when I try to load an idefix-vtk file: when yt_idefix wants to calculate the cell_widths, a shape error occurs:
ValueError: could not broadcast input array from shape (N-1,) into shape (N,)
where N is the number of cells in the first dimension of my grid, at this line

yt_idefix/src/yt_idefix/data_structures.py

Line 209 in 89d9471

cell_widths[idir][:] = np.ediff1d(edges)

Looking at the code, I noticed that the read_grid_coordinates() method, which is supposed to return the coordinates of the cell edges (if I understood correctly), actually returns the coordinates of the cell centers when the native coordinates are in the metadata :

yt_idefix/src/yt_idefix/_io/vtk_io.py

Lines 234 to 241 in 89d9471

	if "native_coordinates" in md:
	nc = md["native_coordinates"]
	return Coordinates(
	nc["X1C_NATIVE_COORDINATES"],
	nc["X2C_NATIVE_COORDINATES"],
	nc["X3C_NATIVE_COORDINATES"],
	array_shape,
	)

I tried replacing XkC_NATIVE_COORDINATES with XkL_NATIVE_COORDINATES and this avoided the problem in my case. Does this work, or would another solution be more appropriate?

Thanks in advance for the help.

[neutrinoceros]

Hi @Anto6453, thanks for your report. I should be able to reproduce this myself, but just in case, could you confirm that you're running the latest version of yt_idefix, and perhaps provide an example setup that you can share ? (Maybe the ini file and definitions.cpp is enough)
I'll take a deeper look later today.

[aborderies]

Hi @neutrinoceros, thanks for the answer.
I'm running on yt_idefix 2.4.1 (installed with PyPI).
Here is the setup I use: Collapse_Setup.zip

[neutrinoceros]

Thanks ! anything special I need to do to compile your setup ? Doing the naive thing doesn't seem to work (I did check out idefix 2.2.00)

/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup/resistivities.cpp:40:24: error: no member named 'dt_hydro' in 'DataBlock'
   40 |   real dt_hydro = data.dt_hydro;
      |                   ~~~~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup/resistivities.cpp:141:24: error: no member named 'dt_hydro' in 'DataBlock'
  141 |   real dt_hydro = data.dt_hydro;
      |                   ~~~~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:275:27: error: no member named 'centralMoment' in 'DumpImage'
  275 |     LCentralGlob = image->centralMoment;
      |                    ~~~~~  ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:277:74: error: no member named 'centralMoment' in 'DumpImage'
  277 |     idfx::cout << "InnitFlow: Old Central Angular Momentum = " << image->centralMoment << std::endl;
      |                                                                   ~~~~~  ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:420:17: error: no member named 'selfPotential' in 'DataBlockHost'
  420 |               d.selfPotential(k+d.loffset[KDIR],j+d.loffset[JDIR],i+d.loffset[IDIR]) = InterpolateFromDump2D(image->arrays["GravPotential"],
      |               ~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:420:35: error: no member named 'loffset' in 'DataBlockHost'
  420 |               d.selfPotential(k+d.loffset[KDIR],j+d.loffset[JDIR],i+d.loffset[IDIR]) = InterpolateFromDump2D(image->arrays["GravPotential"],
      |                                 ~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:420:53: error: no member named 'loffset' in 'DataBlockHost'
  420 |               d.selfPotential(k+d.loffset[KDIR],j+d.loffset[JDIR],i+d.loffset[IDIR]) = InterpolateFromDump2D(image->arrays["GravPotential"],
      |                                                   ~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:420:71: error: no member named 'loffset' in 'DataBlockHost'
  420 |               d.selfPotential(k+d.loffset[KDIR],j+d.loffset[JDIR],i+d.loffset[IDIR]) = InterpolateFromDump2D(image->arrays["GravPotential"],
      |                                                                     ~ ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:432:50: error: no member named 'centralMoment' in 'DumpImage'
  432 |                 real centralOmega = 5./2.*image->centralMoment/(image->centralMass*rIntOld*rIntOld);
      |                                           ~~~~~  ^
/Users/clm/dev/yt-project/yt_idefix/issues/366/Collapse_Setup/setup.cpp:551:88: error: no member named 'centralMoment' in 'DumpImage'
  551 |     idfx::cout << "InnitFlow: Central Angular Momentum if only reduction = " << image->centralMoment*rInt*rInt*rInt*rInt*rInt/(rIntOld*rIntOld*rIntOld*rIntOld*rIntOld) << std::endl;
      |                                                                                 ~~~~~  ^
2 errors generated.

[neutrinoceros]

In the mean time, I tried an even more naive reproduction using a test setup (HD/KHI) from idefix and didn't see any issue with yt.load("data.0000.vtk"), so I'm not even sure I know how the error is supposed to be triggered. Could you also provide a minimal reproducing Python script ?

[neutrinoceros]

Nevermind, using a setup more similar to yours (spherical + 2D) as a basis, and trying to create a yt.SlicePlot suffices to reprod. Let me try to find the correct fix now.

[neutrinoceros]

Hum, the way vtk_io.read_grid_coordinates is used throughout the code base, it seems evident that cell edges are the expected output, so now I'm wondering why I chose to use cell centers in #341... 🤔
Maybe I was just mixing up various needs as I worked on https://github.com/volodia99/nonos/pull/371 in parallel, and nonos usually prefers working with cell centers. Anyway, I'm now leaning towards your solution; I initially thought it was a hack, but now I think it's actually the correct thing to do ! Would you like to open a PR ?

[aborderies]

Ok I'll open one ! Thanks for your help !

[neutrinoceros]

closed by #368
note for future reference: you can link an issue in a PR by including "close #366" or "fix #366" in the PR description, so the issue is automatically closed when the PR is merged. Anyway, your patch is included in yt_idefix 2.4.2, so you should be able to upgrade immediately. Thanks again for taking the time to report and fix this issue !