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

[Håkon Strandenes]

Thanks. Unfortuneatly, I am not in my office until april 7th or 8th, and 
am not able to test this out until then.

I have also seen the other replies to this post. I might have some ideas 
on how to solve this in a bit more flexible way, but again, I cannot 
work on this until I am back in office.

I will follow this up as soon as I am back. Thanks for your help so far.

Regards,
Håkon Strandenes


Den 2015-03-29 04:59, skrev Barry Smith:
> 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