Nek5000 reader

[ashwinvis]

I glanced through your tutorial yesterday. I tried to output a .vtp file after loading the Nek5000 field file in paraview via the visnek utility. It does not seem to directly work with turbulucid's Case class.

I have a rough sketch for a VTK interface generated with help of Paraview here: https://github.com/exabl/sandbox/blob/master/paraview/nekio.py. It might be interesting to have it here.

[timofeymukha]

Hi Ashwin! I did some testing. If you use case = tbl.Case(vtp_path, pointData=True) things work. However, turbulucid is meant to work with cell-centred data, so under the hood the point values at the gll points will be interpolated to the "cell centers", as vtk sees them. This, however gives a rather false representation of how the solution actually looks. So my feeling is that turbulucid can be suitable if you have a very dense mesh and you want to make an overview-type of a 2d plot of some field. In this case the interpolation artifacts can perhaps be ignored. But for extracting data, e.g. linear profiles, the inaccuracy due to interpolation is probably too large to accept. But of course you can use that to get a quick look at your data.

I also did some testing with sts* files generated by the KTH stats toolbox. Here the situation is a bit more difficult, if you save the data from paraview, it will be a vtu file, not vtp. And actually it is a stitched dataset, with separate geometry for each element stitched together. I added some code (now on develop) that will transform this to vtkPolyData, but the cuts along the elements are still there of course. This introduces artifacts with plotting boundaries and also with data extraction. I guess this is fixable through some more advanced manipulation of the data prior to loading it to turbulucid.

[timofeymukha]

Example of a channel flow slice. It's a very coarse mesh, so you clearly see the cells.

And here is how a plot of the sts* file looks, the elements are visible

You can fix the above by passing plotBoudnaries=False to tbl.plot_field but there is still problems with extracting data etc.

[ashwinvis]

Very cool. I tried with the periodic hill example's initial condition, sliced using Paraview. It has funky boundaries:

So there are some "rough edges", but for simple cases it might work.

[timofeymukha]

@ashwinvis In Paraview, untick the "Triangulate the slice" to get rectangular cells.

[ashwinvis]

Thank you so much. Do you think adding a Paraview Python interface to slice and generate VTK files / arrays in memory would be useful for turbulucid?

[timofeymukha]

I'm not sure if the actual slicing should be included, but one could at least add a "from memory" reader, which accepts a vtkPolyData object with a slice somehow obtained elsewhere and passes it to the Case object. The latter will need to be modified a bit, since the way the reader is selected now is based on a file extension, but it should be pretty easy to fix.

One could consider the "from memory" reader to have two constructors, one from a vtkPolyData and one from a filename and some kind of definition of the cutting plane (point and normal). The latter would perform the cutting. However, I would stick to just vtk and not use paraview for the slice. No need for an extra dependency, I think.