TL: generic FieldsIO class for n-dimensional rectilinear grids + VTK support

[tlunet]

pySDC.helpers.fieldsIO reformating

Cart1D and Cart2D have been replace by one Rectilinear class that does the same for all dimensions. It's not backward compatible, but since I was the only one using the old classes, it should be fine ...

⚠️ The new class uses some python 3.11+ features, allowing to extract slices of a ndarray by unpacking a list of slices. For instance :

a = np.arange(2*3*4).reshape((2, 3, 4))
s = [0, slice(None), slice(2)]
a[*s]  # equivalent to a[0, :, :2]

While this is pretty neat to handle n-dimensional data without making the code too complicated, it's still seen as a syntax error in python 3.10 and lower. Hence the module requires now python 3.11 and higher, but it does not prevent to use pySDC with lower version (just do not import pySDC.helpers.fieldsIO ...)

VTK support

Two utility functions added in pySDC.helpers.vtkIO, namely writeToVTR and readFromVTR that allow to go from numpy array to VTR files and vice-versa. Also, a toVTR method has been added to the pySDC.helpers.fieldsIO.Rectilinear class such that conversion of binary file storing fields can be easily done like this :

# Suppose the FieldsIO object is already writen into outputs.pysdc
import os
from pySDC.utils.fieldsIO import Rectilinear
os.makedirs("vtrFiles")  # to store all VTR files into a subfolder
Rectilinear.fromFile("outputs.pysdc").toVTR(baseName="vtrFiles/field", varNames=["u", "v", "w", "T", "p"])

It allows to separate IO during simulation, that can be easilly done in parallel with the FieldsIO features, from the generation of VTR files, that would have been way trickier in space parallel ... Also, this VTR IO feature is now only enabled for 3D fields only, although it could be adapted for 1D or 2D fields (but use rather matplotlib for that, not VTK 😮‍💨)

So thanks to this PR, @brownbaerchen will be able to generate nice 3D movies that could make @danielru happy 🥳

[brownbaerchen]

You can make your code compatible with older Python versions with a slightly less clean syntax:

a = np.arange(2*3*4).reshape((2, 3, 4))
s = [0, slice(None), slice(2)]
a[*s]  # only from 3.11
a[(*s,)]  # works also on older Python versions 

[tlunet]

Oh thanks, I didn't find this yesterday ... but do we need that ? 3.10 will be in end of life next year 😅

[brownbaerchen]

I didn't read too carefully, but looks good to me.
I personally would use the older syntax compatible with 3.10. But I also wouldn't insist on it.
About adding vtk as a dependency: Could that become a problem on some laptops?

[brownbaerchen]

Is this really the number of points in y direction or is it in all directions?

[tlunet]

it's in all directions indeed ... gridSizes instead ?

[brownbaerchen]

:D

[tlunet]

Accidental reverting of a commit that corrected it ... I'll change it if I do another commit

[tlunet]

I didn't read too carefully, but looks good to me. I personally would use the older syntax compatible with 3.10. But I also wouldn't insist on it. About adding vtk as a dependency: Could that become a problem on some laptops?

vtk seems to be packed for most distribution with conda forge, see https://anaconda.org/conda-forge/vtk ... so I would suppose it'll be fine ? We already have dill as a base dependency, I figured that vtk was not too much more

[tlunet]

If no one has any issue anymore, then it's ready to be merged.

[pancetta]

You said "I was the only one using the old classes", but then you need to add vtk as additional base dependency? Plus the code is not backward compatible? To me this sounds more like a project-based addon, so that the impact of this is limited to this project.

[tlunet]

You said "I was the only one using the old classes", but then you need to add vtk as additional base dependency? Plus the code is not backward compatible? To me this sounds more like a project-based addon, so that the impact of this is limited to this project.

Those are two different aspects.

First : there were two classes before, Cart1D and Cart2D classes (wrongly named) for rectilinear 1D and 2D grids. Now the Rectilinear class can do any dimension (and avoid the need of those old classes), in particular for 3D fields so Thomas can use it for parallel solutions output if he manages to adapt the Rayleigh Benard for 3D problems. The only issue is that output files generated by the Cart2D classes cannot be read with the new Rectilinear class, but since I'm the only one that generated such output files before, I can manage on my own the transition to the more generic format. And now, the Rectilinear class can be used in combination to more problems implemented in pySDC, since it allows 3D grids (or even 4D grids, etc ...)

Second : the vtr dependency is only an add-on that can be used independently from the FieldsIO classes. I just added a link between the Rectilinear class because it was quite simple, and makes the generation of vtr files for videos easier. And this can be used for any project that generate 3D rectilinear grid solutions.

[brownbaerchen]

I am generally in favor of having streamlined output in pySDC. To me, this is part of the feature set it should have for prototyping. If you want to test out a new method on the existing problems and then have to work out how to plot the results, even though other people have done that before, that is just a waste of time.
For instance, when I was working with Gusto, I was very appreciative of having the plotting scripts and even reference plots be part of the repository so I could run the test and see if the solution looks as expected without having to deal with learning ParaView. While we don't need reference plots in pySDC because we don't have this kind of test case, having the plotting facilities is useful. I already started adding some plot functions to some problem classes, but this is very rudimentary at this point.
I do like what Thibaut is doing, which is to define a common data structure for pySDC that makes it easy to deal with data. How should you write to disk? Use the pySDC class for that. How to visualize the result? Use whatever ParaView, pyvista, ... you like. No more thousands of pickle files with later untangling the grid from the data and crying because ParaView won't take my numpy.

[pancetta]

So, the idea is to establish this as the default way to dump files for visualization out of pySDC, yes?

[tlunet]

If you are fine with the current minimalist structure, yes 😉

[brownbaerchen]

We should write a hook for doing this as well. But in a follow up PR..

[tlunet]

We should write a hook for doing this as well.

That would be your area of expertise then, not mine 😅

[pancetta]

Good idea. I like this approach and I'm happy to have this as the default way of doing things. Can we add a hook and a demo to this PR, then?

[tlunet]

I can add a more extended documentation + demo of the data IO structure, but concerning the hook class I would prefer letting @brownbaerchen do that (he'll probably do that better than me)

[brownbaerchen]

I can do a hook, but I would prefer doing that in a second PR. This diff is already quite large and this makes sense by itself already.

[tlunet]

So, demo + extended doc + hooks in a next PR, I'm also fine with that. What about you @pancetta ?