save.c

#include <math.h>
#include <stdlib.h>

#include "allvars.h"
#include "proto.h"
#include "power.h"
#include "save.h"

void save_local_data(void)__attribute__((visibility ("hidden")));
int compare_type(const void *a, const void *b)__attribute__((visibility ("hidden")));

void write_particle_data(void)
{
  int nprocgroup, groupTask, masterTask;

  if(ThisTask == 0)
    printf("\nwriting initial conditions... \n");


  if(NTask < NumFilesWrittenInParallel)
    {
      printf
		("Fatal error.\nNumber of processors must be a smaller or equal than `NumFilesWrittenInParallel'.\n");
      FatalError(24131);
    }


  nprocgroup = NTask / NumFilesWrittenInParallel;

  if((NTask % NumFilesWrittenInParallel))
    nprocgroup++;

  masterTask = (ThisTask / nprocgroup) * nprocgroup;


  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))	/* ok, it's this processor's turn */
		save_local_data();
	  
      /* wait inside the group */
      MPI_Barrier(MPI_COMM_WORLD);
    }

  if(ThisTask == 0)
    printf("done with writing initial conditions.\n");
}




void save_local_data(void)
{
#define BUFFER 10
  size_t bytes;
  float *block;
#ifdef USE_64BITID
  long long *blocklongid;
#else
  int *blockid;
#endif
  
  int blockmaxlen, maxlongidlen;
  int32_t dummy;
  FILE *fd;
  char buf[300];
  int i, k, pc;
#ifdef PRODUCEGAS
  double meanspacing, shift_gas, shift_dm;
#endif

  if(NumPart == 0)
    return;

  if(NTaskWithN > 1)
    sprintf(buf, "%s/%s.%d", OutputDir, FileBase, ThisTaskFileNumber);
  else
    sprintf(buf, "%s/%s", OutputDir, FileBase);

  if(!(fd = fopen(buf, "w")))
    {
      printf("Error. Can't write in file '%s'\n", buf);
      FatalError(10);
    }

  for(i = 0; i < 6; i++)
    {
      header.npart[i] = 0;
      header.npartTotal[i] = 0;
      header.mass[i] = 0;
	  header.npartTotalHighWord[i] = 0;
    }


#ifdef MULTICOMPONENTGLASSFILE
  qsort(P, NumPart, sizeof(struct part_data), compare_type);  /* sort particles by type, because that's how they should be stored in a gadget binary file */

#if defined(PRODUCE_CONSISTENT_IDS)
  /* Code from Greg Poole for generating consistent particle IDs across different resolution simulations */
  for(i = 0; i < 3; i++)
    header.npartTotal[i] = header1.npartTotal[i + 1] * (GlassTileFac/GlassTileFacSample) * (GlassTileFac/GlassTileFacSample) * (GlassTileFac/GlassTileFacSample);
#else
  /* Original code from 2LPTic*/
  for(i = 0; i < 3; i++)
    header.npartTotal[i] = header1.npartTotal[i + 1] * GlassTileFac * GlassTileFac * GlassTileFac;
#endif

  for(i = 0; i < NumPart; i++)
    header.npart[P[i].Type]++;

  if(header.npartTotal[0])
    header.mass[0] =
      (OmegaBaryon) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box, 3) / (header.npartTotal[0]);

  if(header.npartTotal[1])
    header.mass[1] =
      (Omega - OmegaBaryon - OmegaDM_2ndSpecies) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box,
																							3) /
      (header.npartTotal[1]);

  if(header.npartTotal[2])
    header.mass[2] =
      (OmegaDM_2ndSpecies) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box, 3) / (header.npartTotal[2]);


#else

  header.npart[1] = NumPart;
  header.npartTotal[1] = TotNumPart;
  //MS - commented out - this is for older versions of Gadget. The new vanilla Gadget versions expect to see the HighWord to be set
  /* header.npartTotal[2] = (TotNumPart >> 32); */
  header.npartTotalHighWord[1] = (TotNumPart >> 32);
  header.mass[1] = (Omega) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box, 3) / TotNumPart;

#ifdef  PRODUCEGAS
  header.npart[0] = NumPart;
  header.npartTotal[0] = TotNumPart;
  header.mass[0] = (OmegaBaryon) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box, 3) / TotNumPart;
  header.mass[1] = (Omega - OmegaBaryon) * 3 * Hubble * Hubble / (8 * PI * G) * pow(Box, 3) / TotNumPart;
#endif
#endif

#ifndef USE_64BITID
  if(TotNumPart > INT_MAX) {
	printf("On task %d: Can not represent particles ID = %lld with 32 bit integers. Enabled the compile time options USE_64BITID\n", ThisTask, TotNumPart);
	FatalError(10000);
  }
#endif

  header.time = InitTime;
  header.redshift = 1.0 / InitTime - 1;

  header.flag_sfr = 0;
  header.flag_feedback = 0;
  header.flag_cooling = 0;
  header.flag_stellarage = 0;
  header.flag_metals = 0;

  header.num_files = NTaskWithN;

  header.BoxSize = Box;
  header.Omega0 = Omega;
  header.OmegaLambda = OmegaLambda;
  header.HubbleParam = HubbleParam;

  header.flag_stellarage = 0;
  header.flag_metals = 0;
  /* header.hashtabsize = 0; */

  dummy = sizeof(header);
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
  my_fwrite(&header, sizeof(header), 1, fd);
  my_fwrite(&dummy, sizeof(dummy), 1, fd);

#ifdef PRODUCEGAS
  meanspacing = Box / pow(TotNumPart, 1.0 / 3);
  shift_gas = -0.5 * (Omega - OmegaBaryon) / (Omega) * meanspacing;
  shift_dm = +0.5 * OmegaBaryon / (Omega) * meanspacing;
#endif

  
  bytes = BUFFER * 1024 * 1024;
  block = malloc(bytes);
  if(block == NULL)
    {
      printf("failed to allocate memory for `block' (%g bytes).\n", (double)bytes);
      FatalError(24);
    }

  blockmaxlen = bytes / (3 * sizeof(float));
  maxlongidlen = bytes / sizeof(long long);
#ifdef USE_64BITID
  blocklongid = (long long *) block;
#else
  blockid = (int *) block;
#endif


  /* write coordinates */
  dummy = sizeof(float) * 3 * NumPart;
#ifdef  PRODUCEGAS
  dummy *= 2;
#endif
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
  for(i = 0, pc = 0; i < NumPart; i++)
    {
      for(k = 0; k < 3; k++)
		{
		  block[3 * pc + k] = P[i].Pos[k];
#ifdef  PRODUCEGAS
		  block[3 * pc + k] = periodic_wrap(P[i].Pos[k] + shift_gas);
#endif
		}

      pc++;

      if(pc == blockmaxlen)
		{
		  my_fwrite(block, sizeof(float), 3 * pc, fd);
		  pc = 0;
		}
    }
  if(pc > 0)
    my_fwrite(block, sizeof(float), 3 * pc, fd);

#ifdef  PRODUCEGAS
  for(i = 0, pc = 0; i < NumPart; i++)
    {
      for(k = 0; k < 3; k++)
		{
		  block[3 * pc + k] = periodic_wrap(P[i].Pos[k] + shift_dm);
		}

      pc++;

      if(pc == blockmaxlen)
		{
		  my_fwrite(block, sizeof(float), 3 * pc, fd);
		  pc = 0;
		}
    }
  if(pc > 0)
    my_fwrite(block, sizeof(float), 3 * pc, fd);
#endif

  my_fwrite(&dummy, sizeof(dummy), 1, fd);



  /* write velocities */
  dummy = sizeof(float) * 3 * NumPart;
#ifdef  PRODUCEGAS
  dummy *= 2;
#endif
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
  for(i = 0, pc = 0; i < NumPart; i++)
    {
      for(k = 0; k < 3; k++)
		block[3 * pc + k] = P[i].Vel[k];

#ifdef MULTICOMPONENTGLASSFILE
      if(WDM_On == 1 && WDM_Vtherm_On == 1 && P[i].Type == 1)
		add_WDM_thermal_speeds(&block[3 * pc]);
#else
#ifndef PRODUCEGAS
      if(WDM_On == 1 && WDM_Vtherm_On == 1)
		add_WDM_thermal_speeds(&block[3 * pc]);
#endif
#endif

      pc++;

      if(pc == blockmaxlen)
		{
		  my_fwrite(block, sizeof(float), 3 * pc, fd);
		  pc = 0;
		}
    }
  if(pc > 0)
    my_fwrite(block, sizeof(float), 3 * pc, fd);
#ifdef PRODUCEGAS
  for(i = 0, pc = 0; i < NumPart; i++)
    {
      for(k = 0; k < 3; k++)
		block[3 * pc + k] = P[i].Vel[k];

      if(WDM_On == 1 && WDM_Vtherm_On == 1)
		add_WDM_thermal_speeds(&block[3 * pc]);

      pc++;

      if(pc == blockmaxlen)
		{
		  my_fwrite(block, sizeof(float), 3 * pc, fd);
		  pc = 0;
		}
    }
  if(pc > 0)
    my_fwrite(block, sizeof(float), 3 * pc, fd);
#endif
  my_fwrite(&dummy, sizeof(dummy), 1, fd);


  /* write particle ID */
#ifdef USE_64BITID
  dummy = sizeof(long long) * NumPart;
#else
  dummy = sizeof(int) * NumPart;
#endif


#ifdef  PRODUCEGAS
  dummy *= 2;
#endif
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
  for(i = 0, pc = 0; i < NumPart; i++)
    {
#ifdef USE_64BITID
      blocklongid[pc] = P[i].ID;
#else
      blockid[pc] = P[i].ID;
#endif

      pc++;

      if(pc == maxlongidlen)
		{
#ifdef USE_64BITID
		  my_fwrite(blocklongid, sizeof(long long), pc, fd);
#else
		  my_fwrite(blockid, sizeof(int), pc, fd);
#endif
		  pc = 0;
		}
    }
  if(pc > 0)
    {
#ifdef USE_64BITID
      my_fwrite(blocklongid, sizeof(long long), pc, fd);
#else
	  my_fwrite(blockid, sizeof(int), pc, fd);
#endif
    }

#ifdef PRODUCEGAS
  for(i = 0, pc = 0; i < NumPart; i++)
    {
#ifdef USE_64BITID
      blocklongid[pc] = P[i].ID + TotNumPart;
#else
	  /*Check that IDs can fit. Both P[i].ID and TotNumPart are 64 bit integers
		so they should not overflow. */
	  if((P[i].ID + TotNumPart) > INT_MAX) {
		printf("On task %d: Can not represent %d particle ID = %lld with 32 bit integers. Enabled the compile time options USE_64BITID\n", 
			   ThisTask, i, P[i].ID);
		FatalError(10000);
	  }
	  
      blockid[pc] = P[i].ID + TotNumPart;
#endif

      pc++;

      if(pc == maxlongidlen)
		{
#ifdef USE_64BITID
		  my_fwrite(blocklongid, sizeof(long long), pc, fd);
#else
		  my_fwrite(blockid, sizeof(int), pc, fd);
#endif
		  pc = 0;
		}
    }
  if(pc > 0)
    {
#ifdef USE_64BITID
      my_fwrite(blocklongid, sizeof(long long), pc, fd);
#else
      my_fwrite(blockid, sizeof(int), pc, fd);
#endif
    }
#endif

  my_fwrite(&dummy, sizeof(dummy), 1, fd);





  /* write zero temperatures if needed */
#ifdef  PRODUCEGAS
  dummy = sizeof(float) * NumPart;
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
  for(i = 0, pc = 0; i < NumPart; i++)
    {
      block[pc] = 0;

      pc++;

      if(pc == blockmaxlen)
		{
		  my_fwrite(block, sizeof(float), pc, fd);
		  pc = 0;
		}
    }
  if(pc > 0)
    my_fwrite(block, sizeof(float), pc, fd);
  my_fwrite(&dummy, sizeof(dummy), 1, fd);
#endif


  /* write zero temperatures if needed */
#ifdef  MULTICOMPONENTGLASSFILE
  if(header.npart[0])
    {
      dummy = sizeof(float) * header.npart[0];
      my_fwrite(&dummy, sizeof(dummy), 1, fd);

      for(i = 0, pc = 0; i < header.npart[0]; i++)
		{
		  block[pc] = 0;

		  pc++;

		  if(pc == blockmaxlen)
			{
			  my_fwrite(block, sizeof(float), pc, fd);
			  pc = 0;
			}
		}
      if(pc > 0)
		my_fwrite(block, sizeof(float), pc, fd);
      my_fwrite(&dummy, sizeof(dummy), 1, fd);
    }
#endif



  free(block);

  fclose(fd);
}


/* This catches I/O errors occuring for my_fwrite(). In this case we better stop.
 */
size_t my_fwrite(void *ptr, size_t size, size_t nmemb, FILE * stream)
{
  size_t nwritten;

  if((nwritten = fwrite(ptr, size, nmemb, stream)) != nmemb)
    {
      printf("I/O error (fwrite) on task=%d has occured.\n", ThisTask);
      fflush(stdout);
      FatalError(777);
    }
  return nwritten;
}


/* This catches I/O errors occuring for fread(). In this case we better stop.
 */
size_t my_fread(void *ptr, size_t size, size_t nmemb, FILE * stream)
{
  size_t nread;

  if((nread = fread(ptr, size, nmemb, stream)) != nmemb)
    {
      printf("I/O error (fread) on task=%d has occured.\n", ThisTask);
      fflush(stdout);
      FatalError(778);
    }
  return nread;
}


#ifdef MULTICOMPONENTGLASSFILE
int compare_type(const void *a, const void *b)
{
  if(((struct part_data *) a)->Type < (((struct part_data *) b)->Type))
    return -1;

  if(((struct part_data *) a)->Type > (((struct part_data *) b)->Type))
    return +1;

  return 0;
}
#endif