path: root/src/Write1D.c

/*@@
   @file      Write1D.c
   @date      March 1997
   @author    Paul Walker
   @desc 
              Output one-dimensional lines in ASCII xgraph/gnuplot format.
   @enddesc 

   @history
   @hauthor   Thomas Radke
   @hdate     30 May 2000
   @hdesc
              Get rid of all the PUGH stuff by using thorn Hyperslab.
   @endhdesc
   @hendhistory
   @version $Id$
@@*/

#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>  /* stat(2) */

#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 char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusBase_IOASCII_Write1D_c)

/*#define DEBUG_IOASCII 1*/

/* macro to output a 1D line (with coordinates for GFs) as typed data */
#define OUTPUT_TYPED_DATA(grouptype, hsize, coord_data, stagger_offset,       \
                          cctk_type, cctk_extract_fn, c_type, data,           \
                          format, file)                                       \
          {                                                                   \
            int h;                                                            \
            cctk_type *typed_data = (cctk_type *) data;                       \
                                                                              \
                                                                              \
            if (grouptype == CCTK_GF)                                         \
            {                                                                 \
              for (h = 0; h < hsize; h++)                                     \
              {                                                               \
                fprintf (file, format,                                        \
                         (double) (coord_data[h] + stagger_offset),           \
                         (c_type) cctk_extract_fn (typed_data[h]));           \
              }                                                               \
            }                                                                 \
            else                                                              \
            {                                                                 \
              for (h = 0; h < hsize; h++)                                     \
              {                                                               \
                fprintf (file, format,                                        \
                         (double) h,                                          \
                         (c_type) cctk_extract_fn (typed_data[h]));           \
              }                                                               \
            }                                                                 \
          }


/* this is needed by some preprocessors to pass into OUTPUT_TYPED_DATA
   as a dummy macro */
#define NOTHING


/*@@
   @routine IOASCII_Write1D
   @date    March 1999
   @author  Gabrielle Allen
   @desc 
            This routine does 1D line output along the orthogonals
            and the diagonal (in case of a cubed grid).
            <p>
            It writes to ASCII files suitable for gnuplot and xgraph.
            A header telling the physical time prefixes the output data.
   @enddesc 
   @calls   IOUtil_RestartFromRecovery
            IOUtil_AdvertiseFile
            Hyperslab_GetHyperslab

   @var     GH
   @vdesc   Pointer to CCTK GH
   @vtype   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 
@@*/

void IOASCII_Write1D (cGH *GH, int vindex, const char *alias)
{ 
  DECLARE_CCTK_PARAMETERS
  asciiioGH *myGH;                  /* IOASCII extension handle */
  int Do_it[4];                     /* flags indicating actual work */
  int coord_index[3];               /* x,y,z coordinate variable indices */
  int is_cubic;                     /* true if variable lives on cubic grid */
  int myproc;                       /* identify processor */
  int i, dir;                       /* Loopers */
  int timelevel;                    /* timelevel of variable to output */
  int groupindex;                   /* variable's group index */
  int have_coords;                  /* boolean for existance of coordinates */
  char coord_system[20];            /* name of the coordinate system */
  char slicename[20];               /* name of current output slice */
  cGroup group_static_data;         /* variable's group static data */
  cGroupDynamicData group_dynamic_data;/* variable's group dynamic data */
  char *fullname;                   /* variable's full name */
  const char *header_fmt;           /* header format string */
  const char *data_fmt_int;         /* data format string for int types */
  const char *data_fmt_real;        /* data format string for float types */
  const char *file_extension;       /* filename extension */
  FILE *file[8];                    /* fds for x,y,z,d output (also COMPLEX)*/
  int num_files;
  int upper, lower;
  struct stat fileinfo;
  const char *openmode;
  static char *extensions[] = {"xl", "yl", "zl", "dl"};
  char *filename, *type_extension;
  ioAdvertisedFileDesc advertisedFile;


  /* get the variable's group index */
  groupindex = CCTK_GroupIndexFromVarI (vindex);

  /* check if variable has storage assigned */
  if (! CCTK_QueryGroupStorageI (GH, groupindex)) 
  {
    fullname = CCTK_FullName (vindex);
    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOASCII_Write1D: No IOASCII_1D output for '%s' (no storage)",
                fullname);
    free (fullname);
    return;
  }

  /* get the handle for IOASCII extensions */
  myGH = (asciiioGH *) CCTK_GHExtension (GH, "IOASCII");

  /* get the variable's group information */
  CCTK_GroupData (groupindex, &group_static_data);
  CCTK_GroupDynamicData (GH, groupindex, &group_dynamic_data);

  /* see what slices should be output */
  Do_it[0] = out1D_x && group_static_data.dim >= 1;
  Do_it[1] = out1D_y && group_static_data.dim >= 2;
  Do_it[2] = out1D_z && group_static_data.dim >= 3;
  /* diagonal slice is done only if grid is cubic and non-staggered. */
  is_cubic = 1;
  for (i = 1; i < group_dynamic_data.dim; i++)
  {
    is_cubic &= group_dynamic_data.gsh[i] == group_dynamic_data.gsh[0];
  }
  if (out1D_d && ! is_cubic)
  {
    fullname = CCTK_FullName (vindex);
    CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOASCII_Write1D: No IOASCII_1D diagonal output for '%s' "
                "(only implemented for cubic grids)",
                fullname);
    free (fullname);
  }
  Do_it[3] = out1D_d &&
             group_static_data.dim == 3 &&
             is_cubic &&
             group_static_data.stagtype == 0;

  /* return if nothing to do */
  if (! (Do_it[0] || Do_it[1] || Do_it[2] || Do_it[3]))
  {
    return;
  }

  /* set header and data format strings */
  i = CCTK_Equals (out_style, "gnuplot");
  if (CCTK_Equals (out_format, "f"))
  {
    header_fmt    = i ? "\n\n#Time = %f\n" : "\n\n\"Time = %f\n";
    data_fmt_int  = "%f\t\t%d\n";
    data_fmt_real = "%f\t\t%.13f\n";
  }
  else
  {
    header_fmt    = i ? "\n\n#Time = %e\n" : "\n\n\"Time = %e\n";
    data_fmt_int  = "%e\t\t%d\n";
    data_fmt_real = "%e\t\t%.13e\n";
  }
  file_extension = i ? ".asc" : ".xg";

#ifdef DEBUG_IOASCII
  printf ("\nIn IOASCII Write1D\n------------------\n");
  printf ("  Variable index is %d\n",vindex);
  printf ("  Alias is -%s-\n",alias);
  fflush (stdout);
#endif

  /* What processor are we on? */
  myproc = CCTK_MyProc (GH);

  /* Processor 0 opens the files with the appropriate name */
  if (myproc == 0) 
  {
    /* 20 extra characters should be enough for '/',
       the type extension, the file extension, and the trailing '\0' */
    i = strlen (myGH->outdir1D) + strlen (alias) + sizeof (slicename) + 20 +
#if 0
    /* FIXME: pughGH->identity_string should be moved into cGH */
          strlen (pughGH->identity_string);
#else
          0;
#endif
    filename = (char *) malloc (i);

    num_files = group_static_data.vartype == CCTK_VARIABLE_COMPLEX ? 8 : 4;

    for (i = 0; i < num_files; i++)
    {
      /* skip empty slices */
      if (! Do_it[i % 4])
      {
        continue;
      }

      /* add suffix for complex type variables */
      if (group_static_data.vartype == CCTK_VARIABLE_COMPLEX)
      {
        if (new_filename_scheme)
        {
          type_extension = i < 4 ? "re_" : "im_";
        }
        else
        {
          type_extension = i < 4 ? "_re" : "_im";
        }
      }
      else
      {
        type_extension = "";
      }
#if 0
      /* FIXME: pughGH->identity_string should be moved into cGH */
      sprintf (filename, "%s/%s%s%s.%s", myGH->outdir1D, alias,
               type_extension, pughGH->identity_string, extensions[i % 4]);
#else
      /* FIXME: this can go after the old filename scheme has gone */
      if (new_filename_scheme)
      {
        if ((i + 1) % 4)
        {
          /* get the indices into spxyz[] */
          lower = ((i % 4) + 1) % 3;
          upper = ((i % 4) + 2) % 3;
          if (upper < lower)
          {
            upper = lower;
            lower = 0;
          }
          /* give the slice origin as range [1 .. n] */
          sprintf (slicename, "%s%c_[%d][%d]", type_extension, 'x' + i % 4,
                   myGH->spxyz[group_static_data.dim-1][i % 4][lower] + 1,
                   myGH->spxyz[group_static_data.dim-1][i % 4][upper] + 1);
        }
        else
        {
          sprintf (slicename, "%s%dD_diagonal", type_extension,
                   group_static_data.dim);
        }
        sprintf (filename, "%s/%s_%s%s", myGH->outdir1D, alias,
                 slicename, file_extension);
      }
      else
      {
        sprintf (filename, "%s/%s%s.%s", myGH->outdir1D, alias,
                 type_extension, extensions[i % 4]);
      }
#endif

      /* see if output file was already created */
      if (GetNamedData (myGH->filenameList1D, filename) == NULL) 
      {
        /* if restart from recovery, all existing files are opened
           in append mode */
        if (IOUtil_RestartFromRecovery (GH))
        {
          openmode = stat (filename, &fileinfo) == 0 ? "a" : "w";
        }
        else
        {
          openmode = "w";
        }

        /* just store a non-NULL pointer in database */
        StoreNamedData (&myGH->filenameList1D, filename, (void *) 1);
      }
      else
      {
        openmode = "a";
      }

      if (! (file[i] = fopen (filename, openmode)))
      {
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "IOASCII_Write1D: Cannot open 1D output file '%s'", 
                    filename);
      } 
      else
      {
        /* advertise new files for downloading */
        if (*openmode == 'w')
        {
          /* FIXME: this can go after the old filename scheme has gone */
          advertisedFile.slice = new_filename_scheme ?
                                 slicename : extensions[i % 4];
          advertisedFile.thorn = CCTK_THORNSTRING;
          advertisedFile.varname = CCTK_FullName (vindex);
          advertisedFile.description = "One-dimensional line plots";
          advertisedFile.mimetype = CCTK_Equals (out_style, "gnuplot") ?
                               "application/gnuplot" : "application/x-graph";

          IOUtil_AdvertiseFile (GH, filename, &advertisedFile);

          free (advertisedFile.varname);
        }
      }
    }
    free (filename);
  }

  /* get the coordinate indices for CCTK_GF variables */
  if (group_static_data.grouptype == CCTK_GF)
  {
    sprintf (coord_system, "cart%dd", group_static_data.dim);
    have_coords = 1;
    for (i = 0; i < group_static_data.dim; i++)
    {
      coord_index[i] = CCTK_CoordIndex (i + 1, NULL, coord_system);
      have_coords &= coord_index[i] >= 0;
    }

    if (! have_coords)
    {
      CCTK_VWarn (8, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOASCII_Write1D: No coordinate system '%s' found for '%s'",
                  coord_system, CCTK_VarName (vindex));
    }
  }
  else
  {
    /* CCTK_ARRAY variables never have coordinates associated */
    have_coords = 0;
  }

  /* get the current time level */
  timelevel = CCTK_NumTimeLevelsFromVarI (vindex) - 1;
  if (timelevel > 0)
  {
    timelevel--;
  }

  /* OK so actually do the I/O in each direction */
  for (dir = 0; dir < 4; dir++) 
  {
    const int length = -1;
    const int downsample = 1;
    const int *origin;
    const int zero_point[3] = {0, 0, 0};
    int directions[3];
    int hsize;
    int coord_timelevel;
    void *data;
    CCTK_REAL *coord_data;


    /* skip empty slices */
    if (! Do_it[dir])
    {
      continue;
    }

    if (group_static_data.grouptype == CCTK_GF)
    {
      /* for GFs: get the coordinate's 1D data (in xyz direction only) */
      if (dir < 3)
      {
        /* set the origin of the line */
        origin = myGH->spxyz[group_static_data.dim-1][dir];

        /* set the direction vector for the 1D line */
        memset (directions, 0, sizeof (directions));
        directions[dir] = 1;

        if (have_coords)
        {
          /* get the current time level for the coordinates */
          coord_timelevel = CCTK_NumTimeLevelsFromVarI (coord_index[dir]) - 1;
          if (coord_timelevel > 0)
          {
            coord_timelevel--;
          }

          if (Hyperslab_GetHyperslab (GH, 0, coord_index[dir], 
                                      coord_timelevel, 1, origin, directions, 
                                      &length, &downsample,
                                      (void **) &coord_data, &hsize) < 0)
          {
            fullname = CCTK_FullName (coord_index[dir]);
            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                        "IOASCII_Write1D: Failed to extract hyperslab for"
                        " variable '%s'",
                        fullname);
            free (fullname);
            continue;
          }
        }
        else
        {
          coord_data = NULL;
        }
      }
      else
      {
        /* set the origin of the line */
        origin = zero_point;

        /* set the direction vector for the diagonal 1D line */
        for (i = 0; i < 3; i++)
        {
          directions[i] = 1;
        }

        /* coordinates are calculated by output processor */
        coord_data = NULL;
      }
    }
    else
    {
      /* set the origin of the line */
      /* FIXME: Should we define some slice center for arrays too ? */
      origin = zero_point;

      /* set the direction vector for the 1D line */
      for (i = 0; i < 3; i++)
      {
        directions[i] = (dir == i || dir == 3) ? 1 : 0;
      }

      /* no coordinates are needed for arrays */
      coord_data = NULL;
    }

    /* get the variable's 1D data */
    if (Hyperslab_GetHyperslab (GH, 0, vindex, timelevel, 1, origin,
                                directions, &length, &downsample,
                                &data, &hsize) < 0)
    {
      fullname = CCTK_FullName (vindex);
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOASCII_Write1D: Failed to extract hyperslab for "
                  "variable '%s'", fullname);
      free (fullname);
      if (coord_data)
      {
        free (coord_data);
      }
      continue;
    }

    /* And write it out on processor 0 */
    if (myproc == 0 && file[dir] != NULL) 
    {
      CCTK_REAL coord_lower, offset;


      if (group_static_data.grouptype == CCTK_GF)
      {
        if (dir < 3)
        {
          /* get the staggering offset for the xyz coordinates */
          offset = CCTK_StaggerDirIndex (dir, group_static_data.stagtype) *
                   0.5 * GH->cctk_delta_space[dir];
          if (! have_coords)
          {
            coord_data = (CCTK_REAL *) malloc (hsize * sizeof (CCTK_REAL));
            for (i = 0; i < hsize; i++)
            {
              coord_data[i] = i;
            }
          }
        }
        else
        {
          /* get the diagonal coordinates */
          offset = have_coords ? GH->cctk_delta_space[0] * sqrt (3) : sqrt(3);
          coord_data = (CCTK_REAL *) malloc (hsize * sizeof (CCTK_REAL));
          for (i = 0; i < hsize; i++)
          {
            coord_data[i] = i * offset;
          }
          if (have_coords)
          {
            CCTK_CoordRange (GH, &coord_lower, &offset, 1, NULL, "cart3d");
            offset = coord_lower * sqrt (3);
          }
          else
          {
            offset = 0.0;
          }
        }
      }

      /* print out header */
      fprintf (file[dir], header_fmt, GH->cctk_time);
      if (group_static_data.vartype == CCTK_VARIABLE_COMPLEX)
      {
        fprintf (file[dir + 4], header_fmt, GH->cctk_time);
      }

      /* and then loop through the line points */
      switch (group_static_data.vartype)
      {
        case CCTK_VARIABLE_CHAR:
          OUTPUT_TYPED_DATA (group_static_data.grouptype, hsize, coord_data,
                             offset, CCTK_BYTE, NOTHING, int, data,
                             data_fmt_int, file[dir]);
          break;

        case CCTK_VARIABLE_INT:
          OUTPUT_TYPED_DATA (group_static_data.grouptype, hsize, coord_data,
                             offset, CCTK_INT, NOTHING, int, data,
                             data_fmt_int, file[dir]);
          break;

        case CCTK_VARIABLE_REAL:
          OUTPUT_TYPED_DATA (group_static_data.grouptype, hsize, coord_data,
                             offset, CCTK_REAL, NOTHING, double, data,
                             data_fmt_real, file[dir]);
          break;

        case CCTK_VARIABLE_COMPLEX:
          OUTPUT_TYPED_DATA (group_static_data.grouptype, hsize, coord_data,
                             offset, CCTK_COMPLEX, CCTK_CmplxReal, double, data,
                             data_fmt_real, file[dir]);
          OUTPUT_TYPED_DATA (group_static_data.grouptype, hsize, coord_data,
                             offset, CCTK_COMPLEX, CCTK_CmplxImag, double, data,
                             data_fmt_real, file[dir + 4]);
          break;

        default:
          CCTK_WARN (1, "Unsupported variable type");
          break;
      }

      /* close the output file(s) */
      fclose (file[dir]);
      if (group_static_data.vartype == CCTK_VARIABLE_COMPLEX)
      {
        fclose (file[dir + 4]);
      }

      /* clean up */
      free (data);
      if (coord_data)
      {
        free (coord_data);
      }
    }
    
  } /* end of loop through all directions */

}