[petsc-users] Load Vec from 1D HDF5 dataset

[Håkon Strandenes]

I have now had the opportunity to look into this again, and have made a 
pull request with a suggested change: 
https://bitbucket.org/petsc/petsc/pull-request/296/vecload_hdf5-can-without-setting-any/diff

In short terms, for loading Vecs, I suggest to do a "best effort guess" 
on the contents of the dataset in the HDF5 file, by doing the following:

1) Assume no extra dimension for bs = 1
2) If the dataset happens to have one dimension more than expected, 
check the size of the dimension corresponding to the bs.
3) If the size of the dimension corresponding to  the bs is == 1, accept 
the dataset with one extra dimension.

Since it can accept Vecs of "both kind" (with and without one extra 
dimension) it is 100% backwards compatible for *loading* Vecs.

The Vec output is not touched, and Barry's "basedimension2"-property is 
left for the output/write part.

The same logic could/should(!) be applied when reading DMDA's with dof=1 
as well, I just thought you might come with some critics on this first try.

Regards,
Håkon


On 29. mars 2015 04:59, Barry Smith wrote:
>
>    Hakon,
>
>      I have pushed to the branch barry/feature-hdf5-flexible-initial-dimension and next the change so that Vecs and Vecs obtained from DMCreateGlobalVector() with a DMDA will NOT have the extra dimension if BS == 1. To add that extra dimension even when bs == 1 with VecView() or to handle a file with that extra dimension with VecLoad() one must call PetscViewerHDF5SetBaseDimension2() or -viewer_hdf5_base_dimension2 true
>
>     Please try it out and let me know if you have any trouble
>
>     Thanks
>
>    Barry
>
>    Matt,
>
>     Scream as much as you want but adding that extra dimension automatically is not intuitive so it cannot be the default.
>
>> On Mar 25, 2015, at 10:36 AM, Håkon Strandenes <haakon at hakostra.net> wrote:
>>
>> Did you come to any conclusion on this issue?
>>
>> Regards,
>> Håkon
>>
>> On 20. mars 2015 22:02, Håkon Strandenes wrote:
>>> 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
>
>