path: root/src/Write3D.c

/*@@
   @file      Write3D.c
   @date      Fri Feb 21 15:27:03 1997
   @author    
   @desc 
   A horrible blocking write and not-quite-so-horrible
   slicer for xgraph output. Basically MPI based IO
   bites the big one, and this is MPI based IO written
   quickly and in the obvious inefficient fashion.
   <p>
   Oh, not it also contains zero-d I/O and a closer also.
   <p>
   This routine can write DFSD style HDF (not the default)
   or Johns IEEEIO 3D files (the default). For 0 and 1D
   I/O it writes xgraph files.
   @enddesc 
   @history
   @hauthor Gabrielle Allen @hdate 26 Sep 1998
   @hdesc Changed routine name (WriteGF->Write3D), Lots of tidying up and 
          renaming of parameters
   @hauthor Gabrielle Allen @hdate 17 Oct 1998
   @hdesc Implemented new parameters IO_verbose,IO_parameter,IO_structures,IO_scalars
   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
   @hdesc Changed names ready for thorn_IO
   @hauthor Thomas Radke @hdate 16 Mar 1999
   @hdesc Converted to Cactus 4.0
   @hauthor Thomas Radke @hdate 17 Mar 1999
   @hdesc Changed naming: IEEEIO -> IOFlexIO
   @hendhistory
   @version $Id$
 @@*/


#include <stdio.h>
#include <stdlib.h>
#include <assert.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 char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusPUGHIO_IOFlexIO_Write3D_c)


/* Function prototypes from other files in this thorn */
void IOFlexIOi_DumpGHExtensions (cGH *GH, IOFile iof);
void IOFlexIOi_DumpParameters (cGH *GH, IOFile iof);



/* local function prototypes */
static IEEEfile_3D_t *IOFlexIO_Get3Dfile (cGH *GH, const char *alias,
                                          int *isNewFile);

/*@@
  @routine IOFlexIO_Write3D
  @author  Paul Walker 
  @date    Feb 1997
  @desc
  This routine does 3D IO for a grid function which has been
  distributed with MPI.
  <p>
  3D MPI IO! Wow, you may say, that sounds pretty clever. Well it isn't
  because <b><i>The method I use sucks rocks</i></b>. That is.
  I spam the data from all processors onto processor 0. Luckily
  it doesn't just put it all in one memory chunk, but recieves it
  bit by bit and writes it out chunkwise, using IEEEIO. Alas, 
  while this happens the other processors <i>sit in a barrier 
  with their metaphoric thumbs up their asses</i>. 
  <p>
  <b>Even Worse</b> it is a blocking send which must come in order
  from the other processors.
  <p> 
  This is inefficient to say the least. I/O costs can go up to 70-80%
  using this technique on the T3E. So there are also a couple of ways to
  get around it.
  <p>
  The easiest option is one-file-per-processor which is
  good on some machines, (eg T3E) and bad on others (eg O2K, I think).
  So you can activate this mode (or de-activate it on the T3E)
  by setting the parameter onefileperproc to "yes" or "no".
  This will create a bunch of files in their own subdirectory,
  and they can be smooshed together using cactus/lib/etc/Recombiner.C
  <p>
  Also, you can downsample output. You can also output one file per
  slice. Various combinations of these work or don't, depending.
  <p>
  Finally notice that @seeroutine main guarantees that this is never
  called on a 1- or 2-D grid.
  @enddesc
  @calledby IOFlexIO_Output3DVarAs, IOFlexIO_TriggerOutput3D
  @history 
  @hdate Wed Sep  2 10:11:41 1998 @hauthor Tom Goodale
  @hdesc Merged in Szu-Wen's Panda calls.
  @hdate Sep 25 1998 1998 @hauthor Gabrielle Allen
  @hdesc Split into subroutines
  @endhistory 
  @var     GH
  @vdesc   Pointer to CCTK GH
  @vtype   cGH
  @vio     in
  @vcomment
  @endvar
  @var     index
  @vdesc   index of variable to output
  @vtype   int
  @vio     in
  @vcomment
  @endvar
  @var     alias
  @vdesc   alias name of variable to output
  @vtype   const char *
  @vio     in
  @vcomment
  @endvar
  @@*/


void IOFlexIO_Write3D (cGH *GH, int index, const char *alias)
{
  DECLARE_CCTK_PARAMETERS
  int timelevel;
  ioGH *ioUtilGH;
  int isNewFile;
  IEEEfile_3D_t *IEEEfile_3D;


  /* first, check if variable has storage assigned */
  if (! CCTK_QueryGroupStorageI (GH, CCTK_GroupIndexFromVarI (index))) {
    char *fullname = CCTK_FullName (index);

    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                "No IOFlexIO output for '%s' (no storage)", fullname);
    free (fullname);
    return;
  }

/*** FIXME: IEEEIO doesn't provide a COMPLEX datatype
            so somehow CCTK_COMPLEX has to be mapped onto REALs or so.
            We have to check this already here because if an IEEEIO file
            is created and nothing is written to, it crashes at re-opening. ***/
  if (CCTK_VarTypeI (index) == CCTK_VARIABLE_COMPLEX) {
    char *fullname = CCTK_FullName (index);

    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                "IOFlexIO output for complex variable '%s' not yet "
                "supported", fullname);
    free (fullname);
    return;
  }

  /* Get the handle for IOUtil and IOFlexIO extensions */
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];

  if (verbose)
    CCTK_VInfo (CCTK_THORNSTRING, "IOFlexIO output of variable '%s', "
                "alias '%s'", CCTK_VarName (index), alias);

  /* Open the file on IO processors */
  if (CCTK_MyProc (GH) != ioUtilGH->ioproc)
  {
    IEEEfile_3D = NULL;
  }
  else
  {
    /* Get the filename and descriptor for output */
    IEEEfile_3D = IOFlexIO_Get3Dfile (GH, alias, &isNewFile);

    /* nothing to do for an already opened file
       which is kept open all the time */
    if (isNewFile || reuse_filehandles)
    {

      if (verbose)
        CCTK_VInfo (CCTK_THORNSTRING, "%s IEEEIO output file '%s'",
                    isNewFile ? "Opening" : "Resuming", IEEEfile_3D->filename);

      /* First time through or one file per slice; open anew */
      if (isNewFile) {
        IEEEfile_3D->iofile = IEEEopen (IEEEfile_3D->filename, "w");
        if (! IOisValid (IEEEfile_3D->iofile))
          CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, "Couldn't "
                      "create IEEEIO output file '%s'", IEEEfile_3D->filename);
      }
      else
      {
        /* resume the file (reopen descriptor)
           This allows descriptor reuse without requiring expensive destruction
           and reallocation of the associated datastructures */
        CACTUS_IEEEIO_ERROR (IOresume (IEEEfile_3D->iofile));
      }

      /* Turn on buffer cache (with default size chosen by IEEEIO lib) */
      if (IEEEfile_3D->iofile > 0)
        IEEEbufferOn (IEEEfile_3D->iofile, 0);
    }
    if (IEEEfile_3D->iofile <= 0)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                  "Invalid file descriptor for IEEEIO ouput file '%s'",
                  IEEEfile_3D->filename);
      return;
    }
  }

  /* get the current timelevel */
  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
  if (timelevel > 0)
    timelevel--;

  /* output the data */
  IOFlexIO_DumpVar (GH, index, timelevel,
                    IEEEfile_3D ? IEEEfile_3D->iofile : (IOFile) -1);

  /* output the GH extensions and parameters */
  /* Don't do this if the output is being byte compared, as the
   * thorn lists might be different.
   */
  if (IEEEfile_3D) {
    /* output parameters necessary for filereader datafiles */
    if (isNewFile) {
      IOFlexIOi_DumpGHExtensions (GH, IEEEfile_3D->iofile);
      if (out3D_parameters)
        IOFlexIOi_DumpParameters (GH, IEEEfile_3D->iofile);
    }

    /* close the file */
    if (out3D_septimefiles || reuse_filehandles)
    {
      if (verbose)
        CCTK_VInfo (CCTK_THORNSTRING, "%s IEEEIO output file '%s'",
                    out3D_septimefiles ? "Closing" : "Pausing",
                    IEEEfile_3D->filename);

      IEEEbufferOff (IEEEfile_3D->iofile); 
      if (out3D_septimefiles)
      {
        CACTUS_IEEEIO_ERROR (IOclose (IEEEfile_3D->iofile));
        free (IEEEfile_3D);
      }
      else
      {
        CACTUS_IEEEIO_ERROR (IOpause (IEEEfile_3D->iofile));
      }
    }
  }
}


/**************************** local routines ****************************/


/*@@
  @routine IOFlexIO_Get3Dfile
  @author  Paul Walker
  @date    Feb 1997
  @desc
  Generates the filename for output
  @enddesc
  @history 
  @hauthor Gabrielle Allen
  @hdate Sep 1998
  @hdesc Split into subroutine
  @endhistory 
  @@*/

static IEEEfile_3D_t *IOFlexIO_Get3Dfile (cGH *GH, const char *alias,
                                          int *isNewFile)
{
  DECLARE_CCTK_PARAMETERS
  ioGH *ioUtilGH;       /* handle for IOUtil extensions */
  pGH *pughGH;          /* handle for PUGH extensions */
  flexioGH *myGH;       /* handle for IOFlexIO extensions */
  int nprocs;
  int myproc;
  char extra [256];     /* Extra stuff in fname based on mode */
  char extradir [256];  /* Extra stuff for an output dir */
  char *outputdir;      /* the output directory to be created */
  IEEEfile_3D_t *IEEEfile_3D;


  /* get the GH extensions for PUGH, IOUtil, and IOFlexIO */
  pughGH = PUGH_pGH (GH);
  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  myGH = (flexioGH *) GH->extensions [CCTK_GHExtensionHandle ("IOFlexIO")];

  IEEEfile_3D = (IEEEfile_3D_t *) GetNamedData (myGH->fileList_3D, alias);
  if (IEEEfile_3D != NULL) {
    /* set flag to indicate that file should be opened in append mode */
    *isNewFile = 0;
    return (IEEEfile_3D);
  }

  extra [0] = '\0';
  extradir [0] = '\0';
  
  nprocs = CCTK_nProcs (GH);
  myproc = CCTK_MyProc (GH);

#if        0
  /* Say if we have a chunked data in the output file(s) */
  if (! ioUtilGH->unchunked && nprocs > 1)
    sprintf (extra, "%s_chunked", extra);
#endif
  
  if (out3D_septimefiles) {
    char *tmp = extra;

    sprintf (extra, "%s.time_%7.3f", extra, GH->cctk_time);
  
    /* And be sure to replace any spaces in the filename with an _ */
    do
      if (*tmp == ' ')
        *tmp = '_';
    while (*++tmp);
  }
 
  /* OUTPUT ONE FILE FOR EACH N PROCESSORS 
   *  ------------------------------------- 
   *  
   *  If only one output file, the single file for each GF is written 
   *  in the normal output dir (that is extradir = ""). Otherwise
   *  a directory is created for each output GF to hold the multiple
   *  file 
   */

  if (ioUtilGH->ioproc_every < nprocs) {
    
    /* Add the output processor number to the extra string */
    sprintf (extra, "%s.file_%d", extra, myproc / ioUtilGH->ioproc_every);
    
    /* If necessary create the output directory */
    if (myproc == 0) {
      int result;
      
      outputdir = (char *) malloc (strlen (myGH->outdir3D) +
                                   strlen (alias) + 5);
      sprintf (outputdir, "%s/%s_3d", myGH->outdir3D, alias);
      
      result = CCTK_CreateDirectory (0755, outputdir);
      if (result < 0)
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Problem creating 3D output directory '%s'", outputdir);
      if (result > 0)
        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "3D output directory '%s' already exists", outputdir);

      free (outputdir);
    }
    
#ifdef CCTK_MPI
    /* Wait for all processors to catch up */
    CCTK_Barrier (GH);
#endif
    
    /* extradir is the relative output directory */
    sprintf (extradir, "%s_3d/", alias);
    
  }
  
  
  /* CREATE THE COMPLETE OUTPUT FILENAME 
     ----------------------------------- */
  
  IEEEfile_3D = (IEEEfile_3D_t *) malloc (sizeof (IEEEfile_3D_t));

  IEEEfile_3D->filename = (char *) malloc (strlen (myGH->outdir3D) +
                                           strlen (extradir) +
                                           strlen (alias) +
                                           strlen (extra) +
                                           strlen (pughGH->identity_string) +
                                           10);
  sprintf (IEEEfile_3D->filename, "%s/%s%s_3d%s%s.ieee", myGH->outdir3D,
           extradir, alias, extra, pughGH->identity_string);

  /* no need to store file info if used only once */
  if (! out3D_septimefiles) {
    if (ioUtilGH->recovered) {
      IEEEfile_3D->iofile = IEEEopen (IEEEfile_3D->filename, "a");
      *isNewFile = ! IOisValid (IEEEfile_3D->iofile);
    } else
      *isNewFile = 1;
    StoreNamedData (&myGH->fileList_3D, alias, IEEEfile_3D);
  } else
    *isNewFile = 1;

  return (IEEEfile_3D);

}