[petsc-users] Load Vec from 1D HDF5 dataset MATT READ THIS EMAIL!

[Håkon Strandenes]

On 20. mars 2015 20:48, Barry Smith wrote:
> Why is 1 dimension a special case that is not worthy of its own format? The same thing would hold for 2d and 3d. One could then argue that we should have a single six dimensional format for the files for all vectors that PETSc produces. Then a 1d problem has five of the dimensions being 1.

This is a very good point, and support my view.

Let me come with two very simple example cases:


Case 1:
Create a list of grid points in an external preprocessor for the purpose 
of loading this into a Vec later:

   x = np.linspace(0.0, 1.0, num=21)
   f.create_dataset('/MESH/nodes/x', data=x)

vs.

   x = np.linspace(0.0, 1.0, num=21)
   x = x.reshape((21,1))
   f.create_dataset('/MESH/nodes/x', data=x)


Case 2:
Read three Vecs written to disk by PETSc, and calculate total "bounding 
box volume" of the grid:

   g = h5py.File('grid.h5', 'r')
   x = g['/MESH/nodes/x']
   y = g['/MESH/nodes/y']
   z = g['/MESH/nodes/z']
   Vol = (xp[-1] - xp[0])*(yp[-1] - yp[0])*(zp[-1] - zp[0])

vs.

   g = h5py.File('grid.h5', 'r')
   x = g['/MESH/nodes/x'][:,0]
   y = g['/MESH/nodes/y'][:,0]
   z = g['/MESH/nodes/z'][:,0]
   Vol = (x[-1] - x[0])*(y[-1] - y[0])*(z[-1] - z[0])


In both cases I think handling this extra, unnecessary dimension makes 
the code less attractive. It's not that either way is difficult, 
problematic or impossible, but it's just that 1D Vecs should intuitively 
be 1D datasets, and not 2D, 3D or 6D. This seriously confused me for 
quite a while until I figured this out, even after having written an 
entire Navier-Stokes DNS solver using the PETSc library for everything 
except time integration and filling these simple 1D coordinate arrays!

Regards,
Håkon