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

[Barry Smith]

> 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?

  No. Different people think a different default is best for different reasons.

  To get the conversation going again I'll propose an approach that allows both; PetscViewerHDF5SetSomething(viewer,PetscBool) where 
Something could be BaseDimension2. And the Load would handle files with either 1 or 2 dimensions as the base as well.

  Barry

> 
> 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