cosmotool/loadRamses.cpp

#include <cmath>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <fstream>
#include "loadRamses.hpp"
#include "load_data.h"

using namespace std;

#if 0
#define read_buf(b, n) \
{ \
  int k; \
  control_size -= n; \
  for (k = (n-1); k >= 0; k--) \
    infile.read(&b[k],1); \
}
#else
#define read_buf(b, n) \
{ \
  int k; \
  control_size -= n; \
  for (k = 0; k < n; k++) \
    infile.read(&b[k], 1); \
}
#endif

#define read_int(i) \
{ \
  char *o = (char*)&(i); \
  read_buf(o, 4); \
}

#define read_real(f) \
{ \
  char *o = (char*)&(f); \
  read_buf(o, 4); \
}

#define read_characters(c, n) { \
  int k; \
  control_size -= n; \
  infile.read(c, n); \
}

#define push_dummy_control(id) \
{ int control_size = 0;

#define pop_dummy_control() }

#define push_control(id) \
{ \
   int control_size = 0; \
   int control_size2 = 0; \
  char *intbuf = (char*)&control_size; \
  read_int(control_size);

#define pop_control(id) \
  assert(control_size == 0); \
  read_int(control_size2); \
}

GadgetData *CosmoTool::loadRamses(const char *name)
{
  GadgetData *gd = (GadgetData *)malloc(sizeof(GadgetData));
  int id = 1;
  uint32_t totPart = 0;

  cout << "Detecting number of files and particles..." << endl;
  while (1)
    {
      ostringstream ss_fname;  
      ss_fname << name << setfill('0') << setw(5) << id;
      string fname = ss_fname.str();

      cout <<  " ... " << fname << endl;
      
      ifstream infile(fname.c_str());
      if (!infile)
	break;

      int nCpu, ndim, nPar;

      push_control(CPU);
      read_int(nCpu);
      pop_control(CPU);

      push_control(DIM);
      read_int(ndim);
      pop_control(DIM);

      push_control(NPAR);
      read_int(nPar);
      pop_control(NPAR);

      cout << "    NUMCPU=" << nCpu <<  " NUMDIM=" << ndim << " NumPart=" << nPar << endl;

      totPart += nPar;
      id++;
    }

  assert (totPart <= ((~(size_t)0)/sizeof(ParticleState)));
  size_t memSize = sizeof(ParticleState)*(size_t)totPart;
  cout << " Needing " << memSize / 1024 << " kBytes" << endl;
  gd->particles = (ParticleState *)malloc(memSize);
  assert(gd->particles != 0);
  gd->NumPart = totPart;
  gd->ntot_withmasses = totPart;
  gd->header.npart[0] = totPart;
  for (int k = 1; k < 6; k++)
    {
      gd->header.npart[k] = 0;
      gd->header.mass[k] = 0;
    }
  gd->header.num_files = id;
  gd->header.BoxSize = 200.0 * 1000; // kPc
  gd->header.Omega0 = 0.30; // ????
  gd->header.OmegaLambda = 0.70; // ????
  gd->header.HubbleParam = 0.70; // ????

  cout << " Total number part=" << totPart << endl
       << "Loading particles ..." << endl;

  uint32_t curPos = 0;
  id = 1;
  while (1)
    {
      ostringstream ss_fname;  
      ss_fname << name << setfill('0') << setw(5) << id;
      string fname = ss_fname.str();

      cout <<  " ... " << id;
      cout.flush();
      
      ifstream infile(fname.c_str());
      if (!infile)
	break;

      int nCpu, ndim, nPar;

      push_control(CPU);
      read_int(nCpu);
      pop_control(CPU);

      push_control(DIM);
      read_int(ndim);
      pop_control(DIM);

      push_control(NPAR);
      read_int(nPar);
      pop_control(NPAR);

      ParticleState *s = &gd->particles[curPos];

      for (uint32_t i = nPar; i > 0; i--,s++)
	s->Mass = 1.0;

      s = &gd->particles[curPos];
      push_control(POSX);
      for (uint32_t i = nPar; i > 0; i--)
	{
	  read_real(s->Pos[0]);
	  s->Pos[0] *= gd->header.BoxSize;
	  s++;
	}
      pop_control(POSX);

      s = &gd->particles[curPos];
      push_control(POSY);
      for (uint32_t i = nPar; i > 0; i--)
	{
	  read_real(s->Pos[1]);
	  s->Pos[1] *= gd->header.BoxSize;
	  s++;
	}
      pop_control(POSY);

      s = &gd->particles[curPos];
      push_control(POSZ);
      for (uint32_t i = nPar; i > 0; i--)
	{
	  read_real(s->Pos[2]);
	  s->Pos[2] *= gd->header.BoxSize;
	  s++;
	}
      pop_control(POSZ);

      push_control(VELX);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(gd->particles[curPos+i].Vel[0]);
	}
      pop_control(VELX);

      push_control(VELY);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(gd->particles[curPos+i].Vel[1]);
	}
      pop_control(VELY);

      push_control(VELZ);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(gd->particles[curPos+i].Vel[2]);
	}
      pop_control(VELZ);

      curPos += nPar;

      id++;
    }

  cout << endl;

  return gd;
}

PurePositionData *CosmoTool::loadRamsesPosition(const char *name)
{
  PurePositionData *gd = (PurePositionData *)malloc(sizeof(PurePositionData));
  int id = 1;
  uint32_t totPart = 0;
  int nCpu = 0;

  cout << "Detecting number of files and particles..." << endl;
  while (1)
    {
      ostringstream ss_fname;  
      ss_fname << name << setfill('0') << setw(5) << id;
      string fname = ss_fname.str();

      cout <<  " ... " << fname << endl;
      
      ifstream infile(fname.c_str());
      if (!infile)
	break;

      int ndim, nPar;

      push_control(CPU);
      read_int(nCpu);
      pop_control(CPU);

      push_control(DIM);
      read_int(ndim);
      pop_control(DIM);

      push_control(NPAR);
      read_int(nPar);
      pop_control(NPAR);

      cout << "    NUMCPU=" << nCpu <<  " NUMDIM=" << ndim << " NumPart=" << nPar << endl;

      totPart += nPar;
      id++;
    }

  assert (totPart <= ((~(size_t)0)/sizeof(FCoordinates)));
  size_t memSize = sizeof(FCoordinates)*(size_t)(totPart+totPart/nCpu);
  cout << " Needing " << memSize / 1024 << " kBytes" << endl;
  gd->pos = (FCoordinates *)malloc(sizeof(FCoordinates)*totPart);
  assert(gd->pos != 0);
  gd->NumPart = totPart;
  gd->BoxSize = 1000.0; 

  cout << " Total number part=" << totPart << endl
       << "Loading particles ..." << endl;

  uint32_t curPos = 0;
  id = 1;
  while (1)
    {
      ostringstream ss_fname;  
      ss_fname << name << setfill('0') << setw(5) << id;
      string fname = ss_fname.str();
      int *idP;

      cout <<  " ... " << id;
      cout.flush();
      
      ifstream infile(fname.c_str());
      if (!infile)
	break;

      int nCpu, ndim, nPar;

      push_control(CPU);
      read_int(nCpu);
      pop_control(CPU);

      push_control(DIM);
      read_int(ndim);
      pop_control(DIM);

      push_control(NPAR);
      read_int(nPar);
      pop_control(NPAR);

      FCoordinates *s = new FCoordinates[nPar];

      push_control(POSX);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(s[i][0]);
	  s[i][0] *= gd->BoxSize;
	}
      pop_control(POSX);

      push_control(POSY);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(s[i][1]);
	  s[i][1] *= gd->BoxSize;
	}
      pop_control(POSY);

      push_control(POSZ);
      for (uint32_t i = 0; i < nPar; i++)
	{
	  read_real(s[i][2]);
	  s[i][2] *= gd->BoxSize;
	}
      pop_control(POSZ);

      // SKIP VELOCITIES
      for (int k = 0; k < 3; k++) {	
	push_control(V);
	for (uint32_t i = nPar; i > 0; i--)
	  {
	    float dummyF;
	    read_real(dummyF);
	  }
	pop_control(V);
      }
      
      float minMass = INFINITY;
      push_control(MASS);
      for (uint32_t i = nPar; i > 0; i--)
	{
	  float dummyF;
	  read_real(dummyF);
	  if (dummyF < minMass) minMass = dummyF;
	}
      pop_control(MASS);
      
      push_control(ID); 
      for (uint32_t i = 0; i < nPar; i++)
	{
	  int id;
	  read_int(id);
	  id--;
	  assert(id >= 0);
	  assert(id < totPart);
	  memcpy(&gd->pos[id], &s[i], sizeof(FCoordinates));
	}
      pop_control(ID);
      
      delete[] s;
      
      curPos += nPar;

      id++;
    }

  cout << endl;

  return gd;
}

#if 0
int main(int argc, char **argv)
{
  GadgetData *gd = loadRamses(argv[1]);
}

#endif

#if 0
int main(int argc, char **argv)
{
  GadgetData *gd = loadRamses(argv[1]);
}

#endif