diff options
Diffstat (limited to 'src/Write1D.c')
| -rw-r--r-- | src/Write1D.c | 436 |
1 files changed, 436 insertions, 0 deletions
diff --git a/src/Write1D.c b/src/Write1D.c new file mode 100644 index 0000000..debf166 --- /dev/null +++ b/src/Write1D.c @@ -0,0 +1,436 @@ +/*@@ + @file Write1D.c + @date Sat 12 June 2004 + @author Thomas Radke + @desc + Output one-dimensional lines in SDF file format. + @enddesc + @version $Id$ +@@*/ + +#include <math.h> /* sqrt(3) */ +#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 = "$Header$"; +CCTK_FILEVERSION(CactusIO_IOSDF_Write1D_c) + + +/******************************************************************** + ******************** Internal Routines ************************ + ********************************************************************/ +static char *OpenFile (const cGH *GH, + const char *fullname, + const char *alias, + const cGroup *gdata, + int dir); + + +/*@@ + @routine IOSDF_Write1D + @date Sat 12 June 2004 + @author Thomas Radke + @desc + This routine does 1D line output along the orthogonals + and the diagonal (in case of a cubed grid). + <p> + It writes to SDF files suitable for gnuplot and xgraph. + A header telling the physical time prefixes the output data. + @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 global index of variable to output + @vtype int + @vio in + @endvar + @var alias + @vdesc alias name (used for creating the output filename) + @vtype const char * + @vio in + @endvar + + @returntype int + @returndesc + 0 for success, or<BR> + -1 if variable has no storage assigned + @endreturndesc +@@*/ +int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias) +{ + ioSDFGH *myGH; + int do_dir[4], coord_index[3]; + 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_INT downsample[4]; + CCTK_INT *origin, *direction; + CCTK_INT hsize, extent; + CCTK_INT vindices[2]; + double *hdata[2]; + const double dtime = GH->cctk_time; + const char *coordnames[] = {"x", "y", "z", "d"}; + DECLARE_CCTK_PARAMETERS + + + /* get the variable's group index and its full name */ + gindex = CCTK_GroupIndexFromVarI (vindex); + fullname = CCTK_FullName (vindex); + + /* check if variable has storage assigned */ + if (! CCTK_QueryGroupStorageI (GH, gindex)) + { + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write1D: No IOSDF_1D output for '%s' (no storage)", + fullname); + free (fullname); + return (-1); + } + + /* get the handle for IOSDF extensions */ + myGH = CCTK_GHExtension (GH, "IOSDF"); + + /* get the variable's group information */ + CCTK_GroupData (gindex, &gdata); + + /* see what slices should be output */ + do_dir[0] = out1D_x && gdata.dim >= 1; + do_dir[1] = out1D_y && gdata.dim >= 2; + do_dir[2] = out1D_z && gdata.dim >= 3; + /* diagonal slice is done only if variable is non-staggered and 3D */ + do_dir[3] = out1D_d && gdata.dim == 3 && gdata.stagtype == 0; + if (out1D_d && ! do_dir[3] && myGH->out1D_last[vindex] < 0) + { + CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write1D: No IOSDF_1D diagonal output for '%s' " + "(only implemented for non-staggered 3D variables)", + fullname); + } + + /* return if nothing to do */ + if (! (do_dir[0] || do_dir[1] || do_dir[2] || do_dir[3])) + { + free (fullname); + return (0); + } + + /* 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; + 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)); + direction = origin + gdata.dim; + extent_int = malloc ((gdata.dim + 1) * sizeof (int)); + + /* set downsampling vector from I/O parameters */ + downsample[0] = out_downsample_x; + downsample[1] = out_downsample_y; + downsample[2] = out_downsample_z; + downsample[3] = 1; + + /* 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); + if (gdata.dim == 3) + { + extent_int[3] = extent_int[0] < extent_int[1] ? + extent_int[0] : extent_int[1]; + if (extent_int[2] < extent_int[3]) + { + extent_int[3] = extent_int[2]; + } + } + /* get the total number of grid points to check for zero-sized variables */ + for (dir = 0, hsize = 1; dir < gdata.dim; dir++) + { + hsize *= extent_int[dir]; + } + + /* 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]) + { + continue; + } + + /* processor 0 opens the files with the appropriate name */ + filename = NULL; + if (myproc == 0) + { + 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++) + { + direction[i] = (dir == i || dir == 3) ? 1 : 0; + } + + /* set the extent */ + extent = extent_int[dir]; + + /* set the origin of the line */ + if (gdata.grouptype == CCTK_GF && dir < 3) + { + for (i = 0; i < gdata.dim; i++) + { + origin[i] = myGH->spxyz[gdata.dim-1][dir][i]; + } + extent -= origin[dir]; + + /* correct extent in the case of staggered grids */ + if (CCTK_StaggerDirIndex(dir,gdata.stagtype)==1) + { + --extent; + } + } + else /* origin for CCTK_ARRAYS is always (0, 0, 0) */ + { + memset (origin, 0, gdata.dim * sizeof (CCTK_INT)); + } + + mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 1, + direction, origin, &extent, + &downsample[dir], + -1, /* table handle */ + NULL /* conversion fn */, + &hsize); + if (mapping < 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write1D: Failed to define hyperslab mapping for " + "variable '%s'", fullname); + continue; + } + if (hsize <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write1D: selected hyperslab has zero size for " + "variable '%s' direction %d", fullname, dir); + Hyperslab_FreeMapping (mapping); + 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; + } + + /* 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); + + /* release the mapping structure */ + Hyperslab_FreeMapping (mapping); + + /* And dump the data to file */ + if (num_returned_hslabs != num_requested_hslabs) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "IOSDF_Write1D: Failed to extract hyperslab for " + "variable '%s'", fullname); + } + else if (filename) + { + 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; + } + } + } + + 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) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error writing to 1D IOSDF output file '%s'", filename); + } + } + + /* clean up */ + free (hdata[0]); + free (hdata[1]); + } /* end of loop through all directions */ + + /* free allocated resources */ + free (origin); + free (fullname); + free (extent_int); + + 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 + @@*/ +static char *OpenFile (const cGH *GH, + const char *fullname, + const char *alias, + const cGroup *gdata, + int dir) +{ + ioSDFGH *myGH; + int upper, lower; + char *filename, *nameddata; + char slicename[40]; + ioAdvertisedFileDesc advertised_file; + DECLARE_CCTK_PARAMETERS + + + /* get handle for and IOSDF GH extensions */ + myGH = CCTK_GHExtension (GH, "IOSDF"); + + nameddata = malloc (strlen (alias) + 3); + sprintf (nameddata, "%s_%d", alias, dir); + filename = GetNamedData (myGH->fileList_1D, nameddata); + if (filename) + { + free (nameddata); + return (filename); + } + + /* get the indices into spxyz[] */ + lower = (dir + 1) % 3; + upper = (dir + 2) % 3; + if (upper < lower) + { + upper = lower; + lower = 0; + } + + if (dir < 3) + { + if (gdata->dim == 1) + { + strcpy (slicename, "1D"); + } + else if (gdata->dim == 2) + { + /* give the slice origin as range [1 .. n] */ + sprintf (slicename, "%c_%d", 'x' + dir, + gdata->grouptype == CCTK_GF ? + myGH->spxyz[gdata->dim-1][dir][lower] : 0); + } + else + { + /* give the slice origin as range [1 .. n] */ + sprintf (slicename, "%c_%d_%d", 'x' + dir, + gdata->grouptype == CCTK_GF ? + myGH->spxyz[gdata->dim-1][dir][lower] : 0, + gdata->grouptype == CCTK_GF ? + myGH->spxyz[gdata->dim-1][dir][upper] : 0); + } + } + else + { + sprintf (slicename, "%dD_diagonal", gdata->dim); + } + + filename = malloc (strlen (myGH->out1D_dir) + strlen (alias) + + sizeof (slicename) + 6); + sprintf (filename, "%s%s_%s.sdf", myGH->out1D_dir, alias, slicename); + StoreNamedData (&myGH->fileList_1D, nameddata, filename); + free (nameddata); + + /* advertise the output file */ + advertised_file.slice = slicename; + advertised_file.thorn = CCTK_THORNSTRING; + advertised_file.varname = fullname; + advertised_file.description = "One-dimensional line plots"; + advertised_file.mimetype = "application/sdf"; + + IOUtil_AdvertiseFile (GH, filename, &advertised_file); + + return (filename); +} |