path: root/src/Write2D.c

/*@@
   @file      Write2D.c
   @date      Thu May 11 2000
   @author    Thomas Radke
   @desc 
              Output two-dimensional slices in ASCII gnuplot format.
   @enddesc 
   @version   $Id$
 @@*/


#include <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "cctk_Parameters.h"
#include "CactusPUGH/PUGHSlab/src/PUGHSlab.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 = "$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_hsizes, coord_i, coord_j, stagger_offset_i,   \
                          stagger_offset_j, data_hsizes, cctk_type, c_type,   \
                          data, grouptype, fmt_string, file)                  \
          {                                                                   \
            int i, j;                                                         \
            cctk_type *typed_data = (cctk_type *) (data);                     \
                                                                              \
                                                                              \
            /* grid functions are output along with their coordinates         \
               grid array points are just indexed */                          \
            if (grouptype == CCTK_GF)                                         \
            {                                                                 \
              for (j = 0; j < data_hsizes[1]; j++)                            \
              {                                                               \
                for (i = 0; i < data_hsizes[0]; i++)                          \
                {                                                             \
                  fprintf (file, fmt_string,                                  \
                           (double) (*coord_i++ + stagger_offset_i),          \
                           (double) (*coord_j++ + stagger_offset_j),          \
                           (c_type) *typed_data++);                           \
                }                                                             \
                fprintf (file, "\n");                                         \
                coord_i += coord_hsizes[0] - data_hsizes[0];                  \
                coord_j += coord_hsizes[0] - data_hsizes[0];                  \
              }                                                               \
              coord_i -= coord_hsizes[0] * data_hsizes[1];                    \
              coord_j -= coord_hsizes[0] * data_hsizes[1];                    \
            }                                                                 \
            else                                                              \
            {                                                                 \
              for (j = 0; j < data_hsizes[1]; j++)                            \
              {                                                               \
                for (i = 0; i < data_hsizes[0]; i++)                          \
                {                                                             \
                  fprintf (file, fmt_string,                                  \
                           (double) i,                                        \
                           (double) j,                                        \
                           (c_type) *typed_data++);                           \
                }                                                             \
                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     cGH *
   @vio       in
   @endvar
   @var       index
   @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 (cGH *GH, int index, const char *alias)
{
  DECLARE_CCTK_PARAMETERS
  int myproc;
  asciiioGH *myGH;
  const char *header_fmt_string;      /* header format string */
  const char *data_fmt_string_int;    /* data format string for int types */
  const char *data_fmt_string_real;   /* data format string for float types */
  int dir;
  int timelevel;
  int groupindex;
  cGroup groupinfo;
  FILE **fdset_2D;                    /* array of output file pointers */
  int coord_index[3];                 /* variable indices for xyz coordinates */
  int coord_timelevel[3];             /* coordinates' current timelevels */
  int origin[3];                      /* the slice origin */
  char *filename;
  char *fullname;
  char slicename[20];
  ioAdvertisedFileDesc advertised_file; 
  const char *extensions[3] = {"yz", "xz", "xy"};


  /* to make the compiler happy */
  fdset_2D = NULL;

  /* get the variable group indormation */
  groupindex = CCTK_GroupIndexFromVarI (index);
  CCTK_GroupData (groupindex, &groupinfo);

  /* check if variable has storage assigned */
  if (! CCTK_QueryGroupStorageI (GH, groupindex)) 
  {
    fullname = CCTK_FullName (index);
    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");

  /* set header and data format strings */
  if (CCTK_Equals (out_format, "f"))
  {
    header_fmt_string    = "\n\n#Time = %f\n";
    data_fmt_string_int  = "%.13f\t\t%.13f\t\t%d\n";
    data_fmt_string_real = "%.13f\t\t%.13f\t\t%.13f\n";
  }
  else
  {
    header_fmt_string    = "\n\n#Time = %e\n";
    data_fmt_string_int  = "%.13e\t\t%.13e\t\t%d\n";
    data_fmt_string_real = "%.13e\t\t%.13e\t\t%.13e\n";
  }

  /* 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 +
#if 0
FIXME: get rid of PUGH here
                               strlen (pughGH->identity_string) + 20);
#else
                               0);
#endif
 
      /* Open/Create files */
      for (dir = 0; dir < 3; dir++)
      {
#if 0
        /* FIXME: move identity_string into cGH ?? */
        sprintf (filename, "%s/%s_2d_%s%s.gnuplot", myGH->outdir2D, alias,
                 extensions[dir], pughGH->identity_string);
#else
        /* FIXME: this can go when we permanently switch the the
                  new filename scheme */
        if (new_filename_scheme)
        {
          /* give the slice origin as range [1 .. n] */
          sprintf (slicename, "%s_[%d]", extensions[dir],
                   myGH->sp2xyz[groupinfo.dim - 1][dir] + 1);

          sprintf (filename, "%s/%s_%s.asc", myGH->outdir2D, alias, slicename);
        }
        else
        {
          sprintf (filename, "%s/%s_2d_%s.gnuplot", myGH->outdir2D, alias,
                   extensions[dir]);
        }
#endif

        /* 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 */
        /* FIXME: this can go when we permanently switch the the
                  new filename scheme */
        advertised_file.slice = new_filename_scheme ?
                                slicename : (char *) extensions[dir];
        advertised_file.thorn = CCTK_THORNSTRING;
        advertised_file.varname = CCTK_FullName (index);
        advertised_file.description = "Two-dimensional slice plots";
        advertised_file.mimetype = "application/gnuplot";

        IOUtil_AdvertiseFile (GH, filename, &advertised_file);

        free (advertised_file.varname);
      }

      /* store file desriptors in database */
      StoreNamedData (&myGH->fileList_2D, alias, fdset_2D);

      free (filename);
    }
  }

  /* get the coordinate indices if we output a grid function */
  if (groupinfo.grouptype == CCTK_GF)
  {
    switch (groupinfo.dim)
    {
      case 1:
        coord_index[0] = CCTK_CoordIndex (1, NULL, "cart1d");
        break;
      case 2:
        coord_index[0] = CCTK_CoordIndex (1, NULL, "cart2d");
        coord_index[1] = CCTK_CoordIndex (2, NULL, "cart2d");
      case 3:
        coord_index[0] = CCTK_CoordIndex (1, NULL, "cart3d");
        coord_index[1] = CCTK_CoordIndex (2, NULL, "cart3d");
        coord_index[2] = CCTK_CoordIndex (3, NULL, "cart3d");
        break;
      default:
        CCTK_WARN (4, "Cannot find appropriate coordinate system");
        break;
    }

    /* get the coordinate timelevels */
    for (dir = 0; dir < 3; dir++)
    {
      coord_timelevel[dir] = CCTK_NumTimeLevelsFromVarI (coord_index[dir]) - 1;
      if (coord_timelevel[dir] > 0)
      {
        coord_timelevel[dir]--;
      }
    }
  }

  /* get the timelevel for the variable to output */
  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
  if (timelevel > 0)
  {
    timelevel--;
  }

  for (dir = 0; dir < 3; dir++) 
  {
    int dir_i, dir_j;
    int directions[3];
    const int lengths[2] = {-1, -1};
    const int downsamples[2] = {1, 1};
    int coord_hsizes[2], data_hsizes[2];
    CCTK_REAL *coord_data_i, *coord_data_j;
    void *data;


    /* get the directions to span the hyperslab */
    if (dir == 0)
    {
      dir_i = 1; dir_j = 2;   /* yz */
    }
    else if (dir == 1)
    {
      dir_i = 0; dir_j = 2;   /* xz */
    }
    else
    {
      dir_i = 0; dir_j = 1;   /* xy */
    }

    /* set the origin using the slice center from IOUtil */
    memset (origin, 0, GH->cctk_dim * sizeof (int));
    origin[dir] = myGH->sp2xyz[groupinfo.dim - 1][dir];

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

    /* get the coordinates for grid function output */
    if (groupinfo.grouptype == CCTK_GF)
    {
      /* get the i-coordinate slice */
      if (Hyperslab_GetHyperslab (GH, 0, coord_index[dir_i],
                                  coord_timelevel[dir_i], 2, origin, directions,
                                  lengths, downsamples,
                                  (void **) &coord_data_i, coord_hsizes) < 0)
      {
        fullname = CCTK_FullName (coord_index[dir_i]);
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Failed to extract 2D hyperslab for variable '%s'",
                    fullname);
        free (fullname);
        return;
      }

      /* get the j-coordinate slice */
      if (Hyperslab_GetHyperslab (GH, 0, coord_index[dir_j],
                                  coord_timelevel[dir_j], 2, origin, directions,
                                  lengths, downsamples,
                                  (void **) &coord_data_j, coord_hsizes) < 0)
      {
        fullname = CCTK_FullName (coord_index[dir_j]);
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Failed to extract 2D hyperslab for variable '%s'",
                    fullname);
        free (fullname);
        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, index, timelevel, 2, origin, directions,
                                lengths, downsamples, &data, data_hsizes) < 0)
    {
      fullname = CCTK_FullName (index);
      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_hsizes, coord_data_i, coord_data_j,
                             stagger_offset_i, stagger_offset_j,
                             data_hsizes, CCTK_CHAR, int, data,
                             groupinfo.grouptype,
                             data_fmt_string_int, fdset_2D[dir]);
          break;

        case CCTK_VARIABLE_INT:
          OUTPUT_TYPED_DATA (coord_hsizes, coord_data_i, coord_data_j,
                             stagger_offset_i, stagger_offset_j,
                             data_hsizes, CCTK_INT, int, data,
                             groupinfo.grouptype,
                             data_fmt_string_int, fdset_2D[dir]);
          break;

        case CCTK_VARIABLE_REAL:
          OUTPUT_TYPED_DATA (coord_hsizes, coord_data_i, coord_data_j,
                             stagger_offset_i, stagger_offset_j,
                             data_hsizes, CCTK_REAL, double, data,
                             groupinfo.grouptype,
                             data_fmt_string_real, fdset_2D[dir]);
          break;

        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 */

}