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

[Barry Smith]

> On Mar 20, 2015, at 1:17 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Fri, Mar 20, 2015 at 12:08 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> > On Mar 20, 2015, at 11:31 AM, Håkon Strandenes <haakon at hakostra.net> wrote:
> >
> > Yes, in the HDF5 file. Currently a dof=1 DMDA Vec is saved as a NZ x NY x NX dataset, wile a dof > 1 is saved as NZ x NY x NX x dof in the HDF5 file.
> 
>   Actually a 1d DMDA is stored as a NX dataset, a 2d DMDA is stored as a NY x NX data set and a 3d DMDA is stored as a NZ x NY x NX data set and an extra dimension is added if dof > 1.   And I think this is great and intuitively obvious for users.
> 
>    I agree with Hakon and advocate NOT sticking the extra dimension in when bs is 1. Here is Matt's entire commit
> 
> 
>    maxDims[dim] = dims[dim];
>    chunkDims[dim] = dims[dim];
>    ++dim;
> -  if (bs > 1) {
> +  if (bs >= 1) {
>      dims[dim]    = bs;
>      maxDims[dim] = dims[dim];
>      chunkDims[dim] = dims[dim];
>    }
>    count[dim] = PetscHDF5IntCast(xin->map->n)/bs;
>    ++dim;
> -  if (bs > 1) {
> +  if (bs >= 1) {
>      count[dim] = bs;
>      ++dim;
>    }
>    }
>    offset[dim] = PetscHDF5IntCast(low/bs);
>    ++dim;
> -  if (bs > 1) {
> +  if (bs >= 1) {
>      offset[dim] = 0;
>      ++dim;
>    }
> 
>   (Then there is Karl's one line fix for reading HDF5 files.)
> 
>   Arguing that handling bs = 1 as a special case greatly (or even modestly) increases the complexity of the PETSc source code is nonsense.
> 
>   So I want a stronger case presented then
> 
> >> My rationale was to simplify our code path, and have one way to handle fields with components. I did not
> >> want to special case a field with one component.
> 
> If we special case here, we must special case in everything we do, such as automatic generation of Xdmf.
> I would rather say that the PETSc format always has the bs in it, which makes it easier to process and
> consistent.

   From what I see the "easier" is just trivially easier; having a tiny bit more complicated code that is intuitively clear for users seems a better choice.

> However, we allowing reading of of a set without block size.

  Hmm, according to the bug report which I confirmed, VecLoad() refuses to load in a 1d dataset
> 
> We would have the ability to write without blocksize by setting the blocksize to -1.

  Ok, so now you are still supporting that "special case" but with some totally non-intuitive value for a block size that will crash if used with the vector in any other usage of the vector. How is this either simpler code or more intuitive code.

 As Hakon notes also, VecView from a DMDA is handled differently with regard to the bs dimension then VecView from a simple "vector". Which is terrible.

Barry

> 
>    Matt
>  
>   otherwise I think we should undo Matt's change.
> 
>   Barry
> 
> 
> 
> >
> > In my humble opinion, it is best to leave the DMDA's as it is, and fix up the code related to the plain Vecs. That code is already full of "if (bs >= 1) dim++" etc. which is useless (as Barry points out).
> >
> > My final point is that one should handle this in a CONSISTENT manner for both "plain Vecs" and "DMDA Vecs", i.e. decide wither one want to always add one dimension for bs and dof or leave that dimension out for the special case with bs=1/dof=1.
> >
> > Håkon
> >
> >
> > On 20. mars 2015 16:57, Jed Brown wrote:
> >> Håkon Strandenes <haakon at hakostra.net> writes:
> >>> I mean that if you decide to stick with the current way of
> >>> reading/writing a "plain 1D" Vec with N elements as a N x 1
> >>> array/dataset, you should treat a dof=1 DMDA Vec as an NZ x NY x NX x 1
> >>> array/dataset.
> >>
> >> In the HDF5 file?  I agree.
> >>
> 
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener