From 0f8fc2db53346e3a533362fc0033013de9bba202 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Fri, 25 Feb 2011 20:42:52 -0500 Subject: CarpetIOHDF5: Store cell centering offset with grid function attributes --- Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc | 4 ++++ Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh | 2 ++ Carpet/CarpetIOHDF5/src/Input.cc | 36 +++++++++++++++++++++++++++------ Carpet/CarpetIOHDF5/src/Output.cc | 12 ++++++++++- Carpet/CarpetIOHDF5/src/OutputSlice.cc | 9 +++++++-- 5 files changed, 54 insertions(+), 9 deletions(-) diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc index d4abcccbd..1305a49d1 100644 --- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc +++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc @@ -299,6 +299,8 @@ int AddSliceAttributes(const cGH* const cctkGH, const vector& origin, const vector& delta, const vector& iorigin, + const vector& ioffset, + const vector& ioffsetdenom, const vector& bbox, const vector& nghostzones, hid_t& dataset) @@ -314,6 +316,8 @@ int AddSliceAttributes(const cGH* const cctkGH, error_count += WriteAttribute(dataset, "origin", &origin[0], origin.size()); error_count += WriteAttribute(dataset, "delta", &delta[0], delta.size()); error_count += WriteAttribute(dataset, "iorigin", &iorigin[0], iorigin.size()); + error_count += WriteAttribute(dataset, "ioffset", &ioffset[0], ioffset.size()); + error_count += WriteAttribute(dataset, "ioffsetdenom", &ioffsetdenom[0], ioffsetdenom.size()); // bbox and nghostzones are only used for grid functions and grid arrays if (bbox.size() > 0) { error_count += WriteAttribute(dataset, "cctk_bbox", &bbox[0], bbox.size()); diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh index e8f44e9b9..dfa6a276e 100644 --- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh +++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh @@ -140,6 +140,8 @@ namespace CarpetIOHDF5 const vector& origin, const vector& delta, const vector& iorigin, + const vector& ioffset, + const vector& ioffsetdenom, const vector& bbox, const vector& nghostzones, hid_t& dataset); diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc index 21c61cea8..cb0e57040 100644 --- a/Carpet/CarpetIOHDF5/src/Input.cc +++ b/Carpet/CarpetIOHDF5/src/Input.cc @@ -35,6 +35,8 @@ typedef struct { int component; int rank; ivect iorigin; + ivect ioffset; + ivect ioffsetdenom; vector shape; } patch_t; @@ -1058,9 +1060,6 @@ static herr_t BrowseDatasets (hid_t group, const char *objectname, void *arg) HDF5_ERROR (attr = H5Aopen_name (dataset, "group_timelevel")); HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &patch.timelevel)); HDF5_ERROR (H5Aclose (attr)); - // in old days (before February 2005) - // timelevels were stored as negative numbers - patch.timelevel = abs (patch.timelevel); HDF5_ERROR (attr = H5Aopen_name (dataset, "iorigin")); HDF5_ERROR (dataspace = H5Aget_space (attr)); assert (H5Sget_simple_extent_npoints (dataspace) == patch.rank); @@ -1068,6 +1067,25 @@ static herr_t BrowseDatasets (hid_t group, const char *objectname, void *arg) patch.iorigin = 0; HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &patch.iorigin[0])); HDF5_ERROR (H5Aclose (attr)); + patch.ioffset = 0; + attr = H5Aopen_name (dataset, "ioffset"); + // ioffset and ioffsetdenom may not be present; if so, use a default + if (attr >= 0) { + HDF5_ERROR (dataspace = H5Aget_space (attr)); + assert (H5Sget_simple_extent_npoints (dataspace) == patch.rank); + HDF5_ERROR (H5Sclose (dataspace)); + HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &patch.ioffset[0])); + HDF5_ERROR (H5Aclose (attr)); + } + patch.ioffsetdenom = 1; + attr = H5Aopen_name (dataset, "ioffsetdenom"); + if (attr >= 0) { + HDF5_ERROR (dataspace = H5Aget_space (attr)); + assert (H5Sget_simple_extent_npoints (dataspace) == patch.rank); + HDF5_ERROR (H5Sclose (dataspace)); + HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &patch.ioffsetdenom[0])); + HDF5_ERROR (H5Aclose (attr)); + } HDF5_ERROR (H5Dclose (dataset)); // try to obtain the map and component number from the patch's name @@ -1141,9 +1159,15 @@ static int ReadVar (const cGH* const cctkGH, } const hid_t datatype = CCTKtoHDF5_Datatype (cctkGH, group.vartype, 0); - const ivect stride = group.grouptype == CCTK_GF ? - maxspacereflevelfact/spacereflevelfact : 1; - ivect lower = patch->iorigin * stride; + centering const refcent = vhh.at(0)->refcent; + int const reffactdenom = + group.grouptype == CCTK_GF ? + (refcent == vertex_centered ? 1 : 2) : 1; + const ivect stride = + group.grouptype == CCTK_GF ? + reffactdenom * maxspacereflevelfact / spacereflevelfact : 1; + assert (all (stride % patch->ioffsetdenom == 0)); + ivect lower = patch->iorigin * stride + patch->ioffset * stride / patch->ioffsetdenom; ivect upper = lower + (shape - 1) * stride; // Traverse all local components on all maps diff --git a/Carpet/CarpetIOHDF5/src/Output.cc b/Carpet/CarpetIOHDF5/src/Output.cc index 3c44d87ac..1b5bf8eec 100644 --- a/Carpet/CarpetIOHDF5/src/Output.cc +++ b/Carpet/CarpetIOHDF5/src/Output.cc @@ -809,7 +809,17 @@ static int AddAttributes (const cGH *const cctkGH, const char *fullname, HDF5_ERROR (H5Aclose (attr)); } - ivect iorigin = bbox.lower() / bbox.stride(); + ivect const ioffsetdenom = bbox.stride(); + HDF5_ERROR (attr = H5Acreate (dataset, "ioffsetdenom", H5T_NATIVE_INT, + dataspace, H5P_DEFAULT)); + HDF5_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &ioffsetdenom[0])); + HDF5_ERROR (H5Aclose (attr)); + ivect const ioffset = ((bbox.lower() % bbox.stride()) + bbox.stride()) % bbox.stride(); + HDF5_ERROR (attr = H5Acreate (dataset, "ioffset", H5T_NATIVE_INT, + dataspace, H5P_DEFAULT)); + HDF5_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &ioffset[0])); + HDF5_ERROR (H5Aclose (attr)); + ivect const iorigin = (bbox.lower() - ioffset) / bbox.stride(); HDF5_ERROR (attr = H5Acreate (dataset, "iorigin", H5T_NATIVE_INT, dataspace, H5P_DEFAULT)); HDF5_ERROR (H5Awrite (attr, H5T_NATIVE_INT, &iorigin[0])); diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc index a87fe0a79..90e9539bc 100644 --- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc +++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc @@ -1381,6 +1381,8 @@ namespace CarpetIOHDF5 { slice_start, NULL, slice_count, NULL)); vector iorigin(rank, 0); + vector ioffset(rank, 0); + vector ioffsetdenom(rank, 1); vector delta(rank, 0), origin(rank, 0); vector bbox(2*rank, 0), nghostzones(rank, 0); for (int d = 0; d < outdim; d++) { @@ -1392,7 +1394,10 @@ namespace CarpetIOHDF5 { delta[d] = (coord_upper[dirs[d]] - coord_lower[dirs[d]]) / (gfext.upper()[dirs[d]] - gfext.lower()[dirs[d]]) * gfext.stride()[dirs[d]]; origin[d] += (org1[dirs[d]] - gfext.lower()[dirs[d]]) * delta[d]; - iorigin[d] /= gfext.stride()[dirs[d]]; + ioffsetdenom[d] = gfext.stride()[dirs[d]]; + ioffset[d] = (iorigin[d] % gfext.stride()[dirs[d]] + gfext.stride()[dirs[d]]) % gfext.stride()[dirs[d]]; + assert ((iorigin[d] - ioffset[d]) % gfext.stride()[dirs[d]] == 0); + iorigin[d] = (iorigin[d] - ioffset[d]) / gfext.stride()[dirs[d]]; } } // store cctk_bbox and cctk_nghostzones (for grid arrays only) @@ -1432,7 +1437,7 @@ namespace CarpetIOHDF5 { H5P_DEFAULT, gfdatas[n]->storage())); error_count += AddSliceAttributes (cctkGH, fullname, rl, ml, m, tl, origin, delta, - iorigin, bbox, nghostzones, dataset); + iorigin, ioffset, ioffsetdenom, bbox, nghostzones, dataset); HDF5_ERROR(H5Dclose (dataset)); free (fullname); -- cgit v1.2.3