/*@@ @file Write2D.c @date Thu May 11 2000 @author Thomas Radke @desc Output two-dimensional slices in ASCII gnuplot format. @enddesc @version $Id$ @@*/ #include #include #include "cctk.h" #include "cctk_Parameters.h" #include "Hyperslab.h" #include "CactusBase/IOUtil/src/ioutil_AdvertisedFiles.h" #include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h" #include "ioASCIIGH.h" /* the rcs ID and its dummy function to use it */ static const char *rcsid = "$Id$"; CCTK_FILEVERSION(CactusBase_IOASCII_Write2D_c) /* enable debug output */ /*#define IOASCII_DEBUG 1*/ /* macro to output a time slice as typed data */ #define OUTPUT_TYPED_DATA(coord_i, coord_j, stagger_offset_i,stagger_offset_j,\ cctk_type, c_type, data, hsize, is_cmplx_type, \ have_coords, fmt_string, file) \ { \ int i, j; \ const cctk_type *typed_data = (const cctk_type *) (data); \ \ \ /* output coordinates if available otherwise just the indices */ \ if (have_coords) \ { \ for (j = 0; j < hsize[1]; j++) \ { \ for (i = 0; i < hsize[0]; i++) \ { \ fprintf (file, fmt_string, \ (double) (*coord_i++ + stagger_offset_i), \ (double) (*coord_j++ + stagger_offset_j), \ (c_type) *typed_data++); \ if (is_cmplx_type) \ { \ fprintf (file, "\t\t"); \ fprintf (file, out_real_format, (double) *typed_data++); \ } \ fprintf (file, "\n"); \ } \ fprintf (file, "\n"); \ } \ coord_i -= hsize[0] * hsize[1]; \ coord_j -= hsize[0] * hsize[1]; \ } \ else \ { \ for (j = 0; j < hsize[1]; j++) \ { \ for (i = 0; i < hsize[0]; i++) \ { \ fprintf (file, fmt_string, \ (double) i, \ (double) j, \ (c_type) *typed_data++); \ if (is_cmplx_type) \ { \ fprintf (file, "\t\t"); \ fprintf (file, out_real_format, (double) *typed_data++); \ } \ fprintf (file, "\n"); \ } \ fprintf (file, "\n"); \ } \ } \ } /*@@ @routine IOASCII_Write2D @date Thu May 11 2000 @author Thomas Radke @desc Writes the 2D slices of a variable into separate ASCII files. @enddesc @calls IOUtil_RestartFromRecovery IOUtil_AdvertiseFile Hyperslab_GetHyperslab @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 @@*/ void IOASCII_Write2D (const cGH *GH, int vindex, const char *alias) { DECLARE_CCTK_PARAMETERS int myproc; asciiioGH *myGH; char header_fmt_string[30]; /* header format string */ char zlabel_fmt_string[30]; /* z-label format string */ char out_real_format[30]; /* data format string for a real type */ char time_fmt_string[30]; /* time format string */ char data_fmt_string_int[30]; /* data format string for int types */ char data_fmt_string_real[30]; /* data format string for real types */ int dir, maxdir; int groupindex; int have_coords; cGroup groupinfo; FILE **fdset_2D; /* array of output file pointers */ int coord_index[3]; /* variable indices for xyz coordinates */ CCTK_REAL coord_lower[3]; /* coordinates' minima */ char coord_system[20]; /* name of the coordinate system */ int origin[3]; /* the slice origin */ CCTK_REAL dummy; char *filename; char *fullname; char slicename[20]; ioAdvertisedFileDesc advertised_file; char buffer[128]; static char *extensions[] = {"xy", "xz", "yz"}; /* to make the compiler happy */ fdset_2D = NULL; /* get the variable group information */ groupindex = CCTK_GroupIndexFromVarI (vindex); CCTK_GroupData (groupindex, &groupinfo); /* check if variable has storage assigned */ if (! CCTK_QueryGroupStorageI (GH, groupindex)) { fullname = CCTK_FullName (vindex); CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, "No IOASCII 2D output for '%s' (no storage)", fullname); free (fullname); return; } /* Get the handle for IOASCII extensions */ myGH = (asciiioGH *) CCTK_GHExtension (GH, "IOASCII"); /* get the number of slices to output */ /* in general: maxdir = groupinfo.dim * (groupinfo.dim - 1) / 2; */ maxdir = groupinfo.dim == 2 ? 1 : 3; /* set header and data format strings */ sprintf (header_fmt_string, "\n\n#Time = %%%s\n", out_format); sprintf (zlabel_fmt_string, " (%%c = %%%s),", out_format); sprintf (out_real_format, "%%%s", out_format); dir = 0; if (CCTK_Equals (out2D_style, "gnuplot f(t,x,y)")) { sprintf (time_fmt_string, "%%%s\t\t", out_format); sprintf (data_fmt_string_int, time_fmt_string, (double) GH->cctk_time); sprintf (data_fmt_string_real, time_fmt_string, (double) GH->cctk_time); dir = strlen (data_fmt_string_int); } sprintf (data_fmt_string_int + dir, "%%%s\t\t%%%s\t\t%%d", out_format, out_format); sprintf (data_fmt_string_real + dir, "%%%s\t\t%%%s\t\t%%%s", out_format, out_format, out_format); /* get the coordinate indices if we output a grid function */ if (groupinfo.grouptype == CCTK_GF) { sprintf (coord_system, "cart%dd", groupinfo.dim); have_coords = 1; for (dir = 0; dir < groupinfo.dim && dir < 3; dir++) { CCTK_CoordRange (GH, &coord_lower[dir], &dummy, dir+1, NULL,coord_system); coord_index[dir] = CCTK_CoordIndex (dir + 1, NULL, coord_system); have_coords &= coord_index[dir] >= 0; } if (! have_coords) { CCTK_VWarn (8, __LINE__, __FILE__, CCTK_THORNSTRING, "IOASCII_Write2D: No coordinate ranges found for '%s'", coord_system); } } else { /* CCTK_ARRAY variables never have coordinates associated */ have_coords = 0; } /* What processor are we on? */ myproc = CCTK_MyProc (GH); /* Open the files on the first trip through if we are proc. 0 */ if (myproc == 0) { /* see if output file for this alias name was already created */ fdset_2D = (FILE **) GetNamedData (myGH->fileList_2D, alias); if (fdset_2D == NULL) { fdset_2D = (FILE **) malloc (3 * sizeof (FILE *)); filename = (char *) malloc (strlen (myGH->outdir2D) + strlen (alias) + sizeof (slicename) + 20); /* get the variable's full name */ fullname = CCTK_FullName (vindex); /* Open/Create files */ for (dir = 0; dir < maxdir; dir++) { /* FIXME: this can go when we permanently switch the the new filename scheme */ if (new_filename_scheme) { if (groupinfo.dim == 2) { strcpy (slicename, "2D"); } else { /* give the slice origin as range [1 .. n] */ sprintf (slicename, "%s_[%d]", extensions[dir], myGH->sp2xyz[groupinfo.dim - 1][dir]); } /* skip pathname if output goes into current directory */ if (strcmp (myGH->outdir2D, ".")) { sprintf (filename, "%s/%s_%s.asc", myGH->outdir2D, alias,slicename); } else { sprintf (filename, "%s_%s.asc", alias, slicename); } } else { /* skip pathname if output goes into current directory */ if (strcmp (myGH->outdir2D, ".")) { if (groupinfo.dim == 2) { sprintf (filename, "%s/%s_2D.gnuplot", myGH->outdir2D, alias); } else { sprintf (filename, "%s/%s_2d_%s.gnuplot", myGH->outdir2D, alias, extensions[dir]); } } else { if (groupinfo.dim == 2) { sprintf (filename, "%s_2D.gnuplot", alias); } else { sprintf (filename, "%s_2d_%s.gnuplot", alias, extensions[dir]); } } } /* if restart from recovery, try to open an existing file ... */ fdset_2D[dir] = NULL; if (IOUtil_RestartFromRecovery (GH)) { fdset_2D[dir] = fopen (filename, "a"); } /* otherwise or if that failed, create a new one */ if (! fdset_2D[dir]) { fdset_2D[dir] = fopen (filename, "w"); } if (! fdset_2D[dir]) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Cannot open 2D output file '%s'", filename); return; } /* advertise the file for downloading and write file info */ /* FIXME: this can go when we permanently switch the the new filename scheme */ advertised_file.slice = new_filename_scheme ? slicename : extensions[dir]; advertised_file.thorn = CCTK_THORNSTRING; advertised_file.varname = fullname; advertised_file.description = "Two-dimensional slice plots"; advertised_file.mimetype = "application/gnuplot"; IOUtil_AdvertiseFile (GH, filename, &advertised_file); if (CCTK_Equals (out_fileinfo, "parameter filename") || CCTK_Equals (out_fileinfo, "all")) { buffer[0] = 0; CCTK_ParameterFilename (sizeof (buffer), buffer); fprintf (fdset_2D[dir], "#Parameter file %s\n", buffer); } if (CCTK_Equals (out_fileinfo, "creation date") || CCTK_Equals (out_fileinfo, "all")) { buffer[0] = 0; Util_CurrentDate (sizeof (buffer), buffer); fprintf (fdset_2D[dir], "#Created %s ", buffer); Util_CurrentTime (sizeof (buffer), buffer); fprintf (fdset_2D[dir], "%s\n", buffer); } if (CCTK_Equals (out_fileinfo, "axis labels") || CCTK_Equals (out_fileinfo, "all")) { fprintf (fdset_2D[dir], "#x-label %c\n", extensions[dir][0]); fprintf (fdset_2D[dir], "#y-label %c\n", extensions[dir][1]); fprintf (fdset_2D[dir], "#z-label %s", advertised_file.varname); if (groupinfo.dim != 2) { if (have_coords) { fprintf (fdset_2D[dir], zlabel_fmt_string, 'x' + (maxdir-dir-1), coord_lower[maxdir-dir-1] + GH->cctk_delta_space[maxdir-dir-1] * myGH->sp2xyz[groupinfo.dim-1][dir]); } fprintf (fdset_2D[dir], " (%ci = %d)", 'x' + (maxdir-dir-1), myGH->sp2xyz[groupinfo.dim-1][dir]); } fputc ('\n', fdset_2D[dir]); } } /* store file desriptors in database */ StoreNamedData (&myGH->fileList_2D, alias, fdset_2D); free (filename); free (fullname); } } for (dir = 0; dir < maxdir; dir++) { int dir_i, dir_j; int directions[3]; const int lengths[2] = {-1, -1}; const int downsamples[2] = {1, 1}; int hsize[2]; CCTK_REAL *coord_data_i, *coord_data_j; void *data; /* get the directions to span the hyperslab */ if (dir == 0) { dir_i = 0; dir_j = 1; /* xy */ } else if (dir == 1) { dir_i = 0; dir_j = 2; /* xz */ } else { dir_i = 1; dir_j = 2; /* yz */ } /* set the origin using the slice center from IOUtil */ memset (origin, 0, GH->cctk_dim * sizeof (int)); origin[maxdir-dir-1] = myGH->sp2xyz[groupinfo.dim - 1][dir]; /* set the directions vector */ memset (directions, 0, sizeof (directions)); directions[maxdir-dir-1] = 1; /* get the coordinates for grid function output */ if (have_coords) { /* get the i-coordinate slice */ if (Hyperslab_GetHyperslab (GH, 0, coord_index[dir_i], 0, 2, origin, directions, lengths, downsamples, (void **) &coord_data_i, hsize) < 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Failed to extract 2D hyperslab for %c-coordinate", 'x' + dir_i); return; } /* get the j-coordinate slice */ if (Hyperslab_GetHyperslab (GH, 0, coord_index[dir_j], 0, 2, origin, directions, lengths, downsamples, (void **) &coord_data_j, hsize) < 0) { CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Failed to extract 2D hyperslab for %c-coordinate", 'x' + dir_j); free (coord_data_i); return; } } else { /* grid arrays don't have coordinates assigned */ coord_data_i = NULL; coord_data_j = NULL; } /* get the variable slice */ if (Hyperslab_GetHyperslab (GH, 0, vindex, 0, 2, origin, directions, lengths, downsamples, &data, hsize) < 0) { fullname = CCTK_FullName (vindex); CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Failed to extract 2D hyperslab for variable '%s'", fullname); free (fullname); if (coord_data_i) { free (coord_data_i); } if (coord_data_j) { free (coord_data_j); } return; } /* proc 0 writes */ if (myproc == 0) { CCTK_REAL stagger_offset_i, stagger_offset_j; /* get the staggering offset for the coordinates */ stagger_offset_i = CCTK_StaggerDirIndex (dir_i, groupinfo.stagtype) * 0.5 * GH->cctk_delta_space[dir_i]; stagger_offset_j = CCTK_StaggerDirIndex (dir_j, groupinfo.stagtype) * 0.5 * GH->cctk_delta_space[dir_j]; /* print out header */ fprintf (fdset_2D[dir], header_fmt_string, GH->cctk_time); switch (groupinfo.vartype) { case CCTK_VARIABLE_CHAR: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_BYTE, int, data, hsize, 0, have_coords, data_fmt_string_int,fdset_2D[dir]); break; case CCTK_VARIABLE_INT: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_INT, int, data, hsize, 0, have_coords, data_fmt_string_int,fdset_2D[dir]); break; case CCTK_VARIABLE_REAL: case CCTK_VARIABLE_COMPLEX: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_REAL, double, data, hsize, groupinfo.vartype == CCTK_VARIABLE_COMPLEX, have_coords, data_fmt_string_real, fdset_2D[dir]); break; #ifdef CCTK_INT2 case CCTK_VARIABLE_INT2: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_INT2, int, data, hsize, 0, have_coords, data_fmt_string_int,fdset_2D[dir]); break; #endif #ifdef CCTK_INT4 case CCTK_VARIABLE_INT4: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_INT4, int, data, hsize, 0, have_coords, data_fmt_string_int,fdset_2D[dir]); break; #endif #ifdef CCTK_INT8 case CCTK_VARIABLE_INT8: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_INT8, int, data, hsize, 0, have_coords, data_fmt_string_int,fdset_2D[dir]); break; #endif #ifdef CCTK_REAL4 case CCTK_VARIABLE_REAL4: case CCTK_VARIABLE_COMPLEX8: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_REAL4, double, data, hsize, groupinfo.vartype == CCTK_VARIABLE_COMPLEX8, have_coords, data_fmt_string_real, fdset_2D[dir]); break; #endif #ifdef CCTK_REAL8 case CCTK_VARIABLE_REAL8: case CCTK_VARIABLE_COMPLEX16: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_REAL8, double, data, hsize, groupinfo.vartype == CCTK_VARIABLE_COMPLEX16, have_coords, data_fmt_string_real, fdset_2D[dir]); break; #endif #ifdef CCTK_REAL16 case CCTK_VARIABLE_REAL16: case CCTK_VARIABLE_COMPLEX32: OUTPUT_TYPED_DATA (coord_data_i, coord_data_j, stagger_offset_i, stagger_offset_j, CCTK_REAL16, double, data, hsize, groupinfo.vartype == CCTK_VARIABLE_COMPLEX32, have_coords, data_fmt_string_real, fdset_2D[dir]); break; #endif default: CCTK_WARN (1, "Unsupported variable type"); break; } /* keep the file open but flush it */ fflush (fdset_2D[dir]); /* free the hyperslabs */ free (data); if (coord_data_i) { free (coord_data_i); } if (coord_data_j) { free (coord_data_j); } } /* end of outputting the data by processor 0 */ } /* end of looping through xyz directions */ }