diff options
Diffstat (limited to 'src/Write1D.c')
| -rw-r--r-- | src/Write1D.c | 113 |
1 files changed, 36 insertions, 77 deletions
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 */ |