diff options
Diffstat (limited to 'src/Write2D.c')
| -rw-r--r-- | src/Write2D.c | 371 |
1 files changed, 371 insertions, 0 deletions
diff --git a/src/Write2D.c b/src/Write2D.c new file mode 100644 index 0000000..551e0c8 --- /dev/null +++ b/src/Write2D.c @@ -0,0 +1,371 @@ +/*@@ + @file Write2D.c + @date Sat 12 June 2004 + @author Thomas Radke + @desc + Output two-dimensional slices in SDF file format. + @enddesc + @version $Id$ + @@*/ + + +#include <stdlib.h> +#include <string.h> + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "util_Table.h" +#include "CactusBase/IOUtil/src/ioGH.h" +#include "CactusBase/IOUtil/src/ioutil_AdvertisedFiles.h" +#include "ioSDFGH.h" + +/* the rcs ID and its dummy function to use it */ +static const char *rcsid = "$Id$"; +CCTK_FILEVERSION(CactusIO_IOSDF_Write2D_c) + + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +static char **OpenFile (const cGH *GH, + const char *fullname, + const char *alias, + int dim, + int maxdir); + + +/*@@ + @routine IOSDF_Write2D + @date Sat 12 June 2004 + @author Thomas Radke + @desc + Writes the 2D slices of a variable into separate SDF files. + @enddesc + @calls Hyperslab_GlobalMappingByIndex + Hyperslab_FreeMapping + Hyperslab_GetList + OpenFile + WriteData + + @var GH + @vdesc Pointer to CCTK GH + @vtype const cGH * + @vio in + @endvar + @var vindex + @vdesc index of variable to output + @vtype int + @vio in + @endvar + @var alias + @vdesc alias name of variable to output + @vtype const char * + @vio in + @endvar + + @returntype int + @returndesc + 0 for success, or<BR> + -1 if variable has no storage assigned<BR> + -2 if output file couldn't be opened<BR> + -3 if hyperslab for coordinates and/or variable couldn't be + extracted + @endreturndesc +@@*/ +int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias) +{ + ioSDFGH *myGH; + int i, total_hsize, num_requested_hslabs, num_returned_hslabs; + 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]; + char **filenames; + const double dtime = GH->cctk_time; + const char *coordnames[] = {"x|y", "x|z", "y|z"}; + DECLARE_CCTK_PARAMETERS + + + /* get the variable name and group information */ + fullname = CCTK_FullName (vindex); + gindex = CCTK_GroupIndexFromVarI (vindex); + CCTK_GroupData (gindex, &gdata); + + /* check if variable has storage assigned */ + if (! CCTK_QueryGroupStorageI (GH, gindex)) + { + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "No IOSDF 2D output for '%s' (no storage)", fullname); + free (fullname); + return (-1); + } + + /* get the handle for IOSDF extensions */ + myGH = CCTK_GHExtension (GH, "IOSDF"); + + /* get the number of slices to output */ + /* in general: maxdir = gdata.dim * (gdata.dim - 1) / 2; */ + maxdir = gdata.dim == 2 ? 1 : 3; + + /* get the coordinate system associated with this grid variable */ + groupname = CCTK_GroupName (gindex); + 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; + 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; + myproc = CCTK_MyProc (GH); + if (myproc == 0) + { + filenames = OpenFile (GH, fullname, alias, gdata.dim, maxdir); + } + + /* get the extents of the variable */ + CCTK_GroupgshVI (GH, gdata.dim, extent_int, 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]; + } + + /* now do the actual I/O looping over all directions */ + for (dir = 0; dir < maxdir; dir++) + { + if (hsize[0] <= 0) + { + continue; + } + + /* get the directions to span the hyperslab */ + if (dir == 0) + { + dir_i = 0; dir_j = 1; /* xy */ + downsample[0] = out_downsample_x; downsample[1] = out_downsample_y; + } + else if (dir == 1) + { + dir_i = 0; dir_j = 2; /* xz */ + downsample[0] = out_downsample_x; downsample[1] = out_downsample_z; + } + else + { + dir_i = 1; dir_j = 2; /* yz */ + downsample[0] = out_downsample_y; downsample[1] = out_downsample_z; + } + + /* set the extent vector */ + extent[0] = extent_int[dir_i]; + extent[1] = extent_int[dir_j]; + + /* set the origin using the slice center from IOUtil */ + memset (origin, 0, sizeof (origin)); + if (have_coords) + { + origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim-1][dir]; + } + + /* set the direction vector */ + memset (direction, 0, sizeof (direction)); + direction[dir_i] = direction[gdata.dim + dir_j] = 1; + + mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2, + direction, origin, extent, + downsample, + -1, /* table handle */ + NULL /* conversion fn */, + hsize); + if (mapping < 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write2D: Failed to define hyperslab mapping for " + "variable '%s'", fullname); + continue; + } + total_hsize = hsize[0] * hsize[1]; + if (total_hsize <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write2D: selected hyperslab has zero size for " + "variable '%s' direction %d", fullname, dir); + Hyperslab_FreeMapping (mapping); + continue; + } + + hdata[0] = hdata[1] = hdata[2]; + if (myproc == 0) + { + /* 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; + } + + /* 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); + + /* release the mapping structure */ + Hyperslab_FreeMapping (mapping); + + /* and dump the data to file */ + if (filenames) + { + if (num_returned_hslabs == num_requested_hslabs) + { + 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]; + + 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) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error writing to 2D IOSDF output file '%s'", + filenames[dir]); + } + } + else + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write2D: Failed to extract hyperslab for " + "variable '%s'", fullname); + } + + /* clean up */ + free (hdata[0]); + free (hdata[1]); + + } /* end of outputting the data by processor 0 */ + + } /* end of looping through xyz directions */ + + free (fullname); + + return (0); +} + + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +/*@@ + @routine OpenFile + @date Sat 12 June 2004 + @author Thomas Radke + @desc + Opens a set of SDF files for a given alias name. + If this is the first time through, it will advertise them + to IOUtil. + @enddesc + + @returntype char ** + @returndesc + the set of full filenames + @endreturndesc + @@*/ +static char **OpenFile (const cGH *GH, + const char *fullname, + const char *alias, + int dim, + int maxdir) +{ + int dir; + char **filenames; + ioSDFGH *myGH; + ioAdvertisedFileDesc advertised_file; + char slicename[30]; + const char *extensions[] = {"xy", "xz", "yz"}; + DECLARE_CCTK_PARAMETERS + + + /* get handle for IOSDF GH extensions */ + myGH = CCTK_GHExtension (GH, "IOSDF"); + + /* see if we are the first time through */ + filenames = GetNamedData (myGH->fileList_2D, alias); + if (filenames) + { + return (filenames); + } + + filenames = malloc (maxdir * sizeof (char **)); + + /* open/create files for each slice */ + for (dir = 0; dir < maxdir; dir++) + { + filenames[dir] = malloc (strlen (myGH->out2D_dir) + strlen (alias) + + sizeof (slicename) + 20); + + if (dim == 2) + { + strcpy (slicename, "2D"); + } + else + { + /* give the slice origin as range [1 .. n] */ + sprintf (slicename, "%s_%d", extensions[dir], + myGH->sp2xyz[dim-1][dir]); + } + + sprintf (filenames[dir], "%s%s_%s.sdf", myGH->out2D_dir, alias, slicename); + + /* advertise the file for downloading and write file info */ + advertised_file.slice = slicename; + advertised_file.thorn = CCTK_THORNSTRING; + advertised_file.varname = fullname; + advertised_file.description = "Two-dimensional slice plots"; + advertised_file.mimetype = "application/sdf"; + + IOUtil_AdvertiseFile (GH, filenames[dir], &advertised_file); + } + + /* store file desriptors in database */ + StoreNamedData (&myGH->fileList_2D, alias, filenames); + + return (filenames); +} |