path: root/src/DumpVar.c

/*@@
   @file      DumpVar.c
   @author    Paul Walker, Gabrielle Allen, Tom Goodale, Thomas Radke
   @desc
              Routines for writing variables into IOFlexIO data or
              checkpoint files.
   @enddesc
   @version   $Id$
 @@*/


#include <stdio.h>
#include <stdlib.h>

#include "cctk.h"
#include "cctk_Parameters.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"
#include "CactusBase/IOUtil/src/ioGH.h"
#include "ioFlexGH.h"


/* the rcs ID and its dummy funtion to use it */
static const char *rcsid = "$Id$";
CCTK_FILEVERSION(CactusPUGHIO_IOFlexIO_DumpVar_c)


/********************************************************************
 ********************    Macro Definitions   ************************
 ********************************************************************/
/* #define IOFLEXIO_DEBUG  1 */

#define IOTAGBASE 20000 /* This may break on more than 2000 processors */


/********************************************************************
 ********************    Internal Typedefs   ************************
 ********************************************************************/
/* info structure describing how to dump a given CCTK variable type */
typedef struct
{
  int flexio_type;              /* the FLEXIO type to use */
  int element_size;             /* size of a single data element */
#ifdef CCTK_MPI
  MPI_Datatype mpi_type;        /* the data type for MPI communication */
#endif
} dumpInfo_t;


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
static int DumpGS (const cGH *GH, int vindex, int timelevel, IOFile file,
                   int flexio_type);
static int DumpGA (const cGH *GH, int vindex, int timelevel, IOFile file,
                   dumpInfo_t *info);
static int GetDumpData (const cGH *GH, int vindex, int timelevel, void **outme,
                        int *free_outme, CCTK_INT4 *geom, int element_size);
static void EachProcDump (const cGH *GH, int vindex, int timelevel, void *outme,
                          CCTK_INT4 *geom, IOFile file, int flexio_type);
static void AddChunkAttributes (const cGH *GH, int vindex, CCTK_INT4 *geom,
                                IOFile file);
static void AddCommonAttributes (const cGH *GH, int vindex, int timelevel,
                                 CCTK_INT4 *gsz, IOFile file);
#ifdef CCTK_MPI
static void ProcDump (IOFile file, const cGH *GH, int vindex, void *outme,
                      CCTK_INT4 *geom, int flexio_type);
#endif


/*@@
   @routine    IOFlexIO_DumpVar
   @date       16 Apr 1999
   @author     Thomas Radke
   @desc
               Generic dump routine, just calls the appropriate dump routine
   @enddesc

   @calls

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      global index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to store
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar

   @returntype int
   @returndesc
               -1 if variable has unknown type
               or return code of either @seeroutine DumpGS or @seeroutine DumpGA
   @endreturndesc
@@*/
int IOFlexIO_DumpVar (const cGH *GH, int vindex, int timelevel, IOFile file)
{
  int vtype, gtype, retval;
  pGH *pughGH;
  ioGH *ioUtilGH;
  dumpInfo_t info;


  /* get the handle for PUGH and IOUtil extensions */
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");
  pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");

  vtype = CCTK_VarTypeI (vindex);

  /* downgrade the output precision if requested */
  if (ioUtilGH->out_single)
  {
    /* Do we only want to downgrade generic CCTK types ? */
#ifdef CCTK_REAL4
    if (vtype == CCTK_VARIABLE_REAL)
    {
      vtype = CCTK_VARIABLE_REAL4;
    }
    else if (vtype == CCTK_VARIABLE_COMPLEX)
    {
      vtype = CCTK_VARIABLE_COMPLEX8;
    }
#endif
#ifdef CCTK_INT4
    if (vtype == CCTK_VARIABLE_INT)
    {
      vtype = CCTK_VARIABLE_INT4;
    }
#endif
  }
  info.element_size = CCTK_VarTypeSize (vtype);
#ifdef CCTK_MPI
  info.mpi_type = PUGH_MPIDataType (pughGH, vtype);
#endif
  info.flexio_type = IOFlexIO_DataType (vtype);

  if (info.element_size <= 0 ||
#ifdef CCTK_MPI
      info.mpi_type == MPI_DATATYPE_NULL ||
#endif
      info.flexio_type < 0)
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Unsupported variable type %d", vtype);
    return (-1);
  }

  /* now branch to the appropriate dump routine */
  gtype = CCTK_GroupTypeFromVarI (vindex);
  if (gtype == CCTK_SCALAR)
  {
    retval = DumpGS (GH, vindex, timelevel, file, info.flexio_type);
  }
  else if (gtype == CCTK_ARRAY || gtype == CCTK_GF)
  {
    retval = DumpGA (GH, vindex, timelevel, file, &info);
  }
  else
  {
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Invalid group type %d", gtype);
    retval = -1;
  }

  return (retval);
}


/*@@
   @routine    IOFlexIO_DataType
   @author     Thomas Radke
   @date       Mon 11 June 2001
   @desc
               Returns the FlexIO datatype for a given CCTK datatype
   @enddesc

   @var        cctk_type
   @vdesc      CCTK datatype
   @vtype      int
   @vio        in
   @endvar

   @returntype int
   @returndesc
               the appropriate FlexIO datatype for success, or -1 otherwise
   @endreturndesc
@@*/
int IOFlexIO_DataType (int cctk_type)
{
  int retval;


  switch (cctk_type)
  {
    case CCTK_VARIABLE_CHAR:   retval = FLEXIO_CHAR; break;
    case CCTK_VARIABLE_INT:    retval = FLEXIO_INT; break;
    case CCTK_VARIABLE_REAL:   retval = FLEXIO_REAL; break;
#ifdef CCTK_INT2
    case CCTK_VARIABLE_INT2:   retval = FLEXIO_INT2; break;
#endif
#ifdef CCTK_INT4
    case CCTK_VARIABLE_INT4:   retval = FLEXIO_INT4; break;
#endif
#ifdef CCTK_INT8
    case CCTK_VARIABLE_INT8:   retval = FLEXIO_INT8; break;
#endif
#ifdef CCTK_REAL4
    case CCTK_VARIABLE_REAL4:  retval = FLEXIO_REAL4; break;
#endif
#ifdef CCTK_REAL8
    case CCTK_VARIABLE_REAL8:  retval = FLEXIO_REAL8; break;
#endif
#ifdef CCTK_REAL16
    case CCTK_VARIABLE_REAL16: retval = FLEXIO_REAL16; break;
#endif

    default: CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                         "Unsupported CCTK variable datatype %d", cctk_type);
             retval = -1;
             break;
  }

  return (retval);
}


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
/*@@
   @routine    DumpGS
   @date       16 Apr 1999
   @author     Thomas Radke
   @desc
               Dumps a SCALAR type variable into a IEEEIO file
   @enddesc

   @calls      AddCommonAttributes

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      global index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to store
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar
   @var        flexio_type
   @vdesc      the IOFlexIO datatype to store
   @vtype      int
   @vio        in
   @endvar

   @returntype int
   @returndesc
               0 for success, or -1 if file handle is invalid
   @endreturndesc
@@*/
static int DumpGS (const cGH *GH, int vindex, int timelevel, IOFile file,
                   int flexio_type)
{
  ioGH *ioUtilGH;
  CCTK_INT gsz[] = {0, 0, 0};    /* not needed here ? */
  int dim[] = {1};               /* size of scalar */


  /* immediately return if file handle is invalid */
  if (! IOisValid (file))
  {
    return (-1);
  }

  /* get the handle for IOUtil extensions */
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  if (CCTK_MyProc (GH) != ioUtilGH->ioproc)
  {
    return (0);
  }

  /* first dump the data then add the attributes */
  FLEXIO_ERROR (IOwrite (file, flexio_type, 1, dim,
                         CCTK_VarDataPtrI (GH, timelevel, vindex)));
  AddCommonAttributes (GH, vindex, timelevel, gsz, file);

  return (0);
}


/*@@
   @routine    DumpGA
   @date       July 1998
   @author     Paul Walker
   @desc
               Dumps a grid array variable into a IEEEIO file.
   @enddesc

   @calls      GetDumpData

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        vindex
   @vdesc      global index of the variable to be dumped
   @vtype      int
   @vio        in
   @endvar
   @var        timelevel
   @vdesc      the timelevel to store
   @vtype      int
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar
   @var        info
   @vdesc      info structure describing
                 - the IOFlexIO datatype to store
                 - the size of an element of variable
                 - the MPI datatype to communicate
   @vtype      dumpInfo_t
   @vio        in
   @endvar

   @returntype int
   @returndesc
               0 for success
               or returncode of @seeroutine GetDumpData
   @endreturndesc
@@*/
static int DumpGA (const cGH *GH, int vindex, int timelevel, IOFile file,
                   dumpInfo_t *info)
{
  int myproc, nprocs, dim, free_outme, retval;
  pGH *pughGH;
  ioGH *ioUtilGH;
  CCTK_INT4 *geom;
  void *outme;
#ifdef CCTK_MPI
  int i, j, incoming, outgoing;
  int *local_size;
  MPI_Status ms;
  void *tmpd;
#endif
  DECLARE_CCTK_PARAMETERS


  /* get the handles for PUGH and IO extensions */
  pughGH = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  myproc = CCTK_MyProc (GH);
  nprocs = CCTK_nProcs (GH);

  /* get the dimension of the variable */
  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));

  /* allocate the geometry buffer */
  geom = (CCTK_INT4 *) malloc (3*dim * sizeof (CCTK_INT4));

  /* get the pointer to the data we want to output (includes downsampling) */
  retval = GetDumpData (GH, vindex, timelevel, &outme, &free_outme, geom,
                        info->element_size);
  if (retval < 0)
  {
    return (retval);
  }

  if (ioUtilGH->ioproc_every == 1)
  {
    EachProcDump (GH, vindex, timelevel, outme, geom, file, info->flexio_type);
  }
#ifdef CCTK_MPI
  else if (myproc != ioUtilGH->ioproc)
  {
    /*
     * Send the geometry and data from non-IO processors to the IO processors
     */
    outgoing = 1;
    for (i = 0; i < dim; i++)
    {
      outgoing *= geom[dim + i];
    }

    /* send geometry */
    CACTUS_MPI_ERROR (MPI_Send (geom, 3*dim, PUGH_MPI_INT4, ioUtilGH->ioproc,
                                2*myproc+IOTAGBASE+1, pughGH->PUGH_COMM_WORLD));
    /* send data */
    if (outgoing > 0)
    {
      CACTUS_MPI_ERROR (MPI_Send (outme, outgoing, info->mpi_type,
                                  ioUtilGH->ioproc, 2*myproc+IOTAGBASE,
                                  pughGH->PUGH_COMM_WORLD));
#ifdef IOFLEXIO_DEBUG
      printf ("Processor %d sent %d data points\n", myproc, outgoing);
#endif
    }
  }
  else if (myproc == ioUtilGH->ioproc)
  {
    if (ioUtilGH->unchunked)
    {
      local_size = (int *) malloc (dim * sizeof (int));

      /* first calculate the local size and reserve the chunk */
      if (ioUtilGH->ioproc_every >= nprocs)
      {
        /* one output file ... the local chunk size is the global size */
        for (i = 0; i < dim; i++)
        {
          local_size[i] = geom[2*dim + i];
        }
      }

      IOreserveChunk (file, info->flexio_type, dim, local_size);

      free (local_size);
    }

    /* write the data first then the attributes */
    ProcDump (file, GH, vindex, outme, geom, info->flexio_type);
    /* delay adding attributes for unchunked files
       until all chunks were written (otherwise attributes get lost !) */
    if (! ioUtilGH->unchunked)
    {
      AddCommonAttributes (GH, vindex, timelevel, &geom[2*dim], file);
    }

    /* dump data from all other processors */
    for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++)
    {
      /* receive geometry */
      CACTUS_MPI_ERROR (MPI_Recv (geom, 3*dim, PUGH_MPI_INT4, i,2*i+IOTAGBASE+1,
                                  pughGH->PUGH_COMM_WORLD, &ms));

      incoming = 1;
      for (j = 0; j < dim; j++)
      {
        incoming *= geom[dim + j];
      }

      /* receive data */
      if (incoming > 0)
      {
        tmpd = malloc (incoming * info->element_size);
        CACTUS_MPI_ERROR (MPI_Recv (tmpd, incoming, info->mpi_type, i,
                                    2*i+IOTAGBASE,pughGH->PUGH_COMM_WORLD,&ms));
      }
      else
      {
        tmpd = NULL;
      }

      /* write data */
      ProcDump (file, GH, vindex, tmpd, geom, info->flexio_type);

      if (tmpd)
      {
        free (tmpd);
      }
    } /* end loop over processors */

    /* now add the attributes for unchunked files */
    if (ioUtilGH->unchunked)
    {
      AddCommonAttributes (GH, vindex, timelevel, &geom[2*dim], file);
    }
  }

  /* wait for every processor to catch up */
  CCTK_Barrier (GH);
#endif /* MPI */

  if (free_outme)
  {
    free (outme);
  }

  free (geom);

  return (retval);
}


/*@@
   @routine    AddCommonAttributes
   @date       July 1998
   @author     Paul Walker
   @desc
               Add "Common" attributes, these are the GF name, the current date,
               simulation time, origin, bounding box, and gridspacings (both
               downsampled and evolution).
   @enddesc

   @calls

   @history
   @hdate      Wed Sep  2 10:15:22 1998
   @hauthor    Tom Goodale
   @hdesc      Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda changes.
   @hdate      Apr 16 1999
   @hauthor    Thomas Radke
   @hdesc      Added attributes groupname, grouptype, ntimelevels,
               and current timelevel to be stored
   @hdate      May 02 1999
   @hauthor    Thomas Radke
   @hdesc      Made chunked attribute nioprocs common so that it can be found by
               the recombiner even for unchunked datasets (eg. SCALARs)
   @hdate      May 05 1999
   @hauthor    Thomas Radke
   @hdesc      Added "unchunked" attribute to distinguish between
               chunked/unchunked output files
   @endhistory
@@*/
static void AddCommonAttributes (const cGH *GH, int vindex, int timelevel,
                                 CCTK_INT4 *gsz, IOFile file)
{
  int i, vdim;
  CCTK_REAL *dtmp;
  CCTK_INT4 itmp;
  char *name;
  ioGH *ioUtilGH;
  char coord_system_name[20];
  DECLARE_CCTK_PARAMETERS


  /* get the handle for IO extensions */
  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  /* get the dimension of the variable */
  vdim = CCTK_GroupDimFromVarI (vindex);

  name = CCTK_FullName (vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "name", FLEXIO_CHAR,
                                  strlen (name) + 1, name));
  free (name);

  name = CCTK_GroupNameFromVarI (vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "groupname", FLEXIO_CHAR,
                                  strlen (name) + 1, name));
  free (name);

  itmp = CCTK_GroupTypeFromVarI (vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "grouptype", FLEXIO_INT4, 1, &itmp));

  itmp = CCTK_NumTimeLevelsFromVarI (vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "ntimelevels", FLEXIO_INT4, 1, &itmp));

  itmp = timelevel;
  FLEXIO_ERROR (IOwriteAttribute (file, "timelevel", FLEXIO_INT4, 1, &itmp));

  FLEXIO_ERROR (IOwriteAttribute (file, "time", FLEXIO_REAL, 1,&GH->cctk_time));

  /* attributes describing the underlying grid
     These are only stored for CCTK_GF variables if they are associated
     with coordinates. */
  /* FIXME: This is hardcoded for cartesian coordinate systems.
            A better solution would be to be able to query the coordinate
            system which is associated with the variable. */
  sprintf (coord_system_name, "cart%dd", vdim);
  if (CCTK_GroupTypeFromVarI (vindex) == CCTK_GF &&
      CCTK_CoordSystemHandle (coord_system_name) >= 0)
  {
    dtmp = (CCTK_REAL *) malloc (3 * vdim * sizeof (CCTK_REAL));
    for (i = 0; i < vdim; i++)
    {
      CCTK_CoordRange (GH, &dtmp[i], &dtmp[i + vdim], i + 1,
                       NULL, coord_system_name);
      dtmp[i + 2*vdim] = GH->cctk_delta_space[i] * ioUtilGH->downsample[i];
    }

    FLEXIO_ERROR (IOwriteAttribute (file, "origin", FLEXIO_REAL, vdim, dtmp));
    FLEXIO_ERROR (IOwriteAttribute (file, "min_ext", FLEXIO_REAL, vdim, dtmp));
    FLEXIO_ERROR (IOwriteAttribute (file, "max_ext", FLEXIO_REAL, vdim,
                                    dtmp + vdim));
    FLEXIO_ERROR (IOwriteAttribute (file, "delta", FLEXIO_REAL, vdim,
                                    dtmp + 2*vdim));

    if (ioUtilGH->downsample[0] > 1 ||
        ioUtilGH->downsample[1] > 1 ||
        ioUtilGH->downsample[2] > 1)
    {
      FLEXIO_ERROR (IOwriteAttribute (file, "evolution_delta", FLEXIO_REAL,
                                      vdim, GH->cctk_delta_space));
    }

    free (dtmp);
  }

  FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, vdim, gsz));

  itmp = CCTK_nProcs (GH);
  FLEXIO_ERROR (IOwriteAttribute (file, "nprocs", FLEXIO_INT4, 1, &itmp));

  itmp = ioUtilGH->ioproc_every;
  FLEXIO_ERROR (IOwriteAttribute (file, "ioproc_every", FLEXIO_INT4, 1, &itmp));

  itmp = ioUtilGH->unchunked;
  FLEXIO_ERROR (IOwriteAttribute (file, "unchunked", FLEXIO_INT4, 1, &itmp));

  itmp = GH->cctk_iteration;
  FLEXIO_ERROR (IOwriteAttribute (file, "iteration", FLEXIO_INT4, 1, &itmp));
}


/*@@
  @routine     AddChunkAttributes
  @author      Paul Walker
  @date        Feb 1997
  @desc
               This routine adds chunk attributes to a data set. That is,
               is adds attributes to the chunks of data passed to the
               I/O processors.
  @enddesc
  @history
  @hauthor     Gabrielle Allen
  @hdate       Jul 10 1998
  @hdesc       Added the name of the grid function
  @hdate       Wed Sep  2 10:08:31 1998
  @hauthor     Tom Goodale
  @hdesc       Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda calls.
  @endhistory
@@*/
static void AddChunkAttributes (const cGH *GH, int vindex, CCTK_INT4 *geom,
                                IOFile file)
{
  int vdim;
  char *name;
  CCTK_INT4 itmp;


  /* there is nothing to do for a serial run */
  if (CCTK_nProcs (GH) == 1)
  {
    return;
  }

  /* get the dimension of the variable */
  vdim = CCTK_GroupDimFromVarI (vindex);

  FLEXIO_ERROR (IOwriteAttribute (file, "chunk_origin", FLEXIO_INT4, vdim,
                                  &geom[0]));
  FLEXIO_ERROR (IOwriteAttribute (file, "subchunk_lb", FLEXIO_INT4, vdim,
                                  &geom[0]));
  FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, vdim,
                                  &geom[2*vdim]));

  itmp = GH->cctk_iteration;
  FLEXIO_ERROR (IOwriteAttribute (file, "chunk_dataset", FLEXIO_INT4, 1,&itmp));

  name = CCTK_FullName (vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "name", FLEXIO_CHAR,
                                  strlen (name)+1, name));
  free (name);
}


/*@@
  @routine GetDumpData
  @author  Paul Walker
  @date    Feb 1997
  @desc
  Bounds and data to be output, takes into account downsampling
  @enddesc
  @history
  @hauthor Gabrielle Allen @hdate Oct 5 1998
  @hdesc Made into subroutine
  @endhistory
  @var     GH
  @vdesc   Identifies grid hierachy
  @vtype   const cGH *
  @vio     in
  @vcomment
  @endvar
  @var     vindex
  @vdesc   index of the variable to dump
  @vtype   int
  @vio     in
  @vcomment
  @endvar
  @var     timelevel
  @vdesc   timelevel of the variable to dump
  @vtype   int
  @vio     in
  @vcomment
  @endvar
  @var     outme
  @vdesc   data segment to output
  @vtype   void **
  @vio     out
  @vcomment This is just GF->data if there is no downsampling
  @endvar
  @var     free_outme
  @vdesc   Specifies whether or not to free the storage associated with outme
  @vtype   int *
  @vio     out
  @vcomment 1: free storage; 0: don't free storage
  @endvar
  @var     geom
  @vdesc   bounds, local size and global shape of the output
  @vtype   CCTK_INT4[3*dim]
  @vio     out
  @vcomment geom[0*dim..1*dim-1] lower bounds and geom[3],geom[4],geom[5]
            geom[1*dim..2*dim-1] local size of data to send
            geom[2*dim..3*dim-1] global shape
  @endvar
  @var     element_size
  @vdesc   the size of an element of variable
  @vtype   int
  @vio     in
  @endvar
@@*/
static int GetDumpData (const cGH *GH, int vindex, int timelevel, void **outme,
                        int *free_outme, CCTK_INT4 *geom, int element_size)
{
  DECLARE_CCTK_PARAMETERS
  int i, myproc, do_downsample, dim;
  ioGH *ioUtilGH;
  pGExtras *extras;
  CCTK_REAL4 *single_ptr;
  CCTK_REAL *real_ptr;
  CCTK_BYTE *char_ptr;
  CCTK_INT *int_ptr;
  void *data = CCTK_VarDataPtrI (GH, timelevel, vindex);


  /* to make the compiler happy */
  single_ptr = NULL;
  real_ptr = NULL;
  char_ptr = NULL;
  int_ptr = NULL;

  myproc = CCTK_MyProc (GH);

  /* get GH extensions for IO */
  ioUtilGH = (ioGH *) GH->extensions[CCTK_GHExtensionHandle ("IO")];

  /* get the pGExtras pointer as a shortcut */
  extras = ((pGA ***)PUGH_pGH (GH)->variables)[vindex][timelevel]->extras;

  /* get the dimension of the variable */
  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));

  do_downsample = 0;
  for (i = 0; i < dim; i++)
    do_downsample |= ioUtilGH->downsample[i] > 1;

  /* All the downsampling code is hard-coded for 3D data.
     For other-dimensional variables we print a warning
     and return the non-downsampled data. */
  if (do_downsample && dim != 3)
  {
    char *fullname = CCTK_FullName (vindex);

    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "No downsampling for %dD variable '%s' ! Downsampling is only "
                "supported for 3D variables.", dim, fullname);
    free (fullname);
    do_downsample = 0;
  }

  /* Simple case if no downsampling */
  if (! do_downsample)
  {

    if (ioUtilGH->out_single)
    {
      single_ptr = (CCTK_REAL4 *) malloc (extras->npoints*sizeof (CCTK_REAL4));

      for (i = 0; i < extras->npoints; i++)
      {
        single_ptr[i] = (CCTK_REAL4) ((CCTK_REAL *) data)[i];
      }

      *outme = single_ptr;
      *free_outme = 1;
    }
    else
    {
      *outme = data;
      *free_outme = 0;
    }

    for (i = 0; i < dim; i++)
    {
      geom[i + 0*dim] = extras->lb[myproc][i];;     /* the bounds */
      geom[i + 1*dim] = extras->lnsize[i];          /* the sizes */
      geom[i + 2*dim] = extras->nsize[i];           /* the global space */
    }

  } else {

    /* NOTE: the following downsampling code is hard-coded for 3D data */

    int start[3], end[3];
    int j, k, l;

    /* Downsampling code ... */
    for (i = 0; i < 3; i++) {
      geom[i + 6] = extras->nsize[i] / ioUtilGH->downsample[i];
      if (extras->nsize[i] % ioUtilGH->downsample[i])
        geom[i + 6]++;
    }

    if (verbose)
      CCTK_VInfo (CCTK_THORNSTRING, "Downsampled sizes (%d, %d, %d) -> "
                  "(%d, %d, %d)",
                  extras->nsize[0], extras->nsize[1], extras->nsize[2],
                  (int) geom[6], (int) geom[7], (int) geom[8]);

    /* Now figure out the local downsampling */
    /* The local starts are the lb modded into the downsample */
    for (i = 0; i < 3; i++) {
      geom[i]   = extras->lb[myproc][i] / ioUtilGH->downsample[i];
      start[i] = geom[i] * ioUtilGH->downsample[i];
      if (start[i] <
          extras->lb[myproc][i] + extras->ownership[PUGH_NO_STAGGER][0][i]) {
        start[i] += ioUtilGH->downsample[i];
        geom[i] ++;
      }
      end[i]   = ((extras->lb[myproc][i] +
                    extras->ownership[PUGH_NO_STAGGER][1][i] - 1) /
                   ioUtilGH->downsample[i]) * ioUtilGH->downsample[i];
      geom[i+3] = (end[i] - start[i]) / ioUtilGH->downsample[i] + 1;
    }

    if (verbose) {
      CCTK_VInfo (CCTK_THORNSTRING, "Downsample ranges (%d, %d, %d) -> "
                  "(%d, %d, %d)", start[0], start[1], start[2],
                  end[0], end[1], end[2]);
      CCTK_VInfo (CCTK_THORNSTRING, "Local size/bound  (%d, %d, %d)  "
                  "(%d, %d, %d)",
                  (int) geom[3], (int) geom[4], (int) geom[5],
                  (int) geom[0], (int) geom[1], (int) geom[2]);
    }

    /* compute local ranges */
    for (i = 0; i < 3; i++) {
      start[i] -= extras->lb[myproc][i];
      end[i] -= extras->lb[myproc][i];
    }

    *outme = malloc (geom[3] * geom[4] * geom[5] * element_size);
    *free_outme = 1;

    /* I hate it to repeat the loops for each case label
       but that way produces much more efficient code */
    l = 0;
    switch (CCTK_VarTypeI (vindex)) {
      case CCTK_VARIABLE_CHAR:
        char_ptr = (CCTK_BYTE *) *outme;
        for (k = start[2]; k <= end[2]; k += ioUtilGH->downsample[2])
          for (j = start[1]; j <= end[1]; j += ioUtilGH->downsample[1])
            for (i = start[0]; i <= end[0]; i += ioUtilGH->downsample[0])
              char_ptr[l++] = ((CCTK_BYTE *) data)[DATINDEX (extras, i, j, k)];
        break;

      case CCTK_VARIABLE_INT:
        int_ptr = (CCTK_INT *) *outme;
        for (k = start[2]; k <= end[2]; k += ioUtilGH->downsample[2])
          for (j = start[1]; j <= end[1]; j += ioUtilGH->downsample[1])
            for (i = start[0]; i <= end[0]; i += ioUtilGH->downsample[0])
              int_ptr[l++] = ((CCTK_INT *) data)[DATINDEX (extras, i, j, k)];
        break;

      case CCTK_VARIABLE_REAL:
        if (ioUtilGH->out_single)
          single_ptr = (CCTK_REAL4 *) *outme;
        else
          real_ptr = (CCTK_REAL *) *outme;
        for (k = start[2]; k <= end[2]; k += ioUtilGH->downsample[2])
          for (j = start[1]; j <= end[1]; j += ioUtilGH->downsample[1])
            for (i = start[0]; i <= end[0]; i += ioUtilGH->downsample[0])
              if (ioUtilGH->out_single)
                single_ptr[l++] = (CCTK_REAL4)
                      (((CCTK_REAL *) data)[DATINDEX (extras, i, j, k)]);
              else
                real_ptr[l++] =
                      ((CCTK_REAL *) data)[DATINDEX (extras, i, j, k)];
        break;

#if 0
/* FIXME: Don't know how to support COMPLEX type too !! */
      case CCTK_VARIABLE_COMPLEX:
        if (ioUtilGH->out_single)
          single_ptr = (CCTK_REAL4 *) *outme;
        else
          cmplx_ptr = (CCTK_COMPLEX *) *outme;
        for (k = start[2]; k <= end[2]; k += ioUtilGH->downsample[2])
          for (j = start[1]; j <= end[1]; j += ioUtilGH->downsample[1])
            for (i = start[0]; i <= end[0]; i += ioUtilGH->downsample[0])
              if (ioUtilGH->out_single) {
                single_ptr[l++] = (CCTK_REAL4)
                       (((CCTK_REAL *) data)[0 + DATINDEX (extras, i, j, k)]);
                single_ptr[l++] = (CCTK_REAL4)
                       (((CCTK_REAL *) data)[1 + DATINDEX (extras, i, j, k)]);
              } else
                cmplx_ptr[l++] =
                       ((CCTK_COMPLEX *) data)[DATINDEX (extras, i, j, k)];
        break;
#endif

      default:
        CCTK_WARN (1, "Unsupported variable type in GetDumpData");
        return (-1);
    }
  }

#ifdef IOFLEXIO_DEBUG
  printf ("Lower bound: %d", (int) geom[0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom[i]);
  printf ("\n");
  printf ("Chunk size : %d", (int) geom[1*dim + 0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom[1*dim + i]);
  printf ("\n");
  printf ("Global size: %d", (int) geom[2*dim + 0]);
  for (i = 1; i < dim; i++)
    printf (" %d", (int) geom[2*dim + i]);
  printf ("\n");
#endif

  return (0);
}


/*@@
  @routine     EachProcDump
  @author      Paul Walker
  @date        Feb 1997
  @desc
               Dump data from each processor
  @enddesc
  @history
  @hauthor     Gabrielle Allen
  @hdate       Oct 5 1998
  @hdesc       Made into subroutine
  @endhistory
@@*/
static void EachProcDump (const cGH *GH, int vindex, int timelevel, void *outme,
                          CCTK_INT4 *geom, IOFile file, int flexio_type)
{
  int i, dim, *chunk_dims;


  /* get the dimension of the variable */
  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));

  chunk_dims = (int *) malloc (dim * sizeof (int));

  /* set up the local shape */
  for (i = 0; i < dim; i++)
  {
    chunk_dims[i] = geom[dim + i];
  }

  /* dump the data */
  FLEXIO_ERROR (IOwrite (file, flexio_type, dim, chunk_dims, outme));

  /* add attributes for global space */
  AddCommonAttributes (GH, vindex, timelevel, &geom[2*dim], file);

  /* add chunk attributes */
  AddChunkAttributes (GH, vindex, geom, file);

  free (chunk_dims);
}


/*@@
  @routine     ProcDump
  @author      Paul Walker
  @date        Feb 1997
  @desc
               Dump data
  @enddesc
  @history
  @hauthor     Gabrielle Allen
  @hdate       Oct 5 1998
  @hdesc       Made into subroutine
  @endhistory
@@*/
#ifdef CCTK_MPI
static void ProcDump (IOFile file, const cGH *GH, int vindex, void *outme,
                      CCTK_INT4 *geom, int flexio_type)
{
  ioGH *ioUtilGH;
  int i, dim, *chunk_dims, *chunk_origin;


  ioUtilGH = (ioGH *) CCTK_GHExtension (GH, "IO");

  /* get the dimension of the variable */
  dim = CCTK_GroupDimI (CCTK_GroupIndexFromVarI (vindex));

  chunk_origin = (int *) malloc (2*dim * sizeof (int));
  chunk_dims = chunk_origin + dim;
  for (i = 0; i < dim; i++)
  {
    chunk_origin[i] = geom[i];
    chunk_dims[i] = geom[dim + i];
  }

  if (ioUtilGH->unchunked)
  {
    /* unchunked data */
    FLEXIO_ERROR (IOwriteChunk (file, chunk_dims, chunk_origin, outme));
  }
  else
  {
    /* uhunked data */
    FLEXIO_ERROR (IOwrite (file, flexio_type, dim, chunk_dims, outme));

    /* write chunk attributes */
    AddChunkAttributes (GH, vindex, geom, file);
  }

  free (chunk_origin);
}
#endif