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$
 @@*/

static char *rcsid = "$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"


/* info structure describing IEEEIO files 
   For each opened file there is one entry in the file database. */
typedef struct {
  IOFile iofile;
  char *filename;
} IEEEfile_3D_t;


/* local function prototypes */
IEEEfile_3D_t *IOFlexIO_Get3Dfile (cGH *GH, const char *alias, int *isNewFile);
void IOFlexIO_Write3D_openFile (cGH *GH, IEEEfile_3D_t *IEEEfile_3D,
                                int isNewFile);
void IOFlexIO_Write3D_closeFile (cGH *GH, IEEEfile_3D_t *IEEEfile_3D);


/*@@
  @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
  @calls    CCTK_GHExtensionHandle
            CCTK_MyProc
            CCTK_nProcs
            CCTK_WARN
  @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 myproc;
  int timelevel;
  ioGH *ioUtilGH;
  int isNewFile;
  IEEEfile_3D_t *IEEEfile_3D;


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

  if (verbose) {
    printf ("-------------------------------------------------------\n");
    printf ("3D I/O on Grid Function %s\n", alias);
    printf ("  downsample (x,y,z): (%d,%d,%d)  ioproc : %d\n", 
            ioUtilGH->downsample_x, ioUtilGH->downsample_y,
            ioUtilGH->downsample_z, ioUtilGH->ioproc);
    fflush (stdout);
  }

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


  /* Get the filename and descriptor for output */
  /* The output directory is also created if needed. */
  if (myproc == ioUtilGH->ioproc)
    IEEEfile_3D = IOFlexIO_Get3Dfile (GH, alias, &isNewFile);
  else
    IEEEfile_3D = NULL;

#if 0
  /* Close the file if it is already open (if IEEEfile_3D is set). 
   *  This shouldn't be needed, since we do it below, 
   *  but better safe than sorry.
   *  Also I have some vague recollection of needing it in some
   *  very strange case ... So leave it here for now (PW 11.5.98)
   */
  if (IEEEfile_3D->iofile >= 0 && out3D_septimefiles) {
    if (verbose)
      printf ("Closing IEEEfile from previous iteration.\n");
    CACTUS_IEEEIO_ERROR (IOclose (IEEEfile_3D->iofile));
    IEEEfile_3D->iofile = -1; /* prevent attempts to close it twice! */
  }
#endif

  /* open the output file (only if we are an output processor) */
  if (IEEEfile_3D)
    IOFlexIO_Write3D_openFile (GH, IEEEfile_3D, isNewFile);
  
  /* 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 scalars 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)
      IOFlexIO_IEEEIOStructDump (GH, IEEEfile_3D->iofile);
#if        0
    if (ioUtilGH->parameters)
      IOFlexIO_IEEEIOparamDump (IEEEfile_3D->iofile);
#endif
    /* close the file */
    IOFlexIO_Write3D_closeFile (GH, IEEEfile_3D);
  }

  if (verbose)
    printf ("-------------------\n");
}


/*@@
  @routine IOFlexIO_Write3D_closeFile
  @author  Paul Walker
  @date    Feb 1997
  @desc
  Closes or pauses the IEEEIO output file.
  @enddesc
  @history 
  @hauthor Gabrielle Allen
  @hdate Sep 1998
  @hdesc Split into subroutine
  @endhistory 
  @@*/


void IOFlexIO_Write3D_closeFile (cGH *GH, IEEEfile_3D_t *IEEEfile_3D)
{
  DECLARE_CCTK_PARAMETERS
  flexioGH *myGH;

  if (IEEEfile_3D == NULL)
    return;

  myGH = (flexioGH *) GH->extensions [CCTK_GHExtensionHandle ("IOFlexIO")];

  if (myGH->out3d_reuse_filehandles) {
    if (verbose)
      printf ("Pausing file %s\n", IEEEfile_3D->filename);
    IEEEbufferOff (IEEEfile_3D->iofile); 
    CACTUS_IEEEIO_ERROR (IOpause (IEEEfile_3D->iofile));
  } else {
    if (verbose)
      printf ("Closing file %s\n", IEEEfile_3D->filename);
    CACTUS_IEEEIO_ERROR (IOclose (IEEEfile_3D->iofile));
  }

  if (out3D_septimefiles)
    free (IEEEfile_3D);
}


/*@@
  @routine IOFlexIO_Write3D_openFile
  @author  Paul Walker
  @date    Feb 1997
  @desc
  Opens or resumes the IEEEIO output file
  @enddesc
  @history 
  @hauthor Gabrielle Allen
  @hdate Sep 1998
  @hdesc Split into subroutine
  @endhistory 
  @@*/

void IOFlexIO_Write3D_openFile (cGH *GH, IEEEfile_3D_t *IEEEfile_3D,
                                int isNewFile)
{
  DECLARE_CCTK_PARAMETERS
  flexioGH *myGH;

  myGH = (flexioGH *) GH->extensions [CCTK_GHExtensionHandle ("IOFlexIO")];

  /* If we are the first time through, or writing one file per slice, we need
     to open a file in write mode. Or, if we are reusing file handles we will
     need to resume the file */
  
  if (isNewFile) {

    /* First time through or one file per slice; open anew */
    if (verbose)
      printf ("Opening file '%s'\n", IEEEfile_3D->filename);

    IEEEfile_3D->iofile = IEEEopen (IEEEfile_3D->filename, "w");
    if (! IOisValid (IEEEfile_3D->iofile)) {
      char *msg = (char *) malloc (strlen (IEEEfile_3D->filename) + 80);

      sprintf (msg, "The file '%s' could not be opened for writing.",
               IEEEfile_3D->filename);
      CCTK_WARN (1, msg);
      free (msg);
    }

  } else if (myGH->out3d_reuse_filehandles) {

    /* resume the file (reopen descriptor)
       This allows descriptor reuse without requiring expensive destruction
       and reallocation of the associated datastructures */
    
    if (verbose)
      printf ("Resuming file '%s'\n", IEEEfile_3D->filename);

    CACTUS_IEEEIO_ERROR (IOresume (IEEEfile_3D->iofile));

  } else {
    /* File is closed, not out3D_septimefiles and not first time.
       Therefore we must be conserving file handles and failed
       to use the pause()/resume methods(), So append.
    */
    if (verbose)
      printf ("Re-opening file '%s' in append mode\n", IEEEfile_3D->filename);
  
    IEEEfile_3D->iofile = IEEEopen (IEEEfile_3D->filename, "a");
    if (! IOisValid (IEEEfile_3D->iofile)) {
      char *msg = (char *) malloc (strlen (IEEEfile_3D->filename) + 80);

      sprintf (msg, "The file %s could not be opened for appending.\n",
               IEEEfile_3D->filename);
      CCTK_WARN (1, msg);
      free (msg);
    }
  }

  /* Turn on buffer cache (with default size chosen by IEEEIO lib) */
  IEEEbufferOn (IEEEfile_3D->iofile, 0);
}


/*@@
  @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 
  @@*/

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;

  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);
  }

  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  pughGH = (pGH *) GH->extensions [CCTK_GHExtensionHandle ("PUGH")];

  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) {
      
      outputdir = (char *) malloc (strlen (myGH->outdir3D) +
                                   strlen (alias) + 5);
      sprintf (outputdir, "%s/%s_3d", myGH->outdir3D, alias);
      
      if (CCTK_mkdir (outputdir) != 0)
        CCTK_WARN (1, "Problem creating 3D output directory");

      free (outputdir);
    }
    
#ifdef MPI
    /* Wait for all processors to catch up */
    CACTUS_MPI_ERROR (MPI_Barrier (pughGH->PUGH_COMM_WORLD));
#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) +
                                           10);
  sprintf (IEEEfile_3D->filename, "%s/%s%s_3d%s.ieee", myGH->outdir3D,
           extradir, alias, extra);

  /* no need to store file info if used only once */
  if (! out3D_septimefiles)
    StoreNamedData (&myGH->fileList_3D, alias, IEEEfile_3D);

  /* set flag to indicate that file should be opened in create mode */
  *isNewFile = 1;

  return (IEEEfile_3D);
}