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 <stdlib.h>
#include <string.h>

#include "cctk.h"
#include "util_Table.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 */

/* tag base for MPI messages (this may break on more than 2000 processors) */
#define MPITAGBASE 20000


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
static int WriteGS (const cGH *GH, const ioRequest *request, IOFile file);
static int WriteGA (const cGH *GH, const ioRequest *request, IOFile file);
static void EachProcDump (const cGH *GH, const ioRequest *request, void *hdata,
                          IOFile file);
static void AddChunkAttributes (const cGH *GH, const ioRequest *request,
                                IOFile file);
static void AddCommonAttributes (const cGH *GH, const ioRequest *request,
                                 IOFile file);
#ifdef CCTK_MPI
static void ProcDump (const cGH *GH, const ioRequest *request, void *hdata,
                      IOFile file);
#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        request
   @vdesc      reference to the I/O request description
   @vtype      const ioRequest *
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar

   @returntype int
   @returndesc
               return code of either @seeroutine WriteGS or @seeroutine WriteGA
               or -1 if variable has unknown type
   @endreturndesc
@@*/
int IOFlexIO_DumpVar (const cGH *GH, const ioRequest *request, IOFile file)
{
  int gtype, retval;


  /* branch to the appropriate dump routine */
  gtype = CCTK_GroupTypeFromVarI (request->vindex);
  if (gtype == CCTK_SCALAR)
  {
    retval = WriteGS (GH, request, file);
  }
  else if (gtype == CCTK_ARRAY || gtype == CCTK_GF)
  {
    retval = WriteGA (GH, request, file);
  }
  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_INT1
    case CCTK_VARIABLE_INT1:   retval = FLEXIO_INT1; break;
#endif
#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 (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                         "Unsupported CCTK variable datatype %d", cctk_type);
             retval = -1;
             break;
  }

  return (retval);
}


/********************************************************************
 ********************    Internal Routines   ************************
 ********************************************************************/
/*@@
   @routine    WriteGS
   @date       16 Apr 1999
   @author     Thomas Radke
   @desc
               Writes a CCTK_SCALAR type variable into a IEEEIO file.
               Only the value from processor 0 is written here.
               If other processors have different values for this scalar
               variable a warning will be issued.
   @enddesc

   @calls      AddCommonAttributes

   @var        GH
   @vdesc      Pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        request
   @vdesc      reference to the I/O request description
   @vtype      const ioRequest *
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar

   @returntype int
   @returndesc
               0 for success, or -1 if file handle is invalid
   @endreturndesc
@@*/
static int WriteGS (const cGH *GH, const ioRequest *request, IOFile file)
{
  int i, myproc, nprocs, hdatatypesize, retval;
  char *buffer, *fullname;
  /* const */ int dim = 1;
  const ioGH *ioUtilGH;


  ioUtilGH = CCTK_GHExtension (GH, "IO");
  myproc = CCTK_MyProc (GH);
  nprocs = CCTK_nProcs (GH);
  fullname = CCTK_FullName (request->vindex);

  hdatatypesize = CCTK_VarTypeSize (request->hdatatype);
  buffer = calloc (nprocs, hdatatypesize);
  memcpy (buffer + myproc*hdatatypesize,
          CCTK_VarDataPtrI (GH, request->timelevel, request->vindex),
          hdatatypesize);

  if (nprocs > 1)
  {
    i = CCTK_ReductionHandle ("sum");
    if (i >= 0)
    {
      i = CCTK_ReduceArray (GH, -1, i, nprocs, request->hdatatype,
                            buffer, 1, 1, request->hdatatype, nprocs, buffer);
    }
    if (i < 0)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "WriteGS: Cannot check whether values on different "
                  "processors are the same for grid scalar '%s'", fullname);

      /* copy this processor's value to the start of buffer */
      memcpy (buffer, buffer + myproc*hdatatypesize, hdatatypesize);
    }
    else
    {
      retval = 0;
      for (i = 1; i < nprocs; i++)
      {
        retval |= memcmp (buffer, buffer + i*hdatatypesize, hdatatypesize);
      }
      if (retval)
      {
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "WriteGS: value of grid scalar variable '%s' (timelevel %d)"
                    " differs between processors, only value from processor 0 "
                    "will be written", fullname, request->timelevel);
      }
    }
  }

  /* only I/O processors write data */
  if (myproc != ioUtilGH->ioproc)
  {
    retval = 0;
  }
  else if (! IOisValid (file))
  {
    retval = -1;
  }
  else
  {
    /* first dump the data then add the attributes */
    FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype), dim,
                           &dim, buffer));
    AddCommonAttributes (GH, request, file);
    retval = 0;
  }

  free (buffer);
  free (fullname);

  return (retval);
}


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

   @calls      Hyperslab_DefineLocalMappingByIndex
               Hyperslab_FreeMapping
               Hyperslab_Get
               WriteData

   @var        GH
   @vdesc      pointer to CCTK grid hierarchy
   @vtype      const cGH *
   @vio        in
   @endvar
   @var        request
   @vdesc      reference to the I/O request description
   @vtype      const ioRequest *
   @vio        in
   @endvar
   @var        file
   @vdesc      the IEEEIO file to dump to
   @vtype      IOFile
   @vio        in
   @endvar

   @returntype int
   @returndesc
                0 for success, or<BR>
               -1 if hyperslab mapping couldn't be defined, or<BR>
               -2 if hyperslab couldn't be extracted
   @endreturndesc
@@*/
static int WriteGA (const cGH *GH, const ioRequest *request, IOFile file)
{
  int i, mapping, hdatasize, table, retval;
  void *hdata;
  char *fullname;
  const ioGH *ioUtilGH;
#ifdef CCTK_MPI
  int myproc, nprocs;
  int *hsize_chunk;
  void *tmpd;
  MPI_Comm comm;
  MPI_Status ms;
  MPI_Datatype mpitype;
#endif
  DECLARE_CCTK_PARAMETERS


  /* define the hyperslab mapping */
  table = -1;
  if (request->with_ghostzones)
  {
    table = Util_TableCreateFromString ("with_ghostzones = 1");
    if (table < 0)
    {
      CCTK_WARN (1, "Failed to hyperslab parameter create table from string");
    }
  }
  mapping = Hyperslab_DefineLocalMappingByIndex (GH, request->vindex,
                                                 request->hdim,
                                                 request->direction,
                                                 request->origin,
                                                 request->extent,
                                                 request->downsample,
                                                 table, NULL,
                                                 request->hsize_chunk,
                                                 request->hsize,
                                                 request->hoffset);
  if (table >= 0)
  {
    Util_TableDestroy (table);
  }

  if (mapping < 0)
  {
    fullname = CCTK_FullName (request->vindex);
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Failed to define hyperslab mapping for variable '%s'",
                fullname);
    free (fullname);
    return (-1);
  }

  /* calculate the size of the hyperslab */
  request->hsize_chunk[request->hdim] = 1;
  for (i = 0; i < request->hdim; i++)
  {
    request->hsize_chunk[request->hdim] *= request->hsize_chunk[i];
  }

  /* get the hyperslab */
  hdatasize = CCTK_VarTypeSize (request->hdatatype);
  hdata = request->hsize_chunk[request->hdim] > 0 ?
          malloc (request->hsize_chunk[request->hdim] * hdatasize) : NULL;
  retval = Hyperslab_Get (GH, mapping, -1, request->vindex, request->timelevel,
                          request->hdatatype, hdata);

  /* release the mapping structure */
  Hyperslab_FreeMapping (mapping);

  if (retval)
  {
    fullname = CCTK_FullName (request->vindex);
    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "Failed to extract hyperslab for variable '%s'", fullname);
    free (fullname);
    if (hdata)
    {
      free (hdata);
    }
    return (-2);
  }

  /* dump data held on I/O processor */
  ioUtilGH = CCTK_GHExtension (GH, "IO");
  if (ioUtilGH->ioproc_every == 1)
  {
    EachProcDump (GH, request, hdata, file);
  }
#ifdef CCTK_MPI
  else
  {
    myproc = CCTK_MyProc (GH);
    nprocs = CCTK_nProcs (GH);
    comm = PUGH_pGH (GH)->PUGH_COMM_WORLD;
    mpitype = PUGH_MPIDataType (PUGH_pGH (GH), request->hdatatype);

    if (myproc == ioUtilGH->ioproc)
    {
      if (ioUtilGH->unchunked)
      {
        /* need to copy CCTK_INT[] to int[] */
        hsize_chunk = malloc (request->hdim * sizeof (int));
        for (i = 0; i < request->hdim; i++)
        {
          hsize_chunk[i] = request->hsize[i];
        }

        IOreserveChunk (file, IOFlexIO_DataType (request->hdatatype),
                        request->hdim, hsize_chunk);

        free (hsize_chunk);
      }

      /* write the data first then the attributes */
      ProcDump (GH, request, hdata, file);
      /* delay adding attributes for unchunked files
         until all chunks were written (otherwise attributes get lost !) */
      if (! ioUtilGH->unchunked)
      {
        AddCommonAttributes (GH, request, file);
      }

      /* dump data from all other processors */
      for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++)
      {
        /* receive geometry (this assumes the geometry arrays
           to be contiguous starting at hoffset */
        CACTUS_MPI_ERROR (MPI_Recv (request->hoffset, 3*request->hdim + 1,
                                    PUGH_MPI_INT, i, 2*i + MPITAGBASE + 1,
                                    comm, &ms));

        /* receive data */
        tmpd = NULL;
        if (request->hsize_chunk[request->hdim] > 0)
        {
          tmpd = malloc (request->hsize_chunk[request->hdim] * hdatasize);
          CACTUS_MPI_ERROR (MPI_Recv (tmpd, request->hsize_chunk[request->hdim],
                                      mpitype, i, 2*i + MPITAGBASE, comm, &ms));
        }

        /* write data */
        ProcDump (GH, request, tmpd, file);

        if (tmpd)
        {
          free (tmpd);
        }

      } /* end loop over processors */

      /* now add the attributes for unchunked files */
      if (ioUtilGH->unchunked)
      {
        AddCommonAttributes (GH, request, file);
      }
    }
    else
    {
      /* send geometry (this assumes the geometry arrays to be contiguous
         starting at hoffset) */
      CACTUS_MPI_ERROR (MPI_Send (request->hoffset, 3*request->hdim + 1,
                                  PUGH_MPI_INT, ioUtilGH->ioproc,
                                  2*myproc + MPITAGBASE + 1, comm));
      /* send data */
      if (request->hsize_chunk[request->hdim] > 0)
      {
        CACTUS_MPI_ERROR (MPI_Send (hdata, request->hsize_chunk[request->hdim],
                                    mpitype, ioUtilGH->ioproc,
                                    2*myproc + MPITAGBASE, comm));
#ifdef IOFLEXIO_DEBUG
        printf ("Processor %d sent %d data points\n",
                myproc, request->hsize_chunk[request->hdim]);
#endif
      }
    }
  }

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

  /* free allocated resources */
  if (hdata)
  {
    free (hdata);
  }

  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, const ioRequest *request,
                                 IOFile file)
{
  int i;
  CCTK_REAL *dtmp;
  CCTK_INT4 *itmp;
  char *name;
  ioGH *ioUtilGH;
  char coord_system_name[20];
  DECLARE_CCTK_PARAMETERS


  /* allocate at least one CCTK_INT4 if hdim is 0 */
  itmp = malloc ((request->hdim+1) * sizeof (CCTK_INT4));

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

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

  itmp[0] = CCTK_GroupTypeFromVarI (request->vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "grouptype", FLEXIO_INT4, 1, itmp));

  itmp[0] = CCTK_MaxTimeLevelsVI (request->vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "ntimelevels", FLEXIO_INT4, 1, itmp));

  itmp[0] = request->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. */
  ioUtilGH = CCTK_GHExtension (GH, "IO");
  sprintf (coord_system_name, "cart%dd", request->hdim);
  if (CCTK_GroupTypeFromVarI (request->vindex) == CCTK_GF &&
      CCTK_CoordSystemHandle (coord_system_name) >= 0)
  {
    dtmp = malloc (3 * request->hdim * sizeof (CCTK_REAL));
    for (i = 0; i < request->hdim; i++)
    {
      CCTK_CoordRange (GH, &dtmp[i], &dtmp[i + request->hdim], i + 1,
                       NULL, coord_system_name);
      dtmp[i + 2*request->hdim] = GH->cctk_delta_space[i] *
                                  ioUtilGH->downsample[i];
    }

    FLEXIO_ERROR (IOwriteAttribute (file, "origin", FLEXIO_REAL, request->hdim,
                                    dtmp));
    FLEXIO_ERROR (IOwriteAttribute (file, "min_ext", FLEXIO_REAL, request->hdim,
                                    dtmp));
    FLEXIO_ERROR (IOwriteAttribute (file, "max_ext", FLEXIO_REAL, request->hdim,
                                    dtmp + request->hdim));
    FLEXIO_ERROR (IOwriteAttribute (file, "delta", FLEXIO_REAL, request->hdim,
                                    dtmp + 2*request->hdim));

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

    free (dtmp);
  }

  if (request->hdim)
  {
    for (i = 0; i < request->hdim; i++)
    {
      itmp[i] = request->hsize[i];
    }
    FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4,
                                    request->hdim, itmp));
  }
  else
  {
    itmp[0] = 0;
    FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4, 1, itmp));
  }

  itmp[0] = CCTK_nProcs (GH);
  FLEXIO_ERROR (IOwriteAttribute (file, "nprocs", FLEXIO_INT4, 1, itmp));

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

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

  itmp[0] = GH->cctk_iteration;
  FLEXIO_ERROR (IOwriteAttribute (file, "iteration", FLEXIO_INT4, 1, itmp));

  free (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, const ioRequest *request,
                                IOFile file)
{
  int i;
  char *fullname;
  CCTK_INT4 *itmp;


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

  itmp = malloc (2 * request->hdim * sizeof (CCTK_INT4));

  /* copy from CCTK_INT[] to CCTK_INT4[] */
  for (i = 0; i < request->hdim; i++)
  {
    itmp[0*request->hdim + i] = request->hoffset[i];
    itmp[1*request->hdim + i] = request->hsize[i];
  }

  FLEXIO_ERROR (IOwriteAttribute (file, "chunk_origin", FLEXIO_INT4,
                                  request->hdim, &itmp[0]));
  FLEXIO_ERROR (IOwriteAttribute (file, "subchunk_lb", FLEXIO_INT4,
                                  request->hdim, &itmp[0]));
  FLEXIO_ERROR (IOwriteAttribute (file, "global_size", FLEXIO_INT4,
                                  request->hdim, &itmp[request->hdim]));

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

  fullname = CCTK_FullName (request->vindex);
  FLEXIO_ERROR (IOwriteAttribute (file, "name", FLEXIO_CHAR,
                                  strlen (fullname)+1, fullname));
  free (fullname);
  free (itmp);
}


/*@@
  @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, const ioRequest *request, void *hdata,
                          IOFile file)
{
  int i, *chunk_dims;


  /* copy from CCTK_INT[] to int[] */
  chunk_dims = malloc (request->hdim * sizeof (int));
  for (i = 0; i < request->hdim; i++)
  {
    chunk_dims[i] = request->hsize_chunk[i];
  }

  /* dump the data */
  FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype),
                         request->hdim, chunk_dims, hdata));

  /* add attributes for global space */
  AddCommonAttributes (GH, request, file);

  /* add chunk attributes */
  AddChunkAttributes (GH, request, 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 (const cGH *GH, const ioRequest *request, void *hdata,
                      IOFile file)
{
  ioGH *ioUtilGH;
  int i, *chunk_dims, *chunk_origin;


  /* copy from CCTK_INT[] to int[] */
  chunk_origin = malloc (2*request->hdim * sizeof (int));
  chunk_dims = chunk_origin + request->hdim;
  for (i = 0; i < request->hdim; i++)
  {
    chunk_origin[i] = request->hoffset[i];
    chunk_dims[i] = request->hsize_chunk[i];
  }

  ioUtilGH = CCTK_GHExtension (GH, "IO");
  if (ioUtilGH->unchunked)
  {
    /* unchunked data */
    FLEXIO_ERROR (IOwriteChunk (file, chunk_dims, chunk_origin, hdata));
  }
  else
  {
    /* chunked data */
    FLEXIO_ERROR (IOwrite (file, IOFlexIO_DataType (request->hdatatype),
                           request->hdim, chunk_dims, hdata));

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

  free (chunk_origin);
}
#endif