Commits

John Wise committed daaa523 Merge

Merged in ngoldbaum/active-particles (pull request #4)

Comments (0)

Files changed (11)

src/enzo/ActiveParticle.C

     return;
 }
 
+int ActiveParticleType::ReadDataset(int ndims, hsize_t *dims, char *name, hid_t group,
+				    hid_t data_type, void *read_to)
+{
+  hid_t file_dsp_id;
+  hid_t dset_id;
+  hid_t h5_status;
+  herr_t      h5_error = -1;
+  int i, j, k, dim;
+  /* get data into temporary array */
+
+  file_dsp_id = H5Screate_simple((Eint32) ndims, dims, NULL);
+  if( file_dsp_id == h5_error ){ENZO_FAIL("Error creating file dataspace");}
+
+  dset_id =  H5Dopen(group, name);
+  if( dset_id == h5_error )ENZO_VFAIL("Error opening %s", name)
+
+  h5_status = H5Dread(dset_id, data_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, (VOIDP) read_to);
+  if( dset_id == h5_error )ENZO_VFAIL("Error reading %s", name)
+
+  h5_status = H5Sclose(file_dsp_id);
+  if( dset_id == h5_error )ENZO_VFAIL("Error closing dataspace %s", name)
+
+  h5_status = H5Dclose(dset_id);
+  if( dset_id == h5_error )ENZO_VFAIL("Error closing %s", name)
+			     
+  return SUCCESS;
+}
+
 int ActiveParticleType::WriteDataset(int ndims, hsize_t *dims, char *name, hid_t group,
-		  hid_t data_type, void *data)
+				     hid_t data_type, void *data)
 {
 
   hid_t file_dsp_id;

src/enzo/ActiveParticle.h

 			  ActiveParticleFormationData &data);
   int static WriteDataset(int ndims, hsize_t *dims, char *name, hid_t group,
 			  hid_t data_type, void *data);
+  int static ReadDataset(int ndims, hsize_t *dims, char *name, hid_t group,
+			 hid_t data_type, void *read_to);
   /* Several pure virtual functions */
   
   /* This should return the number of new star particles created, and should
    int (*ifunc)(),
    int (*feedfunc)(grid *thisgrid_orig, ActiveParticleFormationData &data),
    int (*writefunc)(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id),
+   int (*readfunc)(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id),
    ActiveParticleType *particle
    ){
     this->formation_function = ffunc;
     this->initialize = ifunc;
     this->feedback_function = feedfunc;
     this->write_function = writefunc;
+    this->read_function = readfunc;
     get_active_particle_types()[this_name] = this;
   }
 
   int (*formation_function)(grid *thisgrid_orig, ActiveParticleFormationData &data);
   int (*feedback_function)(grid *thisgrid_orig, ActiveParticleFormationData &data);
   int (*write_function)(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  int (*read_function)(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   void (*describe_data_flags)(ActiveParticleFormationDataFlags &flags);
   ParticleBufferHandler* (*allocate_buffers)(int NumberOfParticles);
   ActiveParticleType* particle_instance;
      (&active_particle_class::InitializeParticleType),
      (&active_particle_class::EvaluateFeedback),
      (&active_particle_class::WriteToOutput),
+     (&active_particle_class::ReadFromOutput),
      pp);
   return pinfo;
 }

src/enzo/ActiveParticle_CenOstriker.C

 #include "TopGridData.h"
 #include "EventHooks.h"
 #include "ActiveParticle.h"
+#include "h5utilities.h"
 
 #ifdef NEW_CONFIG
 
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int InitializeParticleType();
   ENABLED_PARTICLE_ID_ACCESSOR
   int i, j, k, dim, index, offset_y, offset_z;
   int NumberOfNewParticles = 0;
 
-  /* Define some pointers for readability */
 
   float *density = thisGrid->BaryonField[data.DensNum];
   float *velx = thisGrid->BaryonField[data.Vel1Num];
 
 
 	np->pos[0] = thisGrid->CellLeftEdge[0][i] + 0.5*thisGrid->CellWidth[0][i];
-	np->pos[0] = thisGrid->CellLeftEdge[1][j] + 0.5*thisGrid->CellWidth[1][j];
-	np->pos[0] = thisGrid->CellLeftEdge[2][k] + 0.5*thisGrid->CellWidth[2][k];
+	np->pos[1] = thisGrid->CellLeftEdge[1][j] + 0.5*thisGrid->CellWidth[1][j];
+	np->pos[2] = thisGrid->CellLeftEdge[2][k] + 0.5*thisGrid->CellWidth[2][k];
 	
 	if (UnigridVelocities == false) {
 	  float *tvel = thisGrid->AveragedVelocityAtCell(index,data.DensNum,data.Vel1Num);
   ActiveParticleType_CenOstriker *particles = 
     static_cast<ActiveParticleType_CenOstriker*>(*thisGrid->ActiveParticles);
   
-  /* Define some pointers for readability  */
   float *density = thisGrid->BaryonField[data.DensNum];
   float *velx = thisGrid->BaryonField[data.Vel1Num];
   float *vely = thisGrid->BaryonField[data.Vel2Num];
   flags.MetalField = true;
 }
 
+int ActiveParticleType_CenOstriker::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+  /* Open the subgroup within the active particle for active particles of type CenOstriker */
+
+  hid_t CenOstrikerGroupID = H5Gopen(group_id,"CenOstriker");
+
+  readAttribute(CenOstrikerGroupID,HDF5_INT,"number_of_active_particles_of_this_type",n);
+
+  char *ParticlePositionLabel[] =
+     {"position_x", "position_y", "position_z"};
+  char *ParticleVelocityLabel[] =
+     {"velocity_x", "velocity_y", "velocity_z"};
+
+  FLOAT Position[GridRank][*n];
+  float Velocity[GridRank][*n]; 
+  double Mass[*n];
+  float BirthTime[*n];
+  float DynamicalTime[*n];
+  float Metallicity[*n];
+  
+  int i,dim;
+
+  hsize_t TempInt;
+  TempInt = *n;
+  
+  for (dim = 0; dim < GridRank; dim++) {
+    ReadDataset(1,&TempInt,ParticlePositionLabel[dim],
+		  CenOstrikerGroupID, HDF5_FILE_PREC, (VOIDP) Position[dim]);
+  }
+
+  for (dim = 0; dim < GridRank; dim++) {
+    ReadDataset(1,&TempInt,ParticleVelocityLabel[dim],
+		  CenOstrikerGroupID, HDF5_FILE_REAL, (VOIDP) Velocity[dim]);
+  }
+  ReadDataset(1,&TempInt,"mass",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Mass);
+  ReadDataset(1,&TempInt,"creation_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) BirthTime);
+  ReadDataset(1,&TempInt,"dynamical_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) DynamicalTime);
+  ReadDataset(1,&TempInt,"metallicity_fraction",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Metallicity);
+
+  for (i = 0; i < *n; i++) {
+    ActiveParticleType_CenOstriker *np = new ActiveParticleType_CenOstriker;
+    particles_to_read[i] = np;
+    np = new ActiveParticleType_CenOstriker();
+    np->Mass = Mass[i];
+    np->type = CenOstriker;
+    np->BirthTime = DynamicalTime[i];
+    np->Metallicity = Metallicity[i];
+    for (dim = 0; dim < GridRank; dim++){
+      np->pos[dim] = Position[dim][i];
+      np->vel[dim] = Velocity[dim][i];
+    }
+  }
+  return SUCCESS;
+}
+
 int ActiveParticleType_CenOstriker::WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id)
 {
   /* Create a new subgroup within the active particle group for active particles of type CenOstriker */
 
   hid_t CenOstrikerGroupID = H5Gcreate(group_id,"CenOstriker",0);
 
+  writeScalarAttribute(CenOstrikerGroupID,HDF5_INT,"number_of_active_particles_of_this_type",&n);
+
   ActiveParticleType_CenOstriker *ParticlesToWrite = static_cast<ActiveParticleType_CenOstriker*>(these_particles);
 
   char *ParticlePositionLabel[] =
-     {"CenOstrikerParticle_position_x", "CenOstrikerParticle_position_y", "CenOstrikerParticle_position_z"};
+     {"position_x", "position_y", "position_z"};
   char *ParticleVelocityLabel[] =
-     {"CenOstrikerParticle_velocity_x", "CenOstrikerParticle_velocity_y", "CenOstrikerParticle_velocity_z"};
+     {"velocity_x", "velocity_y", "velocity_z"};
 
   /* Create temporary buffers to store particle data */
 
 		  CenOstrikerGroupID, HDF5_FILE_REAL, (VOIDP) Velocity[dim]);
   }
   
-  WriteDataset(1,&TempInt,"CenOstriker_mass",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Mass);
-  WriteDataset(1,&TempInt,"CenOstriker_creation_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) BirthTime);
-  WriteDataset(1,&TempInt,"CenOstriker_dynamical_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) DynamicalTime);
-  WriteDataset(1,&TempInt,"CenOstriker_metallicity_fraction",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Metallicity);
+  WriteDataset(1,&TempInt,"mass",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Mass);
+  WriteDataset(1,&TempInt,"creation_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) BirthTime);
+  WriteDataset(1,&TempInt,"dynamical_time",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) DynamicalTime);
+  WriteDataset(1,&TempInt,"metallicity_fraction",CenOstrikerGroupID,HDF5_FILE_REAL,(VOIDP) Metallicity);
+
+  return SUCCESS;
 }
 
 class CenOstrikerBufferHandler : public ParticleBufferHandler

src/enzo/ActiveParticle_Kravtsov.C

   static int EvaluateFormation(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int InitializeParticleType();
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   return SUCCESS;
 }
 
+int ActiveParticleType_Kravtsov::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+
+
+  return SUCCESS;
+}
+
 class KravtsovParticleBufferHandler : public ParticleBufferHandler
 {
   public:

src/enzo/ActiveParticle_PopIII.C

   static int EvaluateFormation(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int InitializeParticleType();
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   return SUCCESS;
 }
 
+int ActiveParticleType_PopIII::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+
+
+  return SUCCESS;
+}
+
 class PopIIIParticleBufferHandler : public ParticleBufferHandler
 {
   public:

src/enzo/ActiveParticle_SampleParticle.C

   static int EvaluateFormation(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static int InitializeParticleType(void);
   return SUCCESS;
 }
 
+int ActiveParticleType_SampleParticle::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+
+  return SUCCESS;
+}
+
 class SampleParticleBufferHandler : public ParticleBufferHandler
 {
   public:

src/enzo/ActiveParticle_SinkParticle.C

 public:
   static int EvaluateFormation(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   return SUCCESS;
 }
 
+int ActiveParticleType_SinkParticle::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+
+
+  return SUCCESS;
+}
+
+
 class SinkParticleBufferHandler : public ParticleBufferHandler
 {
   public:

src/enzo/ActiveParticle_SpringelHernquist.C

   static int EvaluateFormation(grid *thisgrid_orig, ActiveParticleFormationData &data);
   static void DescribeSupplementalData(ActiveParticleFormationDataFlags &flags);
   static int WriteToOutput(ActiveParticleType *these_particles, int n, int GridRank, hid_t group_id);
+  static int ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id);
   static ParticleBufferHandler *AllocateBuffers(int NumberOfParticles);
   static int InitializeParticleType();
   static int EvaluateFeedback(grid *thisgrid_orig, ActiveParticleFormationData &data);
   return SUCCESS;
 }
 
+int ActiveParticleType_SpringelHernquist::ReadFromOutput(ActiveParticleType **particles_to_read, int *n, int GridRank, hid_t group_id)
+{
+
+
+  return SUCCESS;
+}
+
 class SpringelHernquistParticleBufferHandler : public ParticleBufferHandler
 {
 public:

src/enzo/Make.mach.darwin

File contents unchanged.

src/enzo/New_Grid_ReadGrid.C

 //  Input a grid from file pointer fpt
  
 #include <hdf5.h>
+#include <map>
+#include <iostream>
+#include <stdexcept>
 #include <string.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <unistd.h>
 #include <math.h>
 #include <assert.h>
+#include <vector>
 #include "h5utilities.h"
  
 #include "ErrorExceptions.h"
 #include "GridList.h"
 #include "ExternalBoundary.h"
 #include "Grid.h"
+#include "fortran.def"
+#include "CosmologyParameters.h"
+#include "ActiveParticle.h"
+
 void my_exit(int status);
  
 #ifdef PROTO /* Remove troublesome HDF PROTO declaration. */
 		ENZO_FAIL("Error reading GravityBoundaryType.");
       }
 
+    /* 5) Read active particle info */
+    
+    if (fscanf(fptr, "NumberOfActiveParticles = %"ISYM"\n", &NumberOfActiveParticles) != 1) {
+      ENZO_FAIL("error reading NumberOfActiveParticles.");
+    }
+    
+    if (NumberOfActiveParticles > 0) {
+      
+      /* Read particle file name. */
+      
+      if (fscanf(fptr, "ParticleFileName = %s\n", procfilename) != 1) {
+	ENZO_FAIL("Error reading ParticleFileName.");
+      }
+    }
+    
     // If HierarchyFile has different Ghostzones (which should be a parameter not a macro ...)
     // (useful in a restart with different hydro/mhd solvers) 
     int ghosts =DEFAULT_GHOST_ZONES;
 
   } // end: if (NumberOfParticles > 0) && ReadData && (MyProcessorNumber == ProcessorNumber)
  
+  if (NumberOfActiveParticles > 0 && ReadData && (MyProcessorNumber == ProcessorNumber)) {
+
+    /* Open file if not already done (note: particle name must = grid name). */
+
+    if ((NumberOfBaryonFields == 0 || ReadParticlesOnly) && NumberOfParticles == 0) {
+
+#ifndef SINGLE_HDF5_OPEN_ON_INPUT 
+      file_id = H5Fopen(procfilename, H5F_ACC_RDONLY, H5P_DEFAULT);
+      if( file_id == h5_error )ENZO_VFAIL("Error opening file %s", name)
+#endif
+ 
+      group_id = H5Gopen(file_id, name);
+      if( group_id == h5_error )ENZO_VFAIL("Error opening group %s", name)
+
+    } 
+
+    /* Open the active particles group */
+
+    hid_t ActiveParticleGroupID = H5Gopen(group_id, "ActiveParticles");
+
+    /* Loop over enabled active particle types */
+
+    int NumberOfActiveParticlesOnDisk;
+    
+    ActiveParticleType **ActiveParticlesOnDisk;
+
+    for (i = 0; i < EnabledActiveParticlesCount; i++)
+      {
+    
+	/* Instantiate an active particle helper of this type
+	   This class contains the function that allows us to read from disk */
+
+	ActiveParticleType_info *ActiveParticleTypeToEvaluate = EnabledActiveParticles[i];
+	
+	/* Create an empty buffer that particles in the data file will be read into */
+
+	ActiveParticleType **ActiveParticlesOfThisTypeOnDisk;
+	int NumberOfActiveParticlesOfThisType;
+
+	/* Read the active particles from disk */
+
+	ActiveParticleTypeToEvaluate->read_function(ActiveParticlesOfThisTypeOnDisk,
+						    &NumberOfActiveParticlesOfThisType,
+						    GridRank,
+						    ActiveParticleGroupID);
+
+	/* Add the active particles read from disk to the active particle buffer */
+
+	int OldNumberOfActiveParticles = NumberOfActiveParticlesOnDisk;
+
+	ActiveParticleType **OldActiveParticles = ActiveParticlesOnDisk;
+
+	NumberOfActiveParticlesOnDisk += NumberOfActiveParticlesOfThisType;
+	
+	ActiveParticlesOnDisk = new ActiveParticleType*[NumberOfActiveParticlesOnDisk];
+
+	for (i = 0; i < OldNumberOfActiveParticles; i++) {
+	  ActiveParticlesOnDisk[i] = OldActiveParticles[i];
+	}
+	for (i = 0; i < NumberOfActiveParticlesOfThisType; i++) {
+	  ActiveParticlesOnDisk[OldNumberOfActiveParticles + i] = ActiveParticlesOfThisTypeOnDisk[i];
+	}
+
+	/* clean up */
+
+	delete [] OldActiveParticles;
+	delete [] ActiveParticlesOfThisTypeOnDisk;
+	
+      } // end: for EnabledActiveParticlesCount
+
+    /* Assign the active particle buffer to this grid */
+
+    this->ActiveParticles = ActiveParticlesOnDisk;
+
+  } // end: if (NumberOfActiveParticles > 0) && ReadData && (MyProcessorNumber == ProcessorNumber)
+
 
   /* Close file. */
  

src/enzo/New_Grid_WriteGrid.C

 	
 	/* Write them to disk */
 
-	ActiveParticleTypeToEvaluate->write_function(*ActiveParticlesOfThisType,NumberOfActiveParticlesOfThisType,
-						    GridRank,ActiveParticleGroupID);
-
+	ActiveParticleTypeToEvaluate->write_function(*ActiveParticlesOfThisType,
+						     NumberOfActiveParticlesOfThisType,
+						     GridRank,
+						     ActiveParticleGroupID);
+						     
 	/* Clean up */
 
 	delete [] ActiveParticlesOfThisType;