[Users] (possibly) inconsistent handling of ioffsetdenom in CarpetIOHDF5

[Erik Schnetter]

On Sun, May 4, 2014 at 4:39 PM, Roland Haas <rhaas at tapir.caltech.edu> wrote:

> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> Hello all,
>
> due to some issues we have with recovering checkpoints I find myself
> reading CarpetIOHDF5's Output.cc and Input.cc files.
>
> In particular lines 849--860 of Output.cc:
>  849   ivect const ioffsetdenom = bbox.stride();
>  850   HDF5_ERROR (attr = H5Acreate (dataset, "ioffsetdenom", \
>                                      H5T_NATIVE_INT,
>  851                                 dataspace, H5P_DEFAULT));
>  852   HDF5_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &ioffsetdenom[0]));
>  853   HDF5_ERROR (H5Aclose (attr));
>  854   ivect const ioffset =
>  855     ((bbox.lower() % bbox.stride()) + bbox.stride()) % \
>          bbox.stride();
>  856   HDF5_ERROR (attr = H5Acreate (dataset, "ioffset", H5T_NATIVE_INT,
>  857                                 dataspace, H5P_DEFAULT));
>  858   HDF5_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &ioffset[0]));
>  859   HDF5_ERROR (H5Aclose (attr));
>  860   ivect const iorigin = (bbox.lower() - ioffset) / bbox.stride();
> which as far as I can tell always just sets ioffsetdenom to the stride.
>
> On the other hand, Input.cc in line 1399-1400:
> 1399 assert (all (stride % patch->ioffsetdenom == 0));
> 1400 ivect lower = patch->iorigin * stride + patch->ioffset * stride / \
>                    patch->ioffsetdenom;
> only requires stride to be a multiple of ioffsetdenom (rather than
> being identical).
>
> Finally if I read the code correctly then if actually stride !=
> ioffsetdenom then the code produces incorrect results since Input.cc
> would compute lower as:
>
> lower = bbox.lower() - ioffset + patch->ioffset * stride / \
>                                  patch->ioffsetdenom
>
> (substituting for patch->iorigin) which only ever evaluates to the
> (presumably correct) bbox.lower() if stride == ioffsetdenom.
>
> Any thoughts?
>

With vertex centring, the offset denominator is always 1, and this is how
CarpetIOHDF5 was originally developed. There may be a bug in handling the
denominator, which is only relevant for cell centring. I would also expect
that most of the region descriptions calculated here are later ignored, in
favour of those calculated by CarpetRegrid2 and CarpetLib. The region
descriptions that survive should be the superregions as stored in the grid
structure.

Unfortunately, the attributes in CarpetIOHDF5's data format are not well
documented. Since these attributes determine whether simulation data can be
read back in months or years later, there is a strong pressure to remain
backward and bug compatible, and rather to introduce a new attribute than
to correct one that is set wrong; presumably, code that reads these
attributes (if any) will already deal with these bugs.

In this case, it seems that ioffset is the offset between the actual origin
of the grid and the location of the first data point, apparently measured
in units of the level's grid spacing. If so, the calculation would be
correct: the stride is the distance between two grid points of this level
(as measured in CarpetLib's internal integer coordinate system), and
ioffset is calculated as the region's lower boundary modulo this stride.
(This modulo calculation is slightly complicated by C++'s rules for the
module operator for negative numbers.) The expression ioffset /
ioffsetdenom, interpreted as real number, is thus non-negative, strictly
less than one, and the particular grid spacing used by CarpetLib's internal
integer coordinate system is irrelevant since it cancels.

The expression ioffset / ioffsetdenom should be interpreted as real number,
i.e. this fraction could be arbitrarily reduced, although CarpetIOHDF5
doesn't do that. When reading, the actual value of ioffsetdenom does not
matter, since only the fraction ioffset / ioffsetdenom should be evaluated;
however, since integer arithmetic is used, the newly used denominator needs
to be sufficiently large to hold the (reduced) fraction.

It seems that the assert in line 1399 is too strict if the fraction ioffset
/ iofsetdenom can reduced.

I do not understand what you want to say in the paragraph starting with
"Finally". bbox.lower is the location of the first grid point of a region;
the actual domain covered by the region would be 1/2 grid spacing larger if
cell centring is used. However, I don't know whether this is actually the
calculation performed here.

-erik

-- 
Erik Schnetter <schnetter at cct.lsu.edu>
http://www.perimeterinstitute.ca/personal/eschnetter/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.einsteintoolkit.org/pipermail/users/attachments/20140510/6e898c3f/attachment.html