path: root/src/Write3D.c

/*@@
   @file      Write3D.c
   @date      Sat 12 June 2004
   @author    Thomas Radke
   @desc
              Three-dimensional output of variables 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 = "$Header$";
CCTK_FILEVERSION(CactusIO_IOSDF_Write3D_c)


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
static char *OpenFile (const cGH *GH, const char *fullname, const char *alias);

/*@@
   @routine   IOSDF_Write3D
   @date      Sat 12 June 2004
   @author    Thomas Radke
   @desc
              Writes the 3D volume of a variable into a gnuplot SDF file.
   @enddesc
   @calls     Hyperslab_GlobalMappingByIndex
              Hyperslab_FreeMapping
              Hyperslab_GetList
              OpenFile

   @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_Write3D (const cGH *GH, int vindex, const char *alias)
{
  int i, total_hsize;
  int myproc, gindex, have_coords;
  int num_requested_hslabs, num_returned_hslabs;
  cGroup gdata;
  CCTK_INT coord_system_handle, coord_handles[3];
  char *fullname, *groupname, *filename;
  double *hdata[4];
  int extent_int[3];
  double offset[3];
  int mapping;
  CCTK_INT vindices[4], extent[3], downsample[3], hsize[3];
  const double dtime = GH->cctk_time;
  const CCTK_INT origin[] = {0, 0, 0},
                 direction[] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
  DECLARE_CCTK_PARAMETERS


  /* get the variable 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 3D output for '%s' (no storage)", fullname);
    free (fullname);
    return (-1);
  }

  /* get the coordinate system associated with this grid variable */
  vindices[0] = vindex;
  groupname = CCTK_GroupName (gindex);
  coord_system_handle = Coord_GroupSystem (GH, groupname);
  free (groupname);

  have_coords = coord_system_handle >= 0 &&
                Util_TableGetIntArray (coord_system_handle, 3,
                                       coord_handles, "COORDINATES") >= 0;
  if (have_coords)
  {
    /* get the coordinate functions and coordinate physical minimum */
    for (i = 1; i <= 3; i++)
    {
      vindices[i] = -1;
      Util_TableGetInt (coord_handles[i-1], &vindices[i], "GAINDEX");
      have_coords &= vindices[i] >= 0;
    }
  }
  num_requested_hslabs = have_coords ? 4 : 1;

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

  /* Open the file on processor 0 */
  filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL;

  /* get the total number of grid points to check for zero-sized variables */
  /* set the extent vector (copy from 'int' to 'CCTK_INT') */
  CCTK_GroupgshVI (GH, 3, extent_int, vindex);
  for (i = 0, total_hsize = 1; i < 3; i++)
  {
    total_hsize *= extent_int[i];
    extent[i] = extent_int[i];
  }
  if (total_hsize <= 0)
  {
    free (fullname);
    return (0);
  }

  downsample[0] = out_downsample_x;
  downsample[1] = out_downsample_y;
  downsample[2] = out_downsample_z;

  /* get the hyperslab mapping */
  mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 3,
                                            direction, origin, extent,
                                            downsample,
                                            -1,   /* table handle */
                                            NULL  /* conversion fn */,
                                            hsize);
  if (mapping < 0)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOSDF_Write3D: Failed to define hyperslab mapping for "
                "variable '%s'", fullname);
    free (fullname);
    return (-1);
  }
  for (i = 0, total_hsize = 1; i < 3; i++)
  {
    extent_int[i] = hsize[i];
    total_hsize *= hsize[i];
  }
  if (total_hsize <= 0)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOSDF_Write3D: selected hyperslab has zero size for "
                "variable '%s'", fullname);
    Hyperslab_FreeMapping (mapping);
    free (fullname);
    return (-1);
  }

  hdata[0] = hdata[1] = hdata[2] = hdata[3];
  if (myproc == 0)
  {
    /* allocate hyperslab buffers */
    hdata[0] = malloc (total_hsize * sizeof (double));
    hdata[1] = have_coords ? malloc (3 * total_hsize * sizeof (double)) : NULL;
    hdata[2] = hdata[1] + 1*total_hsize;
    hdata[3] = hdata[1] + 2*total_hsize;
  }
 
  /* get the hyperslabs */
  num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs,
                                           NULL, vindices, NULL, NULL,
                                           (void **) hdata, NULL);

  /* And dump the data to file */
  if (myproc == 0)
  {
    if (num_returned_hslabs == num_requested_hslabs)
    {
      if (have_coords)
      {
        /* get the staggering offset for the coordinates */
        for (i = 0; i < 3; i++)
        {
          offset[i] = 0.5 * GH->cctk_delta_space[i] *
                      CCTK_StaggerDirIndex (i, gdata.stagtype);
        }
        for (i = 0; i < total_hsize; i++)
        {
          hdata[1][i] += offset[0];
          hdata[2][i] += offset[1];
          hdata[3][i] += offset[2];
        }
      }

      i = have_coords ? gft_out_full (filename, dtime, extent_int, "x|y|z", 3,
                                      hdata[1], hdata[0]) :
                        gft_out_brief (filename, dtime, extent_int, 3,hdata[0]);
      if (i <= 0)
      {
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Error writing to 3D IOSDF output file '%s'", filename);
      }
    }
    else
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "IOSDF_Write3D: Failed to extract hyperslab for "
                  "variable '%s'", fullname);
    }

    /* clean up */
    free (hdata[0]);
    free (hdata[1]);
  } /* end of outputting the data by processor 0 */

  free (fullname);

  return (0);
}


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
/*@@
   @routine    OpenFile
   @date       Sat 12 June 2004
   @author     Thomas Radke
   @desc
               Opens an SDF file for a given alias name.
               If this is the first time through, it will advertise it
               to IOUtil.
   @enddesc

   @returntype char *
   @returndesc
               the full filename of the SDF output file
   @endreturndesc
 @@*/
static char *OpenFile (const cGH *GH, const char *fullname, const char *alias)
{
  char *filename;
  ioSDFGH *myGH;
  ioAdvertisedFileDesc advertised_file;
  DECLARE_CCTK_PARAMETERS


  /* get handle for IOSDF GH extensions */
  myGH = CCTK_GHExtension (GH, "IOSDF");

  /* see if we are the first time through */
  filename = GetNamedData (myGH->fileList_3D, alias);
  if (! filename)
  {
    filename = malloc (strlen (myGH->out3D_dir) + strlen (alias) + 9);

    /* open/create the file */
    sprintf (filename, "%s%s_3D.sdf", myGH->out3D_dir, alias);

    /* advertise the file for downloading and write file info */
    advertised_file.slice = "";
    advertised_file.thorn = CCTK_THORNSTRING;
    advertised_file.varname = fullname;
    advertised_file.description = "Full-dimensional variable contents";
    advertised_file.mimetype = "application/sdf";

    IOUtil_AdvertiseFile (GH, filename, &advertised_file);

    /* store filename in database */
    StoreNamedData (&myGH->fileList_3D, alias, filename);
  }

  return (filename);
}