path: root/src/Write2D.c

/*@@
   @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);
}