The entropy_in_IC flag is now set and read correctly. Restart from snapshot is possible.

examples/main.c

@@ -322,18 +322,21 @@ int main(int argc, char *argv[]) {
   size_t Ngas = 0, Ngpart = 0;
   double dim[3] = {0., 0., 0.};
   int periodic = 0;
+  int flag_entropy_ICs = 0;
   if (myrank == 0) clocks_gettime(&tic);
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
   read_ic_parallel(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+                   &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
+                   MPI_INFO_NULL, dry_run);
 #else
   read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+                 &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
+                 MPI_INFO_NULL, dry_run);
 #endif
 #else
   read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                 dry_run);
+                 &flag_entropy_ICs, dry_run);
 #endif
   if (myrank == 0) {
     clocks_gettime(&toc);
@@ -466,7 +469,7 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* Initialise the particles */
-  engine_init_particles(&e);
+  engine_init_particles(&e, flag_entropy_ICs);
 
   /* Legend */
   if (myrank == 0)

src/engine.c

@@ -2140,8 +2140,10 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
  *forward in time.
  *
  * @param e The #engine
+ * @param flag_entropy_ICs Did the 'Internal Energy' of the particles actually
+ *contain entropy ?
  */
-void engine_init_particles(struct engine *e) {
+void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
 
   struct space *s = e->s;
 
@@ -2209,7 +2211,7 @@ void engine_init_particles(struct engine *e) {
   TIMER_TOC(timer_runners);
 
   /* Apply some conversions (e.g. internal energy -> entropy) */
-  space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
+  if (!flag_entropy_ICs) space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
 
   clocks_gettime(&time2);
 

src/engine.h

@@ -214,7 +214,7 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
                    unsigned int submask);
 void engine_prepare(struct engine *e);
 void engine_print(struct engine *e);
-void engine_init_particles(struct engine *e);
+void engine_init_particles(struct engine *e, int flag_entropy_ICs);
 void engine_step(struct engine *e);
 void engine_maketasks(struct engine *e);
 void engine_split(struct engine *e, struct partition *initial_partition);

src/hydro/Default/hydro_io.h

@@ -122,3 +122,10 @@ void writeSPHflavour(hid_t h_grpsph) {
   writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
                    const_max_u_change);
 }
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }

src/hydro/Gadget2/hydro_io.h

@@ -83,8 +83,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
   writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "SmoothingLength",
              FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, us,
              UNIT_CONV_LENGTH);
-  writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy", FLOAT, N,
-             1, parts, N_total, mpi_rank, offset, entropy, us,
+  writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "InternalEnergy",
+             FLOAT, N, 1, parts, N_total, mpi_rank, offset, entropy, us,
              UNIT_CONV_ENTROPY_PER_UNIT_MASS);
   writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
              ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id, us,
@@ -113,3 +113,10 @@ void writeSPHflavour(hid_t h_grpsph) {
   writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
   writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
 }
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 1; }

src/hydro/Minimal/hydro_io.h

@@ -115,3 +115,10 @@ void writeSPHflavour(hid_t h_grpsph) {
   writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
                    const_max_u_change);
 }
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }

src/parallel_io.c

@@ -367,8 +367,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  */
 void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
                       struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                      int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                      MPI_Info info, int dry_run) {
+                      int* periodic, int* flag_entropy, int mpi_rank,
+                      int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -404,6 +404,7 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
   if (h_grp < 0) error("Error while opening file header\n");
 
   /* Read the relevant information and print status */
+  readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
   readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
   readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
   readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -626,6 +627,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
   double MassTable[6] = {0., 0., 0., 0., 0., 0.};
   writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
   unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+  flagEntropy[0] = writeEntropyFlag();
   writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                  NUM_PARTICLE_TYPES);
   writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);

src/parallel_io.h

@@ -36,8 +36,8 @@
 
 void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
                       struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                      int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                      MPI_Info info, int dry_run);
+                      int* periodic, int* flag_entropy, int mpi_rank,
+                      int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run);
 
 void write_output_parallel(struct engine* e, const char* baseName,
                            struct UnitSystem* us, int mpi_rank, int mpi_size,

src/serial_io.c

@@ -424,8 +424,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  */
 void read_ic_serial(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                    int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                    MPI_Info info, int dry_run) {
+                    int* periodic, int* flag_entropy, int mpi_rank,
+                    int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -462,6 +462,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
     if (h_grp < 0) error("Error while opening file header\n");
 
     /* Read the relevant information and print status */
+    readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
     readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
     readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
     readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -699,6 +700,7 @@ void write_output_serial(struct engine* e, const char* baseName,
     double MassTable[6] = {0., 0., 0., 0., 0., 0.};
     writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
     unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+    flagEntropy[0] = writeEntropyFlag();
     writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                    NUM_PARTICLE_TYPES);
     writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);

src/serial_io.h

@@ -36,8 +36,8 @@
 
 void read_ic_serial(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                    int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                    MPI_Info info, int dry_run);
+                    int* periodic, int* flag_entropy, int mpi_rank,
+                    int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run);
 
 void write_output_serial(struct engine* e, const char* baseName,
                          struct UnitSystem* us, int mpi_rank, int mpi_size,

src/single_io.c

@@ -322,6 +322,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @param Ngas (output) number of Gas particles read.
  * @param Ngparts (output) The number of #gpart read.
  * @param periodic (output) 1 if the volume is periodic, 0 if not.
+ * @param flag_entropy 1 if the ICs contained Entropy in the InternalEnergy
+ * field
  * @param dry_run If 1, don't read the particle. Only allocates the arrays.
  *
  * Opens the HDF5 file fileName and reads the particles contained
@@ -334,7 +336,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  */
 void read_ic_single(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                    int* periodic, int dry_run) {
+                    int* periodic, int* flag_entropy, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -368,6 +370,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
   if (h_grp < 0) error("Error while opening file header\n");
 
   /* Read the relevant information and print status */
+  readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
   readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
   readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
   readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -549,6 +552,7 @@ void write_output_single(struct engine* e, const char* baseName,
   double MassTable[NUM_PARTICLE_TYPES] = {0};
   writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
   unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+  flagEntropy[0] = writeEntropyFlag();
   writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                  NUM_PARTICLE_TYPES);
   writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);

src/single_io.h

@@ -31,7 +31,7 @@
 
 void read_ic_single(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ndm,
-                    int* periodic, int dry_run);
+                    int* periodic, int* flag_entropy, int dry_run);
 
 void write_output_single(struct engine* e, const char* baseName,
                          struct UnitSystem* us);