From 2837c74e58aaa313073b0baaa24a6af491018d7a Mon Sep 17 00:00:00 2001 From: tradke Date: Thu, 26 Aug 2004 10:23:44 +0000 Subject: Don't get hyperslabs from coordinate grid variables. Get the bounding box information directly from coordinate tables. This should fix wrong bboxes in the output. git-svn-id: http://svn.cactuscode.org/arrangements/CactusIO/IOSDF/trunk@5 3af55ef0-e5e4-43b4-b020-ca5761ff09b2 --- interface.ccl | 16 ++++----- src/Write1D.c | 113 +++++++++++++++++++--------------------------------------- src/Write2D.c | 98 +++++++++++++++++++------------------------------- src/Write3D.c | 101 ++++++++++++++++++++------------------------------- src/ioSDFGH.h | 1 + 5 files changed, 120 insertions(+), 209 deletions(-) diff --git a/interface.ccl b/interface.ccl index 2c84453..7f0609b 100644 --- a/interface.ccl +++ b/interface.ccl @@ -5,15 +5,13 @@ implements: IOSDF inherits: IO CCTK_INT FUNCTION \ - Hyperslab_GetList (CCTK_POINTER_TO_CONST IN cctkGH, \ - CCTK_INT IN mapping_handle, \ - CCTK_INT IN num_arrays, \ - CCTK_INT ARRAY IN procs, \ - CCTK_INT ARRAY IN vindices, \ - CCTK_INT ARRAY IN timelevels, \ - CCTK_INT ARRAY IN hdatatypes, \ - CCTK_POINTER ARRAY IN hdata, \ - CCTK_INT ARRAY OUT retvals) + Hyperslab_Get (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN mapping_handle, \ + CCTK_INT IN proc, \ + CCTK_INT IN vindex, \ + CCTK_INT IN timelevel, \ + CCTK_INT IN hdatatype, \ + CCTK_POINTER IN hdata) CCTK_INT FUNCTION \ Hyperslab_GlobalMappingByIndex (CCTK_POINTER_TO_CONST IN cctkGH, \ diff --git a/src/Write1D.c b/src/Write1D.c index debf166..31073b4 100644 --- a/src/Write1D.c +++ b/src/Write1D.c @@ -76,22 +76,19 @@ static char *OpenFile (const cGH *GH, int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) { ioSDFGH *myGH; - int do_dir[4], coord_index[3]; + int do_dir[4]; int i, dir, myproc, gindex, have_coords, mapping; - int num_requested_hslabs, num_returned_hslabs; int *extent_int; cGroup gdata; char *fullname, *groupname, *filename; CCTK_INT coord_system_handle, coord_handles[3]; - double offset; - CCTK_REAL coord_lower[3]; + CCTK_REAL minext[3], delta[3]; CCTK_INT downsample[4]; CCTK_INT *origin, *direction; CCTK_INT hsize, extent; - CCTK_INT vindices[2]; - double *hdata[2]; + double *hdata; + double bbox[2]; const double dtime = GH->cctk_time; - const char *coordnames[] = {"x", "y", "z", "d"}; DECLARE_CCTK_PARAMETERS @@ -141,24 +138,6 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) coord_system_handle = Coord_GroupSystem (GH, groupname); free (groupname); - dir = gdata.dim < 3 ? gdata.dim : 3; - - have_coords = coord_system_handle >= 0 && - Util_TableGetIntArray (coord_system_handle, dir, - coord_handles, "COORDINATES") >= 0; - if (have_coords) - { - /* get the coordinate functions and coordinate physical minimum */ - for (i = 0; i < dir; i++) - { - coord_index[i] = -1; - coord_lower[i] = 0; - Util_TableGetInt (coord_handles[i], &coord_index[i], "GAINDEX"); - Util_TableGetReal (coord_handles[i], &coord_lower[i], "COMPMIN"); - have_coords &= coord_index[i] >= 0; - } - } - myproc = CCTK_MyProc (GH); origin = calloc (2*gdata.dim, sizeof (CCTK_INT)); @@ -171,6 +150,22 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) downsample[2] = out_downsample_z; downsample[3] = 1; + dir = gdata.dim < 3 ? gdata.dim : 3; + + have_coords = coord_system_handle >= 0 && + Util_TableGetIntArray (coord_system_handle, dir, + coord_handles, "COORDINATES") >= 0; + if (have_coords) + { + /* get the coordinates ranges */ + for (i = 0; i < dir; i++) + { + Util_TableGetReal (coord_handles[i], &minext[i], "COMPMIN"); + Util_TableGetReal (coord_handles[i], &delta[i], "DELTA"); + delta[i] *= downsample[i]; + } + } + /* get the variable's extents, compute the extent for 3D-diagonals as the minimum of grid points in each direction */ CCTK_GroupgshVI (GH, gdata.dim, extent_int, vindex); @@ -192,13 +187,8 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) /* now do the actual I/O looping over all directions */ for (dir = 0; dir < 4; dir++) { - if (hsize <= 0) - { - continue; - } - /* skip empty slices */ - if (! do_dir[dir]) + if (hsize <= 0 || ! do_dir[dir]) { continue; } @@ -210,10 +200,6 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) filename = OpenFile (GH, fullname, alias, &gdata, dir); } - /* get the number of hyperslabs to extract - (ie. whether to get a coordinate hyperslab too or not) */ - num_requested_hslabs = have_coords && dir < 3 ? 2 : 1; - /* set the direction vector */ for (i = 0; i < gdata.dim; i++) { @@ -233,9 +219,9 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) extent -= origin[dir]; /* correct extent in the case of staggered grids */ - if (CCTK_StaggerDirIndex(dir,gdata.stagtype)==1) + if (CCTK_StaggerDirIndex (dir, gdata.stagtype) == 1) { - --extent; + extent--; } } else /* origin for CCTK_ARRAYS is always (0, 0, 0) */ @@ -265,26 +251,17 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) continue; } - /* allocate hyperslab buffers on I/O processor */ - hdata[0] = hdata[1] = NULL; - if (myproc == 0) - { - hdata[0] = malloc (hsize * sizeof (double)); - hdata[1] = have_coords ? malloc (hsize * sizeof (double)) : NULL; - } + /* allocate hyperslab buffer on I/O processor */ + hdata = myproc == 0 ? malloc (hsize * sizeof (double)) : NULL; - /* get the hyperslabs */ - vindices[0] = vindex; - vindices[1] = coord_index[dir]; - num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs, - NULL, vindices, NULL, NULL, - (void **) hdata, NULL); + /* get the hyperslab */ + i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata); /* release the mapping structure */ Hyperslab_FreeMapping (mapping); /* And dump the data to file */ - if (num_returned_hslabs != num_requested_hslabs) + if (i != 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "IOSDF_Write1D: Failed to extract hyperslab for " @@ -292,35 +269,17 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) } else if (filename) { + bbox[0] = -1; bbox[1] = +1; if (have_coords) { - if (dir < 3) - { - /* get the staggering offset for the xyz coordinates */ - offset = 0.5 * GH->cctk_delta_space[dir] * - CCTK_StaggerDirIndex (dir, gdata.stagtype); - for (i = 0; i < hsize; i++) - { - hdata[1][i] += offset; - } - } - else - { - /* calculate the diagonal coordinates */ - offset = GH->cctk_delta_space[0] * sqrt (3); - for (i = 0; i < hsize; i++) - { - hdata[1][i] = coord_lower[0]*sqrt (3) + i*offset; - } - } + bbox[0] = minext[dir]; + bbox[1] = dir < 3 ? bbox[0] + delta[dir]*(hsize-1) : + (minext[0] + hsize*GH->cctk_delta_space[0]) * sqrt (3); } + gft_out_set_bbox (bbox, 1); extent_int[0] = hsize; - i = have_coords ? gft_out_full (filename, dtime, extent_int, - coordnames[dir], 1, hdata[1], hdata[0]): - gft_out_brief (filename, dtime, extent_int, - 1, hdata[0]); - if (i <= 0) + if (gft_out_brief (filename, dtime, extent_int, 1, hdata) <= 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Error writing to 1D IOSDF output file '%s'", filename); @@ -328,8 +287,8 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) } /* clean up */ - free (hdata[0]); - free (hdata[1]); + free (hdata); + } /* end of loop through all directions */ /* free allocated resources */ diff --git a/src/Write2D.c b/src/Write2D.c index 551e0c8..2906c79 100644 --- a/src/Write2D.c +++ b/src/Write2D.c @@ -43,7 +43,7 @@ static char **OpenFile (const cGH *GH, @enddesc @calls Hyperslab_GlobalMappingByIndex Hyperslab_FreeMapping - Hyperslab_GetList + Hyperslab_Get OpenFile WriteData @@ -75,21 +75,19 @@ static char **OpenFile (const cGH *GH, int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) { ioSDFGH *myGH; - int i, total_hsize, num_requested_hslabs, num_returned_hslabs; + int i, total_hsize; int dir, dir_i, dir_j, maxdir, myproc, gindex, have_coords; int mapping; cGroup gdata; - int coord_index[3]; CCTK_INT coord_system_handle, coord_handles[3]; char *fullname, *groupname; - int extent_int[3]; - double offset[2]; - CCTK_INT vindices[3], origin[3], extent[2], direction[6], downsample[2], - hsize[2]; - double *hdata[3]; + int full_extent[3], slice_extent[2]; + double bbox[2*2]; + CCTK_INT origin[3], extent[2], direction[6], downsample[2], hsize[2]; + CCTK_REAL minext[2], delta[2]; + double *hdata; char **filenames; const double dtime = GH->cctk_time; - const char *coordnames[] = {"x|y", "x|z", "y|z"}; DECLARE_CCTK_PARAMETERS @@ -124,17 +122,6 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) have_coords = coord_system_handle >= 0 && Util_TableGetIntArray (coord_system_handle, dir, coord_handles, "COORDINATES") >= 0; - if (have_coords) - { - /* get the coordinate functions and coordinate physical minimum */ - for (i = 0; i < dir; i++) - { - coord_index[i] = -1; - Util_TableGetInt (coord_handles[i], &coord_index[i], "GAINDEX"); - have_coords &= coord_index[i] >= 0; - } - } - num_requested_hslabs = have_coords ? 3 : 1; /* processor 0 opens the files on the first trip through */ filenames = NULL; @@ -145,12 +132,12 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) } /* get the extents of the variable */ - CCTK_GroupgshVI (GH, gdata.dim, extent_int, vindex); + CCTK_GroupgshVI (GH, gdata.dim, full_extent, vindex); /* get the total number of grid points to check for zero-sized variables */ for (dir = 0, hsize[0] = 1; dir < gdata.dim; dir++) { - hsize[0] *= extent_int[dir]; + hsize[0] *= full_extent[dir]; } /* now do the actual I/O looping over all directions */ @@ -179,8 +166,8 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) } /* set the extent vector */ - extent[0] = extent_int[dir_i]; - extent[1] = extent_int[dir_j]; + extent[0] = full_extent[dir_i]; + extent[1] = full_extent[dir_j]; /* set the origin using the slice center from IOUtil */ memset (origin, 0, sizeof (origin)); @@ -216,22 +203,30 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) continue; } - hdata[0] = hdata[1] = hdata[2]; - if (myproc == 0) + /* get the bounding box information */ + bbox[0] = bbox[2] = -1; + bbox[1] = bbox[3] = +1; + if (have_coords) { - /* allocate hyperslab buffers */ - hdata[0] = malloc (total_hsize * sizeof (double)); - hdata[1] = have_coords ? malloc (2 * total_hsize * sizeof(double)) : NULL; - hdata[2] = hdata[1] + total_hsize; + Util_TableGetReal (coord_handles[dir_i], &minext[0], "COMPMIN"); + Util_TableGetReal (coord_handles[dir_i], &delta[0], "DELTA"); + Util_TableGetReal (coord_handles[dir_j], &minext[1], "COMPMIN"); + Util_TableGetReal (coord_handles[dir_j], &delta[1], "DELTA"); + + delta[0] *= downsample[0]; + delta[1] *= downsample[1]; + bbox[0] = minext[0]; + bbox[1] = bbox[0] + delta[0]*(hsize[0]-1); + bbox[2] = minext[1]; + bbox[3] = bbox[2] + delta[1]*(hsize[1]-1); } + gft_out_set_bbox (bbox, 2); - /* get the hyperslabs */ - vindices[0] = vindex; - vindices[1] = coord_index[dir_i]; - vindices[2] = coord_index[dir_j]; - num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs, - NULL, vindices, NULL, NULL, - (void **) hdata, NULL); + /* allocate hyperslab buffer */ + hdata = myproc == 0 ? malloc (total_hsize * sizeof (double)) : NULL; + + /* get the hyperslab */ + i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata); /* release the mapping structure */ Hyperslab_FreeMapping (mapping); @@ -239,30 +234,12 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) /* and dump the data to file */ if (filenames) { - if (num_returned_hslabs == num_requested_hslabs) + if (i == 0) { - if (have_coords) - { - /* get the staggering offset for the coordinates */ - offset[0] = 0.5 * GH->cctk_delta_space[dir_i] * - CCTK_StaggerDirIndex (dir_i, gdata.stagtype); - offset[1] = 0.5 * GH->cctk_delta_space[dir_j] * - CCTK_StaggerDirIndex (dir_j, gdata.stagtype); - for (i = 0; i < total_hsize; i++) - { - hdata[1][i] += offset[0]; - hdata[2][i] += offset[1]; - } - } - - extent[0] = hsize[0]; - extent[1] = hsize[1]; + slice_extent[0] = hsize[0]; + slice_extent[1] = hsize[1]; - i = have_coords ? gft_out_full (filenames[dir], dtime, extent_int, - coordnames[dir], 2, hdata[1], hdata[0]): - gft_out_brief (filenames[dir], dtime, extent_int, - 2, hdata[0]); - if (i <= 0) + if (gft_out_brief (filenames[dir], dtime, slice_extent, 2, hdata) <= 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Error writing to 2D IOSDF output file '%s'", @@ -277,8 +254,7 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) } /* clean up */ - free (hdata[0]); - free (hdata[1]); + free (hdata); } /* end of outputting the data by processor 0 */ diff --git a/src/Write3D.c b/src/Write3D.c index f07ef55..314146b 100644 --- a/src/Write3D.c +++ b/src/Write3D.c @@ -38,7 +38,7 @@ static char *OpenFile (const cGH *GH, const char *fullname, const char *alias); @enddesc @calls Hyperslab_GlobalMappingByIndex Hyperslab_FreeMapping - Hyperslab_GetList + Hyperslab_Get OpenFile @var GH @@ -69,16 +69,16 @@ static char *OpenFile (const cGH *GH, const char *fullname, const char *alias); int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias) { int i, total_hsize; - int myproc, gindex, have_coords; - int num_requested_hslabs, num_returned_hslabs; + int myproc, gindex; cGroup gdata; CCTK_INT coord_system_handle, coord_handles[3]; char *fullname, *groupname, *filename; - double *hdata[4]; + double *hdata; int extent_int[3]; - double offset[3]; int mapping; - CCTK_INT vindices[4], extent[3], downsample[3], hsize[3]; + CCTK_INT extent[3], downsample[3], hsize[3]; + double bbox_dbl[2*3]; + CCTK_REAL bbox[2*3], delta[3]; const double dtime = GH->cctk_time; const CCTK_INT origin[] = {0, 0, 0}, direction[] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; @@ -100,32 +100,10 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias) } /* get the coordinate system associated with this grid variable */ - vindices[0] = vindex; groupname = CCTK_GroupName (gindex); coord_system_handle = Coord_GroupSystem (GH, groupname); free (groupname); - have_coords = coord_system_handle >= 0 && - Util_TableGetIntArray (coord_system_handle, 3, - coord_handles, "COORDINATES") >= 0; - if (have_coords) - { - /* get the coordinate functions and coordinate physical minimum */ - for (i = 1; i <= 3; i++) - { - vindices[i] = -1; - Util_TableGetInt (coord_handles[i-1], &vindices[i], "GAINDEX"); - have_coords &= vindices[i] >= 0; - } - } - num_requested_hslabs = have_coords ? 4 : 1; - - /* What processor are we on? */ - myproc = CCTK_MyProc (GH); - - /* Open the file on processor 0 */ - filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL; - /* get the total number of grid points to check for zero-sized variables */ /* set the extent vector (copy from 'int' to 'CCTK_INT') */ CCTK_GroupgshVI (GH, 3, extent_int, vindex); @@ -144,6 +122,12 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias) downsample[1] = out_downsample_y; downsample[2] = out_downsample_z; + /* What processor are we on? */ + myproc = CCTK_MyProc (GH); + + /* Open the file on processor 0 */ + filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL; + /* get the hyperslab mapping */ mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 3, direction, origin, extent, @@ -174,46 +158,40 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias) return (-1); } - hdata[0] = hdata[1] = hdata[2] = hdata[3]; - if (myproc == 0) + /* get the bounding box information */ + for (i = 0; i < 3; i++) { - /* allocate hyperslab buffers */ - hdata[0] = malloc (total_hsize * sizeof (double)); - hdata[1] = have_coords ? malloc (3 * total_hsize * sizeof (double)) : NULL; - hdata[2] = hdata[1] + 1*total_hsize; - hdata[3] = hdata[1] + 2*total_hsize; + bbox_dbl[2*i + 0] = -1; + bbox_dbl[2*i + 1] = +1; } + if (coord_system_handle >= 0 && + Util_TableGetIntArray (coord_system_handle, 3, coord_handles, + "COORDINATES") >= 0) + { + for (i = 0; i < 3; i++) + { + Util_TableGetReal (coord_handles[i], &bbox[2*i + 0], "COMPMIN"); + Util_TableGetReal (coord_handles[i], &delta[i], "DELTA"); + + delta[i] *= downsample[i]; + bbox_dbl[2*i + 0] = bbox[2*i + 0]; + bbox_dbl[2*i + 1] = bbox_dbl[2*i + 0] + delta[i]*(hsize[i]-1); + } + } + gft_out_set_bbox (bbox_dbl, 3); + + /* allocate hyperslab buffer */ + hdata = myproc == 0 ? malloc (total_hsize * sizeof (double)) : NULL; - /* get the hyperslabs */ - num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs, - NULL, vindices, NULL, NULL, - (void **) hdata, NULL); + /* get the hyperslab */ + i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata); /* And dump the data to file */ if (myproc == 0) { - if (num_returned_hslabs == num_requested_hslabs) + if (i == 0) { - if (have_coords) - { - /* get the staggering offset for the coordinates */ - for (i = 0; i < 3; i++) - { - offset[i] = 0.5 * GH->cctk_delta_space[i] * - CCTK_StaggerDirIndex (i, gdata.stagtype); - } - for (i = 0; i < total_hsize; i++) - { - hdata[1][i] += offset[0]; - hdata[2][i] += offset[1]; - hdata[3][i] += offset[2]; - } - } - - i = have_coords ? gft_out_full (filename, dtime, extent_int, "x|y|z", 3, - hdata[1], hdata[0]) : - gft_out_brief (filename, dtime, extent_int, 3,hdata[0]); - if (i <= 0) + if (gft_out_brief (filename, dtime, extent_int, 3, hdata) <= 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Error writing to 3D IOSDF output file '%s'", filename); @@ -227,8 +205,7 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias) } /* clean up */ - free (hdata[0]); - free (hdata[1]); + free (hdata); } /* end of outputting the data by processor 0 */ free (fullname); diff --git a/src/ioSDFGH.h b/src/ioSDFGH.h index e31867b..3706d07 100644 --- a/src/ioSDFGH.h +++ b/src/ioSDFGH.h @@ -68,6 +68,7 @@ int IOSDF_Output3DVarAs (const cGH *GH, const char *var, const char *alias); /* other function prototypes */ void IOSDF_Startup (void); +void IOSDF_Terminate (void); int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias); int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias); -- cgit v1.2.3