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

[Roland Haas]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hello Erik, all,

> 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.
These are the bboxes that CarpetIOHDF5 constructs to decide where in
the grid to put data read from the file. This part of the code does
not know or use the superregions or the grid structure dataset in the
metadata group.

> 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.)
(Unrelated): should there ever be negative numbers present in this
expression? Can bbox.lower() ever be negative (for an actual extend)?

> 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.
Do you mean numerator, ie the term "patch->ioffset * stride" in line 1400?

> It seems that the assert in line 1399 is too strict if the fraction
> ioffset / iofsetdenom can reduced.
That it true. The cases where this happens would seem strange though.
I suspect it would only really become interesting if we were to allow
say reflevel 3 of the files to be read into reflevel 0 of the new
simulation (say in the file reader).

> 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.
I think I was just confused by exactly how things are used and how
they have to be interpreted. Let me see if I understand this all
correctly:

* iorigin, ioffset and ioffsetdenom together describe the coordinate
of the first point in the patch
* the coordinate of the first point is given by x0 = CCTK_ORIGIN(0) +
(iorigin * stride + ioffset * stride / ioffsetdenom) *
CCTK_DELTA_SPACE(0) / stride
* iorigin contains the number of whole cells (for the level) from the
origin, ioffset/ioffsetdenom contains the fraction of a cell that
needs to be added due to cell centering
* only the ratio ioffset/ioffsetdenom has any meaning, neither one
separately is a meaningful number

So, if I wanted to use (for the sake of argument) 3*stride as
ioffsetdenom then the values of the attributes in the file would
change as such:

iorigin' = iorigin (no change)
ioffset' = 3*ioffset
ioffsetdenom' = 3*ioffsetdenom

and I would change only the code in Output.cc to read:
 849   ivect const ioffsetdenom = 3*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() * 3;
 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/3) / bbox.stride();
ie change lines 849, 855 and 860 to contains factors of 3

So the rule would be that reader routines may only assume that
ioffset/ioffsetdenom makes sense and writer thorns are free to write
ioffset and ioffsetdenom in any way as long as the (reduced) fraction
is identical (which I admit the assert violates)?

When first reading the code I had assumed that the reader would also
make the assumption that ioffsetdenom == stride but that does not seem
to be the case.

Yours,
Roland

- -- 
My email is as private as my paper mail. I therefore support encrypting
and signing email messages. Get my PGP key from http://keys.gnupg.net.
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1
Comment: Using GnuPG with Icedove - http://www.enigmail.net/

iEYEARECAAYFAlNu/JoACgkQTiFSTN7SboWxXQCggn6hRtcQRQ/acWaUvmsUR6tK
kwoAoJ5zPn9bOQP9lxgMk9XuhwN4WXuM
=6e6J
-----END PGP SIGNATURE-----