cosmotool/src/h5_readFlash.cpp

/*+
This is CosmoTool (./src/h5_readFlash.cpp) -- Copyright (C) Guilhem Lavaux (2007-2014)

guilhem.lavaux@gmail.com

This software is a computer program whose purpose is to provide a toolbox for cosmological
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)

This software is governed by the CeCILL license under French law and
abiding by the rules of distribution of free software.  You can  use,
modify and/ or redistribute the software under the terms of the CeCILL
license as circulated by CEA, CNRS and INRIA at the following URL
"http://www.cecill.info".

As a counterpart to the access to the source code and  rights to copy,
modify and redistribute granted by the license, users are provided only
with a limited warranty  and the software's author,  the holder of the
economic rights,  and the successive licensors  have only  limited
liability.

In this respect, the user's attention is drawn to the risks associated
with loading,  using,  modifying and/or developing or reproducing the
software by the user in light of its specific status of free software,
that may mean  that it is complicated to manipulate,  and  that  also
therefore means  that it is reserved for developers  and  experienced
professionals having in-depth computer knowledge. Users are therefore
encouraged to load and test the software's suitability as regards their
requirements in conditions enabling the security of their systems and/or
data to be ensured and,  more generally, to use and operate it in the
same conditions as regards security.

The fact that you are presently reading this means that you have had
knowledge of the CeCILL license and that you accept its terms.
+*/
/*+
!!

This file has been developped by P. M. Sutter

!!

*/
/* This file contains the functions that read the data from the HDF5 file
 * The functions accept the PARAMESH data through arguments.
 */

#include "h5_readFlash.hpp"

using namespace H5;

/* indices of attributes in fof_particle_type */
int iptag_out = 0;
int ipx_out   = 0;
int ipy_out   = 1;
int ipz_out   = 2;
int ipvx_out   = 3;
int ipvy_out   = 4;
int ipvz_out   = 5;

/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */

/*  n*_runtime_parameters should be set by the caller to
    the maximum number of runtime parameters to read.
*/

void h5_read_runtime_parameters
     (H5File* file,    /* file handle */
      double* LBox,    // box size
      int* numPart,
      double *hubble,
      double *omegam,
      double *omegalambda)       // number of particles
{

 int MAX_PARM = 200;

 int nreal_runtime_parameters, nint_runtime_parameters, nstr_runtime_parameters;

 char real_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
 char int_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
 char str_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];

 double real_runtime_parameter_values[MAX_PARM];
 int    int_runtime_parameter_values[MAX_PARM];
 char   str_runtime_parameter_values[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];

  int rank;
  hsize_t dimens_1d, maxdimens_1d;

  real_runtime_params_t *real_rt_parms;
  int_runtime_params_t *int_rt_parms;
  str_runtime_params_t *str_rt_parms;
  log_runtime_params_t *log_rt_parms;

  double omegarad;

  int i;

  nint_runtime_parameters = MAX_PARM;
  nreal_runtime_parameters = MAX_PARM;

  /* integer runtime parameters */
  int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t));

  rank = 1;
  DataSet dataset = file->openDataSet("integer runtime parameters");

  IntType int_rt_type = dataset.getIntType();
  //int_rt_type = H5Dget_type(dataset);

  DataSpace dataspace = dataset.getSpace();
  //dataspace = H5Dget_space(dataset);

  int ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
  //H5Sget_simple_extent_dims(dataspace, &dimens_1d, &maxdimens_1d);

  /* don't read in more than we can handle */
  if (nint_runtime_parameters < dimens_1d) {
    dimens_1d = nint_runtime_parameters;
  } else {
    nint_runtime_parameters = dimens_1d;
  }
  DataSpace memspace(rank, &dimens_1d);
  //memspace = H5Screate_simple(rank, &dimens_1d, NULL);

  dataset.read(int_rt_parms, int_rt_type, memspace, dataspace);


  for (i = 0; i < nint_runtime_parameters; i++) {
    strncpy(int_runtime_parameter_names[i], int_rt_parms[i].name, RUNTIME_PARAMETER_STRING_SIZE);
    int_runtime_parameter_values[i] = int_rt_parms[i].value;
  }

  free(int_rt_parms);
  memspace.close();
  dataspace.close();
  dataset.close();
  //H5Sclose(dataspace);
  //H5Dclose(dataset);

  /* done with int runtime parameters */

  /* reals */

  real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t));

  rank = 1;
  dataset = file->openDataSet("real runtime parameters");
  //dataset = H5Dopen(*file_identifier, "real runtime parameters");

  dataspace = dataset.getSpace();
  FloatType real_rt_type = dataset.getFloatType();
  ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
  //dataspace = H5Dget_space(dataset);
  //real_rt_type = H5Dget_type(dataset);
  //H5Sget_simple_extent_dims(dataspace, &dimens_1d, &maxdimens_1d);

  /* don't read in more than we can handle */
  if (nreal_runtime_parameters < dimens_1d) {
    dimens_1d = nreal_runtime_parameters;
  } else {
    nreal_runtime_parameters = dimens_1d;
  }
   memspace = DataSpace(rank, &dimens_1d);
  //memspace = H5Screate_simple(rank, &dimens_1d, NULL);

  dataset.read(real_rt_parms, real_rt_type, memspace, dataspace);


  for (i = 0; i < nreal_runtime_parameters; i++) {
    strncpy(real_runtime_parameter_names[i], real_rt_parms[i].name, RUNTIME_PARAMETER_STRING_SIZE);
    real_runtime_parameter_values[i] = real_rt_parms[i].value;
  }

  free(real_rt_parms);
  memspace.close();
  dataspace.close();
  dataset.close();
  //H5Sclose(dataspace);
  //H5Dclose(dataset);

  /* done with reals */

  // grab the data we want
  for (i = 0; i < nreal_runtime_parameters; i++) {
    if (strncmp(real_runtime_parameter_names[i],"xmax",4) == 0 ) {
      *LBox = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"hubbleconstant", 14) == 0 ) {
      *hubble = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"omegamatter", 11) == 0 ) {
      *omegam = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"omegaradiation", 11) == 0 ) {
      omegarad = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"cosmologicalconstant", 20) == 0 ) {
      *omegalambda = real_runtime_parameter_values[i];
    }
  }

  for (i = 0; i < nint_runtime_parameters; i++) {
    if (strncmp(int_runtime_parameter_names[i],"pt_numx",7) == 0 ) {
      *numPart = int_runtime_parameter_values[i];
      *numPart = *numPart * *numPart * *numPart;
    }
  }
}


/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
void h5_read_flash3_particles (H5File* file,
                        int* totalparticles,
                        int* localnp,
                        int* particle_offset,
                        float *pos1,
                        float *pos2,
                        float *pos3,
                        float *vel1,
                        float *vel2,
                        float *vel3,
                        int64_t    *id)
{

  herr_t   status;
  hsize_t  dimens_1d, maxdims_1d;
  hsize_t  start_1d, stride_1d, count_1d;
  int rank;

  int      numProps, i, p;
  int      numPartBuffer = 5000, sizePartBuffer, sizePart;
  int      pstack, poffset, pcount;
  int      iptag, ipx, ipy, ipz, ipvx, ipvy, ipvz;
  char     *propName;
  double   *partBuffer;
  char    part_names[50][OUTPUT_PROP_LENGTH];
  int     string_size;

//  char part_names[NPART_PROPS][OUTPUT_PROP_LENGTH];
  hsize_t dimens_2d[2], maxdimens_2d[2];
  hsize_t start_2d[2], count_2d[2], stride_2d[2];

  /* skip this routine if no particles to read */
  if ((*localnp) == 0) {
    return;
  }

 /* first determine how many particle properties are
     present in the input data file, and determine which of these
     are the properties we are interested in */
  DataSet dataset = file->openDataSet("particle names");
  DataSpace dataspace = dataset.getSpace();
  //dataset = H5Dopen(*file_identifier, "particle names");
  //dataspace = H5Dget_space(dataset);

  int ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
  //H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);

  //total number of particle properties
  numProps = dimens_2d[0];

  string_size = OUTPUT_PROP_LENGTH;
  StrType string_type = H5Tcopy(H5T_C_S1);
  string_type.setSize(string_size);
  //status = H5Tset_size(string_type, string_size);

  rank = 2;

  start_2d[0] = 0;
  start_2d[1] = 0;

  stride_2d[0] = 1;
  stride_2d[1] = 1;

  count_2d[0] = dimens_2d[0];
  count_2d[1] = dimens_2d[1];

  dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
  //status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
  //                      stride_2d, count_2d, NULL);

  DataSpace memspace(rank, dimens_2d);
  //memspace = H5Screate_simple(rank, dimens_2d, NULL);

  dataset.read(part_names, string_type);


  string_type.close();
  memspace.close();
  dataspace.close();
  dataset.close();
  //H5Tclose(string_type);
  //H5Sclose(memspace);
  //H5Sclose(dataspace);
  //H5Dclose(dataset);

  for (i=0;i<numProps;i++) {
    if (strncmp(part_names[i], "tag", 3) == 0)  { iptag = i+1; }
    if (strncmp(part_names[i], "posx", 4) == 0) { ipx = i+1; }
    if (strncmp(part_names[i], "posy", 4) == 0) { ipy = i+1; }
    if (strncmp(part_names[i], "posz", 4) == 0) { ipz = i+1; }
    if (strncmp(part_names[i], "velx", 4) == 0) { ipvx = i+1; }
    if (strncmp(part_names[i], "vely", 4) == 0) { ipvy = i+1; }
    if (strncmp(part_names[i], "velz", 4) == 0) { ipvz = i+1; }
  }

  if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) ||
       (ipvy < 0) || (ipvz < 0) ) {
    printf("One or more required particle attributes not found in file!\n");
    return;
  }

  //printf("iptag = %d, ipx = %d, ipy = %d, ipz = %d\n", iptag, ipx, ipy, ipz);
  //printf("ipvx = %d, ipvy = %d, ipvz = %d\n", ipvx, ipvy, ipvz);

  //read particles
  dataset = file->openDataSet("tracer particles");
  //dataset = H5Dopen(*file_identifier, "tracer particles");

  FloatType datatype = dataset.getFloatType();
  //datatype = H5Dget_type(dataset);

  /* establish read-in particle buffer */
  sizePart = numProps*(sizeof(double));
  sizePartBuffer = numPartBuffer * sizePart;
  partBuffer = (double *)malloc(sizePartBuffer);

  dataspace = dataset.getSpace();
  ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
  //dataspace = H5Dget_space(dataset);
  //H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);

  /*insert particle properties (numPartBuffer) particles at a time*/
  pstack = (*localnp);
  poffset = 0;
  if (pstack > numPartBuffer) {
   pcount = numPartBuffer;
  }
  else {
   pcount = pstack;
  }

  while ( pstack > 0) {
    rank       = 2;
    maxdimens_2d[0] = (hsize_t) (*totalparticles);
    maxdimens_2d[1] = (hsize_t) (numProps);

    start_2d[0]  = (hsize_t) (*particle_offset + poffset);
    start_2d[1]  = (hsize_t) 0;

    stride_2d[0] = 1;
    stride_2d[1] = 1;

    count_2d[0]  = (hsize_t) (pcount);
    count_2d[1]  = (hsize_t) (numProps);

    dimens_2d[0] = (pcount);
    dimens_2d[1] = (numProps);

    dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
    //status     = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
    //                                 stride_2d, count_2d, NULL);

    memspace = DataSpace(rank, dimens_2d);
    //memspace   = H5Screate_simple(rank, dimens_2d, maxdimens_2d);

    /* read data from the dataset */
   dataset.read(partBuffer, datatype, memspace, dataspace);

    /* convert buffer into particle struct */

    if (id) {
    for(p=0; p < (pcount); p++) {
      id[p+poffset] = (int64_t)  *(partBuffer+iptag-1+p*numProps);
    } }

    if (pos1 && pos2 && pos3) {
    for(p=0; p < (pcount); p++) {
      pos1[p+poffset]  = (float) *(partBuffer+ipx-1+p*numProps);
      pos2[p+poffset]  = (float) *(partBuffer+ipy-1+p*numProps);
      pos3[p+poffset]  = (float) *(partBuffer+ipz-1+p*numProps);
    }
    }


    if (vel1 && vel2 && vel3) {
    for(p=0; p < (pcount); p++) {
      vel1[p+poffset]  = (float) *(partBuffer+ipvx-1+p*numProps);
      vel2[p+poffset]  = (float) *(partBuffer+ipvy-1+p*numProps);
      vel3[p+poffset]  = (float) *(partBuffer+ipvz-1+p*numProps);
    }
    }

    memspace.close();
    //status = H5Sclose(memspace);
   /* advance buffer */
    pstack = pstack - pcount;
    poffset = poffset + pcount;
    if (pstack > numPartBuffer) {
      pcount = numPartBuffer;
    }
    else {
      pcount = pstack;
    }

  } /* end while */


  datatype.close();
  dataspace.close();
  dataset.close();
  //status = H5Tclose(datatype);
  //status = H5Sclose(dataspace);
  //status = H5Dclose(dataset);
}


/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

void h5_read_flash3_header_info(H5File* file,
				double* time,                  /* simulation time */
				double* redshift)              /* redshift of checkpoint */
{

  herr_t status;

  int file_version;

  hid_t sp_type, si_type;

  hsize_t dimens_1d, maxdimens_1d;
  hid_t string_type;
  real_list_t *real_list;
  int* num_real, num_int;
  int MAX_SCALARS = 100;
  char real_names[MAX_SCALARS][MAX_STRING_LENGTH];
  double real_values[MAX_SCALARS];
  char int_names[MAX_SCALARS][MAX_STRING_LENGTH];
  int int_values[MAX_SCALARS];
  int i;

  H5std_string DATASET_NAME;

  string_type = H5Tcopy(H5T_C_S1);
  H5Tset_size(string_type, MAX_STRING_LENGTH);

  DataSet dataset = file->openDataSet("real scalars");
  DataSpace dataspace = dataset.getSpace();

  /* read extent of 'dataspace' (i.e. # of name/value pairs) into 'dimens_1d' */
  int ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);

  if (dimens_1d > MAX_SCALARS) {
    printf("Error: reading more than MAX_SCALARS runtime parameters in checkpoint file!\n");
  }

  /* malloc a pointer to a list of real_list_t's */
  real_list = (real_list_t *) malloc(dimens_1d * sizeof(real_list_t));

  // create a new simple dataspace of 1 dimension and size of 'dimens_1d'
  DataSpace memspace(1, &dimens_1d);

  // create an empty vessel sized to hold one real_list_t's worth of data
  CompType real_list_type( sizeof(real_list_t) );

  // subdivide the empty vessel into its component sections (name and value)
  real_list_type.insertMember(
          "name",
          HOFFSET(real_list_t, name),
          string_type);

  real_list_type.insertMember(
          "value",
          HOFFSET(real_list_t, value),
          PredType::NATIVE_DOUBLE);

  // read the data into 'real_list'
  dataset.read( real_list, real_list_type, memspace, dataspace);


  if (status < 0) {
    printf("Error readingruntime parameterss from data file\n");
  }

  for (i = 0; i < dimens_1d; i++) {
    strncpy(real_names[i], real_list[i].name, MAX_STRING_LENGTH);
    real_values[i] = real_list[i].value;

    if (strncmp(real_names[i],"time",4) == 0 ) {
      *time = real_values[i];
    }
    if (strncmp(real_names[i],"redshift",8) == 0 ) {
      *redshift = real_values[i];
    }
  }

  free(real_list);
  real_list_type.close();
  memspace.close();
  dataspace.close();
  dataset.close();
}